Integrated microRNA-mRNA analyses of distinct expression profiles in follicular thyroid tumors.

MicroRNAs (miRNAs/miRs) are small non-coding RNAs identified in plants, animals and certain viruses; they function in RNA silencing and post-transcriptional regulation of gene expression. miRNAs also serve an important role in the pathogenesis, diagnosis and treatment of tumors. However, few studies have investigated the role of miRNAs in thyroid tumors. In the present study, the expression of miRNA and mRNA was compared between follicular thyroid carcinoma (FTC) and follicular thyroid adenoma (FA) samples, and then miRNA-mRNA regulatory network analysis was performed. Microarray datasets (GSE29315 and GSE62054) were downloaded from the Gene Expression Omnibus, and profiling data were processed with R software. Differentially expressed miRNAs (DEMs) and differentially expressed genes (DEGs) were determined, and Gene Ontology enrichment analysis was subsequently performed for DEGs using the Database for Annotation, Visualization and Integrated Discovery. The target genes of the DEMs were identified with miRWalk, miRecords and TarMir databases. Network analysis of the DEMs and DEMs-targeted DEGs was performed using Cytoscape software. In GSE62054, 23 downregulated and 9 upregulated miRNAs were identified. In GSE29315, 42 downregulated and 44 upregulated mRNAs were identified. A total of 36 miRNA-gene pairs were also identified. Network analysis indicated a co-regulatory association between miR-296-5p, miR-10a, miR-139-5p, miR-452, miR-493, miR-7, miR-137, miR-144, miR-145 and corresponding targeted mRNAs, including TNF receptor superfamily member 11b, benzodiazepine receptor (peripheral) -associated protein 1, and transforming growth factor β receptor 2. These results suggest that miRNA-mRNAs networks serve an important role in the pathogenesis, diagnosis and treatment of FTC and FA.

Introduction

Thyroid follicular cells are found in the thyroid gland, specifically in the epithelial monolayer. In total, >95% of thyroid tumors are derived from these follicular cells (1). In 2016, the incidence of thyroid tumors rose globally, largely due to technological and diagnostic advances (2). However, it remains difficult to distinguish whether a thyroid nodule is benign or malignant. Follicular thyroid tumors may be divided into malignant follicular thyroid carcinoma (FTC) and benign follicular thyroid adenoma (FA). Only 5–10% of thyroid nodules are malignant (3). Patients with follicular tumors usually must undergo thyroid lobectomy for diagnosis, which is often an unnecessary surgery, as the disease is usually benign. Fine-needle aspiration cytology is considered the most accurate method for the diagnosis of FTC and FA (4).

Previously, microRNAs (miRNAs/miRs) have been demonstrated to be involved in the pathogenesis of various diseases, like cancer, diabetes and osteoarthritis (5–7). miRNAs are small (18–25 nucleotides) non-coding, single-stranded RNA molecules that bind to targets in a base pair-mediated manner, resulting in the degradation or inhibition of the expression and function of protein-coding mRNAs. miRNAs often bind to the 3′-untranslated region (3′UTR) of target genes (8), although they are usually only partially complementary to the target (9). miRNAs regulate ~30% of the human genes associated with proliferation, apoptosis, metastasis, cell immunity and differentiation (10). Each miRNA is able to regulate several hundred mRNAs, and each mRNA may be the target of several miRNAs. Therefore, a regulatory control network exists between miRNAs and mRNAs (11). Furthermore, miRNAs have been associated with several types of tumors, including non-small cell lung cancer, colon and esophageal cancer, and FTC (12–14). However, there are few studies of specific miRNA and mRNA analyses of follicular thyroid tumors.

Several microarray studies have already described the differentially expressed genes (DEGs) between malignant and benign thyroid nodules. However, these studies have several restrictions, including the fact that the samples are limited, they contain significant false-negatives, and they require external analysis at an offsite company laboratory (15–17). Certain studies have aimed to reveal the potential miRNAs associated with follicular thyroid tumors (18).

