A miRNA-clinicopathological nomogram for the prediction of central lymph node metastasis in papillary thyroid carcinoma-analysis from TCGA database.

It is of significance to evaluate central lymph node status in patients with papillary thyroid carcinoma (PTC), because it can decrease postoperative complications resulting from unnecessary prophylactic central lymph node dissection (CLND). Due to the low sensitivity and specificity of neck ultrasonography in the evaluation of central lymph node metastasis (CLNM), it is urgently required to find alternative biomarkers to predict CLNM in PTC patients, which is the main purpose of this study.RNA-sequencing datasets and clinical data of 506 patients with thyroid carcinoma from the Cancer Genome Atlas (TCGA) database were downloaded and analyzed to identify differentially expressed miRNAs (DEMs), which can independently predict CLNM in PTC. A nomogram predictive of CLNM was developed based on clinical characteristics and the identified miRNAs. Receiver operating characteristics curves were drawn to evaluate the predictive performance of the nomogram. Bioinformatics analyses, including target genes identification, functional enrichment analysis, and protein-protein interaction network, were performed to explore the potential roles of the identified DEMs related to CLNM in PTC.A total of 316 PTC patients were included to identify DEMs. Two hundred thirty-seven (75%) PTC patients were randomly selected from the 316 patients as a training set, while the remaining 79 (25%) patients were regarded as a testing set for validation. Two DEMs, miRNA-146b-3p (HR: 1.327, 95% CI = 1.135-1.551, P = .000) and miRNA-363-3p (HR: 0.714, 95% CI = 0.528-0.966, P = .029), were significantly associated with CLNM. A risk score based on these 2 DEMs and calculating from multivariate logistic regression analysis, was significantly lower in N0 group over N1a group in both training (N0 vs N1a: 2.04 ± 1.01 vs 2.73 ± 0.61, P = .000) and testing (N0 vs N1a: 2.20 ± 0.93 vs 2.79 ± 0.68, P = .003) sets. The nomogram including risk score, age, and extrathyroidal extension (ETE) was constructed in the training set and was then validated in the testing set, which showed better prediction value than the other three predictors (risk score, age, and ETE) in terms of CLNM identification. Bioinformatics analyses revealed that 5 hub genes, SLC6A1, SYT1, COL19A1, RIMS2, and COL1A2, might involve in pathways including extracellular matrix organization, ion transmembrane transporter activity, axon guidance, and ABC transporters.On the basis of this study, the nomogram including risk score, age, and ETE showed good prediction of CLNM in PTC, which has a potential to facilitate individualized decision for surgical plans.

Introduction

