Linc00941 Is a Novel Transforming Growth Factor β Target That Primes Papillary Thyroid Cancer Metastatic Behavior by Regulating the Expression of Cadherin 6

Abstract

Background:

Papillary thyroid cancers (PTCs) are common, usually indolent malignancies. Still, a small but significant percentage of patients have aggressive tumors and develop distant metastases leading to death. Currently, it is not possible to discriminate aggressive lesions due to lack of prognostic markers. Long noncoding RNAs (lncRNAs), which are selectively expressed in a context-dependent manner, are expected to represent a new landscape to search for molecular discriminants. Transforming growth factor β (TGFβ) is a multifunctional cytokine that fosters epithelial-to-mesenchymal transition and metastatic spreading. In PTCs, it triggers the expression of the metastatic marker Cadherin 6 (CDH6). Here, we investigated the TGFβ-dependent lncRNAs that may cooperate to potentiate PTC aggressiveness.

Methods:

We used a genome-wide approach to map enhancer (ENH)-associated lncRNAs under TGFβ control. Linc00941 was selected and validated using functional in vitro assays. A combined approach using bioinformatic analyses of the thyroid cancer (THCA)—the cancer genome atlas (TCGA) dataset and RNA-seq analysis was used to identify the processes in which linc00941 was involved in and the genes under its regulation. Correlation with clinical data was performed to evaluate the potential of this lncRNA and its targets as prognostic markers in THCA.

Results:

Linc00941 was identified as transcribed starting from one of the TGFβ-induced ENHs. Linc00941 expression was significantly higher in aggressive cancer both in the TCGA dataset and in a separate validation cohort from our institution. Loss of function assays for linc00941 showed that it promotes response to stimuli and invasiveness while restraining proliferation in PTC cells, a typical phenotype of metastatic cells. From the integration of TCGA data and linc00941 knockdown RNA-seq profiling, we identified 77 genes under the regulation of this lncRNA. Among these, we found the prometastatic gene CDH6. Linc00941 knockdown partially recapitulates the effects observed upon CDH6 silencing, promoting cell cytoskeleton and membrane adhesions rearrangements and autophagy. The combined expression of CDH6 and linc00941 is a distinctive feature of highly aggressive PTC lesions.

Conclusions:

Our data provide new insights into the biology driving metastasis in PTCs and highlight how lncRNAs cooperate with coding transcripts to sustain these processes.

Introduction

Papillary thyroid cancer (PTC) is the most common thyroid cancer (THCA) subtype, usually characterized by slow proliferation and rare distant metastatic spread, resulting in a favorable prognosis for most patients (1). However, ∼5% of PTCs behave aggressively and are associated with distant metastasis. PTCs are characterized by low mutational load. Although low-risk PTC could be sufficiently treated with an ipsilateral lobectomy, some patients undergo aggressive treatment regimens, including total thyroidectomy and radioactive iodine administration (2). The major clinical challenge in this setting remains to find molecular features that correlate with PTC aggressiveness to better stratify and manage patients. In the last decade, some markers have been proposed. Re-expression of telomerase is a common feature in tumor, somatic mutations in telomerase reverse transcriptase (TERT) coding region are uncommon, but germline and somatic mutations in noncoding TERT promoter region were found in many cancers. The mutations occurring −124 bp (C228T) and −146 bp (C250T), enhancing promoter activity, were described as cancer hotspot mutations and proposed, alone or in combination with other markers, as prognostic markers to discriminate aggressive tumors (3 –9). Dissecting the molecular bases underlying PTCs metastatic spreading is pivotal to implement predictive aggressiveness signatures and individualize patients' treatment.

The dual role of transforming growth factor β (TGFβ) in cancer is widely recognized. In the early phases of tumor development, TGFβ acts as repressor signal, restraining proliferation and favoring apoptosis. During tumor progression, on the contrary, TGFβ promotes metastatic spread by inducing epithelial-to-mesenchymal transition (EMT) (10 –13). We previously demonstrated that TGFβ partakes in the induction of aggressive phenotype in PTCs by promoting EMT and inducing a transcriptional program that support cellular invasiveness and motility as well as changes in the metabolic condition of cancer cells (14).

A major mediator of TGFβ signaling in PTCs is Cadherin 6 (CDH6), a mesenchymal type-II cadherin, which showed to be a marker of metastatic spreading and reduced survival probability in PTCs (13,15). Mechanistically, we demonstrated that CDH6 promotes EMT completion by inhibiting autophagy (15). This cadherin coordinates profound changes in cytoskeleton architecture and cell–cell interaction, which serve cancer cells to escape the tight epithelial organization and overcome the mechanical stress that invasion into adjacent sites requires. CDH6 also blocks autophagy by sequestering GABARAP and restraining the formation of autophagosomes and the negative effect of this catabolic process on cell motility (15).

The improved ability of understanding the genome revolutionized the idea about its functioning and brought to light the existence of an impressive number of functionally relevant noncoding transcripts (16). Among these, long noncoding RNAs (lncRNAs) are a class of noncoding RNAs >200 bp sharing characteristics with coding transcripts, including intron–exon organization and transcription by RNA Polymerase II (17). Enhancer (ENH)-associated lncRNAs are a peculiar class of lncRNAs transcribed from active ENHs. The current hypothesis is that they represent a further and finer layer of regulation of ENH activity in the gene expression processes (18). Indeed, the expression of ENH-associated lncRNAs correlates with transcription of specific protein-coding genes, implicating a role in their expression regulation (19,20).

LncRNAs expression is strictly regulated in space and time, and this specificity is thought to contribute to many physiological and pathological cellular processes (21).

The large number and high expression specificity of lncRNAs in cancer make these molecules promising candidates as diagnostic and prognostic markers. The role of lncRNAs in different tumor settings and the possibility to exploit these molecules, also in signatures composed of coding and noncoding transcripts, are currently under investigation. Despite the increasing number of studies on this topic, the role of lncRNAs in mediating TGFβ signaling during cancer progression is not completely known.

In this study, we explored the genome-wide landscape of lncRNAs under TGFβ regulation to define their role in mediating metastatic phenotype in THCA. Among the noncoding transcripts identified, we functionally characterized linc00941, showing that this molecule partakes in TGFβ signaling by regulating CDH6 expression, a crucial mediator of PTC invasiveness and metastatic spread.

Materials and Methods

Cell culture and transfection

B-CPAP, TPC1, Cal62, and 8505c cell lines were obtained from Dr. Santoro, University of Naples, Nthy-ori3-1 from Dr. Greco, Fondazione IRCCS-INT, Milan, FTC133 was purchased from Sigma-Aldrich (St. Louis, MO), and SW579 was purchased from ATCC (Manassas, VA). Cells were grown at 37°C and 5% CO₂ in Dulbecco's modified Eagle's medium (DMEM), RPMI or DMEM Ham's F12 1:1 (Gibco) with 10% fetal bovine serum (FBS). Twenty nanomolars of either small interference RNA oligonucleotide against linc00941 (exon 5; LOC40019 Silencer Pre-designed siRNA [small interfering RNA], Ambion; Thermo Fisher Scientific, Waltham, MA) and Silencer Select Negative Control (Thermo Fisher Scientific) were used to silence linc00941 and as a negative control, respectively. Lipofectamine RNAiMax and Lipofectamine 2000 (Thermo Fisher Scientific) were used for siRNAs and plasmids transfection, respectively. Plasmids used are listed in Supplementary Table S1. Plasmids used were pEGFP-C3_LC3B, p3XFlagCMV10_GABARAP, and pCDNA-3.1.

TGFβ induction and ChIP-seq

Nthy-ori3-1, TPC1, B-CPAP, and Cal62 cells were starved with 1% FBS growth medium overnight, then treated with TGFβ (Peprotech, Rocky Hill, NJ) for 24 hours. Chromatin immunoprecipitation (ChIP) on Nthy-ori3-1 was performed as previously described (22). In brief, 10⁷ Nthy-ori3-1 cells were crosslinked for 15 minutes with 1% formaldehyde, lysed and sonicated for 10 cycles (30 minutes ON, 30 minutes OFF) using Bioruptor Pico Sonicator (Diagenode, Liege, Belgium), to obtain 100–200 bp chromatin fragments. Chromatin was precipitated overnight using Magna ChIP Protein G Magnetic Beads (Merck, Darmstadt, Germany) and 1 μg of H3K27ac (Abcam, Trumpington, United Kingdom) or Rabbit IgG (Santa Cruz Biotechnology, Dallas, TX) antibodies. One percent of total chromatin was used as input. Samples were quantified at Qubit (Thermo Fisher Scientific), and the quality was evaluated by Bioanalyzer (Agilent, Santa Clara, CA). Library for sequencing was obtained following the TruSeq ChIP Sample Preparation protocol (Illumina, Inc., San Diego, CA) using 5–10 ng ChIP DNA as starting material. Triplicates were sequenced on Illumina NextSeq500 high-output cartridge (single stranded, reads length 75 bp—1 × 75).