In the present study, an integrated analysis of differentially expressed miRNAs (DEMs) and DEGs between FTC and FA was performed. A Gene Ontology (GO) analysis of the DEGs was performed. A total of 36 miRNA-gene pairs were identified between the DEGs and the target genes of the DEMs. A miRNA-mRNA network analysis was then performed to additionally investigate the pathogenesis of FTC.

Materials and methods

Analysis of mRNA and miRNA profiling datasets

Expression profile datasets containing mRNA and miRNA were acquired from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/). The expression profiling data of GSE29315 (Tomas et al, unpublished) are mRNA profiling data, originally obtained from a cohort of 9 FTC and 17 FA samples. The GSE62054 dataset contains miRNA profiling data, which was originally obtained from 17 FTC and 8 FA samples (19). Additionally, GSE29315 was hybridized on the Affymetrix U95 GeneChip platform (Affymetrix; ThermoFisher Scientific, Inc., Waltham, MA, USA) and GSE62054 was performed on the Illumina Human v2 miRNA expression BeadChip (Illumina, Inc., San Diego, CA, USA).

Preprocessing of profiling data

GSE29315 and GSE62054 data were first preprocessed by the Affy package in R language version 3.4.0 and then were processed by log2 transformation, background correction and data normalization using the Robust Multi-array Average algorithm (20).

Identification of DEMs, DEGs and GO enrichment analysis

Identification of DEMs and DEGs were conducted by the Limma package version 3.32.5 in R software (21). The threshold values for different expression were log2 (fold-change)>0.5 or log2 (fold-change)<-0.5 with P<0.05 (22). GO enrichment analysis for DEGs was performed with the Database for Annotation, Visualization and Integrated Discovery (DAVID) (23).

Overlapping genes of DEGs and the predicted target genes of the DEMs

The predicted target mRNAs of the DEMs were generated using the miRWalk (24), miRecords (25) and TarMir databases (26). The overlapping DEGs and the predicted target mRNAs of the DEMs were identified for additional network analysis.

Construction and analysis of miRNA-mRNA regulatory network

To obtain an improved understanding of the biological function of the miRNA-mRNA regulatory network, node-degree analysis was performed, based on the overlapping genes and their upstream miRNAs. The network was visualized using the Cytoscape platform software version 3.0.1 (27).

Results

Identifying DEMs and DEGs between FTC and FA

GSE29315 and GSE62054 were downloaded from GEO and then normalized, and corrected by the quantile normalization method and hierarchical clustering analysis using R software. DEMs and DEGs were identified between the FTC and FA. A total of 86 DEGs and 32 DEMs were obtained when the threshold values were set at P<0.05 and log2(fold-change)>0.5 or log2 (fold-change)<-0.5. The top 5 downregulated DEMs were miR-7, miR-1179, miR-7-2, miR-486-5p and miR-130b. The top 5 upregulated DEMs were miR-663b, miR-137, miR-30c-1, miR-767-5p and miR-603 (Table I). As for the DEGs, the top 5 downregulated genes were fatty acid binding protein 4 (FABP4), cytidine monophospho-N-acetylneruaminic acid (CMAHP), integral membrane protein 2A (ITM2A), carbonic anhydrase 4 (CA4) and family with sequence similarity 189 member A2 (FAM189A2), and the top 5 upregulated genes were erythrocyte membrane protein band 4.1 like 3 (EPB41L3), secretogranin V (SCG5), paired box 1 (PAX1), methylenetetrahydrofolate dehydrogenase (NADP + dependent) 2, methenyltetrahydrofolate cyclohydrolase (MTHFD2) and cadherin 2 (CDH2) (Table II). A volcano plot was constructed to identify the DEGs (Fig. 1).

Table I.

Top 5 differentially expressed miRNAs of malignant follicular thyroid carcinoma compared with benign follicular thyroid adenoma.

miRNA	P-value	log2(fold-change)
Downregulated
miR-7	0.0041392	−1.7320437
miR-1179	0.0081728	−1.3950195
miR-7–2	0.0006626	−1.2525509
miR-486-5p	0.0412501	−1.0502825
miR-130b	0.0028172	−0.9176468
Upregulated
miR-663b	0.0009353	0.9881272
miR-137	0.0088044	0.9341108
miR-30c-1	0.0059237	0.8695624
miR-767-5p	0.0036048	0.7353497
miR-603	0.0392875	0.6646499