In recent years, the incidence of papillary thyroid carcinoma (PTC) has risen remarkably all over the world mainly due to the widespread utility of the high-resolution ultrasonography (US), and US-guided fine-needle aspiration cytology. More than 56,000 new patients were diagnosed with PTC in the United States in 2015.[ In Korea, PTC is the most common endocrine malignancy in females, the incidence of which is the highest globally[; while in China, PTC is the 3rd most common malignancy affecting individuals.[

Generally, PTC grows slowly, differentiates well, and has an indolent clinical course, with a disease-specific mortality rate of less than 1%.[ However, lymph node metastasis (LNM) is not uncommon in PTC, with a prevalence ranging from 37% to 64%.[ LNM often occurs at the early stage and initially appears in the central compartment in PTC, which is associated with locoregional recurrence and is regarded as an unfavorable prognostic factor in PTC patients.[ However, according to the 2015 guidelines of the American Thyroid Association (ATA), lymph node dissection (LND) in the central compartment is not recommended in PTC patients in the absence of clinically detectable cervical nodal metastasis.[ On the other hand, prophylactic LND is an effective procedure to decrease locoregional recurrence, but it is not accepted by all surgeons due to the higher incidence of postoperative complications such as hypoparathyroidism and recurrent laryngeal nerve injury.[

Due to the limited role of neck ultrasonography in the evaluation of central lymph node metastasis (CLNM), with a low sensitivity of 23.0% to 53.2%,[ it is necessary to find new risk factors of CLNM to identify patients at high risk and to individualize surgical plans. MicroRNA (miRNA), a small nonencoding RNA negatively regulating gene expression, plays a vital role in tumor pathogenesis and progression in various carcinomas including PTC.[ Previous studies have revealed the relationships between miRNAs and LNM in PTC, but the clinical values of the related miRNAs were not well addressed.[

The purpose of this study is to mine datasets of PTC from the Cancer Genome Atlas (TCGA) database and to identify differentially expressed miRNAs (DEMs) related to CLNM. Furthermore, we proposed a nomogram based on the identified DEMs and clinical data to predict CLNM in PTC. Bioinformatics analyses were also performed to explore the potential roles of the identified miRNAs in regulating the biological process (BP) of PTC.

Methods

Data acquisition and grouping

We downloaded research data including miRNA and mRNA expression profiles and clinical information of patients with thyroid carcinoma from TCGA database (https://cancergenome.nih.gov/) in April 2020. Five hundred sixty-seven records with RNA-sequencing expressions were available including 56753 mRNAs and 2190 miRNAs, while clinical data were obtained from 506 patients. The patient inclusion criteria in this study were as follows:

patients pathologically diagnosed with PTC;

pathological status of lymph node was N0 or N1a.

Patients diagnosed with the other types of thyroid carcinoma and patients with lymph node status of Nx, N1, or N1b were excluded. Finally, 316 PTC patents with clinicopathological data including age, sex, tumor size, extrathyroidal extension (ETE) status, primary neoplasm focus type (unifocal or multifocal), primary neoplasm location (bilaterality) were included to identify DEMs associated with PTC. In addition, a training set and a testing set, with a ratio of 3:1, were randomly generated from the 316 patients by R package for further analyses (Fig. 1).

Figure 1

Grouping criteria. CLNM = central lymph node metastasis, PTC = papillary thyroid carcinoma.

Data processing, potential risk miRNAs identification and risk score definition

The miRNA expression profiles were analyzed by R/Bioconductor package of edgeR. DEMs between normal tissues and tumor tissues were filtered by |log2FC| ≥ 1 and FDR < 0.05, and hierarchical clustering heatmaps and volcano plots of the RNA-sequencing count datasets from the TCGA database were drawn using R software. Low-abundance DEMs were excluded if the corresponding LogCPM was less than 1.0 to identify top 10 up-regulated and top 10 down-regulated significant DEMs. These miRNAs were put into univariate and then multivariate analyses to find potential risk miRNAs related to CLNM in the training set. The coefficients of the risk miRNAs from multivariate logistic regression analysis were used to calculate miRNA-based risk score of CLNM in both training and testing sets: Risk score = β1 × level1 +β2 × level2 +…+βn × leveln, where β is the individual coefficient and level indicates the expression level of the corresponding miRNA.

Establishment of a nomogram to predict CLNM

To test whether the miRNA-based risk score can predict CLNM independent of clinical characteristics in patients with PTC, the potential clinical risk factors of CLNM (including age, sex, tumor size, ETE, primary neoplasm focus type, primary neoplasm location) and the risk score were analyzed by univariate and multivariate analyses. A nomogram to predict CLNM was established based on the multivariate logistic analysis in the training set and then validated in the testing set.

Functional enrichment analysis

The target genes of the risk miRNAs were predicted using the miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), Targetscan (http://www.targetscan.org) and miRDB (http://www.mirdb.org/) databases. Target genes that are covered in at least 2 of these 3 databases were obtained. Then these target genes were intersected with the differentially expressed genes between N0 PTC and N1a PTC to get the CLNM-specific target genes. Venn diagrams were used to show CLNM-specific target genes and Cytoscape was utilized to visualize the relationships between the significant miRNAs and the CLNM-specific target genes. All of the CLNM-specific target genes were analyzed by Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway using R language “clusterProfiler” package and the “org.Hs.eg.db” package. The Gene Ontology (GO) enrichment analysis was conducted via http://www.webgestalt.org/option.php. Protein–protein interaction (PPI) network was built using https://string-db.org/, and the hub genes were filtered and visualized using Cytoscape.

Ethical approval and informed consent were not necessary, because the present study was based on, a public database, TCGA database.

Statistical analyses

For quantitative data, the results were expressed as mean ± standard deviation. For categorical data, the number and percentage of cases were used. Two quantitative data were transformed to binary variables (age ≥ 55 years, age < 55 years; tumor size >2 cm, tumor size ≤2 cm). The filtered risk miRNAs were firstly evaluated by univariate analysis using Student t test, then the risk miRNAs with P value <.05 were put into multivariate logistic regression analysis with the forward stepwise method and likelihood ratio to confirm the independent risk miRNAs of CLNM. Risk scores were calculated based on the logistic analysis result. The potential clinical risk factors of CLNM were analyzed using univariate and multivariate analyses, and then a nomogram including clinical factors and the risk score was constructed on the basis of the multivariate logistic regression analysis (In order to identify more potential predictors, the level for rejection of the null hypothesis was set at a P value <.10 for nomogram establishment.). Receiver operator characteristic (ROC) curves were drawn to compare the CLNM prediction power of different risk factors in both training and testing sets. SPSS 16.0 (SPSS, Chicago, IL) for windows was used to perform statistical analyses, and the P value <.05 was considered to be statistically significant if not specifically indicated.

Results

Clinical data of 237 PTC patients in the training set and 79 PTC patients in the testing set

Two hundred thirty-seven (75% of the 316 patients) PTC patients in the training set and 79 (25% of the 316 patients) PTC patients in the testing set were randomly generated with a ratio of 3:1, whose clinical characteristics were shown in Table 1. The clinical data were comparable between the 2 sets, suggesting the testing set can be used to evaluate the predictive performance of the predictive nomogram developed from the data in the training set.

Table 1

Characteristics of included patients in training and testing sets.

Selection of risk miRNAs

A total of 2190 miRNAs were obtained from TCGA database. Among the 2190 miRNAs, 127 miRNAs (|log2FC| ≥ 1 and FDR < 0.05) were significantly dysregulated in PTC tissues comparing to normal tissues. Heatmap (Fig. 2) and volcano plots (Fig. 3) showed different miRNA expression patterns between normal and PTC tissues. The correlations between the expression level of DEMs and central lymph node status were firstly determined in the training set. The Univariate and multivariate analyses results showed that 2 DEMs, miRNA-146b-3p (HR: 1.327, 95% CI = 1.135–1.551, P = .000) and miRNA-363–3p (HR: 0.714, 95% CI = 0.528–0.966, P = .029), can independently predict CLNM (Tables 2 and 3) in the training set. The miRNA-based risk score in N0 group was significantly lower than that in N1a group in the training (N0 vs N1a: 2.04 ± 1.01 vs 2.73 ± 0.61, P = .000) and testing (N0 vs N1a: 2.20 ± 0.93 vs 2.79 ± 0.68, P = .003) sets (Fig. 4A and B).

Figure 2

The heatmap of top 40 differentially expressed miRNAs.

Figure 3

The volcano plots of differentially expressed miRNAs.

Table 2

Top 20 DEMs associated with CLNM in the training set.

Table 3

Multivariate analyses to identify independent risk DEMs associated with CLNM in the training set.

Figure 4

A, Risk score in N0 and N1a groups in the training set (N0 vs N1a: 2.04 ± 1.01 vs 2.73 ± 0.61, P = .000); B, Risk score in N0 and N1a groups in the testing set (N0 vs N1a: 2.20 ± 0.93 vs 2.79 ± 0.68, P = .003).

Establishment of nomogram

As is shown in Table 4, in addition to risk score, 2 clinical variables, age and ETE, were significantly associated with CLNM not only in the training set, but also in the testing set. These 2 clinical factors, together with risk score, were put into multivariate logistic analysis (Table 5), which confirmed that these 3 factors remained significantly associated with CLNM in both training and testing sets. Based on the analyses above, a nomogram including age, ETE, and risk score was constructed (Fig. 5). The risk possibility of CLNM can be calculated for an individual patient. A line can be drawn for each risk factor (line 2–4) to obtain point values (line 1). Then the sum of all the 3 points can be plotted on the total point axis (line 5), and the corresponding risk axis can be drawn (line 6).

Table 4

Univariate analyses to identify potential predictors of CLNM.

Table 5

Multivariate analyses of factors to predicting CLNM in training and testing sets.

Figure 5

Nomogram for predicting central lymph node metastasis. The “0” represents age ≥55 years, while the “1” represents age <55 years in line 2. The “0” represents PTC with ETE, while the “1” represents PTC with ETE in line 3. ETE = extrathyroidal extension, PTC = papillary thyroid carcinoma.

Test of nomogram

To test the predictive power of the nomogram, ROC curves were drawn to compare the predictive performance of CLNM between the nomogram and the other 3 independent predictors in both training and testing sets. And the areas under the ROC curves suggested that the nomogram performed better than the other 3 risk factors in terms of CLNM prediction (Fig. 6).

Figure 6

ROC curves to compare the prediction performance among different predictors in both training (A) and testing (B) sets. AUC = area under curve, ETE = extrathyroidal extension, ROC = receiver operating characteristics.

Target genes

The overlapping genes predicted by at least 2 of the 3 databases (miRTarBase, TargetScan, and miRDB) were regarded as target genes of the 2 miRNAs. The Venn diagrams (Fig. 7) showed that 293 and 1068 target genes were identified for miRNA-146b-3p and miRNA-363–3p, respectively. In order to explore whether these genes were involved in the invasiveness of PTC, the target genes of miRNA-146b-3p were intersected with the down-regulated differentially expressed genes in PTC, while the target genes of miRNA-363–3p were intersected with the up-regulated differentially expressed genes in PTC. The results showed that a total of 15 and 73 target genes were identified for miRNA-146b-3p and miRNA-363–3p, respectively. The network between the two miRNAs and the target genes were visualized in Figure 8.

Figure 7

Venn diagrams of target genes.

Figure 8

The network between the miRNAs and the corresponding target genes.

Functional enrichment

The results of GO annotation of the 88 target genes associated with CLNM were shown in Figure 9. Ten biological process (BP) terms, 6 molecular function (MF) terms, and 8 cellular component (CC) terms were identified from GO categories. BP analysis showed that the target genes were enriched in extracellular structure organization, regulation of protein complex assembly, second-messenger-mediated signaling and so on. BP analysis showed that the target genes were enriched in molecular adaptor activity, passive/active transmembrane transporter activity, anion/metal ion transmembrane transporter activity, and monovalent inorganic cation transmembrane transporter activity. CC analysis indicated that these genes were mainly enriched in extracellular matrix (ECM), endoplasmic reticulum lumen, and synaptic membrane. Furthermore, the KEGG analysis (Fig. 10) indicated that these target genes were significantly enriched in pathways, including fatty acid metabolism, protein digestion and absorption, axon guidance, viral protein interaction with cytokine and cytokine receptor, and ABC transporters.

Figure 9

Gene Ontology analysis of the target genes associated with central lymph node metastasis: A, biological process term; B, molecular function term; C, cellular component term.

Figure 10

Kyoto Encyclopedia of Genes and Genomes pathway enrichment.

PPI network

To establish a PPI network, the target genes of the 2 miRNAs were analyzed using STRING database (Fig. 11). All of the 88 target genes were included in the PPI network, containing 88 nodes, and 38 edges. Top 5 hub genes, including SLC6A1, SYT1, COL19A1, RIMS2, and COL1A2, were illustrated using Cytoscape (Fig. 12).

Figure 11

Protein–protein interaction network of target genes associated with central lymph node metastasis.

Figure 12

Five hub genes associated with central lymph node metastasis.

Discussion

PTC, the predominant histologic type of differentiated thyroid cancer, is characterized by cervical LNM. Previous studies have confirmed that LNM is an unfavorable prognosis factor in PTC patients.[ LND can effectively decrease postoperative recurrence rate, but prophylactic central lymph node dissection (CLND) for clinically node-negative PTC remains a matter of debate.[ Preoperative neck ultrasonography is a reliable tool to detect LNM in the lateral neck, however its sensitivity and specificity for identifying central compartment lymphadenopathy are weakened.[ Therefore, alternative biomarkers are urgently required to assess lymph node status in the central compartment. This can decrease postoperative complications due to unnecessary neck dissection. In this study, 2 miRNAs significantly associated with CLNM in PTC were identified from TCGA database and a nomogram combining the risk score with 2 clinical parameters (age and ETE) was built and validated, which showed better prediction value than the other individual predictor in terms of CLNM identification.

MiRNA-146b-3p is derived from MIR146B gene, which is located on chromosome 10 at position q24.32. Yu C reported miRNA-146b-3p promoted PTC tumor metastasis by suppressing its target NF2.[ Han PA suggested miRNA-146b-3p could be predictive of CLNM preoperatively.[ Riesco-Eizaguirre G discovered that a miRNA-146b-3p/PAX8/NIS regulatory circuit might modulate thyroid cell differentiation and iodide uptake for improved treatment of advanced thyroid cancer.[ These results suggested that miRNA-146b-3p was significantly involved in PTC metastasis. In this study, miRNA-146b-3p is an independent predictor of CLNM, and this outcome can be mechanically explained by the above-mentioned basic researches. On the other hand, miRNA-363–3p is a component of the miRNA-106a-363 cluster on chromosome X. As reported, miRNA-363–3p exerted an anti-tumor effect on a variety of carcinomas. Chang found that miRNA-363–3p could inhibited migration, invasion, and epithelial-mesenchymal transition by suppressing NEDD9 and SOX4 in non-small-cell lung cancer.[ Dong reported that tumor growth and metastasis of colorectal cancer could be restrained by miRNA-363–3p targeting SphK2.[ For PTC, miRNA-363–3p was significantly associated with LNM, and over-expression of miRNA-363–3p could prohibit anoikis resistance of PTC cells via targeting ITGA6.[ Also, miRNA-363–3p was found to inhibit PTC progression by targeting PIK3CA.[ Consistently, in the present study, univariate and multivariate analyses showed that miRNA-363–3p served as a protective factor against CLNM in PTC. In addition to these 2 miRNA molecules, potential clinical risk factors of CLNM were also evaluated. In line with previous studies, age and ETE were confirmed to be predictive of CLNM,[ and were included in the final nomogram prediction model, which showed better performance than the other individual risk factor in terms of predicting CLNM.

To gain insight into the potential mechanism of the 2 miRNAs in CLNM of PTC, functional enrichment analysis and PPI network establishment of the identified target genes were performed. The GO annotation results indicated that the target genes of the 2 miRNAs were mainly involved in extracellular matrix (ECM) organization, ion transmembrane transporter activity, synaptic membrane and so on. And the KEGG enrichment analysis results revealed that these genes participated in pathways, such as fatty acid metabolism, protein digestion and absorption, axon guidance, viral protein interaction with cytokine and cytokine receptor, and ABC transporters. It is well known that ECM plays a vital role in tumor progression. Zhai performed weighted gene coexpression network analysis based on TCGA database, reporting that the regulation of the ECM was associated with LNM of PTC.[ Takeyama showed that the expression of an ECM molecule, fibronectin, was positively correlated with LNM of PTC.[ Furthermore, integrins, the principal receptors to bind ECM, were also associated with the progression of thyroid carcinoma.[ Ion transmembrane transporter activity can modulate the cell membrane potential, which is important for cell proliferation, differentiation and migration.[ A growing body of evidence suggested that the level of membrane potential had functional roles in tumor cell progression. Bikas reported that mitotane exerted its anticancer effects in thyroid cancer cells by inhibition of mitochondrial membrane potential and induction of apoptosis.[ Emerging data indicated that metabolite profiling was significantly changed in tumor progression. Previous studies showed that fatty acid metabolism might be associated with the neoplasia, invasion and metastasis of the colorectal tumor, nonsmall cell lung cancer, and cholangiocarcinoma.[ However, the roles of fatty acid metabolism in the CLNM of PTC are not well addressed. Protein digestion and absorption, axon guidance, viral protein interaction with cytokine and cytokine receptor, and ABC transporters pathways are all closely related to cancer initiation and progression. Axon guidance factor SLIT2 can inhibit the proliferation, migration and invasion of thyroid cancer cells by inhibiting transcriptional activity of beta-catenin and regulating Rho GTPase activity.[ Mato indicated that the cell subpopulation identified in TPC-1 expressed ABCG2/BCRP gene (belonging to ABC transporter family), probably contributing to the malignant progression of PTC.[ Nevertheless, the underlying functions of protein digestion and absorption, and viral protein interaction with cytokine and cytokine receptor in the progression of PTC should be further explored. Based on PPI network, 5 hub genes were screened, including SLC6A1, SYT1, COL19A1, RIMS2, and COL1A2. SLC6A1 and COL1A2 are unfavorable prognostic markers in renal cancer.[SYT1 promotes colon cancer cell proliferation, migration, and invasion.[COL19A1 encodes the alpha chain of type XIX collagen, but the function of this collagen is not well known. Presumably this collagen may serve to maintain the integrity of ECM, which is important for tumor invasiveness. RIMS2 encodes a presynaptic protein that interacts with RAB3, a protein important for normal neurotransmitter release. Although RIMS2 was reported to be mutated in melanoma, its function in carcinogenesis is still largely unknown.[ Regretfully, the roles of these genes in CLNM of PTC have not been reported. Future functional experiments are required to verify bioinformatics predictions in this study, which may further elucidate the mechanism of CLNM of PTC.

This study established a nomogram to predict CLNM of PTC combining clinical data and RNA-segueing data from TCGA database, which is the original novelty of the present study. The other strengths of this research include: a training group and a testing group were randomly generated to establish and validate the nomogram; the ROC curves were drawn to investigate the diagnostic performance between the nomogram and the other independent risk factor; functional enrichment analysis of the target genes of the 2 miRNAs was performed to provide ideas for further experimental researches. Admittedly, this study is limited by the fact that there was not an external validation set to further verify the nomogram outcomes; only bioinformatics analyses without experimental data might significantly weaken the conclusions in this study, however, the roles of the identified 2 miRNAs in PTC progression have been confirmed in previous studies, which may in turn enhance the power of our results to some extent. On the other hand, it is meaningful to investigate how this nomogram alters or improves the clinical outcomes in PTC patients. Preoperative quantitative analysis of miRNA-146b-3p and miRNA-363–3p levels using specimen from US-guided fine-needle aspiration are encouraged, which can be used to calculate individual miRNA-based risk score. If a total point of more than 110 was obtained in a PTC patient, the possibility of CLNM is over 50%. And prophylactic LND can be recommended in this patient, which may decrease locoregional recurrence. However, this is another topic and not the main purpose of the current study.

Taken together, on the basis of this study, we found age, ETE, miRNA-146b-3p, and miRNA-363–3p were independent predictors of CLNM in PTC. The nomogram including risk score, age, and ETE showed good prediction of CLNM in PTC, which may help surgeons individualize surgical plans. Future basic studies are needed to elucidate the roles of the 2 miRNAs-target genes networks in the progression of PTC.

Acknowledgments

Mingjun Wang wants to thank the great support from Sisi Liu for this study.

Author contributions

Conceived and designed the study: Mingjun Wang, Xiuhe Zou, Rixiang Gong, Tao Wei, Jingqiang Zhu, and Zhihui Li

Final approval of the version to be published: Jingqiang Zhu and Zhihui Li

Performed the experiments: Mingjun Wang and Rongjing Li

Revised the paper: Mingjun Wang, Jingqiang Zhu, and Zhihui Li

Wrote the paper: Mingjun Wang