RNA extraction, quantitative reverse transcription PCR, and RNA-seq

Matched normal and tumor fresh frozen patients' samples were obtained from the Research Tissue BioBank of the azienda unità sanitaria locale (AUSL)—istituto di ricovero e cura a carattere scientifico (IRCCS) of Reggio Emilia after written informed consent was obtained. The project was approved by the local Ethical Committee (Protocol No. 2014/0014425 of 06/05/2014).

RNA from patients' samples and from cell lines was extracted using TRIzol (Thermo Fisher Scientific) or by automatic extractor MaxwellRSC (simply RNA Cells kit; Promega, Madison, WI). Bio-Rad (Hercules, CA) iScript cDNA (complementary DNA) kit was used for reverse transcription, and Bio-Rad Sso Fast EvaGreen SuperMix was used to perform quantitative reverse transcription PCR (q-RT-PCR), read at Bio-Rad CFX96 Real Time PCR Detection System. Relative expression of target genes was calculated by the 2^−ΔΔCt method normalizing on the geometric mean of reference genes CyclophilinA (CYPA) or glucuronidase beta (GUSB). Primer sequences are listed in Supplementary Table S2. RNA-seq libraries were obtained starting from 500 ng of total RNA following Illumina TruSeq Stranded TotalRNA preparation protocol. Sequencing was performed using Illumina NextSeq high-output cartridge (double stranded, reads length 75 bp—2 × 75).

Next-generation sequencing analysis

Next-generation sequencing (NGS) was performed using a customized panel of genomic regions and sequenced by MiSeq Sequencer (Illumina, Inc.), as previously described (23).

Western blot and immunofluorescence

Cells (5 × 10⁵ cells) were lysed in Passive Lysis Buffer (Promega). Sodium dodecyl sulphate–polyacrylamide gel electrophoresis gels using 15–30 μg of lysates were used. Proteins were transferred by wet transfer method or by semidry method using Transblot Turbo (Bio-Rad). The membranes were blocked in 5% milk (Bio-Rad) and incubated with primary antibodies for 1 hour overnight. Primary antibodies used for the study were as follows: 1:5000 anti-β-actin (Sigma-Aldrich), 1:1000 antiphosphorylated polymerase II (Ser5), 1:1000 anti-p62/SQSTM1 (Abcam), 1:1000 anti-c-Jun, 1:1000 anti-LC3A/B, 1:1000 anti-GABARAP (Cell Signaling Technology, Danvers, MA), 1:1000 anti-β-catenin (BD Biosciences, San José, CA) (for more information, see Supplementary Table S3). Secondary antibodies were HRP-conjugated donkey antirabbit and sheep antimouse (NA934V and NXA931V; Used 1:5000; GE Healthcare, Chicago, IL). Membranes were incubated with secondary HRP-conjugated antibodies for 1 hour. Bio-Rad Clarity ECL was used for detection using ChemiDoc Imager (Bio-Rad) following the manufacturer's instructions. ImageJ software was used to normalize signals and quantify band intensities.

For immunofluorescence, 7 × 10⁴ cells were seeded in Chamber slides (Nunc; Thermo Fisher Scientific) 24 hours upon linc00941 silencing. After 24 hours, cells were fixed with 4% PFA, permeabilized with 0.1% Triton (Sigma-Aldrich), blocked with FBS-bovine serum albumin (BSA; Gibco, Thermo Fisher Scientific), and incubated with primary antibodies for 1 hour. The following primary antibodies were used: 1:100 anti-β-catenin (BD Biosciences), 1:100 anti-Flag M2 (Sigma-Aldrich). Antibody binding was revealed with secondary antimouse Alexa 555 or antirabbit Alexa 594-conjugated antibodies diluted 1:1000 (Invitrogen, Thermo Fisher Scientific) (for more information see Supplementary Table S3). Phalloidin conjugated with green fluorophore (Thermo Fisher Scientific) was used to mark actin filaments. DAPI (Thermo Fisher Scientific) was used to stain cell nuclei. ECLIPSE Ni fluorescent microscope was used for detection and analysis (Nikon Instruments, Melville, NY).

All antibodies used are listed in Supplementary Table S3.

Cell proliferation assay

Thirty hours after siRNA or control oligos transfection, 1.5 × 10³ TPC1 and 2 × 10³ B-CPAP were seeded into 96-well plates. IncuCyte NucLight Rapid Red diluted 1:500 for TPC1 and 1:250 for B-CPAP in complete medium was added to cells for staining. Cells were left to settle for 30 minutes at room temperature, then placed in IncuCyte S3 Live-Cell Analysis System (Essen BioScience, Ann Arbor, MI). Cell proliferation was monitored every 2 hours using the phase contrast confluence metric, or the labeling efficiency using the red fluorescence channel. To confirm IncuCyte data, trypan blue cell counting was performed. Thirty hours after siRNA transfection, 5 × 10³ cells were seeded into 96-well plates. Viable cells were counted after 24, 48, and 72 hours using trypan blue (Sigma-Aldrich) staining and the Countess Automated Cell Counter (Thermo Fisher Scientific).

Chemotaxis and invasion assays

Chemotaxis assay was performed using IncuCyte S3 Live-Cell Analysis System (Essen BioScience) following the manufacturer's protocol. TPC1 cells were transfected with siCT or silinc00941, seeded on IncuCyte ClearView inserts under different serum concentrations (no FBS, FBS 10%, FBS 20%), and counted on the top and on the bottom of invasion membrane overtime.

Thirty hours upon transfection, cells were plated on the IncuCyte ClearView inserts (Essen BioScience). Chemotaxis plate membranes were coated with 50 μg/mL matrigel diluted in DMEM without FBS. After 30 minutes of incubation at 37°C and 30 minutes at room temperature, excess matrigel was aspired. Cells (1 × 10³ cells/well) were seeded in DMEM without FBS. Cells were left to settle at room temperature for 30 minutes, then DMEM complemented with 10% or 20% FBS or without FBS was added. The ClearView Cell Migration plate was placed into the IncuCyte system and allowed to warm at 37°C for 15 minutes before scanning every 2 hours for 62 hours.

Invasion assay was performed as previously described (24) using BioCoat Matrigel Invasion Chambers (Corning, Tewksbury, MA). Thirty hours after transfection, 3 × 10⁴ cells were seeded into invasion chambers inserts in DMEM without FBS. The wells in which the inserts were placed were filled with DMEM without or with 10% FBS. After 24 hours, the excess matrigel and the noninvading cells were removed from the inserts. Invaded cells on the inserts' membranes were fixed with methanol and stained with crystal violet (Sigma-Aldrich). Invading cells were counted using ImageJ and normalized on control inserts following the manufacturer's instructions.

Cell fractionation

Cell fractionation was performed as previously described (25). Cells (15 × 10⁶ cells) were harvested, 15% of the supernatant was used for total RNA extraction (TRIzol Reagent; Thermo Fisher Scientific) and 15% for total protein extraction (PLB; Promega), and 70% of the supernatant was used to isolate cellular fractions. To isolate cytosolic fraction, cells were incubated in 250 μL buffer A (10 mM HEPES pH 7.9, 10 mM KCl, 1.5 mM MgCl₂, 0.5% NP40, RNAse and protease inhibitors), then centrifuged. Incubation time was 8 minutes for TPC1, 4 minutes for B-CPAP, and 2 minutes for Nthy-ori 3-1. The supernatant was equally divided for protein and RNA extraction. Nuclei pellets were washed and resuspended in 250 μL buffer C (20 mM HEPES pH 7.9, 25% glycerol, 0.42 M NaCl, 1.5 mM MgCl₂, 0.2 mM EDTA, RNAse, and protease inhibitors) and incubated 30 minutes on ice, vortexed, then centrifuged. The supernatant was equally divided for protein and RNA extraction. Chromatin pellet was resuspended in TRIzol.

Flow cytometry assays