miR/miRNA, microRNA.

Table II.

Top 5 differentially expressed mRNAs of malignant follicular thyroid carcinoma compared with benign follicular thyroid adenoma.

mRNA	P-value	log2(fold-change)
Downregulated
FABP4	0.001621719	−2.100023748
CMAHP	0.024414059	−1.066127774
ITM2A	0.016922069	−1.060957028
CA4	0.003105875	−1.030691864
FAM189A2	0.001993414	−1.000814956
Upregulated
EPB41L3	0.000527208	1.020917517
SCG5	0.044745635	0.990761798
PAX1	0.040281329	0.936795356
MTHFD2	0.002917107	0.870767506
CDH2	0.038829964	0.862357097

Figure 1.

Distribution of differentially expressed genes between the malignant and benign follicular thyroid tumors in GSE29315. Blue triangles indicate 42 significantly downregulated genes and red crosses indicate 44 significantly upregulated genes.

GO enrichment analysis of DEGs

GO analysis of all the DEGs (Table III) identified 8 associated biological processes: Positive regulation of macromolecule metabolic process, regulation of cell motion, cell proliferation, tube development, regulation of response to external stimulus, regulation of locomotion, response to drug and T cell activation (Fig. 2).

Table III.

Differentially expressed mRNAs of malignant follicular thyroid carcinoma compared with benign follicular thyroid adenoma.

Gene	P-value	log2 (fold-change)
FABP4	0.001621719	−2.100023748
CMAHP	0.024414059	−1.066127776
ITM2A	0.016922069	−1.060957028
CA4	0.003105875	−1.030691864
FAM189A2	0.001993414	−1.000814956
MPPED2	0.007206792	−0.981309728
HGD	0.023615869	−0.852452356
TNFRSF11B	0.036234193	−0.836887107
SLC16A4	0.020876545	−0.826437495
BZRAP1	0.006323889	−0.784389089
PDGFRL	0.019724555	−0.769696713
TFF3	0.004160978	−0.758375163
LAMB1	0.019590536	−0.754903323
UBR2	0.001897906	−0.746515077
TGFBR2	0.016942031	−0.736954014
CPQ	0.040123761	−0.722559518
RDX	0.014625862	−0.700582583
PTPRN2	0.008418827	−0.692563257
SALL1	0.039302461	−0.687132893
LRRN3	0.006382367	−0.685656353
TRAM2	0.016638751	−0.664343953
SLITRK5	0.040721395	−0.652918137
RGS16	0.027819411	−0.648828034
STARD13	0.013643606	−0.642603503
THBD	0.010984131	−0.635676793
GJA1	0.035287842	−0.633044868
TNFSF10	0.009577025	−0.620052341
PKIA	0.022943634	−0.617442203
CLDN8	0.014099121	−0.606720582
IL11RA	0.016027339	−0.589253369
CDC27	0.012718038	−0.584208837
SLC35D2	0.010299633	−0.579905832
IFI44L	0.015727697	−0.576791484
NR2F2	0.013187137	−0.575425266
CHD9	0.048840737	−0.573993055
FCGRT	0.011662768	−0.566086321
AGTR1	0.013773676	−0.538397573
CXCL12	0.003888448	−0.538189987
HSD17B8	0.019222754	−0.533271237
LOC728093	0.035060908	−0.524392385
SELE	0.017742851	−0.507241011
AKAP12	0.010628133	−0.506141506
PVALB	0.045712444	0.509393237
NNT	0.027397598	0.513526356
ADORA1	0.036372247	0.513841883
GPI	0.018981531	0.514219746
AREG	0.042790083	0.515091119
ELF4	0.001468499	0.516939945
NAB2	0.011359365	0.527142113
GPC1	0.030827537	0.529427807
SLC25A5	0.029858761	0.534298394
RYR1	0.014981235	0.538436907
PEG10	0.018255692	0.549699307
SLIT2	0.005701378	0.551929669
ESYT1	0.003009396	0.561757654
CRLF1	0.029356956	0.562850973
DPP4	0.018431425	0.567301902
CES2	0.003842477	0.571035541
CYCS	0.030715945	0.585863898
BASP1	0.026758782	0.594217929
KLHL21	0.022505203	0.599053201
LCN2	0.017669691	0.601124995
EEF1A2	0.018392179	0.602560272
TGIF1	0.030147532	0.605396529
TSPYL2	0.046736047	0.617648983
CKS2	0.026759687	0.623716633
SPAG5	0.027923235	0.624868361
CKMT2	0.043578063	0.645621487
ALDH1A3	0.010848438	0.653685247
ZNF148	0.031227695	0.661077991
UCHL1	0.003566593	0.663217388
ASNS	0.006539605	0.667304137
NPTXR	0.009448284	0.674187256
FKBP5	0.027541727	0.694255151
GOT1	0.001282557	0.730143285
IGFBP3	0.009116606	0.771327366
SERPINE2	0.020069598	0.785032654
SOX4	0.006616428	0.808876177
BSG	0.001981905	0.809699733
SCNN1A	0.027983691	0.825462627
NPC2	0.005654992	0.837641446
CDH2	0.038829964	0.862357097
MTHFD2	0.002917107	0.870767506
PAX1	0.040281329	0.936795356
SCG5	0.044745635	0.990761798
EPB41L3	0.000527208	1.020917517

