This article is mentioned in:

Introduction

Medullary thyroid carcinoma (MTC) arises from the calcitonin-producing parafollicular C cells of the thyroid and was first described by Hazard in 1959 (1). MTC comprises 5–10% of all primary thyroid malignancies (2,3). MTC is mainly sporadic (termed SMTC) and only 20–30% of cases are hereditary MTC (HMTC) (4–6). Activating mutations in the RET proto-oncogene are responsible for medullary thyroid carcinoma, and these mutations are believed to be associated with a more persistent disease and a lower overall survival (7–9). However, to date, the molecular mechanisms of MTC carcinogenesis remain unclear.

MicroRNAs (miRNAs or miRs) are a large subgroup of non-coding RNAs, 18–25 nucleotides in length, that are evolutionary conserved. These molecules control post-transcriptional gene expression through the inhibition of mRNA translation or induction of its degradation (10) miRNAs can act as oncogenes, or tumor suppressor genes and can be used as diagnostic and predictive biomarkers, and even for the treatment of diseases, including MTC (11,12). Shabani et al reported that a high expression of hsa-miR-144 and hsa-miR-183miR-34a could be considered as biomarkers of MTC (13). Furthermore, in a previous study, hsa-miR-375 was upregulated and hsa-miR-9* was downregulated in SMTC vs. HMTC, and the overexpression of hsa-miR-183 and hsa-miR-375 in MTC predicted lateral lymph node metastases (14).

Although advances have been made in understanding the mechanisms of MTC, studies on these are limited and further confirmation is required. This study aimed to identify biomarkers by analyzing differentially expressed miRNAs (DEMs) and differentially expressed genes (DEGs) between MTC and normal thyroid tissues. Subsequently, Gene Ontology (GO) and pathway enrichment analysis, protein-protein interaction (PPI) and transcription factor (TF)-DEMs-genes networks were conducted to identify key miRNAs and elucidate the potential molecular mechanisms of MTC.

Materials and methods

Data source

The original datasets comparing the gene expression profiles between MTC and normal thyroid tissue were downloaded from the NCBI GEO databases. The accession numbers were GSE97070 (15), GSE40807 (16) and GSE27155 (17,18). The microarray data of GSE97070 were based on GPL18402 (Agilent-046064 Unrestricted_Human_miRNA_V19.0_Microarray), and included 8 MTC samples, 9 lymph node metastasis samples and 3 normal samples. GSE40807 was based on GPL18227 (Agilent-019118 Human miRNA Microarray 2.0 G4470B), which consisted of 40 MTC samples and 40 normal samples. The platform of GSE27155 was GPL96, [HG-U133A] Affymetrix Human Genome U133A Array, including 96 samples with 2 MTC and 4 normal samples. Platform and series matrix file(s) were downloaded as TXT files.

Data pre-processing and DEM analysis

A comparison between the 2 sample groups, MTC and lymph node metastasis vs. normal thyroid tissues, was performed in each GEO dataset to identify DEMs and DEGs. The R software package was used to process the downloaded files and to convert and reject the unqualified data, and the limma R package was then used to identify DEMs and DEGs. Before we decided to conduct this study, we carried out research on bioinformatics analysis, and found that the ‘false discovery rate (FDR) <0.05, |log2FC|>1’ was a common criteria for screening DEMs or DEGs (19–21). Thus, in this study, samples with the criteria mentioned above were considered DEMs or DEGs. The TXT results were preserved for subsequent analysis. Heat maps of DEMs and DEGs were generated using FunRich 3.1.3 software.

Identification of miRNA targets