For cell cycle analysis, 5 × 10⁵ to 1 × 10⁶ cells were transfected for 24 hours, then resuspended and incubated in Nicoletti solution (50 μg/mL PI, Triton-X 0.1%, sodium citrate 0.1%), then read on the BD FACS Canto II (BD Biosciences). Blank cells were also read. TPC1 cells siCT or silinc00941 were incubated with 1 μM 5-(and-6)-carboxyfluorescein diacetate, succinimidyl ester (CSFE), washed, suspended in culture medium, and seeded. At different time points, cells were washed and suspended in phosphate buffered saline 1 × , 0.5% BSA, 2 mM EDTA for reading on the BD FACS CantoII (BD Biosciences).

Statistical and bioinformatic analyses

Statistical significance was determined by two-tailed t test in all experiments with a minimum of three biological replicates. For RNA-seq, raw paired-end reads were assessed for their sequencing quality using FastQC software; low-quality sequences were discarded. Filtered reads were aligned to the human reference transcriptome (GRCh38.p12) using STAR 2.7 (26). Mapped reads were determined, and the expression at gene and transcript level was computed by RSEM (27). Differential expression analysis was performed with R package DESeq2 (28), using 0.05 FDR (false discovery rate) threshold. Genes significantly differentially expressed underwent enrichment analysis, performed on GO biological processes and KEGG pathways through enrichR package (29). Enriched pathways with Benjamini–Hochberg adjusted p-value ≤0.05 were selected for network-based modeling using networkx package (30) within python 2.7 environment and Cytoscape software.

For ChIP-seq analysis, reads were mapped on the human reference genome hg38, and filtered using BWA and Samtools. Peak calling was performed using Macs2 after Picard duplicates removal.

NGS mutational analysis

The sequences obtained were analyzed using VariantStudio Software (Illumina, Inc.) and the Integrative Genomics Viewer 2.3 (IGV) tool. Only mutations present in ∼5% of the total number of reads analyzed and observed in both strands were considered for mutational calls (23).

For the THCA—the cancer genome atlas (TCGA) analysis, we used R TCGAbiolinks package to download and analyze RNAseq data (workflow.type “HTSeq - FPKM”) and mi-RNA seq data (workflow.type “BCGSC miRNA Profiling”) from THCA-TCGA project, and expression differences were evaluated applying the Kruskal–Wallis test. For correlation analysis of linc00941 with BRAF and TERT promoter mutations, we used TCGA mutational data from TCGA, Cell 2014 (31).

Correlations and survival analysis were performed using R Corrplot and Survival package, respectively. For microRNAs (miRNAs) correlation analysis, Linc00941-correlated miRNAs were analyzed with DIANA-microT-CDS algorithm to predict miRNA targets and perform enrichment analysis (32).

Results

Linc00941 is a novel TGFβ target, and it is highly expressed in PTCs

To identify ENH-associated lncRNAs under the direct control of TGFβ, we used a ChIP-seq approach to map TGFβ-responsive transcriptional elements. Changes in the enrichment profile of H3K27ac, histone marker of active transcription, were investigated in Nthy-ori 3-1 upon TGFβ stimulation. Based on our previous data, Nthy-ori 3-1 is the most suitable model to study TGFβ-mediated EMT in thyroid system (13). Analysis of treated and untreated cells identified 22,775 peaks with significantly different levels of H3K27ac. Peak-to-target assignment identified 9677 transcriptional units potentially under the regulation of TGFβ. Of these, 1025 presented at least one peak enriched upon TGFβ treatment (Fig. 1A). Categorization showed that the vast majority of TGFβ-associated transcriptional units (70.5%, 723/1025) were protein-coding genes, 28% (287/1025) were noncoding RNAs, and 1.5% (15/1025) were pseudogenes (Fig. 1B). Of the noncoding transcripts, 38.7% (106/287) were noncoding or antisense RNAs, 28.8% (79/287) were lncRNAs, 28.1% (77/287) were miRNAs, and 4.4% (12/287) were small-nucleolar RNAs (snoRNAs) (Fig. 1C).

FIG. 1.

TGFβ-induced linc00941 expression in thyroid cell lines and PTC patients. (A) Venn diagram showing all genes correlated to enriched H3K27ac peaks resulted from ChIP-seq analysis in Nthy-ori 3-1 NT, TGFβ treated or both conditions. (B) Pie chart representing the different typologies of transcripts distribution among the genes correlated to enriched H3K27ac peaks in TGFβ condition. (C) Pie chart displaying the distribution of noncoding RNA typologies among the noncoding transcripts associated with H3K27ac enrichment in TGFβ condition. (D) q-RT-PCR of seven lncRNAs selected for validation in Nthy-ori 3-1 NT or TGFβ treated. (E) q-RT-PCR of linc00941 expression in normal thyrocyte-derived and thyroid carcinoma cell lines. (F) linc00941 expression in matched normal and tumor samples from PTC patients of THCA-TCGA dataset. (G) linc00941 expression in matched normal and tumor patient tissues from our institutional biobank. (H) Differential expression of linc00941 in normal thyroids, primary PTCs, and metastases (THCA-TCGA). (I) linc00941 expression in matched tumor and metastatic tissues (THCA-TCGA). *p < 0.05, **p < 0.01, ***p < 0.001. ChIP, chromatin immunoprecipitation; lncRNA, long noncoding RNA; NT, not treated; PTC, papillary thyroid cancer; q-RT-PCR, quantitative reverse transcription PCR; TCGA, the cancer genome atlas; TGFβ, transforming growth factor β; THCA, thyroid cancer.

To functionally prove that changes in the chromatin transcriptional status reflected different levels of expression of the identified targets, we selected a subset of the top scoring lncRNAs with enriched H3K27ac peaks in their locus upon TGFβ treatment (Supplementary Fig. S1A). q-RT-PCR confirmed that all the analyzed transcripts underwent expression changes, consistent with chromatin status measured by ChIP-seq in response to TGFβ, confirming the validity of our approach (Fig. 1D).

Among the lncRNAs analyzed, linc00941 has been associated with cancer progression and was chosen for further investigations (33,34). First, we aimed to verify that linc00941 was overexpressed in THCA cells. Analysis of linc00941 levels in a panel of THCA cell lines of differentiated PTCs (TPC1 and B-CPAP), anaplastic thyroid cancer (ATC; 8505c and Cal62, SW579), and follicular THCA (FTC133), the only cell line derived from metastasis, showed that this transcript was dramatically upregulated compared with normal thyrocytes (Fig. 1E). Interestingly, its expression in PTC and ATC cell lines was not further enhanced by TGFβ (Supplementary Fig. S1B). To confirm these observations, we took advantage of the THCA gene expression dataset of TCGA project. Analysis of 58 matched PTCs and normal thyroid samples showed that linc00941 is overexpressed in tumor tissue as compared with normal (Fig. 1F). We validated these data by q-RT-PCR analysis in a retrospective cohort of 12 matched PTCs and normal thyroid samples derived from the biobank at our institution (Fig. 1G). Next, to evaluate the role of linc00941 in THCA pathobiology, we interrogated the THCA-TCGA dataset to explore changes in the expression levels of this transcript and disease aggressiveness. Expression data of 502 PTCs, 58 normal tissues, and 8 metastasis samples were available in the THCA-TCGA dataset. Figure 1H, I shows linc00941 expression trend in these samples. Linc00941 expression in tumor was higher compared with normal thyroid tissues; this analysis showed a significant upregulation of this transcript in metastasis as compared with primary lesions (Fig. 1H, I). Analysis of matched primary and metastatic samples (n = 8) showed a higher expression of linc00941 in metastatic tissue in 75% of patients, confirming the correlation of this lncRNA with metastatic behavior of PTC, although the number of cases was too small for statistical significance (Fig. 1I).

The expression of linc00941 correlates with aggressive features in PTC patients

We considered the association of linc00941 with clinical features of PTC (Fig. 2 and Table 1). Linc00941 expression was not associated with PTC stage (Fig. 2A) and sex (Table 1), but associated with the follicular variant of PTC (0.083 ± 0.109) compared with classical (0.197 ± 0.236) and tall cell variants (0.167 ± 0.122) (Fig. 2B). Mutation data were available for 396 samples, with 5 (1.3%) having BRAF mutations different from V600E, which were not considered in further analysis. Interestingly, linc00941 expression was significantly higher in BRAF^V600E PTCs (0.210 ± 0.200 vs. 0.124 ± 0.254 in wild type samples, p < 0.0001), which was associated with extrathyroidal extension. We found no association between linc00941 expression and TERT promoter mutation status (Fig. 2C and Supplementary Fig. S1C). Samples with both mutations showed higher levels of linc00941 expression (Table 1) (Fig. 2C). The association of linc00941 expression with BRAF^V600E mutation may account for the observed higher expression of linc00941 in tall cell variants as compared with follicular and classical variants of PTCs (Fig. 2B). We also evaluated the correlation of linc00941 with BRAF^V600E and TERT promoter in 11 samples, 7 with BRAF^V600E mutation (63.6%) (Supplementary Table S4). In this samples set, we found significantly higher linc00941 expression in BRAF^V600E mutant samples (Supplementary Fig. S1D).