Figure 2.

GO analysis of differentially expressed genes between malignant and benign follicular thyroid tumors in GSE29315. GO, Gene Ontology.

Integrated network analysis of miRNA-mRNA interaction

From the miRWalk, miRecords and TarMir databases, target genes of the DEMs were identified. A total of 24 overlapping genes were identified between the targets genes and DEGs (Table IV). Furthermore, 36 miRNA-gene pairs were obtained among the 24 overlapping genes and 9 DEMs (Table V). Node-degree analysis is summarized in Table VI. The regulation network between those overlapping genes and their upstream miRNAs is presented in Fig. 3.

Table IV.

Significant regulation of mRNAs in the specific miRNA-mRNA interacting regulatory network.

Gene	P-value	log2 (fold-change)
TNFRSF11B	0.036234	−0.836892
BZRAP1	0.006324	−0.784391
TGFBR2	0.016942	−0.736951
SALL1	0.039302	−0.687134
LRRN3	0.006382	−0.685662
TRAM2	0.016639	−0.664343
RGS16	0.027819	−0.648836
STARD13	0.013644	−0.642635
THBD	0.010984	−0.635682
GJA1	0.035288	−0.633043
PKIA	0.022944	−0.617443
CDC27	0.012718	−0.584214
NR2F2	0.013187	−0.575433
CHD9	0.048841	−0.573992
CXCL12	0.003888	−0.538192
NAB2	0.011359	0.527142
SLC25A5	0.029859	0.534298
PEG10	0.018256	0.549699
SLIT2	0.005701	0.551934
BASP1	0.026759	0.594218
TGIF1	0.030148	0.605397
ALDH1A3	0.010848	0.653685
ZNF148	0.031228	0.661078
CDH2	0.038832	0.862357

miRNA, microRNA.

Table V.

Significant regulation of miRNAs in the specific miRNA-mRNA interacting regulatory network.

miRNA	P-value	log2 (fold-change)
miR-7	0.0041392	−1.7320437
miR-137	0.0088043	0.9341108
miR-144	0.0059723	−0.8403724
miR-139-5p	0.0147834	−0.5824534
miR-145	0.0087553	−0.5871226
miR-296-5p	0.0124989	−0.5545413
miR-10a	0.0061396	−0.7961598
miR-452	0.0325921	−0.5399825
miR-493	0.0186862	0.6062014

miR/miRNA, microRNA.

Table VI.

Node-degree analysis of miRNA-mRNA interactions.

Node	Degree
miR-144	10
miR-137	9
miR-145	5
miR-7	5
ZNF148	4
miR-139-5p	3
PKIA	3
TGFBR2	3
TRAM2	3
TGIF1	2

miR/miRNA, microRNA.

Figure 3.

Regulatory network between miRNAs and their target mRNAs. miRNAs/miRs, microRNAs.

Discussion

