Expression profiles of pivotal microRNAs and targets in thyroid papillary carcinoma: an analysis of The Cancer Genome Atlas.

In the present study, we analyzed microRNA (miRNA) and gene expression profiles using 499 papillary thyroid carcinoma (PTC) samples and 58 normal thyroid tissues obtained from The Cancer Genome Atlas database. A pivotal regulatory network of 18 miRNA and 16 targets was identified. Upregulated miRNAs (miR-222, miR-221, miR-146b, miR-181a/b/d, miR-34a, and miR-424) and downregulated miRNAs (miR-9-1, miR-138, miR-363, miR-20b, miR-195, and miR-152) were identified. Among them, the upregulation of miR-424 and downregulation of miR-363, miR-195, and miR-152 were not previously identified. The genes CCNE2 (also known as cyclin E2), E2F1, RARA, CCND1 (cyclin D1), RUNX1, ITGA2, MET, CDKN1A (p21), and COL4A1 were overexpressed, and AXIN2, TRAF6, BCL2, RARB, HSP90B1, FGF7, and PDGFRA were downregulated. Among them, CCNE2, COL4A1, TRAF6, and HSP90B1 were newly identified. Based on receiver operating characteristic curves, several miRNAs (miR-222, miR-221, and miR-34a) and genes (CCND1 and MET) were ideal diagnostic indicators, with sensitivities and specificities greater than 90%. The combination of inversely expressed miRNAs and targets improved diagnostic accuracy. In a clinical feature analysis, several miRNAs (miR-34a, miR-424, miR-20b, and miR-152) and genes (CCNE2, COL4A1, TRAF6, and HSP90B1) were associated with aggressive clinical features, which have not previously been reported. Our study not only identified a pivotal miRNA regulatory network associated with PTC but also provided evidence that miRNAs and target genes can be used as biomarkers in PTC diagnosis and clinical risk evaluation.

Introduction

Thyroid carcinoma is the most common endocrine-related malignancy worldwide. According to cell origin and pathological pattern, thyroid carcinomas are normally divided into four types: papillary, follicular, medullary, and anaplastic. The most common type of thyroid carcinoma is papillary thyroid carcinoma (PTC), accounting for 80% of all cases. Thyroid carcinoma is considered a curable disease, but an accurate diagnosis and appropriate risk evaluation are still necessary. Currently, the gold standard for thyroid nodule diagnosis is ultrasound-guided fine-needle aspiration biopsy, achieving correct final diagnoses in 70%–80% of cases; the remaining 20%–30% of cases are considered indeterminate for malignancy.1 Operative strategy and prognosis are related to risk evaluation. High-risk tumors are often characterized by progressive features such as extrathyroid extension, lymph node metastasis, and distant metastasis.2 In these cases, an extended surgery strategy should be performed.