FIG. 2.

Linc00941 expression and clinical features. (A) Distribution of linc00941 expression among different PTC stages. (B) Distribution of linc00941 expression among different PTC histotypes. (C) Association of linc00941 expression with BRAF^V600E and TERT promoter mutations. (D) Linc00941 expression distribution in alive and dead patients. (E) Overall survival in PTC patients from THCA-TCGA subdivided into quartiles based on linc00941 levels: I quartile: linc00941 expression <25%, IV quartile: linc00941 expression >75%. Follow-up 24 months. *p < 0.05, **p < 0.01, ***p < 0.001. Where not specified, the correlations are not significant. ns, not significant.

Table 1.

Association Between linc00941 Expression and Clinical Features in Papillary Thyroid Cancer Patients

linc00941 expression	Mean ± SD		p
Sample
Normal [58]	0.058 ± 0.028	Normal vs. primary tumor	<0.0001
Primary tumor [502]	0.171 ± 0.213	Normal vs. metastasis	<0.0001
Metastasis [8]	0.255 ± 0.12	Primary tumor vs. metastasis	0.0137
Sex			0.8063
Male [135]	23.4 ± 26.018
Female [367]	22.686 ± 29.877
Tumor stage			0.09548
Stage I [281]	0.175 ± 0.22
Stage II [52]	0.152 ± 0.288
Stage III [112]	0.162 ± 0.158
Stage IV [55]	0.188 ± 0.203
Status at 2 years of follow-up			0.09548
Alive [117]	0.187 ± 0.222
Dead [4]	0.286 ± 0.156
Histology
Papillary classical [356]	0.197 ± 0.236	Classical vs. follicular	<0.0001
Papillary follicular variant [101]	0.083 ± 0.109	Classical vs. tall cell	0.349
Papillary tall cell variant [36]	0.167 ± 0.122	Follicular vs. tall cell	<0.0001
Other [9]	0.146 ± 0.196
BRAF and TERT promoter mutations
WT [144]	0.124 ± 0.254	BRAF^V600E vs. BRAF^V600E +TERT promoter	0.358
BRAF^V600E [211]	0.210 ± 0.200	BRAF^V600E vs. TERT promoter	0.009
TERT promoter [8]	0.075 ± 0.085	TERT promoter vs. BRAF^V600E +TERT promoter	0.545
BRAF^V600E +TERT promoter [28]	0.271 ± 0.279	WT vs. BRAF^V600E	<0.0001
		WT vs. BRAF^V600E +TERT promoter	<0.0001
		WT vs. TERT promoter	0.562

Significant values are reported in bold.

SD, standard deviation; WT, wild type.

Two years follow-up was available for 121 of the 502 THCA-TCGA PTCs analyzed. Of these, 4 patients were dead while 117 alive. Linc00941 expression was higher in deceased (0.286 ± 0.156) compared with alive (0.187 ± 0.222) patients (Fig. 2D), and patients with high expression of linc00941 (IV quartile) displayed worst overall survival than the low expressing ones (I quartile) (Fig. 2E). Together, these data indicate a potential association of linc00941 with aggressive features of PTCs.

Linc00941 regulates cellular proliferation and invasion

To provide functional validation of the association between linc00941 and PTC aggressiveness in the THCA-TCGA dataset, we used cellular models of thyrocyte-derived cell line to perform in vitro experiments. First, since lncRNAs function is determined by their cellular localization, we performed fractionation assay on TPC1, B-CPAP (Fig. 3A), and Nthy-ori 3-1 (Supplementary Fig. S1E) cell lines. q-RT-PCR of linc00941 levels showed that it was ubiquitously present in different cellular compartments, implying its possible activity in both nuclear and cytoplasmic compartments. To define the biological relevance of linc00941 in PTC, we targeted it by siRNAs (Supplementary Fig. S1F). Linc00941 silencing (silinc00941) led to a significant increase in cell proliferation compared with control (siCT) (Fig. 3B), which was observed to a greater extent in TPC1 than in B-CPAP cells. Cell cycle analysis by PI staining showed that cell cycle phases were not affected by linc00941 knockdown (Fig. 3C). CSFE labeling confirmed that linc00941 silencing leads to increased rate of cell division (Fig. 3B, D), suggesting a restraining function of this lncRNA on cellular proliferation.

FIG. 3.

Linc00941 silencing effects on cellular proliferation and invasion of PTC cell lines. (A) linc00941 localization in TPC1 and B-CPAP cell lines by fractionation. Transcription levels of linc00941 and Xist, used as positive control for chromatin fraction, were normalized to Glucuronidase Beta (GUSB). Data shown as FC on total (±SD). q-RT-PCR graphs show mean and SD of three independent experiments. A representative Western blot shows the distribution of PolII and β-actin as control for fraction purity. (B) Proliferation assay in TPC1 and B-CPAP cells. Cell count at each time point was normalized to initial cell number. (C) Cell cycle distribution by PI of TPC1 siCT and silinc00941 cells. (D) Representative FACS analysis of CSFE labeling of TPC1 siCT and silinc00941 cells. (E) Chemotaxis assay on TPC1 siCT and silinc00941 with No, 10% or 20% FBS. Graph shows bottom cell count normalized to initial top cell count. (F) Images and cell count of invading siCT and silinc00941 TPC1 cells. The assays were performed at least in three independent experiments. *p < 0.05, **p < 0.01, ***p < 0.001. CSFE, 5-(and-6)-carboxyfluorescein diacetate, succinimidyl ester; FBS, fetal bovine serum; FC, fold change; SD, standard deviation.

The ability to degrade extracellular matrix and invade surrounding tissues is a major feature of metastatic cells. Thus, we measured the effect of linc00941 silencing on cellular chemotaxis and invasiveness. Chemotaxis analysis using live imaging in TPC1 cells showed that silencing linc00941 resulted in the inability of cells to respond to chemical stimuli, independently from the concentration of chemoattractant used, compared with siCT cells (Fig. 3E). Moreover, invasion assay showed that silencing linc00941 expression significantly reduced the invasive potential of TPC1 cells (Fig. 3F). These features could not be evaluated in B-CPAP, which do not respond to chemical stimuli (Supplementary Fig. S1G).

Linc00941 expression correlates with coding genes and miRNAs implicated in tumor aggressiveness in PTC

To get a deeper understanding of the molecular mechanisms linked to linc00941 function in PTCs, we used an integrated approach, combining gene expression data from human samples and functional readouts from cell lines. First, we searched the TCGA-THCA dataset for genes and miRNAs expression, which correlated with linc00941 expression in human samples. We found 77 genes, 11 negatively and 66 positively correlated with linc00941 expression in PTC samples (Fig. 4A and Supplementary Table S5). Next, we performed RNA-seq analysis on TPC1 cell line transfected with scramble or linc00941 siRNAs to define genes potentially affected by this lncRNA. Comparative analysis identified 2452 differentially expressed genes upon linc00941 knockdown. Among these, 1080 genes were upregulated, and 1372 genes were downregulated with silencing of linc00941 expression compared with siCT (Fig. 4B). A subset of 37 genes (reported in Supplementary Fig. S2A) was validated in an independent set of experiments. All the genes tested (30 down- and 7 upregulated) were confirmed to be affected by linc00941 silencing (Supplementary Fig. S2A).

FIG. 4.

Linc00941 correlation patterns. (A) Correlation matrix of the 77 genes correlated to linc00941 expression in PTC patients (THCA-TCGA). In red inversely correlated, in blue directly correlated genes. *indicates significance, dot size, and color (legend on the right) are proportional to the correlation coefficient. (B) Volcano plot of the deregulated genes in TPC1 silinc00941 compared with siCT ranked by p-value and FC. Red dots are genes with significant p-value and Log₂ FC ≤−1 or ≥+1. (C) Deregulation level of the 77 genes correlated to linc00941 expression in PTC patients (THCA-TCGA) in our RNA-seq analysis in TPC1 cells upon silinc00941. (D) Visualization of linc00941-correlated genes (blue circles) and the pathways in which they are enriched (red circles). (E) Visualization of linc00941 inversely (I) and directly (D) correlated miRNAs clusters (yellow diamonds) and the pathways in which they are enriched (colored circles). Pathways falling in the same macrocategory are highlighted with the same colors. Arrows indicate CDH6 (red) and linc00941 (blue). CDH6, cadherin 6; miRNA, microRNA. For ease of reading all data, please refer to the online version of the figure.