The important roles of miRNAs in the pathogenesis of FTC have been identified previously (28). miRNAs exhibit different expression patterns within different tumor types, and are closely associated with the diagnosis, treatment and prognosis of tumors (29–31). Ak et al (21) observed that DEMs and differentially expressed mRNAs vary between benign and malignant tumors, which may suggest the different roles of these miRNAs and mRNAs. miR-197 and miR-346 have been indicated to be overexpressed in FTC, resulting in the dysregulation of their target genes (32). However, studies regarding DEMs and DEGs in FTC are rare. In the present study, the difference between miRNA-mRNA regulatory networks from FTC and FA samples were compared in order to investigate the mechanism of FTC. It was identified that miR-7, miR-1179, miR-7-2, miR-486-5p and miR-130b were the top downregulated miRNAs, and that miR-663b, miR-137, miR-30c-1, miR-767-5p and miR-603 were the top upregulated miRNAs. For the DEGs, the top downregulated genes were FABP4, CMAHP, ITM2A, CA4 and FAM189A2, and the top upregulated genes were EPB41L3, SCG5, PAX1, MTHFD2 and CDH2. In addition, miR-7, miR-296-5p, miR-10a, miR-144, miR-139-5p, miR-452 and miR-145 were downregulated, and miR-137 and miR-493 were upregulated in the FTC miRNA-mRNA regulatory network compared with those in FA. The gene arrays identified DEGs, in which leucine rich repeat neuronal 3, chromodomain helicase DNA binding protein 9, PKIA, zinc finger protein 148 (ZNF148), TGFB induced factor homeobox 1, transforming growth factor β receptor 2, gap junction protein α1 and CDH2 were observed to be target genes inversely correlated with miR-7, miR-144, miR-139-5p, miR-145 and miR-137. In other studies, FTC or FA have been compared with normal tissue, and differences in miRNA expression were observed to occur in the range between 1.2- and 2-fold, which was similar to the data of the present study (33–35).

In the present study, it was identified that miR-7, miR-296-5p, miR-10a, miR-144, miR-139-5p, miR-452, miR-145, miR-137 and miR-493 are important miRNAs that are differentially expressed between carcinoma and adenoma samples. Certain studies have suggested that miR-7 is not only a tumor promoter, but also a tumor suppressor. As a tumor suppressor, miR-7 is downregulated in tumors, such as thyroid cancer, breast cancer and castration-resistant prostate cancer, leading to a derepression of the oncogenes epidermal growth factor receptor, insulin receptor substrate 1, Raf-1 proto-oncogene, serine/threonine kinase, tyrosine kinase non-receptor 2, phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit δ, mechanistic target of rapamycin kinase, Ribosomal protein S6 kinase β-1 and phosphatidylinositol-4,5-bisphosphate 3-kinase (36–38). miR-296-5p has been revealed to be significantly inversely correlated with post-contrast T1 values for diffuse myocardial fibrosis in patients with hypertrophic cardiomyopathy, and is a downstream effector under conditions that promote glioblastoma stem cell stemness, and inhibit glioblastoma cell stemness and their capacity to self-renew as spheres and propagate glioma xenografts in vivo (39,40). miR-10a has been identified as a downregulated miRNA associated with human metastatic medullary thyroid carcinoma, and it may be important for tumor development and/or reflect C-cell lineage (41,42). miR-144 may suppress the invasion and migration capability of thyroid cancer and suppress the expression of zinc finger E-box-binding homeobox (ZEB)1 and ZEB2, the two E-cadherin suppressors, by directly binding their 3′UTRs (43). miR-137 was indicated to participate in hematopoiesis, particularly in the efficacy of warfarin, wherein miR-137 may cause aberrant vitamin K epoxide reductase complex subunit 1 expression (44). miR-139-5p is an oncogenic molecule in the process of tumorigenesis, and has been demonstrated to be a sensitive and specific biomarker for the diagnosis of thyroid tumors and others tumor types (45). Furthermore, it may be of use as a tractable therapeutic target to decrease the mortality rate and increase the survival rate (46). miR-145 has primarily been indicated as being downregulated in colorectal tumors. Previously, certain studies have identified that miR-145 is highly expressed in mesenchymal cells such as fibroblasts and smooth muscle cells (47). The miRNA was demonstrated to directly regulate the expression of thyroid hormone receptor TRβ1 in renal cancer cells and to correlate with intracellular triiodothyronine concentrations in renal tumors (48). miR-493 also promoted the invasion and chemoresistance of gastric cancer cells. However, dickkopf related-protein 1 overexpression reversed its effects on proliferation, invasion and chemo-sensitivity (49). Based on these data, we hypothesize that these miRNAs serve important roles in FTC with different pathways.

