Multi-omics analysis revealed TEK and AXIN2 are potential biomarkers in multifocal papillary thyroid cancer

Papillary thyroid carcinoma (PTC), the most common endocrine cancer, accounts for 80–85% of all malignant thyroid tumors. This study focused on identifying targets that affect the multifocality of PTC. In a previous study, we determined 158 mRNAs related to multifocality in BRAF-mutated PTC using The Cancer Genome Atlas. We used multi-omics data (miRNAs and mRNAs) to identify the regulatory mechanisms of the investigated mRNAs. miRNA inhibitors were used to determine the relationship between mRNAs and miRNAs. We analyzed the target protein levels in patient sera using ELISA and immunohistochemical staining of patients’ tissues. We identified 44 miRNAs that showed a negative correlation with mRNA expression. Using in vitro experiments, we identified four miRNAs that inhibit TEK and/or AXIN2 among the target mRNAs. We also showed that the downregulation of TEK and AXIN2 decreased the proliferation and migration of BRAF ( +) PTC cells. To evaluate the diagnostic ability of multifocal PTC, we examined serum TEK or AXIN2 in unifocal and multifocal PTC patients using ELISA, and showed that the serum TEK in multifocal PTC patients was higher than that in the unifocal PTC patients. The immunohistochemical study showed higher TEK and AXIN2 expression in multifocal PTC than unifocal PTC. Both TEK and AXIN2 play a potential role in the multifocality of PTC, and serum TEK may be a diagnostic marker for multifocal PTC.


Introduction
Thyroid cancer is the most common malignant endocrine tumor and over the last three decades, its incidence has increased continuously worldwide [1,2].
(2) intraglandular dissemination of a single common malignant clone [7,8]. Previously, we identified the overexpression of 158 mRNAs in multifocal BRAF ( +) PTCs, compared to unifocal BRAF ( +) PTCs [9]. However, none of the mRNAs in BRAF ( +) PTCs showed lower expression than that of the mRNAs in unifocal PTCs [9]. Additionally, none of the mRNAs showed significantly different expression between multifocal and unifocal BRAF (−) PTCs [9]. In addition, mRNAs that were overexpressed in multifocal BRAF ( +) PTCs were associated with Wnt-and pluripotency-related pathways, which might account for the difference between multifocal and unifocal BRAF ( +) PTCs [9]. The multifocality of PTC has a clinical impact on the treatment of PTC. Although multifocality was not considered as a risk factor in the American Joint Committee on Cancer/Union for International Cancer Control (AJCC/UICC) staging system [4] or American Thyroid Association (ATA) guidelines [4], multifocal PTC has an increased risk of lymph node metastasis, recurrence, and distant metastasis [10,11]. In addition, multifocality of PTC is associated with the decision making on the optimal extent of surgery for PTCs, if it was identified before the surgery. However, fine needle aspiration cannot be performed for multiple thyroid nodules in clinical settings. The extent of surgical management affects the prognosis of PTC in terms of recurrence [12] and survival [13]. Therefore, adequate information regarding the multifocality of PTC should be provided to improve decision making.
MicroRNAs (miRNAs), a class of non-coding RNAs (15-27 nucleotide RNA molecules), control the expression of mRNA post-transcriptionally by binding to the 3ʹ-untranslated region of mRNA by blocking translation or mRNA degradation [14,15]. As a large number of miRNAs have been discovered with the continuous progress of technology, the importance of miRNA-mRNA regulatory mechanisms in physiological and pathological states is slowly becoming highlighted [14,16]. Recently, studies of miRNAs have focused mostly on malignant neoplasms, and miR-NAs have been shown to play pivotal roles in various cancers by regulating the expression of their target mRNAs [17][18][19].
In this study, we hypothesized that miRNAs interact with overexpressed mRNAs in multifocal BRAF ( +) PTCs, regulating protein expression. Therefore, we investigated target miRNAs by exploring the following: (1) screening miRNAs that interact with mRNAs in multifocal BRAF ( +) PTCs, (2) validation of miR-NAs with functional assays, and (3) protein expression in blood samples from patients with PTCs.