GO and KEGG enrichment analysis was performed, and the pathways affected by linc00941 silencing are displayed and ranked by p-value in Supplementary Figure S2B. We identified several categories, including signal transduction, cell growth and death, cell–cell adhesion and cytoskeleton organization, folding, sorting and degradation, transport and catabolism, besides tumor-specific pathways.

Downregulated genes were involved in the TNF, PI3K-Akt, MAPK, mTOR, NF-Kβ, cell adhesion, apoptosis, and anoikis pathways. Interestingly, several members of the integrin family, playing key roles in TGFβ signaling transduction, and membrane-associated proteins regulating actin cytoskeleton, cell polarity, proliferation, and chemotaxis were found. Many of the downregulated targets upon linc00941 silencing are unfavorable prognostic markers in different cancer types (35).

Interrogating the databases for the genes found upregulated following linc00941 silencing, fewer pathways resulted enriched. Among these were cell cycle, p53 signaling pathway, and endocytosis. We also found CDC25B, required for entry into mitosis and several cyclins, which promote transition through G1/S and G2/M. We found several genes responsible for vesicle trafficking, membrane shaping and remodeling and cytoskeletal organization to be upregulated.

The 77 genes we found significantly associated with linc00941 in the THCA-TCGA dataset were affected by linc00941 silencing in THCA cells in the RNA-seq analysis. Indeed, the 11 genes inversely correlated with linc00941 expression were upregulated on linc00941 silencing versus siCT condition, while the 66 genes positively correlated with linc00941 expression were downregulated upon linc00941 silencing (Fig. 4A, C). To reconstruct the biological role of these genes, exploiting significant categories emerged from the enrichment pathways analysis of our RNA-seq data, we found that these 77 genes belong to pathways that fit well within the linc00941-associated phenotype (Fig. 4D). Three main functional clusters emerged from this analysis: (i) cell metabolism and vesicle trafficking, (ii) cell adhesion and extracellular matrix interaction, and (iii) cancer-related signal transduction. In addition, we identified 504 miRNAs that showed either a positive (n = 251, 49.8%) or negative (n = 253, 50.2%) correlation with linc00941 expression. Several studies reported that lncRNAs bind miRNAs sequestering them and blocking their activity (36 –38). Thus, we mapped the functional pathways in which the top scoring linc00941-associated miRNAs were enriched. Setting an absolute correlation value of 0.3, a total of 58 miRNAs were used to interrogate DIANA-mirPath v.3 (Fig. 4E). Thirty-five miRNAs had a positive correlation with linc00941 (group D: R ≥ 0.3) and 23 were inversely correlated (group I: R ≤ −0.3). The predicted pathways overlapped with the ones in Figure 4D, confirming what we observed with linc00941-associated coding genes. The three major clusters identified were as follows: (i) cell adhesion and extracellular matrix interactions, (ii) cancer-related signal transduction, and (iii) metabolic and catabolic processes.

Linc00941 regulates CDH6 to modulate cytoskeleton architecture and autophagy

We recently reported that CDH6 mediates TGFβ signaling, promoting PTC spreading and modulating the balance between autophagy and cellular architecture (13,15). Interestingly, CDH6 was one of the 77 genes of the linc00941-associated PTCs signature. CDH6 expression was positively correlated with linc00941 in the THCA-TCGA dataset (Fig. 4C) and significantly downregulated upon its silencing in PTC cells (Supplementary Fig. S2A). Thus, we reasoned that modulation of this protein could explain part of the linc00941-associated phenotype. To test this hypothesis, we assessed whether linc00941 silencing recapitulated the effects we observed upon CDH6 loss of function.

We previously observed that CDH6 silencing profoundly affects cytoskeleton structure and cell–cell junctions' organization in PTCs (15). Similarly, linc00941 silencing compared with siCT cells, cell membranes staining with β-catenin highlighted changes in cell–cell interactions, increased intimate connections and roof tile-like contacts among adjacent cells (Fig. 5A). Recently, Yan et al. described the interplay of linc00941 with ANXA2 and miR-34a to allow β-catenin nuclear translocation and EMT promotion in hepatocellular carcinoma (33). β-Catenin distribution analysis by immunofluorescence (Fig. 5A) and cellular fractionation (Fig. 5B) did not show significant alterations in β-catenin localization upon linc00941 silencing. However, ANXA2 was among the 77 linc00941-associated genes in PTC (Fig. 4A).

FIG. 5.

Linc00941 and CDH6 in autophagy. (A) Immunofluorescence staining of β-catenin (red), actin filaments (green), and merged in TPC1 siCT and silinc00941 cells. (B) Western blot analysis of β-catenin cellular localization in TPC1 siCT and silinc00941 cells. c-Jun and β-actin were used as positive controls for nuclear and cytoplasmic fractions, respectively (C) Western blot of autophagy-related proteins in TPC1 at different time points upon linc00941 silencing. (D) Immunofluorescence analysis showing colocalization of Flag-GABARAP (red) and EGFP-LC3 (green) in siCT and si-linc00941-transfected TPC1 cells. Nuclei were stained with DAPI (blue). (E) Quantification of the number of autophagic dots per cell. **p < 0.01. Magnification 400 × . Scale bars IF: 10 μm. All the assays were performed at least in three independent experiments. ImageJ was used for WB quantification. **p < 0.01. IF, immunofluorescence; WB, Western blot. For ease of reading all data, please refer to the online version of the figure.

Consistent with what is reported for CDH6 silencing, actin filaments staining demonstrated a peculiar reorganization of cytoskeleton in cells with linc00941 silencing compared with siCT cells. Linc00941 silencing induced actin assembly into thick fibers (stress fibers) spanning across the cytoplasm along the major cell axis, not visible in control cells (Fig. 5A). These thick and stable actin fibers are distinctive features of nonmotile cells, in line with the reduced invasion capability observed upon linc00941 silencing (39). Next, since we reported that CDH6 restrains autophagy, we evaluated linc00941's role in this process. At early time points (24–48 hours) after linc00941 silencing, autophagy markers LC3, p62/SQSTM1, and GABARAP levels increased compared with siCT. By contrast, at later time points (96 hours) LC3, p62/SQSTM1 showed a significant downregulation (Fig. 5C). This modulation is compatible with autophagosomes formation and subsequent degradation of cargoes. Autophagic machinery activation upon linc00941 silencing was confirmed by immunofluorescence analysis of autophagosomes by EGFP-LC3 and Flag-GABARAP staining (Fig. 5D). In siCT cells, both these markers showed diffuse and homogeneous distribution in the cytoplasm, but dramatically changed on linc00941 silencing, with both markers colocalizing into distinct dots highlighting autophagosome formation. The dots quantification is given in Figure 5E. Overall, these data indicate that linc00941 and CDH6 functionally cooperate along the same axis to potentiate PTC aggressive features.

Linc00941 and CDH6 expression are associated with PTC aggressiveness

A well-established mechanism by which lncRNAs post-transcriptionally control protein function is by interfering with miRNAs, acting as competing endogenous RNAs (ceRNAs).

We showed that linc00941 expression in PTCs is associated with a subset of cancer-related miRNAs, thus we reasoned that physical sequestration of CDH6-targeting miRNAs promotes CDH6 expression and activity. We searched the THCA-TCGA for miRNAs correlated with both transcripts. CDH6 expression was associated with 616 miRNAs, 405 of which (66%) were significantly correlated also with linc00941 (Supplementary Fig. S2C). Restricting the analysis to miRNAs with correlation coefficient R ≤ −0.3 or R ≥ 0.3, 47 miRNAs remained commonly associated to both transcripts (Supplementary Table S6 and Fig. 6A). Nineteen were inversely and 28 directly correlated to both linc00941 and CDH6 (Fig. 6B). Among these, miR-223 has been reported to target CDH6 (40). These data support the hypothesis that linc00941 regulates CDH6 expression.

FIG. 6.