In the present study, several genes that were overlapping were identified between the DEGs and the target genes of the DEMs. These may be upregulated or downregulated. However, they all contributed to the development of FTC. Certain functions of these genes in cancer have been studied. For example, ZNF148 is a member of the human zinc finger Krüppel family and it maps to regions implicated in recurrent chromosomal rearrangements in hematological malignancies (50). The present study identified spalt-like transcription factor 1 (SALL1), which is one of the four human family members of the Spalt family. Members of the Spalt family are highly conserved zinc-finger transcription factors that are conserved from Caenorhabditis elegans to vertebrates, with regulatory functions in organogenesis, limb formation and cell fate assignment during neural development. SALL1 expression has been identified to correlate with the expression of CDH1, which is consistent with its tumor suppressive function and suggests its potential involvement in epithelial-to-mesenchymal transition (51,52). Cell division cycle 27 (CDC27) is a core component of the anaphase-promoting complex and is involved in the regulation of mitotic checkpoints to ensure chromosomal integrity (53). CDC27 may significantly affect the function of the polymeric protein complex and is also a target of certain anticancer drugs (54,55). Nuclear receptor subfamily 2 group F member 2 (NR2F2), also known as chicken ovalbumin upstream promoter transcription factor, is highly prioritized as a candidate gene associated with hypertension (56). Certain studies have demonstrated that NR2F2 is nuclear receptor transcription factor vital for angiogenesis and heart development (57). These data suggest that several genes have functions in numerous pathways involved in tumorigenesis and progression.

Each miRNA is able to regulate several hundred mRNAs. In addition, each mRNA may be targeted by several miRNAs and each mRNA participates in several biological functions in the human body. Therefore, each miRNA may affect different biological processes and pathways through a miRNA-mRNA network (58,59). It is important to understand the pathogenesis and treatment of tumors by investigating the specific miRNA-mRNA co-regulation effects. In the present study, mRNAs and their functions were described with GO enrichment analysis. There were 86 mRNAs and 8 biological functions involved. In total, ~80% of follicular carcinomas contain Ras gene mutations or a paired box gene 8/peroxisome proliferator-activated receptor γ gene rearrangement, which leads to uncontrolled proliferation. Mutations in the phosphatase and tensin homologue suppressor gene and the phosphatidylinositol 3-kinase pathway may be an important factor in the development of more aggressive thyroid cancer types and may be more common in follicular cancer, which is responsible for cell motility, locomotion and response to external stimulus (60–62). Other factors that have been implicated in the pathogenesis of FTC include gene mutations in p53, c-myc, c-fos and the thyrotropin receptor (63–66). These molecules serve functions in cell proliferation, apoptosis, cytoskeleton rearrangement and responses to drugs. Additionally, FTC, but not adenoma, recruits tumor-associated macrophages by releasing Chemokine (C-C motif) ligand 5; therefore, an abnormal immune response, including T cell activation, may be involved in follicular cancer. Other GO terms may be validated in future studies (3,67).

In conclusion, the present study identified 86 DEGs and 32 DEMs between FTC and FA. A total of 24 overlapping genes were identified between the DEGs and the target genes of the DEMs. Network analysis indicated a co-regulatory association between miR-296-5p, miR-10a, miR-139-5p, miR-452, miR-493, miR-7, miR-137, miR-144, miR-145 and corresponding targeted mRNAs in FTC. However, the present study has limitations, such as the small sample size, although attention was paid to ensure the use of two genetically homogenous populations to avoid population stratification. The mechanism of the miRNA-mRNA network and the roles of these genes in FTC require additional study and validation in vitro and in vivo.