MicroRNAs (miRNAs) are nonprotein-coding, endogenous small RNA molecules. They act as oncogenes as well as tumor suppressors by regulating the expression of target genes at the posttranscriptional level. There is evidence that miRNAs are involved in thyroid tumorigenesis and progression.3–8 In the present study, we analyzed the expression profiles of miRNAs and genes retrieved from The Cancer Genome Atlas (TCGA; http://cancergenome.nih.gov/). A pivotal regulatory network of miRNAs associated with PTC was identified. The diagnostic values of miRNAs and genes in the regulatory network were examined, and the results indicated that diagnostic accuracy could be improved by the combination of miRNAs and targets. Furthermore, an association between miRNA and target expression levels and aggressive clinical features was detected.

Materials and methods

TCGA data access

The miRNA expression microarray data (Level 3), RNA expression microarray data (Level 3), and corresponding clinical data for thyroid cancer patients were obtained from TCGA. In total, 499 PTC samples and 58 control samples (paracancerous tissues) with complete miRNA, RNA, and clinical data were selected for inclusion in this study. The expression levels of 1,046 human miRNAs and 20,531 human genes in each sample were assessed using the Illumina’s HiSeq system.

TCGA background

TCGA was launched by the National Cancer Institute and the National Human Genome Research Institute in 2006. The program grew to include samples from 11,000 patients across 33 tumor types and currently represents the largest tumor collection, including key genomic and molecular characteristics. Findings for 13 cancer types have already been published, and as of late 2014, TCGA scientists have nearly completed the sequencing of protein-coding regions (exomes) for most tumor types and completed whole-genome sequencing of 1,000 tumor samples. The latter data characterize the complete genomic DNA sequence. Results from TCGA analyses have led to more than 2,700 articles in research journals to date.

Data acquisition

From the TCGA homepage, http://cancergenome.nih.gov, the THCA (thyroid carcinoma) page was selected. All miRNASeq, RNASeqV2, and clinical data were chosen, and the Build Archive function was used. The miRNA and RNA expression data are classified as Level 3. In the miRNAseq section, files with the suffixes .miRNA, .quantification, and .txt were selected, and all files in the metadata and clinical sections were examined.

Inclusion and exclusion criteria

The TCGA database consisted of 507 cases of thyroid carcinoma and 59 normal tissue samples. Excluding non-PTC samples and those lacking complete clinical information, 499 PTC samples and 58 normal tissue samples were included in the subsequent analyses.

Microarray analysis

The raw data obtained from miRNA microarrays were normalized to the 75th percentile signal intensity as recommended by the vendor. After normalization, all negative signal values were replaced with 0.01 and the values for multiple replicate spots for each miRNA were summarized as median signals, which were used for statistical analyses.

To identify significantly upregulated and downregulated miRNAs, the pairwise Welch’s t-test was applied to detect differences between c-Raf-transgenic and non-transgenic animal groups. The significance thresholds for the Welch’s t-test were set to P<0.05, fold-change [FC] >2, and a mean signal of all samples >50th percentile. Statistical and hierarchical clustering analyses were implemented in SAM (Significance Analysis of Microarrays) 3.01. In the tables and hierarchical clustering analyses, the FCs gene expression levels are given as signal logarithm ratios with base 2.

Prediction of miRNA target genes

Target mRNAs were predicted using the web tools TargetScan (http://www.targetscan.org/) and PITA (http://genie.weizmann.ac.il/pubs/mir07/index.html). Both miRNA target prediction programs rely on sequence complementarities of the miRNA seed region (nucleotides 2–7) to the 3′-UTR sequences in candidate target genes. The results of the two analyses are summarized in a Venn diagram.

Statistical analysis

Receiver operating characteristic curves for diagnostic value calculations were estimated using MedCalc 12.1.3.0 software (MedCalc Software, Mariakerke, Belgium). The area under the curve (AUC) represents the diagnostic accuracy for a specific criterion. The criterion is the threshold value for a laboratory test. Criteria listed in the tables are automated threshold values corresponding to the highest AUC values. Expression ratios of miRNAs and targets were calculated following the methods of Gordon et al.9 Other statistical analyses were carried out using SPSS 17.0 (SPSS Inc., Chicago, IL, USA). Mann–Whitney U-tests were used to identify possible associations between miRNA and target expression levels and clinical features. A P-value less than 0.05 was considered statistically significant.

Results

Expression profiles of miRNAs and genes

In total, 48 miRNAs and 2,731 genes were differentially expressed between the 499 PTC samples and 58 normal tissue samples. Of the 48 miRNAs, 18 were upregulated and 30 were downregulated (false discovery rate [FDR] <0.01 and FC>2). Of the 2,731 genes, 1,536 were upregulated and 1,195 were downregulated (FDR <0.01 and FC>2).

Regulatory network of miRNAs in PTC

Using starBase v2.0 (Jianhua Yang, Sun Yat-sen University, 2010–2013) to predict target genes, 30 miRNAs and 277 target genes were successfully predicted (data not shown). After Kyoto Encyclopedia of Genes and Genomes pathway mapping, 16 of 277 genes involved in cancer-related pathways were identified. Corresponding miRNAs were 17 of 30. The expression levels of 17 miRNAs and 16 target genes were visualized and are represented in Figure 1. These miRNAs and target genes were pivotal in PTC because all of them were differentially expressed by at least twofold and played roles in cancer signaling pathways. The regulatory network is represented in Figure 2.

Figure 1

Differential expression of 17 miRNAs (A) and 16 target genes (B) in 499 PTC samples and 58 normal tissue samples retrieved from TCGA.

Note: Red to green indicates decreasing expression levels of individual miRNAs or genes.

Abbreviations: miRNAs, microRNAs; PTC, papillary thyroid carcinoma; TCGA, The Cancer Genome Atlas.

Figure 2

Regulatory network of pivotal miRNAs in PTC.

Notes: Red-marked miRNAs and genes were upregulated and blue-marked miRNAs and genes were downregulated. Key regulatory network nodes were identified, including many genes, such as AXIN2, BCL2, CCND1, RUNX1, CCNE2, and miRNAs such as miR-424, miR-181a/b/d, miR-20b, miR-195, miR-152, and miR-144. Regulatory networks formed by these key nodes are shown, including the internal relationships among thyroid cancer miRNAs and genes.

Abbreviations: miRNAs, microRNAs; PTC, papillary thyroid carcinoma.

Diagnostic value of pivotal miRNAs and genes

To determine whether pivotal miRNAs and genes could be used to distinguish PTC from normal samples, receiver operating characteristic curves were applied to calculate diagnostic utility. AUC values, sensitivity, specificity, and the criterion for each miRNA and gene are listed in Tables 1 and 2. The miRNAs miR-222, miR-221, miR-146b, miR-34a, and miR-9-2 showed ideal diagnostic values with respect to AUC, sensitivity, and specificity, each of which was greater than 90%. The CCND1 and MET genes showed ideal diagnostic values with respect to AUC, sensitivity, and specificity, each of which was greater than 90%.

Table 1

Diagnostic value of differentially expressed miRNAs

miRNAs	AUC	SEN	SPE	Criterion
miR-221	0.961	91.4	96.6	>503.6685
miR-222	0.958	91.4	94.8	>125.782
miR-146b	0.95	84.5	96.6	>2,194.6538
miR-34a	0.944	91.4	94.8	>178.0286
miR-9-1	0.927	86.2	94.8	<53.4904
miR-144	0.924	94.8	81	<502.158
miR-9-2	0.924	91.4	91.4	<74.6703
miR-181b-2	0.919	86.2	91.4	>13.2275
miR-363	0.916	82.8	87.9	<7.4676
miR-195	0.913	77.6	96.6	<38.8908
miR-181a-2	0.897	87.9	89.7	>2,924.7844
miR-152	0.892	86.2	91.4	<500.1376
miR-138-1	0.883	81.0	94.8	<45.2887
miR-138-2	0.88	75.9	96.6	<28.8837
miR-20b	0.875	79.3	87.9	<11.7008
miR-181b-1	0.872	84.5	81.0	>710.2916

Abbreviations: miRNAs, microRNAs; AUC, area under the curve; SEN, sensitivity (%); SPE, specificity (%).

Table 2

Diagnostic value of differentially expressed genes

Genes	AUC	SEN	SPE	Criterion
CCND1	0.963	91.4	93.1	>5,178.2373
MET	0.942	93.1	91.4	>2,623.9624
BCL2	0.942	84.5	93.1	<2,647.5432
RARB	0.94	82.8	94.8	<226.8591
E2F1	0.936	98.3	79.3	>62.7434
RARA	0.935	91.4	87.9	>747.9881
TRAF6	0.93	84.5	91.4	<182.3617
RUNX1	0.903	69	100	>864.4198
FGF7	0.894	81	93.1	<115.4518
ITGA2	0.89	82.8	96.6	>405.4988
PDGFRA	0.847	74.1	93.1	<235.9651
HSP90B1	0.841	58.6	100	<30,834.7969
AXIN2	0.836	70.7	91.4	<400.454
CCNE2	0.805	81	77.6	>24.1202
COL4A1	0.801	74.1	82.8	>3,715.8176
CDKN1A	0.793	86.2	63.8	>2,494.7954

Abbreviations: AUC, area under the curve; SEN, sensitivity (%); SPE, specificity (%).

To further improve the diagnostic power of miRNAs and genes, their combinations were examined. Inversely expressed miRNAs and target genes were considered novel markers. The expression ratios were estimated following the methods of Gordon et al.9 Consistent with expectations, several miRNA and target combinations improved PTC diagnosis accuracy, as evidenced by the increased AUC value, sensitivity, and specificity. These combination markers are listed in Table 3. Among them, the miR-34a/BCL2 combination substantially improved diagnostic accuracy (AUC of 0.989, sensitivity of 98.3%, and specificity of 98.4%) compared with miR-34a alone (AUC of 0.944, sensitivity of 0.983, and specificity of 0.984) and BCL2 alone (AUC of 0.942, sensitivity of 84.5%, and specificity of 93.1%).

Table 3

Diagnostic value of combinations of miRNAs and targets

Combinations	AUC	SEN	SPE	Criterion
miR-34a/BCL2	0.989	98.3	98.4	>0.0731
miR-181b-2/BCL2	0.97	94.8	94.8	>0.0036
miR-221/AXIN2	0.966	94.8	91.4	>0.6077
miR-222/AXIN2	0.963	93.1	94.1	>0.2262
miR-146b/TRAF6	0.962	91.4	93.1	>6.8312
miR-181a-2/BCL2	0.962	91.4	96.6	>1.071
miR-181b-2/HSP90B1	0.942	93	95	>0.0003

Abbreviations: miRNAs, microRNAs; AUC, area under the curve; SEN, sensitivity (%) SPE, specificity (%).

Association with high-risk clinical features

In view of the essential roles of miRNAs and genes in the cancer signaling pathway and their use in effective tumor identification, we inferred that they may be related to progressive clinical signatures. We analyzed 499 PTCs and 58 normal tissue samples. Associations between lymph node metastasis (Table 4), miRNA and gene expression levels and extrathyroidal extension (Table 5), and high TNM classification stage (Table 6) were detected. Distant metastasis is commonly considered a high-risk clinical feature. Owing to the limited number of appropriate cases, distant metastasis was not explored in the present study.

Table 4

Associations between miRNAs and genes with lymph node metastasis

miRNA/gene	Lymph node metastasis (+)	Lymph node metastasis (−)	P-value
miR-146b	3.09 (2.78–3.16)*	2.97 (2.47–3.14)	0.024
miR-222	2.75 (2.62–2.90)	2.65 (2.41–2.84)	0.000
miR-221	3.31 (3.14–3.41)	3.23 (2.97–3.42)	0.004
miR-181a-1	3.78 (3.68–3.91)	3.72 (3.59–3.88)	0.005
miR-152	2.43 (2.30–2.53)	2.48 (2.31–2.68)	0.005
AXIN2	2.39 (2.19–2.60)	2.47 (2.26–2.67)	0.044
TRAF6	2.14 (2.05–2.22)	2.17 (2.08–2.25)	0.048
BCL2	2.92 (2.73–3.06)	3.03 (2.87–3.06)	0.000
HSP90B1	4.42 (4.31–4.57)	4.47 (4.36–4.67)	0.029
PDGFRA	2.24 (1.93–2.58)	1.98 (1.54–2.39)	0.000
RARA	3.09 (3.03–3.17)	3.04 (2.93–3.12)	0.000
RUNX1	3.15 (2.99–3.31)	3.03 (2.68–3.21)	0.000
MET	4.11 (3.93–4.22)	4.00 (3.57–4.17)	0.000

Note:

Original expressing data were log transformed.

Abbreviation: miRNAs, microRNAs.

Table 5

Associations between miRNAs and genes with extra-thyroid extension

miRNA/gene	Extrathyroidal extension (+)	Extrathyroidal extension (−)	P-value
miR-146b	4.49 (4.21–4.65)	4.30 (3.45–4.62)	0.000
miR-222	2.72 (2.60–2.86)	2.63 (2.39–2.82)	0.000
miR-221	3.31 (3.12–3.46)	3.21 (2.95–3.39)	0.000
miR-34a	2.68 (2.50–2.84)	2.61 (2.45–2.79)	0.013
miR-424	2.53 (2.29–2.76)	2.64 (2.38–2.84)	0.009
miR-138-1	1.26 (1.10–1.54)	1.38 (1.19–1.60)	0.005
miR-20b	0.72 (0.56–1.06)	0.98 (0.66–1.19)	0.000
miR-152	2.45 (2.32–2.58)	2.49 (2.32–2.72)	0.014
AXIN2	2.36 (2.18–2.55)	2.51 (2.31–2.69)	0.000
TRAF6	2.14 (2.05–2.20)	2.18 (2.11–2.27)	0.000
BCL2	2.89 (2.74–3.06)	3.06 (2.90–3.33)	0.000
RARB	2.06 (1.91–2.21)	2.19 (1.97–2.43)	0.000
HSP90B1	4.39 (4.30–4.52)	4.52 (4.38–4.65)	0.000
PDGFRA	2.24 (1.90–2.58)	1.96 (1.47–2.39)	0.000
RARA	3.09 (3.12–3.06)	3.04 (2.93–3.03)	0.000
RUNX1	3.19 (2.95–3.33)	2.99 (2.58–3.17)	0.000
ITGA2	3.28 (2.85–3.5)	2.92 (2.63–3.27)	0.000
MET	4.13 (3.95–4.27)	3.96 (3.55–4.15)	0.000
CDKN1A	3.80 (3.63–3.96)	3.74 (3.54–3.90)	0.019

Abbreviation: miRNAs, microRNAs.

Table 6

Association between miRNAs and genes with later tumor stage

miRNA/gene	III/IV stage	I/II stage	P-value
miR-146b	4.49 (4.08–4.68)	4.31 (3.58–4.61)	0.001
miR-222	2.70 (2.52–2.85)	2.64 (2.38–2.83)	0.012
miR-34a	2.69 (2.50–2.84)	2.61 (2.45–2.76)	0.004
miR-424	2.53 (2.28–2.77)	2.63 (2.38–2.84)	0.014
miR-138-1	1.27 (1.13–1.57)	1.38 (1.18–1.64)	0.045
miR-20b	0.82 (0.56–1.9)	0.95 (0.67–1.16)	0.01
AXIN2	2.38 (2.20–2.61)	2.49 (2.30–2.67)	0.005
TRAF6	2.14 (2.05–2.20)	2.18 (2.11–2.26)	0.000
BCL2	2.94 (2.75–3.11)	3.06 (2.89–3.30)	0.000
RARB	2.09 (1.92–2.32)	2.18 (1.97–2.39)	0.004
HSP90B1	4.40 (4.30–4.55)	4.52 (4.38–4.66)	0.000
PDGFRA	2.22 (1.77–2.56)	1.96 (1.55–2.39)	0.002
RUNX1	3.14 (2.88–2.28)	2.99 (2.67–3.20)	0.000
MET	4.06 (3.86–4.21)	3.97 (3.54–4.17)	0.001

Abbreviation: miRNAs, microRNAs.

The expression levels of miR-146b, miR-222, miR-221, miR-34a, miR-181a, miR-424, miR-138-1, miR-20b, and miR-152 were associated with PTC invasion or progression (P<0.05). Among them, miR-146b and miR-222 were associated with all high-risk clinical features. miR-181b-1, miR-181a-2, miR-181b-2, miR-181d, miR-9-1, miR-9-2, miR-138-2, and miR-195 did not show a relationship with PTC invasion or progression.

The expression levels of AXIN2, TRAF6, BCL2, RARB, HSP90B1, PDGFRA, RARA, RUNX1, ITGA2, MET, and CDKN1A were associated with high-risk clinical features. Among them, AXIN2, TRAF6, BCL2, HSP90B1, PDGFRA, RUNX1, and MET were correlated with all high-risk clinical features. FGF7, CCNE2, E2F1, CCND1, and COL4A1 showed no relationship with PTC invasion or progression.

Discussion

The high incidence of thyroid disease in the population is problematic. As many as 50% of individuals have microscopic nodules, 3.5% have occult papillary carcinoma, and 15% have palpable goiters.10 Currently, ultrasound-guided fine-needle aspiration biopsy is the most reliable method for detecting thyroid nodules. However, definitive diagnoses still cannot be made for 20%–30% of cases.1 Therefore, additional methods that improve the sensitivity and specificity of diagnosis are highly desirable. Molecular markers such as BRAF, RAS, RET/PTC, PAX8/PPARγ, and galectin-3 may be considered for indeterminate cytology according to the American Thyroid Association guidelines.11 Overcoming the challenges of accurate assessments of the risk for individual patients is important to establish appropriate treatment plans and optimize outcomes. An increasing number of mutations in thyroid tumors from low grade to high grade have been reported, indicating frequent recurrence and death.12 In this context, molecular markers would facilitate tumor stage identification and risk stratification, assisting clinicians in determining appropriate treatment strategies, and in clinical monitoring.

miRNAs play important roles in multiple biological and metabolic processes, such as cell differentiation, proliferation, survival, and malignancy.13,14 Although numerous miRNAs involved in PTC have been identified, the global regulation of miRNAs in PTC remains unclear. In our study, we identified a pivotal regulatory network of miRNAs associated with PTC. All differentially expressed miRNAs and genes were screened using as many as 547 samples obtained from TCGA. Targeting relationships were reliable based on successful predictions by multiple common programs. All target genes played roles in cancer-related pathways. In all, 18 miRNAs and 20 target genes were core and vital in the regulation of PTC. Among the regulatory relationships, miR-146b targeting TRAF6, miR-222 targeting AXIN2, miR-181a/b targeting BCL2, miR-424 targeting BCL2, miR-34a targeting PDGFRA and BCL2, miR-138 targeting RARA, miR-20b targeting CCND1 (also known as cyclin D1) and RUNX1, miR-195 targeting CCND1, and miR-152 targeting MET have been previously reported.15–22 miR-181a/b targeting HSP90B1, miR-424 targeting AXIN2/RARB/FGF7, miR-181d targeting BCL2 and HSP90B1, miR-144 targeting RUNX1 and CCNE2 (also known as cyclin E2), miR-363 targeting CCNE2, miR-195 targeting ITGA2, miR-20 targeting LAMA3/CDKN1A (also known as p21)/E2F1, and miR-152 targeting COL4A1 are newly identified.

Next, we sought to evaluate the diagnostic ability of the miRNAs and targets. The miRNAs miR-221, miR-222, miR-34a, and miR-9-2 displayed ideal diagnostic values based on the AUCs, sensitivities, and specificities greater than 90%. The CCND1 and MET genes displayed ideal diagnostic values. All had potential use in PTC diagnosis for indefinite cases. Some studies have indicated that miRNA–target combinations can improve diagnostic accuracy;9,23 accordingly, we determined combinations based on miRNAs and targets with opposing expression patterns. Several combinations of miRNAs and targets showed higher diagnostic values compared with single factors. Particularly, the miR-34a/BCL2 combination had an AUC value of 0.989, sensitivity of 98.3%, and specificity of 98.4% for certain criterion (expression ratio >0.0731), which indicates that they are suitable tumor markers. Our results provide a novel combination method to improve diagnostic accuracy and can be validated in a wide range of cancer types.

Thyroid tumor markers are often related to high-risk clinical features and indicate poor prognosis.12,24 In the present study, we examined the associations between all pivotal miRNAs and genes and aggressive factors, including lymph node metastasis, extrathyroidal extension, and late-stage tumors. Nine of 17 miRNAs and eleven of 17 target genes were associated with at least one high-risk clinical feature. Clinical associations involving miR-221, miR-222, miR-146b, and miR-138-1 were supported by previous studies.25–27 miR-34a, miR-424, miR-20b, and miR-152, which played roles in PTC progression and invasion, have not been previously reported. Differential expression of CCND1, BCL2, MET, E2F1, RUNX1, ITGA2, RARA, PDGFRA, and CDKNA1 has been demonstrated for aggressive PTC in other studies.28–36 Our study is the first to identify the association between AXIN2, HSP90B1, and CDKNA1 and aggressive clinical features. The aberrant expression levels of miRNAs and genes identified in the present study can be used as diagnostic markers as well as treatment targets.

We identified novel miRNAs for PTC diagnosis and associations with aggression. Biofunction experiments and clinical validation of these miRNAs are our next goals to confirm their significance in PTC.

Conclusion

In summary, we identified a set of pivotal thyroid cancer-related miRNAs and genes. Several novel targeting relationships were discovered. Many of these miRNAs and genes displayed potential value for the early diagnosis of PTC. Combinations of inversely expressed miRNAs and targets further improved the observed diagnostic accuracy. Based on their roles in the cancer pathway and clinical examinations, most of the miRNAs and genes in the present study were related to high-risk clinical features, indicating their potential use for risk stratification and prognosis for PTC patients. Our study indicates potential clinical applications of the miRNAs and genes for diagnosis, prognosis, and targeted treatment in thyroid malignant disease.