Linc00941 and CDH6 coexpression in PTCs. (A) Venn diagram showing the number of miRNAs correlated (R ≤ −0.3 and R ≥ 0.3) to linc0041, CDH6, or both transcripts in the THCA-TCGA dataset. (B) Correlation coefficient of linc00941 and CDH6 with the 47 commonly correlated miRNAs. (C) Percentage of normal, metastasis, and tumor samples from THCA-TCGA distributed for the combined expression of CDH6 and linc00941 in quartiles. (D) Overall survival in PTC patients from THCA-TCGA subdivided into high or low linc00941 and CDH6 combined expression. Follow-up 24 months.

To translate these observations in a clinical setting, we explored whether the combined expression of CDH6 and linc00941 was associated with PTCs aggressiveness in the THCA-TCGA dataset. The expression of CDH6 and linc00941was subdivided into quartiles (I-II low, III-IV high expression). Then, we analyzed the percentage of patients for each subset distributed based on the expression of both (Supplementary Table S7 and Fig. 6C). 93% of normal samples were in the low expression quartiles for both markers, while 100% of metastases fell into the high expression quartile for at least one of the two markers. Tumor samples displayed a heterogeneous distribution; 41.3% of tumor samples expressed high levels and 31.4% low levels of both transcripts (Supplementary Table S7 and Fig. 6C). We performed overall survival analysis considering both CDH6 and linc00941 expression levels, the only three deaths in the THCA-TCGA cohort had high expression of both markers (Fig. 6D).

Because our experiments indicated that CDH6 and linc00941 are both target of TGFβ in thyroid-derived cells, we evaluated whether their expression correlated with TGFβ and receptors' expression in the THCA-TCGA dataset. The expression of both CDH6 and linc00941 was significantly correlated with the levels of TGFβ1, TGFβ2, TGFβ3, and TGFβR1 in PTC samples. In contrast, a poor correlation was observed with TGFβR2 (Supplementary Fig. S2D). Together, these observations, in line with our previous data, suggest a role for CDH6 and linc00941 in mediating tumor aggressiveness and the potential utility of these markers in improving PTCs risk-based stratification.

Discussion

The molecular mechanisms underlining PTC aggressiveness and metastatic potential remain poorly understood, limiting development of molecular-based risk stratification that could inform optimal patient management. One of the reasons for this is the rarity of highly aggressive PTCs and thus the difficulty of collecting significant cohorts to profile. While aggressiveness of some cancer types relies on high mutational burden and specific mutations in coding genes, this seems not true for PTCs (31).

We recently reported that highly metastatic PTCs are not characterized by a different overall mutational burden compared with nonmetastatic tumors (9). No mutations in coding genes were significantly associated with metastatic spreading and PTC patients' survival. By contrast, we and others reported mutations in noncoding elements like TERT promoter as main drivers of PTCs progression and dissemination (3 –9). These observations indicate that noncoding genome, including noncoding transcripts, may represent a promising approach to identify molecular discriminants that explain differential aggressive behavior.

LncRNAs are interesting molecules to investigate. Unbiased genome-wide searches discovered thousands of lncRNAs, and their total number is estimated to be triple the number of protein-coding genes (41). Based on their large number and expression specificity, lncRNAs are thought to contribute to many cellular processes, including oncogenic signaling. More than 8000 lncRNA genes were recently discovered to be highly specific in cancer, representing a potential reservoir for biomarkers and therapeutic targets (42,43). The present work seems to support the idea that lncRNAs are major players in modulating PTCs biology. Here, we reported that linc00941 potentiates PTC aggressiveness enhancing responsiveness to chemical stimuli and invasiveness capability. We also showed that linc00941 expression is associated with disease aggressiveness.

The role of lncRNAs in PTC initiation and progression has been reported by many investigators. Among these, MALAT1 (Metastasis-Associated Lung Adenocarcinoma Transcript 1) expression is higher in PTCs compared with normal tissues, but is significantly lower in ATCs, suggesting spatial-temporal expression regulation (44). LOC100507661 has been reported to promote THCA cell proliferation and to be associated with BRAF^V600E mutation (45). HOTAIR (HOX transcript antisense RNA) has also been reported to promote cancer cells colony formation, and its expression was associated with worse survival in patients with PTC (46).

Many lncRNAs were also found to have tumor-suppressive function in THCA. Low expression of MEG3 has been found to strongly correlate with lymph node metastases and increased invasion in TPC1 cells (36,47). LINC00271 downregulation has been reported to be associated with extrathyroidal extension and advanced PTC stage (48). NONHSAT037832 expression is negatively correlated with lymph node metastases in PTC patients (49). In a recent genome-wide association study, single nucleotide polymorphisms leading to the transcriptional deregulation of the lncRNA PTCSC3 were found significantly associated with increased risk of developing PTC. Indeed, this lncRNA has been described to play tumor suppression function in THCA (50).

Several lncRNAs have been proposed to partake TGFβ signaling in cancer (51,52). Here, we provide evidence that linc00941, previously linked to EMT execution in other cancer types (33,34,53), is a TGFβ target and may mediate TGFβ-dependent phenotype. We also showed that linc00941 plays a dual role in this context. Its silencing led to a significant promotion of cell proliferation, while completely abolishing cell response to chemical stimuli and invasive capability. This effect may seem paradoxical, but it is established that during metastatic spreading inhibition of proliferation represents a necessary condition to allow cancer cells to reorganize cellular architecture to resist mechanical stresses and allow cell movement (54 –56). Thus, by slowing cell division and promoting invasiveness, linc00941 affects both, namely metastatic behavior.

Many functionally characterized lncRNAs act by conditioning the activity or stability of specific proteins, bridging the interplay between coding and noncoding genome. Our data indicate that lin00941 controls or modulates the expression of a set of coding genes whose function is related to cellular organization and transduction of cancer-related signals. Among these, CDH6 is a linc00941-relevant target and a major mediator of its biological effect in PTC. At the interface between cellular structure and autophagy regulation, CDH6 candidates membrane adhesion molecules to be not just structural components but sensors of environmental stimuli and signaling centers for cellular pathways (15,57). We described, for the first time, that linc00941 controls CDH6 activity in PTCs and, through its modulation, restrains autophagy by promoting cytoskeleton rearrangements required for invasiveness. Their functional interplay and relevance for PTC progression are highlighted in the human PTC samples analyzed. We showed that the combined high expression of both these molecules is a distinctive feature of metastatic samples. We also showed that the expression of both CDH6 and linc00941 is positively associated with TGFβ1, 2, and 3 and TGFβR1 expression in PTC samples.

Many studies described the ability of cytoplasmic lncRNAs to sponge miRNAs. Sequestering these small transcripts, lncRNAs indirectly stabilize their target proteins, becoming an additional layer of post-transcriptional regulation (58). For instance, NEAT1 was described as upregulated in PTCs where it has established different functions, among which sponging miRNA-214 to aid THCA progression (37). Liu et al. showed how H19 expression in PTC cells promoted their capacity to proliferate and invade. Its mode of action includes miR-17-5p sponging to regulate YES1 oncogene levels (38). We did not provide experimental evidence that linc00941 effect on CDH6 regulation is mediated by its activity as miRNA sponge. However, we found that CDH6 and linc00941 expression in PTC patients is significantly associated with a specific subset of miRNAs. The biological pathways in which these miRNAs are involved overlapped with the ones predicted to be affected by linc00941 and are compatible with the biological function described for CDH6 in cancer. Experimental validation is required to support the idea that linc00941 may condition CDH6 transcript levels sequestering miRNAs. However, there are studies that support this hypothesis as linc00941 has been shown to act as ceRNA (33) and CDH6 expression has been shown to be inhibited by miR-223, which is predicted to be involved in our model (40).

In summary, our data provide new evidence on the role of linc00941 in PTC and on the possibility of exploiting integrative signature to improve patients' risk-based stratification to optimize their management.

Footnotes

Acknowledgments

We thank Marina Grassi for the great technical support to this study, the Research Tissue BioBank of the Azienda USL-IRCCS of Reggio Emilia for providing patients' samples, and all the colleagues of our laboratory for participating in stimulating discussion on this work.

Author Disclosure Statement

The authors declare that they have no competing financial interests.

Funding Information

This work was supported by the Italian Ministry of Health, RF-2016-02364167. M.G. is supported by Associazione Italiana Ricerca sul Cancro (AIRC) fellowship for Italy, year of submission 2017, project code 20933. A.C. is supported by AIRC Investigator Grant number AIRC IG21772.

Supplementary Material

Supplementary Figure S1

Supplementary Figure S2