starBase (http://starbase.sysu.edu.cn/index.php) provides certain miRNA-target regulatory association pairs, which are verified by experiments and predicted by 7 programs, including microT, miRanda, miRmap, PITA, RNA22, PicTar and TargetScan (22). StarBase v3.0 identifies >1.1 million miRNA-ncRNA, 2.5 million miRNA-mRNA, 2.1 million RBP-RNA and 1.5 million RNA-RNA interactions from multi-dimensional sequencing data. In this study, miRNA-target gene regulatory association pairs were verified according to the following standards: CLIP Data (low stringency), Degradome Data (low stringency), Pan-cancer (1 cancer type) and Program Number (3 programs).

GO and kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses of target genes

GO is a common method for annotating genes, gene products and sequences to underlying biological phenomena; KEGG is an integrated database resource for the biological interpretation of genome sequences and other high-throughput data (23). Metascape (http://metascape.org) is an online program that aims to develop a set of reliable, productive and intuitive tools that help biomedical research community to analysis gene/protein lists and make better data-driven decisions (24). In this study, Metascape was used to perform GO and KEGG pathway analysis on target genes of DEMs. In addition, GO terms consisted of 3 aspects: Biological process (BP), cellular component (CC) and molecular function (MF).

PPI network construction and analysis of modules

NetworkAnalyst (https://www.networkanalyst.ca/) is a series of web-based tools for statistical meta-analysis, visual data mining and data integration, through the rapid generation of biological networks. It supports the meta-analysis of gene lists, and data integration is achieved through robust statistical procedures and subsequently visually examined within PPI networks (25,26). In this study, STRING was selected as a PPI database. To access more objective and reliable results, this study restricted the sources requiring experimental evidence and the cut-off score was set at a high confidence (900). Nodes with more degrees were considered as hub genes and may serve as core proteins or key candidates with important physiological regulatory functions. The pathway enrichment analysis of genes in the modules was performed, and P<0.05 was considered to indicate a statistically significant difference.

TF-DEMs-target gene and DEGs regulatory network construction

TransmiR (http://www.cuilab.cn/transmir) is a database for TF-miR regulations, through which regulatory associations between TFs and miRNAs can be identified. To date, TransmiR v2.0 contains 3,730 literature-curated TF-miRNA regulations from 1,349 publications and 1,785,998 TF-miRNA regulations derived from ChIP-seq evidence (27,28). This study inputted the overlapped DEMs into the database to examine the regulatory association pairs between TFs and DEMs. The inclusion criteria of evidence were supported by high-throughput experiments from literature.

Based on the data this study obtained, TF-DEMs-target gene and DEGs regulatory network was constructed and visualized by Cytoscape 3.6.1 software to show the overlapped TFs, DEMs, target genes and DEGs. Therefore, these TFs, DEGs, DEMs and target genes may play a potential role in the pathogenesis and treatment of MTC.

Results

Microarray data information and identification of DEMs and DEGs

The MTC expression microarray datasets GSE97070, GSE40807 and GSE27155 were standardized. The datasets were subsequently screened using the limma package (FDR<0.05, |log2FC|>1); 57 DEMs were obtained in GSE97070. Among these, 29 upregulated and 28 downregulated DEMs were identified. Overall, 134 DEMs were screened from the GSE40807 dataset, including 70 upregulated genes and 64 downregulated genes. Among these, hsa-miR-375, hsa-miR-127-3p and hsa-miR-429 were significantly upregulated in both the GSE97070 and GSE40807 datasets, whereas the expression of hsa-miR-199b-5p and hsa-miR-199a-3p was downregulated. In addition, 235 DEGs were screened from the GSE27155 dataset, including 135 upregulated genes and 100 downregulated genes (Fig. 1 and Table SI). Heatmaps were generated based on the expression levels, where each column represented a biological sample and each row in the heat map represented a DEM or DEG. The color indicated the relative expression levels of miRNA in tissue specimens (Figs. S1-S3).

Figure 1. Differential expression of data between two sets of samples. The red points represent upregulated miRNAs or genes and green points represent the downregulation of the expression of miRNAs or genes.

Target genes of DEMs

The target genes of DEMs were identified according to the standards described above. For the 5 commonly altered miRNAs, a total of 300 target genes were obtained, including 202 genes of upregulated DEMs and 98 genes of downregulated DEMs (Table SII).

SNTB2 and MED13 were predicted as the potential targets of hsa-miR-375 and hsa-miR-429. A total of 4 genes were the potential targets of hsa-miR-429 and hsa-miR-199a-3p, including SCD, CFL2, IREB2 and DTNA. In addition, 5 genes, including CITED2, DUSP1, FHL1, LDHB and C1orf115, of 100 downregulated DEGs were potentially targeted by DEMs (Figs. 2 and 3).

Figure 2. Venn diagram of overlapping genes between target genes of DEMs and DEGs. DEMs, differentially expressed miRNAs; DEGs, differentially expressed genes.

Figure 3. TF-DEMs-target genes and DEGs regulatory network. Red represents upregulated DEMs, green represents downregulated DEMs, dark blue represents hub genes, light blue represent TFs, yellow represents overlapped genes between DEGs and target genes of DEMs, pink represents overlapped genes between hsa-miR-429 and hsa-miR-199a-3p; grey represents overlapped genes between hsa-miR-429 and hsa-miR-375; red arrows indicate activation and black arrows indicate repression. TF, transcription factor; DEMs, differentially expressed miRNAs; DEGs, differentially expressed genes.

Significant functions and pathway enrichment analysis

Metascape was used to analyze the downstream target genes of the five common DEMs. For each given gene list, a pathway and process enrichment analysis was carried out with the following ontology sources: KEGG pathway, GO BP, GO CC and GO MF. All genes in the genome were used as the enrichment background. Terms with a P-value <0.01, a minimum count of 3, and an enrichment factor >1.5 were collected and grouped into clusters based on their membership similarities. The enrichment factor is the ratio between the observed counts and the counts expected by chance. For GO term enrichment analysis, the target genes of upregulated DEMs were mainly enriched in adherens junction, anchoring junction and cell-substrate junction assembly, while the target genes of downregulated DEMs were mainly enriched in non-canonical Wnt signaling pathway, heat shock protein binding and regulation of small molecule metabolic process. KEGG analysis revealed that the target genes of upregulated DEMs were mostly enriched in focal adhesion, regulation of actin cytoskeleton and adherens junction. The target genes of downregulated DEMs were enriched in RNA transport, the lysosome and mRNA surveillance pathway (Table I and Fig. 4).

Figure 4. Top 20 clusters with their representative enriched terms (colored by cluster ID). where nodes that share the same cluster ID are typically close to each other.

Table I. Results of top 10 GO terms and top 3 KEGG pathway analysis.

Table I.

Results of top 10 GO terms and top 3 KEGG pathway analysis.

Category	Term	Description	Gene counts	LogP
Target genes of upregulated DEMs
GO:0005912	CC	Adherens junction	20	−7.110
GO:0070161	CC	Anchoring junction	20	−6.902
GO:0007044	CC	Cell-substrate junction assembly	9	−6.880
GO:0051347	BP	Positive regulation of transferase activity	21	−6.251
GO:0005925	CC	Focal adhesion	16	−6.123
GO:0005924	CC	Cell-substrate adherens junction	16	−6.080
GO:0030055	CC	Cell-substrate junction	16	−6.024
GO:0019900	MF	Kinase binding	22	−5.909
GO:0034329	BP	Cell junction assembly	12	−5.872
GO:0019904	MF	Protein domain specific binding	21	−5.842
hsa04510	KEGG	Focal adhesion	12	−6.581
hsa04810	KEGG	Regulation of actin cytoskeleton	10	−4.637
hsa04520	KEGG	Adherens junction	6	−4.333
Target genes of downregulated DEMs
GO:0035567	BP	Non-canonical Wnt signaling pathway	5	−4.074
GO:0031072	MF	Heat shock protein binding	5	−3.724
GO:0062012	BP	Regulation of small Molecule metabolic process	7	−3.187
GO:0001738	BP	Morphogenesis of a polarized epithelium	4	−3.151
GO:0016055	BP	Wnt signaling pathway	8	−3.042
GO:0198738	BP	cell-cell signaling by wnt	8	−3.030
GO:0019904	MF	Protein domain specific binding	10	−3.026
GO:0005793	CC	Endoplasmic reticulum-Golgi intermediate compartment	4	−2.684
GO:0060071	BP	Wnt signaling pathway, planar cell polarity pathway	3	−2.599
GO:0090175	BP	Regulation of establishment of planar polarity	3	−2.560
hsa03013	KEGG	RNA transport	5	−3.037
hsa04142	KEGG	Lysosome	4	−2.684
hsa03015	KEGG	mRNA surveillance pathway	3	−2.132

[i] GO, gene ontology; BP, biological process, CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEMs, differentially expressed miRNAs.

PPI network construction and analysis of modules

The PPI network was constructed and visualized by NetworkAnalyst. The minimum network was used to keep seed proteins, as well as essential non-seed proteins, that were suitable for simplifying a dense network to study key associations. A zero-order network was subsequently used to keep only seed proteins that directly interact with each other and to visualize modules. Degree >20 was set as the cut-off criterion. A total of 13 genes were chosen as hub genes, such as UBC, EP300, MAPK8, FYN, RHOA, PPP2CA, SMAD2, UBE2D1, HNRNPD, YES1, YWHAG, NR3C1 and TBP. A total of 9 genes were target genes of hsa-miR-429 (EP300, FYN, RHOA, PPP2CA, UBE2D1, HNRNPD, YWHAG, NR3C1 and TBP). In addition, the top 5 modules were selected, and the KEGG pathway enrichment analysis revealed that target genes in these modules were mostly enriched in the neurotrophin signaling pathway, the MAPK signaling pathway, the pathways in cancer, the Wnt signaling pathway and the focal adhesion (Fig. 5 and Table SIII).

Figure 5. PPI network and the selected module. (A) PPI network from minimum network, red nodes with labels indicated hub genes; (B) top 5 significant module selected from zero-order network, each module was represented by a different color. PPI, protein-protein interaction.

Analysis of TF-DEMs-target genes and DEGs regulatory network

From the data of TransmiR, upregulated DEMs were regulated by 51 TFs, and downregulated DEMs were regulated by 24 TFs. In addition, 17 TFs regulated 2 miRNAs, including upregulated or downregulated DEMs, while CREB1 regulated all upregulated DEMs. Furthermore, EP300 was detected as a TF of hsa-miR-375 and hsa-miR-199b-5p, which was also a target gene of has-miR-429 (Table SIV). In summary, based on the aforementioned results from the bioinformatics analysis, a regulatory network was constructed to indicate the overlapped TFs, DEMs, target genes and DEGs (Fig. 3).

Discussion

MTC is a rare malignancy with poor prognosis, as lymph node metastases are found in 55% patients by the time of diagnosis (29). Surgical resection remains the most effective therapy of this disease, however in advanced cases and patients with distant metastases, this treatment method is not sufficient (30–33). Therefore, it is important to study the molecular mechanisms of the carcinogenesis and development of MTC.

In the present study, a bioinformatics approach was used to identify candidate biomarker and therapeutic targets of MTC. Following the analysis, 191 DEMs, including 99 upregulated DEMs and 92 downregulated DEMs were identified. Among these, hsa-miR-375, hsa-miR-127-3p and hsa-miR-429 were upregulated, and hsa-miR-199a-3p and hsa-miR-199b-5p were downregulated in 2 miRNA profiles, suggesting that they may function as carcinogens or tumor suppressors in MTC. As shown by the OncomiR database, hsa-miR-199a-3p and hsa-miR-199b-5p were upregulated in normal tissues, and hsa-miR-375 was upregulated in thyroid carcinoma, which is in accordance with the analysis of this study. Furthermore, a number of researchers have reported that the overexpression of hsa-miR-375 significantly contributes to the pathophysiology and development of MTC (34–36).

GO and KEGG analysis results showed that target genes of upregulated DEMs were significantly enriched in adherens junction and focal adhesion, which might play important roles in the tumor development and progression. Target genes of downregulated DEMs were mainly involved in Wnt signaling pathway, protein binding and RNA transport. The Wnt signaling pathway is a group of signal transduction pathways, which begins with proteins that pass signals into a cell through cell surface receptors. It is an important signaling pathway in biological development and tumorigenesis (37–39). In the comparison of MTCM918T and MTC634, in addition to MTCM918T and MTCWT, Maliszewska et al reported that many biochemical pathways were involved in the malignant behavior of MTC, including the Wnt pathway (40).

By constructing the PPI network, the present study identified 13 hub genes, and 10 of these, EP300, MAPK8, FYN, RHOA, PPP2CA, SMAD2, UBE2D1, HNRNPD, YES1 and NR3C1 were involved in the top 5 modules. KEGG pathway enrichment analysis of modules showed that Focal adhesion, TGF-β signaling pathway and Wnt signaling pathway contained 4 hub genes respectively. The TGF-β signaling pathway mediates intracellular signaling and participates in embryonic development, tumorigenesis, and physiological processes (41). Furthermore, TGF-β can cause enhanced adhesion and motility of tumor cells (42). Santarpia et al observed the effect of miRNAs in MTC tumorigenesis, migration, proliferation and invasion. The cell lines were treated with miR-200 inhibitor and analysis was performed in accordance to the array data, showing that the members of the miR-200 family regulate the expression of E-cadherin by directly targeting ZEB1,2 and through the enhanced expression of tumor growth factor β-1,2 (43).

The hub genes MAPK8 and RHOA were involved in the 5 of the top 10 pathways, including the neurotrophin signaling pathway, MAPK signaling pathway, colorectal cancer, adherens junction, TGF-beta signaling pathway and Wnt signaling pathway. RHOA encodes a member of the Rho family of small GTPases, which cycle between inactive GDP-bound and active GTP-bound states, and function as a molecular switch in signal transduction cascades. The overexpression of this gene is associated with tumor cell proliferation and metastasis. A number of studies have suggested that RHOA can serve as a biomarker of colorectal cancer, with regards to a therapeutic target (44,45). Takahashi et al conducted an experiment on Rb1(+/-)Nras(+/-) animals. The results of the aforementioned study revealed that distant MTC metastases were associated with the loss of the remaining wild-type Nras allele. In addition the loss of Nras in Rb1-deficient C cells results in an elevated RHOA activity, therefore leading to the malignant behavior of these cells (46). MAPK8 is a member of the MAP kinase family. The activation of this kinase by TNF-α is found to be required for TNF-α induced apoptosis. The importance of the MAPK pathway has been well established in the tumorigenesis of papillary thyroid cancer (47,48). For MTC, Chang et al conducted an exome-wide analysis of the mutational spectrum and indicated that a number of pathways was involved in the variant process, including MAPK pathway (49).

In the constructed TF-DEMs-target genes and DEGs regulatory network, 5 downregulated DEGs (C1orf15, CITED2, DUSP1, FHL1 and LDHB) were also the target genes of 3 upregulated DEMs, namely hsa-miR-375, hsa-miR-127-3p and hsa-miR-429, and the downregulated DEM, hsa-miR-199b-5p. In addition, DUSP1 and FHL1 were involved in the module of PPI network, and DUSP1 participated in the MAPK pathway. The protein encoded by DUSP1 was involved in several cellular processes and resulted in chemotherapy and radiotherapy resistance, which indicated that DUSP1 can serve as a target for cancer therapy. Despite that DUSP1 has been reported as an oncogene, therapeutic target and a biomarker in many different types of tumor (50–52), the study of DUSP1 expression and function in thyroid carcinoma is limited and further investigations are required.

As shown in the regulatory network, hsa-miR-429 and hsa-miR-199a-3p were regulated by TGFB1, a TF which encodes a secreted ligand of the TGF-β superfamily and regulates cell proliferation, differentiation and growth. In addition, TGFB1 has been recognized as an activator of hsa-miR-199a-3p and a repressor of hsa-miR-429, based on the data from TransmiR by the evidence level of literature. According to the results of target prediction, 4 genes, including CFL2, DTNA, IREB2 and SCD, were overlapped between has-miR-429 and has-miR-199a-3p, whose effects have been reported in tumors for their differential expression or potential function in diagnosis and treatment (53–57). However, to the best of our knowledge, CFL2, DTNA, IREB2 and SCD has not been mentioned in MTC. In addition, 10 of the 13 hub genes were the target genes of hsa-miR-429 and hsa-miR-199a-3p, including MAPK8 and RHOA. Taken together, the regulatory network of TGFB1, hsa-miR-429/hsa-miR-199a-3p and target genes may play important roles in the development of MTC and warrant further investigation.

Cancer is a complex disease caused by multiple factors and integrated bioinformatics analysis can help in the investigation and understanding of its molecular mechanism. The aim of present study was to identify key DEMs and genes and to determine potential biomarkers to predict the progression of MTC. However, this study presents with a number of limitations. First, the dataset sample size was limited, due to the difficulty to obtain clinical samples. Second, there were only 5 DEMs that overlapped in different GSE chips, suggesting that some valuable miRNAs may be missing. Third, the incidence of MTC was low; thus, studies of how those genes affect the prognosis of MTC were seldom reported. Our team are collecting the clinical and pathological data of MTC in order to carry out further investigations.

In conclusion, this study identified numerous DEMs that may contribute to the initiation and development of MTC. Furthermore, the present study also identified a series of significant pathways and mechanisms for treatment. The regulatory association between TGFB1, hsa-miR-429 and hsa-miR-199a-3p may provide novel insight for the diagnosis and treatment of MTC. In addition, we aim to perform further experiments to examine the expression of the identified DEMS and genes in different sample types, and subsequently confirm their utility in the diagnosis and molecular therapy of MTC.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

QP, DL and MZ contributed to the design of the study. LZ, ML and QP performed the bioinformatics analysis and wrote the manuscript. LZ and DL were responsible for article revision. LZ and ML contributed to language editing and the revision of this manuscript. All authors have read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References