Data acquisition
Clinical characteristics and gene expression data (mRNAs and miRNAs) for PTC were downloaded from the Genomic Data Commons Data Portal (https:// gdcportal. nci. nih. gov/). The Cancer Genome Atlas (TCGA) data were available without restrictions on publications or presentations according to TCGA publication guidelines. Patients were categorized according to BRAF mutation status. In addition, patients were divided into two groups according to the multifocality of PTC. Of the 237 patients with BRAF ( +) PTCs, 110 had multifocal PTC and 127 had unifocal PTC according to our previous research [9].

Cell culture
BCPAP (PTC cell line) harboring the BRAF mutation was purchased from DSMZ Korea and was maintained in Roswell Park Memorial Institute (RPMI) 1640 medium containing 10% fetal bovine serum (FBS) and 100 µg/mL penicillin-streptomycin. All cultures were incubated at 37 °C in the presence of 5% CO 2 .

siRNA transfection
Negative control (Bioneer, South Korea, #SN-1003) siRNA and target gene siRNA were purchased from Bioneer Corporation (Daejeon, South Korea). Sequences are in Table 1 To create the knockdown cell line, BCPAP cells were seeded at 8 × 10 4 cells per well in a 6-well plate in RPMI containing 10% FBS. The cells were transfected using DharmaFECT 1 (Thermo Fisher Scientific), as per the manufacturer's instructions, with 300 nM siRNA and incubated for 48 h.

Cell proliferation assay
Cells were seeded at a density of 5 × 10 3 cells per well in 96-well plates in RPMI containing 10% FBS. After transfection with siRNA for 72 h, BCPAP cell viability was measured using the Cyto X cell viability assay kit from LPS solution Corporation (South Korea). Cyto X (10 µL) was added to each well and incubated for 1 h in a CO 2 incubator. Optical density (OD) values were quantitatively measured at 450 nm using an enzyme-linked immunosorbent assay reader.

Real-time PCR
Total RNA was extracted using the RNeasy mini kit (Qiagen, Germany), miRNA was extracted using miRNeasy mini kit (Qiagen), and cDNA synthesis was performed using a cDNA synthesis kit (Smart gene, South Korea, # SG-CDNAC100). Real-time monitoring of PCR reactions was performed using the LightCycler TM system (Roche Applied Science) and SYBR Green Q-PCR Master Mix with Low Rox (Smart gene, South Korea, #SG-SYBR-ROXL). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal control and primer sequences for real-time PCR (Table 2).

Wound healing assay
In order to measure the cell migration during wound healing, 48 h after the transfection was carried out (in a 100 mm cell culture dish), BCPAP cells were replated in 24-well plates at a density of 3 × 10 5 cells per well. Twenty-four hours later, BCPAP cells were cultured in RPMI containing 10% FBS and treated with 1 µg/mL mitomycin-C (Sigma-Aldrich, USA) for 3 h and then wounded with a linear scratch using SPLScar ™ Scratcher (SPL Science, South Korea). The average extent of the wound area was evaluated by measuring the width of the wound using ImageJ software.

Three-dimensional (3D) spheroid formation assay
3D spheroid formation was examined by culturing the cells in 200 μL complete medium containing 300 cells in each well and cultured in an ultra-low attachment 96-well plate (Corning, USA, #7007). After 1, 3, and 5 days of incubation, spheroid formation was photographed using phase contrast microscopy (4 × magnification).

Western blot
Cells were washed with PBS, dissolved in radioimmunoprecipitation assay (RIPA) buffer supplemented, and centrifuged at 15,000×g for 10 min at 4 °C. Protein concentration was determined by the BCA protein assay (#23,227, Thermo Fisher Scientific) using bovine serum albumin (BSA) as the standard. The proteins were separated by SDS-PAGE and transferred to hybridization nitrocellulose filter membranes (Merck Millipore, USA). The membranes were blocked for non-specific binding with 5% BSA in Tris-buffered saline containing TBS-T (TBS with 0.1% Tween 20) for 1 h at room temperature and then incubated with specific primary antibodies (diluted 1:1000) in TBS-T at 4 °C overnight. After washing with TBS-T three times, the proteins were identified using appropriate secondary antibodies (diluted 1:2000 with 5% BSA). Chemiluminescence was detected with SuperSignal ™ West Dura Extended Duration Substrate (Thermo Scientific, USA, #34,075) and visualized using an Amersham Imager 680 (GE Healthcare, USA).

Immunohistochemistry
Human tumor tissue paraffin blocks were processed into sections and deparaffinized, followed by incubation with 3% hydrogen peroxide for 20 min to block the endogenous peroxidases. Next, the tissue sections were blocked with 1% bovine serum albumin for 30 min and then incubated with a polyclonal antibody against AXIN2 and TEK (ABclonal Technology; 1:100 dilution) for overnight at 4 °C. Immunoreaction was visualized using the EnVision detection system kit (Dako), and Mayer's hematoxylin solution was used to stain the nuclei. After staining, images were obtained using an Axio Scan Z1 Digital Slide Scanner (Zeiss, Germany) and analyzed using Zen Blue software (Carl Zeiss. Germany) (magnification: × 100).

Statistical analysis
Statistical differences in clinical variables were analyzed using the chi-square test. To determine the relationship between the 145 mRNAs and total miRNAs, we used the Spearman correlation method based on the Hmisc R package (Hmisc version 4.0-3 and R version 3.4.3). From the correlation analysis results, we selected targets with negative correlations of 0.1 or more, and searched the genes related to the selected miRNAs in Tarbase v.8 (http:// carol ina. imis. athena-innov ation. gr/ diana_ tools/ web/ index. php?r= tarba sev8) to confirm the correlation. We sorted the selected miRNAs by 13 pathways identified in previous studies. In vitro experimental data were statistical significance of the differences among groups was determined by a one-way analysis of variance (ANOVA) followed by Tukey's test for multiple comparison using the GraphPad Prism 5 software (GraphPad Software, La Jolla, CA) and differences were considered statistically significant at p < 0.05. A comparison of receiver operating characteristic (ROC) was performed to test the difference between the area under the curve of two ROCs from blood samples using MedCalc 19.8 (MedCalc Software Ltd, Ostend, Belgium).

Screening of potential regulatory miRNAs of 145 mRNAs related to multifocality
In the previous study, we reported the enriched genes and pathways in multifocal PTC compared to unifocal PTC [9]. Since multifocal PTC is a more aggressive cancer, we focused on the functional role and regulatory mechanisms of the genes based on their pathways.  To obtain more convincing potential regulatory mechanisms, we used multi-omics data from TCGA. 145 mRNAs in the enriched pathways and whole miRNAs were included in the correlation analysis. miRNAs that correlated negatively with mRNAs were discovered in each pathway as follows [9]: 13 (Fig. 1). To select highly relevant genes, the relationship between mRNAs and miRNAs was confirmed by using Tarbase (Table 3 and Additional file 3: Table S1). Combining the results of Tarbase and correlation analysis, we selected four miRNAs (miR21, miR34a, miR203, and miR362) as potential regulators of target mRNAs.

Selection of target genes through miRNA inhibitors treatment
To validate the relationship between target mRNAs and miRNAs, BCPAP cells were treated with inhibitors of miR21, miR34a, miR203, and miR362. The expression of each miRNA was suppressed after 24 h ( Fig. 2A). Four upregulated genes (TEK, AXIN2, ADAMTS9, and ADAMTSL2) were shown to be miRNA-regulated genes after treatment with miRNA inhibitors (Fig. 2B). Negative results are presented in Additional file 1: Figure S1.

Functional assays of target mRNAs through siRNA treatment
The knockdown efficiency was assessed by RT-PCR and western blotting, and the results showed that TEK and AXIN2 expression was knocked down effectively by siRNA in BCPAP cells (Fig. 3A, B and C, and Additional file 2). To study whether this downregulation could inhibit the proliferation of BCPAP cells, their proliferation was monitored. As shown in Fig. 3D, the proliferation rate had significantly decreased on day 3 in siRNA-transfected cells compared to the negative control. In the 3D spheroid formation assay, the BCPAP cells were transfected prior to the generation of spheroids, which were then allowed to grow for 3 days. After 3 days of treatment with TEK and AXIN2 siRNA, the diameter of 3D spheroids was smaller than that of the negative control group, indicating that the TEK and AXIN2 siRNA complexes could significantly inhibit cell growth (Fig. 3E).
To evaluate whether blocking TEK and AXIN2 gene expression affects the oncogenic behavior of BCPAP  . 2 Inhibition of miR21, miR34a, miR203, and miR362 in BCPAP cells. A miRNA inhibitor transfection reduced the expression levels of all four miRNAs by more than 40% in BCPAP as analyzed by qPCR. B AXIN2, ADAMTS9 and ADAMTSL2 mRNA expression was up-regulated by miR21 inhibitor (gray), TEK, AXIN2 mRNA expression was up-regulated by miR34a inhibitor (green). And AXIN2, ADAMTS9 and ADAMTSL2 mRNA expression was up-regulated by miR203 inhibitor (orange), AXIN2 mRNA expression was up-regulated by miR362 inhibitor (blue). *P < 0.05, **P < 0.01, ***P < 0.001, compared to the control cells, we performed a cell migration assay. In a woundhealing assay, the artificial wound gap in plates of the negative siRNA-transfected BCPAP cells was significantly narrower than that of the TEK-and AXIN2 siRNAtransfected BCPAP cells at 6 h (Fig. 3F). ADAMTS9 and ADAMTSL2 were excluded from the experiment because of the repeated heterogeneous results.

TEK and AXIN2 from patients with BRAF ( +) PTCs: multifocal vs unifocal
42 patients with BRAF ( +) unifocal PTCs and 48 patients with BRAF ( +) multifocal PTC were included in this analysis. The sizes of the largest tumors were not significantly different (P = 0.5710) between unifocal (1.6 ± 0.7 cm) and multifocal (1.7 ± 0.8 cm) PTCs. The level of TEK in serum from multifocal PTCs was higher than that from unifocal PTCs (P < 0.0001, Fig. 4A), while the level of AXIN2 from multifocal PTCs was lower than that from unifocal PTCs (P < 0.0001, Fig. 4B). To test the performance of TEK and AXIN2, the areas under the curves were compared, which were not significantly different (TEK 0.854; AXIN2 0.779) (Fig. 4C). In the IHC results, the TEK and AXIN2 expression in multifocal PTC are much higher than unifocal PTC (Fig. 4D).

Discussion
The results of the current study are summarized as follows: (1) 13 miRNAs interacted with mRNAs overexpressed in multifocal BRAF ( +) PTCs, (2) after validation with miRNA inhibitors, and functional assays, the mRNA expression of TEK and AXIN2 was associated with the multifocality of BRAF ( +) PTCs, (3) TEK and AXIN2 in blood samples of patients with multifocal PTCs were significantly different from those with unifocal PTCs. PTC is the most common form of thyroid cancer with a favorable prognosis; however, a minority of patients develop locoregional recurrence, which eventually leads to mortality in some of these patients [4]. Risk factors that affect the risk of recurrence in PTC include extrathyroidal extension, lymph node involvement, BRAF mutation status, tumor size, and sex [42]. Multifocal PTCs, defined as the presence of two or more than 2 anatomically separated tumor foci in the thyroid gland [43], have been associated with an increased risk of lymph node and distant metastases, as well as disease recurrence [44]. In clinical settings, physicians often encounter patients with multiple thyroid nodules; however, fine needle aspiration of multiple nodules is rarely performed. The selection of thyroid nodules for fine needle aspiration is mainly determined by the cancer probability from ultrasonographic  3). B, C Western blotting analysis of AXIN2, TEK knockdown efficiency in BCPAP cells. D Cell proliferation assay showed that AXIN2, TEK siRNA suppressed the proliferation of BCPAP cells. E After transfection with AXIN2, TEK siRNA or control siRNA, spheroid formation was evaluated at 3 days post-seeding. The spheroid area was measured using ImageJ software, encompassing the entire spheroid to the outside edge of spheroid. Scale bar indicates 1000 µm. F Wound healing assay showed that AXIN2, TEK knockdown resulted in lower migration capacity in BCPAP cells than in negative control cells. Scale bar indicates 1000 µm. (n = 4) Values are means ± SEM *P < 0.05, **P < 0.01, ***P < 0.001 findings, such as the content of the nodule, echogenicity, shape, margin, calcification, and vascularity [45]. According to the ATA guidelines [4], hemithyroidectomy is sufficient for unifocal PTCs without prior irradiation of the head and neck area, a history of familial thyroid cancer, and known lymph node metastasis. In contrast, total thyroidectomy can be recommended if the nodules are confirmed to be malignant in the bilateral lobes [46]. Therefore, decisions regarding the optimal extent of surgery for patients with multiple thyroid nodules should be made after careful consideration.
Based on the results of the correlation analysis of the TCGA data and the in vitro experiments, we identified the oncogenic role of TEK and AXIN2 through miRNA regulatory mechanisms. Although we did not find the exact role of miRNAs in PTC, we observed that 4 miR-NAs (miR21, miR34a, miR203, and miR363) inhibition increased the expression of AXIN2 while TEK is only affected by mi34a in BRAF ( +) PTC cells. Along with the discovery of the candidate regulatory mechanisms, four miRNAs and two mRNAs were also validated in this study. TEK (i.e., Tie2) is a receptor tyrosine kinase, which is mainly expressed in endothelial cells and controls vascular regeneration and stabilization [47]. The angiopoietin-Tie system is known to be involved in inflammation, metastasis, and lymphangiogenesis; therefore, multiple clinical trials were performed with selective Tie2 inhibitors [48]. However, the role of TEK in cancer cells remains unclear. Knockdown of TEK increased the proliferation and migration of clear cell renal cell carcinoma [49]; however, the overexpression of TEK in glioma cells was associated with tumor malignancy and drug resistance [50,51]. Therefore, TEK in BRAF ( +) PTC cells may affect cell proliferation, invasion, and multifocality.
Since AXIN2 acts as an inhibitor of canonical Wnt signals in normal cells, many studies have focused on elucidating their role as tumor suppressors. However, silencing AXIN2 decreases the invasive and metastatic characteristics of colon cancer [52]. In the current study, AXIN2 expression was higher in multifocal PTC than in unifocal PTC, and knockdown of AXIN2 inhibited cell proliferation in BRAF ( +) PTC cells. In addition, it was found that the expression of Wnt target genes, including AXIN2, differed based on the presence/absence of the BRAF mutation [53]. Our previous results showed that AXIN2 expression did not differ between BRAF (−) multifocal and unifocal PTC [9]. Taken together, these results suggest that the role of AXIN2 may be different for each cancer, and BRAF mutations may affect AXIN2 expression. In this study, the level of TEK in the serum from multifocal PTC patients was higher than that in unifocal PTCs, while the level of AXIN2 from multifocal PTCs was lower than that from unifocal PTCs. On the contrary, a lower level of AXIN2 was observed in multifocal PTC. The biological reason for this phenomenon may be the half-life of the protein [54,55] as well as posttranscriptional and post-translational modifications. According to a recent analysis, protein expression correlates with the corresponding mRNA level by 20-40%, and mRNA expression levels are not completely representative of the corresponding protein concentration [56,57]. Therefore, serum TEK and AXIN2 levels may provide information on multifocality in patients with BRAF ( +) PTC in a separate way.
Summarily, we used big database to obtain potential regulatory mechanism of target genes and highlighted the interaction of miRNAs with the expression of target genes and proteins in multifocal BRAF ( +) PTC. miRNA inhibition increased the mRNA expression of TEK and AXIN2. Serum TEK and AXIN2 levels may provide information on the multifocality of BRAF ( +) PTC. Although there are some limitations, this study provides evidences for the regulatory mechanism of the genes and their contribution to PTC multifocality for the first time.