Supplementary Table S1

Supplementary Table S2

Supplementary Table S3

Supplementary Table S4

Supplementary Table S5

Supplementary Table S6

Supplementary Table S7

References

Cabanillas

, McFadden

, Durante

. 2016. Thyroid cancer. Lancet, 388:2783–2795.

Filetti

, Durante

, Hartl

, Leboulleux

, Locati

, Newbold

, Papotti

, Berruti A

, MO Guidelines

Committee

. 2019. Thyroid cancer: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up. Ann Oncol, 30:1856–1883.

Liu

, Bishop

, Shan

, Pai

, Liu

, Murugan

, Sun

, El-Naggar

, Xing

. 2013. Highly prevalent TERT promoter mutations in aggressive thyroid cancers. Endocr Relat Cancer, 20:603–610.

Landa

, Ganly

, Chan

, Mitsutake

, Matsuse

, Ibrahimpasic

, Ghossein

, Fagin

. 2013. Frequent somatic TERT promoter mutations in thyroid cancer: higher prevalence in advanced forms of the disease. J Clin Endocrinol Metab, 98:E1562–E1566.

Melo

, da Rocha

, Vinagre

, Batista

, Peixoto

, Tavares

, Celestino

, Almeida

, Salgado

, Eloy

, Castro

, Prazeres

, Lima

, Amaro

, Lobo

, Martins

, Moura

, Cavaco

, Leite

, Cameselle-Teijeiro

, Carrilho

, Carvalheiro

, Máximo

, Sobrinho-Simões

, Soares

. 2014. TERT promoter mutations are a major indicator of poor outcome in differentiated thyroid carcinomas. J Clin Endocrinol Metab, 99:E754–E765.

Liu

, Qu

, Liu

, Sheng

, Shi

, Zhu

, Murugan

, Guan

, Yu

, Wang

, Sun

, Shan

, Teng

, Xing

. 2014. TERT promoter mutations and their association with BRAF V600E mutation and aggressive clinicopathological characteristics of thyroid cancer. J Clin Endocrinol Metab, 99:E1130–E1136.

Gandolfi

, Ragazzi

, Frasoldati

, Piana

, Ciarrocchi

, Sancisi

. 2015. TERT promoter mutations are associated with distant metastases in papillary thyroid carcinoma. Eur J Endocrinol, 172:403–413.

de Biase

, Gandolfi

, Ragazzi

, Eszlinger

, Sancisi

, Gugnoni

, Visani

, Pession

, Casadei

, Durante

, Costante

, Bruno

, Torlontano

, Paschke

, Filetti

, Piana

, Frasoldati

, Tallini

, Ciarrocchi

. 2015. TERT promoter mutations in papillary thyroid microcarcinomas. Thyroid, 25:1013–1019.

Gandolfi

, Ragazzi

, de Biase

, Visani

, Zanetti

, Torricelli

, Sancisi

, Gugnoni

, Manzotti

, Braglia

, Cavuto

, Merlo

, Tallini

, Frasoldati

, Piana

, Ciarrocchi

. 2018. Genome-wide profiling identifies the THYT1 signature as a distinctive feature of widely metastatic papillary thyroid carcinomas. Oncotarget, 9:1813–1825.

10.

Dongre

, Weinberg

. 2019. New insights into the mechanisms of epithelial-mesenchymal transition and implications for cancer. Nat Rev Mol Cell Biol, 20:69–84.

11.

Xue

, Restuccia

, Lan

, Hynx

, Dirnhofer

, Hess

, Rüegg

, Hemmings

. 2012. Akt/PKB-mediated phosphorylation of Twist1 promotes tumor metastasis via mediating cross-talk between PI3K/Akt and TGF-β signaling axes. Cancer Discov, 2:248–259.

12.

David

, Huang

Y-H

, Chen

, Su

, Zou

, Bardeesy

, Iacobuzio-Donahue

, Massagué

. 2016. TGF-β tumor suppression through a lethal EMT. Cell, 164:1015–1030.

13.

Sancisi

, Gandolfi

, Ragazzi

, Nicoli

, Tamagnini

, Piana

, Ciarrocchi

. 2013. Cadherin 6 is a new RUNX2 target in TGF-β signalling pathway. PLoS One, 8:e75489.

14.

Sancisi

, Borettini

, Maramotti

, Ragazzi

, Tamagnini

, Nicoli

, Piana

, Ciarrocchi

. 2012. Runx2 isoform I controls a panel of proinvasive genes driving aggressiveness of papillary thyroid carcinomas. J Clin Endocrinol Metab, 97:E2006–E2015.

15.

Gugnoni

, Sancisi

, Gandolfi

, Manzotti

, Ragazzi

, Giordano

, Tamagnini

, Tigano

, Frasoldati

, Piana

, Ciarrocchi

. 2017. Cadherin-6 promotes EMT and cancer metastasis by restraining autophagy. Oncogene, 36:667–677.

16.

Djebali

, Davis

, Merkel

, Dobin

, Lassmann

, Mortazavi

, Tanzer

, Lagarde

, Lin

, Schlesinger

, Xue

, Marinov

, Khatun

, Williams

, Zaleski

, Rozowsky

, Röder

, Kokocinski

, Abdelhamid

, Alioto

, Antoshechkin

, Baer

, Bar

, Batut

, Bell

, Chakrabortty

, Chen

, Chrast

, Curado

, Derrien

, Drenkow

, Dumais

, Duttagupta

, Falconnet

, Fastuca

, Fejes-Toth

, Ferreira

, Foissac

, Fullwood

, Gao

, Gonzalez

, Gordon

, Gunawardena

, Howald

, Jha

, Johnson

, Kapranov

, King

, Kingswood

, Luo

, Park

, Persaud

, Preall

, Ribeca

, Risk

, Robyr

, Sammeth

, Schaffer

, See

L-H

, Shahab

, Skancke

, Suzuki

, Takahashi

, Tilgner

, Trout

, Walters

, Wang

, Wrobel

, Yu

, Ruan

, Hayashizaki

, Harrow

, Gerstein

, Hubbard

, Reymond

, Antonarakis

, Hannon

, Giddings

, Ruan

, Wold

, Carninci

, Guigó

, Gingeras

. 2012. Landscape of transcription in human cells. Nature, 489:101–108.

17.

Ransohoff

, Wei

, Khavari

. 2018. The functions and unique features of long intergenic non-coding RNA. Nat Rev Mol Cell Biol, 19:143–157.

18.

Ørom

, Shiekhattar

. 2013. Long noncoding RNAs usher in a new era in the biology of enhancers. Cell, 154:1190–1193.

19.

Engreitz

, Haines

, Perez

, Munson

, Chen

, Kane

, McDonel

, Guttman

, Lander

. 2016. Local regulation of gene expression by lncRNA promoters, transcription and splicing. Nature, 539:452–455.

20.

Engreitz

, Ollikainen

, Guttman

. 2016. Long non-coding RNAs: spatial amplifiers that control nuclear structure and gene expression. Nat Rev Mol Cell Biol, 17:756–770.

21.

Kopp

, Mendell

. 2018. Functional classification and experimental dissection of long noncoding RNAs. Cell, 172:393–407.

22.

Sancisi

, Manzotti

, Gugnoni

, Rossi

, Gandolfi

, Gobbi

, Torricelli

, Catellani

, Faria do Valle

, Remondini

, Castellani

, Ragazzi

, Piana

, Ciarrocchi

. 2017. RUNX2 expression in thyroid and breast cancer requires the cooperation of three non-redundant enhancers under the control of BRD4 and c-JUN. Nucleic Acids Res, 45:11249–11267.

23.

de Biase

, Torricelli

, Ragazzi

, Donati

, Kuhn

, Visani

, Acquaviva

, Pession

, Tallini

, Piana

, Ciarrocchi

. 2018. Not the same thing: metastatic PTCs have a different background than ATCs. Endocr Connect, 7:1370–1379.

24.

Petrachi

, Romagnani

, Albini

, Longo

, Argenziano

, Grisendi

, Dominici

, Ciarrocchi

, Dallaglio

. 2017. Therapeutic potential of the metabolic modulator phenformin in targeting the stem cell compartment in melanoma. Oncotarget, 8:6914–6928.

25.

Rossi

, Pistoni

, Sancisi

, Gobbi

, Torricelli

, Donati

, Ribisi

, Gugnoni

, Ciarrocchi

. 2019. RAIN is a novel enhancer-associated lncRNA that controls RUNX2 expression and promotes breast and thyroid cancer. Mol Cancer Res, 18:140–152.

26.

Dobin

, Davis

, Schlesinger

, Drenkow

, Zaleski

, Jha

, Batut

, Chaisson

, Gingeras

. 2013. STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29:15–21.

27.

, Dewey

. 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12:323.

28.

Love

, Huber

, Anders

. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15:550.

29.

Kuleshov

, Jones

, Rouillard

, Fernandez

, Duan

, Wang

, Koplev

, Jenkins

, Jagodnik

, Lachmann

, McDermott

, Monteiro

, Gundersen

, Ma'ayan

. 2016. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res, 44:W90–W97.

30.

Hagberg

, Schult

, Swart

. 2008 Exploring network structure, dynamics, and function using NetworkX. Proceedings of the 7th Python in Science Conference (SciPy2008), Gäel Varoquaux, Travis Vaught, and Jarrod Millman.

31.

Cancer Genome Atlas Research Network. 2014. Integrated genomic characterization of papillary thyroid carcinoma. Cell, 159:676–690.

32.

Vlachos

, Zagganas

, Paraskevopoulou

, Georgakilas

, Karagkouni

, Vergoulis

, Dalamagas

, Hatzigeorgiou

. 2015. DIANA-miRPath v3.0: deciphering microRNA function with experimental support. Nucleic Acids Res, 43:W460–W466.

33.

Yan

, Zhang

, Wu

, Qian

, Hao

, Yan

, Zhu

, Wu

, Huang

, Luo

, Liu

, Chen

, Du

, Chen

, Fan

. 2017. Mesenchymal stem cells promote hepatocarcinogenesis via lncRNA-MUF interaction with ANXA2 and miR-34a. Cancer Res, 77:6704–6716.

34.

Liu

, Wu

, Zhang

, Zhong

, Zhang

, Guo

, Nie

, Liu

. 2019. Long non-coding RNA LINC00941 as a potential biomarker promotes the proliferation and metastasis of gastric cancer. Front Genet, 10:5.

35.

Uhlen

, Zhang

, Lee

, Sjöstedt

, Fagerberg

, Bidkhori

, Benfeitas

, Arif

, Liu

, Edfors

, Sanli

, von Feilitzen

, Oksvold

, Lundberg

, Hober

, Nilsson

, Mattsson

, Schwenk

, Brunnström

, Glimelius

, Sjöblom

, Edqvist

P-H

, Djureinovic

, Micke

, Lindskog

, Mardinoglu

, Ponten

. 2017. A pathology atlas of the human cancer transcriptome. Science, 357:eaan2507.

36.

Liu

, Yue

, Zhou

, Zhang

, Wang

, Chen

. 2018. LncRNA MEG3 enhances 131I sensitivity in thyroid carcinoma via sponging miR-182. Biomed Pharmacother, 105:1232–1239.

37.

J-H

, Zhang

S-Q

, Qiu

X-G

, Zhang

S-J

, Zheng

S-H

, Zhang

D-H

. 2017. Long non-coding RNA NEAT1 promotes malignant progression of thyroid carcinoma by regulating miRNA-214. Int J Oncol, 50:708–716.

38.

Liu

, Yang

, Zhu

, Li

, Lv

, Zhang

. 2016. Long noncoding RNA H19 competitively binds miR-17-5p to regulate YES1 expression in thyroid cancer. FEBS J, 283:2326–2339.

39.

Tojkander

, Gateva

, Lappalainen

. 2012. Actin stress fibers—assembly, dynamics and biological roles. J Cell Sci, 125:1855–1864.

40.

, Xu

, Song

, Xu

, Tai

, Goodman

, Bi

, Xu

, Jiao

, Maloney

, Wang

. 2018. miR-223-3p inhibits human osteosarcoma metastasis and progression by directly targeting CDH6. Mol Ther, 26:1299–1312.

41.

Hangauer

, Vaughn

, McManus

. 2013. Pervasive transcription of the human genome produces thousands of previously unidentified long intergenic noncoding RNAs. PLoS Genet, 9:e1003569.

42.

Iyer

, Niknafs

, Malik

, Singhal

, Sahu

, Hosono

, Barrette

, Prensner

, Evans

, Zhao

, Poliakov

, Cao

, Dhanasekaran

, Wu

Y-M

, Robinson

, Beer

, Feng

, Iyer

, Chinnaiyan

. 2015. The landscape of long noncoding RNAs in the human transcriptome. Nat Genet, 47:199–208.

43.

Schmitt

, Chang

. 2016. Long noncoding RNAs in cancer pathways. Cancer Cell, 29:452–463.

44.

Zhang

, Hardin

, Huang

, Chen

, Asioli

, Righi

, Maletta

, Sapino

, Lloyd

. 2017. MALAT1 long non-coding RNA expression in thyroid tissues: analysis by in situ hybridization and real-time PCR. Endocr Pathol, 28:7–12.

45.

Kim

, Lee

, Jeong

, Seol

M-Y

, Kim

K-S

, Lee

, Jo

. 2016. Upregulation of long noncoding RNA LOC100507661 promotes tumor aggressiveness in thyroid cancer. Mol Cell Endocrinol, 431:36–45.

46.

Zhu

, Lv

, An

, Shi

, Pan

, Zhou

, Yang

. 2016. Onco-lncRNA HOTAIR and its functional genetic variants in papillary thyroid carcinoma. Sci Rep, 6:31969.

47.

Wang

, Yan

, Zhang

, Jia

, Bu

. 2015. Long non-coding RNA MEG3 suppresses migration and invasion of thyroid carcinoma by targeting of Rac1. Neoplasma, 62:541–549.

48.

, Liao

, Wen

, Dong

, Zhou

, Yang

, Wang

, Ji

. 2016. Long intergenic non-coding RNA 271 is predictive of a poorer prognosis of papillary thyroid cancer. Sci Rep, 6:36973.

49.

Lan

, Sun

, Zhang

, He

, Dong

, Wang

, Liu

, Zhang

. 2016. Downregulation of long noncoding RNA NONHSAT037832 in papillary thyroid carcinoma and its clinical significance. Tumour Biol, 37:6117–6123.

50.

Jendrzejewski

, He

, Radomska

, Li

, Tomsic

, Liyanarachchi

, Davuluri

, Nagy

, de la Chapelle

. 2012. The polymorphism rs944289 predisposes to papillary thyroid carcinoma through a large intergenic noncoding RNA gene of tumor suppressor type. Proc Natl Acad Sci U S A, 109:8646–8651.

51.

Gugnoni

, Ciarrocchi

. 2019. Long noncoding RNA and epithelial mesenchymal transition in cancer. Int J Mol Sci, 20:1924.

52.

Farooqi

, Attar

, Qureshi

, Fayyaz

, Sohail

, Sabitaliyevich

, Nurmurzayevich

, Yelekenova

, Yaylim

, Alaaeddine

. 2018. Interplay of long non-coding RNAs and TGF/SMAD signaling in different cancers. Cell Mol Biol (Noisy-le-grand), 64:1–6.

53.

Luo

, Tao

, Zhang

, Zhu

, Minyao

, Haleem

, Dong

, Zhang

, Zhao

, Liao

. 2018. Regulatory network analysis of high expressed long non-coding RNA LINC00941 in gastric cancer. Gene, 662:103–109.

54.

Wang

, Goswami

, Lapidus

, Wells

, Wyckoff

, Sahai

, Singer

, Segall

, Condeelis

. 2004. Identification and testing of a gene expression signature of invasive carcinoma cells within primary mammary tumors. Cancer Res, 64:8585–8594.

55.

Hoek

, Eichhoff

, Schlegel

, Döbbeling

, Kobert

, Schaerer

, Hemmi

, Dummer

. 2008. In vivo switching of human melanoma cells between proliferative and invasive states. Cancer Res, 68:650–656.

56.

Schedin

, Eckel-Mahan

, McDaniel

, Prescott

, Brodsky

, Tentler

, Gutierrez-Hartmann

. 2004. ESX induces transformation and functional epithelial to mesenchymal transition in MCF-12A mammary epithelial cells. Oncogene, 23:1766–1779.

57.

, Yang

, Li

, Zhang

. 2019. Cadherin signaling in cancer: its functions and role as a therapeutic target. Front Oncol, 9:989.

58.

Poliseno

, Salmena

, Zhang

, Carver

, Haveman

, Pandolfi

. 2010. A coding-independent function of gene and pseudogene mRNAs regulates tumour biology. Nature, 465:1033–1038.