Identification of R-Spondin Gene Signature Predictive of Metastatic Progression in BRAFV600E-Positive Papillary Thyroid Cancer

Papillary thyroid carcinoma (PTC) is the most common malignancy of the thyroid gland and early stages are curable. However, a subset of PTCs shows an unusually aggressive phenotype with extensive lymph node metastasis and higher incidence of locoregional recurrence. In this study, we investigated a large cohort of PTC cases with an unusual aggressive phenotype using a high-throughput RNA sequencing (RNA-Seq) to identify differentially regulated genes associated with metastatic PTC. All metastatic PTC with mutated BRAF (V600E) but not BRAF wild-type expressed an up-regulation of R-Spondin Protein 4 (RSPO4) concomitant with an upregulation of genes involved in focal adhesion and cell-extracellular matrix signaling. Further immunohistochemistry validation confirmed the upregulation of these target genes in metastatic PTC cases. Preclinical studies using established PTC cell lines support that RSPO4 overexpression is associated with BRAF V600E mutation and is a critical upstream event that promote activation of kinases of focal adhesion signaling known to drive cancer cell locomotion and invasion. This finding opens up the potential of co-targeting B-Raf, RSPO and focal adhesion proteins as a pharmacological approach for aggressive BRAF V600E PTC.


Introduction
Thyroid cancer is the most common endocrine malignancy worldwide [1]. Among the four differentiated histological types, papillary thyroid cancer (PTC) represents 85% of thyroid gland cancers [2]. Epidemiological studies have shown gradual increased in PTC incidence over the past decades, particularly among young patients and women [3,4]. The main risk factors include exposure to radiation [5] and inherent genetic factors related to family history and ethnicity [4]. Although early stages of PTC are curable and the 10-year overall survival rates are >97%, there is a subgroup of high-risk PTC patients presenting frequently with lymph node metastasis [6], distant metastasis with high tumor recurrence rate [7][8][9].
PTC are staged according to the American Joint Committee on Cancer (AJCC) TNM system. Few years ago, the American Thyroid Association (ATA) proposed a new threetiered clinicopathological risk stratification system that classified patients as having a low, intermediate, or high risk of recurrence [10]. However, both systems do not adequately predict individual risk of recurrence or metastasis [11]. Methods to predict and/or stratify the risk of PTC recurrence are crucial to avoid overtreatment of low-risk and undertreatment of high-risk patients. It will also help to identify those patients who are at risk of recurrence over time in order to tailor their treatment plan accordingly.
Multi-omics experiments in thyroid cancer have allowed a comprehensive analysis and understanding of the molecular changes linked to tumor progression [12][13][14][15][16]. Mutations in BRAF and in particular V600E are found in approximately 50% tumors, but there is no consensus on whether this biomarker is an independent predictor for recurrence risk of PTC [17][18][19][20]. Consequently, the 2015 ATA guideline does not routinely recommend BRAF V600E alone for initial postoperative risk stratification [11]. However, an increasing number of studies show that genetic signatures may supersede tumor size and other classical staging criteria in predicting the presence of aggressive pathologic features in PTC [21,22].
PTC recurrence is estimated to occur at a rate of 4% to 28% [9]. The primary goal of this study was to identify transcriptomic differences in a unique cohort of rare PTCs with an unusual aggressive phenotype presenting with extensive lymph node metastasis and higher incidence of locoregional recurrence (metastatic), compared to non-metastatic PTC. All the metastatic PTC had BRAF mutation and up-regulation of R-Spondin protein (RSPO4). Genes specifically deregulated in BRAF-mutated PTC was investigated to determine their implication in the metastatic signaling and further validated in a large independent set of PTC samples. Preclinical studies using a panel of metastatic PTC cell lines supported that RSPO4 overexpression associated with BRAF V600E mutation can promote activation of focal adhesion kinases signaling in metastatic PTC.

Study Population
After Ethics Review Board approval (#13-093), patients diagnosed with PTC were prospectively enrolled in the fresh tissue biobank of Sir Mortimer B. Davis-Jewish General Hospital in Montreal, QC, Canada. Patients were consented at the time of surgery for use of their samples. Full access to their medical records was granted for research purposes. Detailed demographics (clinical stage, lymph nodes involvement, histological subtype and variant, extra-thyroidal extension, multifocality, and surgical margins), treatment, follow-up and recurrence data were obtained retrospectively. Clinicopathological data were handled in a coded fashion according to ethical guidelines. Exclusion criteria were papillary microcarcinoma, pregnancy, age less than 18 years, and missing follow-up. Strengthening the reporting of observational studies (STROBE Statement) was used to ensure appropriate methodological quality (http://www.strobe-statement.org/, accessed on 1 January 2021).

Sample Preparation and RNA Isolation
RNA samples were obtained from a unique cohort of 20 fresh-frozen PTCs (10 metastatic vs. 10 non-metastatic) followed-up for over than 10 years. Fresh-frozen tumor material was collected at the time of surgical resection, snap-frozen in liquid nitrogen, and stored at −80 • C until RNA extraction. Time between tumor resection and storage in liquid nitrogen did not exceed 45 min. Metastatic and non-metastatic samples were matched for age, sex and tumor stage (TNM) to assure comparability of the two groups (Table 1). RNA extraction was done using mRNeasy Kit (Qiagen, Inc., Valencia, CA, USA) according to the manufacturer's instruction. All RNA samples were assessed for quality using the RNA 6000 Nano assay on the 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA) and for quantity by Nanodrop (Thermo Scientific, Wilmington, DE, USA).

Bioinformatic Analysis
Sequence alignment and quantification of exon expression were carried out using an internally developed RNA-Seq analytical pipeline [23] Reads were trimmed from the 3 end to have a phred score of at least 30. Illumina sequencing adapters were removed from the reads, and all reads were required to have a length of at least 32. Trimming and clipping were done with the Trimmomatic software [24]. The filtered reads were aligned to a reference genome (hg19). The alignment was done with the combination of tophat/bowtie software [25]. The Cufflinks program [26] was used to assemble aligned RNA-Seq reads into transcripts and to estimate their abundance (FPKM). Several metrics and exploratory analysis to control the data quality and to verify the biological reliability of the data were used and are presented in the first section of the results. The differential gene expression analysis was done using DESeq [27] and edgeR [28] R Bioconductor package. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses was implemented using the clusterProfiler package [29].

CRISPR-Cas9 for the Generation of RSPO4 and BRAF Knockout Cells
Five target guide sequences for human RSPO4 and BRAF knockdown (Table 2) were cloned into lentiCRISPRv2 vector. Metastatic and non-metastatic thyroid cancer cell lines (8505c, THJ-16T were then transduced by lentiviral particles. The wide-type control was only the lentiCRISPR vector. The transduced cells were selected with puromycin at 1 µg/mL. RSPO4 and BRAF knockout clones were confirmed by qRT-PCR and immunoblotting.

qRT-PCR
Total RNA was isolated from control and knockdown THJ-16T, 8505c, and TPC-1 cell lines. cDNA was synthesized from 500 ng of RNA using Superscript II reverse transcriptase (Invitrogen ® , Carlsbad, CA, USA) and oligo-dT primers (Invitrogen ® , CA, USA ). qRT-PCR was be performed in the ABI PrismTM 7900 (Applied Biosystems ® , Foster City, CA, USA) using SYBR ® Green (Applied Biosystems ® , Foster City, CA, USA) in a 10 µL total volume and quality controls were used as proposed by MIQE Guidelines [30]. Specific primer set sequences was designed ( Table 3). The reactions were carried out in triplicate. GeNorm (https://genorm.cmgg.be/, accessed on 1 January 2022) algorithm were used to determine the most stable genes. The software packages were used as excel add-ons. ACTB, GAPDH and HRPT1 was selected by the algorithm as endogenous control. Fold differences in the relative gene expression were calculated using Pfaffl model [31]. Table 3. Specific primer set sequences to evaluate and validate the gene expression in this study.

Cell Motility and Invasion Assays
Invasion was measured by evaluating the migratory cell rate through a polycarbonate membrane (8-µm pore diameter) coated with BD Matrigel Matrix (BD Biosciences, Bedford, MA, USA) in a modified Boyden chamber (Corning, cat. no. 3422) as previously described [30]. Collective cell migration was evaluated using qualitative wound-healing assay [30]. Three independent experiments were performed for each condition and results were expressed as mean ± SD. Statistical analysis was done using Student's t test.

Immunohistochemistry (IHC) in the PTC Tissue Microarray (TMA)
Immunohistochemistry reaction was carried out on the TMA with 171 samples [32]. In brief, the slides were incubated with primary antibodies diluted in PBS overnight at 4 • C using: anti-RSPO4 (PA525615, Thermo Fisher Scientific Inc, Cleveland, OH, USA; 1:50) and anti-BRAF (20899-1-AP, Proteintech Group Inc, Chicago, IL, USA; 1:200). Sections were incubated with secondary antibodies (Advanced TM HRP Link, DakoCytomation, Glostrup, Denmark) for half-hour followed by the polymer detection system (Advanced TM HRP Link, DakoCytomation, Glostrup, Denmark) for half-hour at room temperature. Reactions were developed using a solution of 0.6 mg/mL of DAB (Sigma, St Louis, MO, USA) and 0.01% H 2 O 2 and then counter-stained with hematoxylin. Positive controls were included in all reactions in accordance with manufacturer's recommendations. Negative control consisted in omitting the primary antibody and replacing the primary antibody by normal serum. IHC reactions were replicated on distinct TMA slides to represent different tissues levels in the same lesion. The second slide was 25-30 sections deeper than the first slide, resulting in a minimum of 300 µm distance between sections representing 4-fold redundancy with different cell populations for each tissue.
Two independent certified pathologists conducted the IHC analysis blindly to the clinical data. Cores were scanned in 10× power field to settle on the foremost to marked area predominant in a minimum of 10% of the neoplasia [32,33]. IHC reaction was considered as positive if of a clearly visible dark brown precipitation occurred. IHC analysis considered the percentage and intensity of staining as: 0 (no detectable reaction or little staining in <10% of cells), 1 (weak but positive IHC expression in >10% of cells) and 2 (strong positivity in >10% of cells) [33].

Statistical Analysis
For in vitro analysis, statistical analyses were performed using the two-tailed Student's t-test for unpaired samples to one independent experiment and one-way analysis of variance (ANOVA) with post hoc comparisons based on the Tukey's multiple comparisons to two independent experiments. The analysis was carried on GraphPad Prism 5 software (GraphPad Software Inc., San Diego, CA, USA). Data were presented as mean ± SD, and the level of significance considered p ≤ 0.05. SD indicates dispersion of the data from mean and it was used to summarized the variability. For dichotomous variables, the two-sided Fischer Exact test was used to compare proportions. Odds ratio (OR) and 95% confidence interval (95% confidence interval; CI) were calculated according to the Mantel-Haenzel method. For discrete variables showing normal distribution, means and standard errors of means (SEM) are given and comparisons were made using the t-test. SEM was used to estimate the dispersion of the mean (average) of the data to be true among the population mean and it was limited to compute (CI) which measures the precision of population estimate. Alternatively, median, first quartile (Q25), and third quartile (Q75) are indicated and the non-parametric Kruskal-Wallis test was used. Statistical analyses were performed using SPSS ® 21.0.0 software (IBM©, Armonk, NY, USA).

Functional Analysis of Genes Associated with Metastatic PTC
Function annotation of DEGs was done using clusterprofiler [29]. Gene set enrichment analysis (GSEA) was done using publicly available knowledge of databases such Gene Ontology (GO) terms [23] or Kyoto Encyclopedia of Genes and Genomes (KEGG) [24]. A p-value of ≤ 0.05 and FDR < 0.05 are considered strongly enriched. (Figure 2A). The

Functional Analysis of Genes Associated with Metastatic PTC
Function annotation of DEGs was done using clusterprofiler [29]. Gene set enrichment analysis (GSEA) was done using publicly available knowledge of databases such Gene Ontology (GO) terms [23] or Kyoto Encyclopedia of Genes and Genomes (KEGG) [24]. A p-value of ≤0.05 and FDR < 0.05 are considered strongly enriched. (Figure 2A). The GO and KEEG gene set enrichment pathways analysis revealed the largest gene ratio differences in molecular function (MF), cellular component (CC), biological process (BP) and signaling pathways are linked to cytokine-cytokine receptor interaction (CCRI) and focal adhesion signaling (FA) (Figure 2A,B). When we combined MF, CC, BP and KEEG analysis, FA genes were overlapped in more than one network, such as the cell adhesion molecules and ECM-receptor interaction pathways (Figure 2A). These genes were clearly able to distinguish between metastatic vs. non-metastatic groups ( Figure 2B), suggesting that FA process had disrupted key genes in metastatic PTC. Hence, we decided to focus on this pathway and analyze the metastatic potential in PTC.

RSPO4 Is the Most Significantly Overexpressed Gene in Metastatic BRAF V600E PTC
In the unsupervisioned analysis ( Figure 1C), a gene called RSPO4 was highly overexpressed in metastatic PTC (<5% FDR, p < 0.05, log2FC > |2|). Based on the lists of pathways and networks potentially enriched in metastatic PTC (Figure 2A). RSPO4 showed to

RSPO4 Is the Most Significantly Overexpressed Gene in Metastatic BRAF V600E PTC
In the unsupervisioned analysis ( Figure 1C), a gene called RSPO4 was highly overexpressed in metastatic PTC (<5% FDR, p < 0.05, log2FC > |2|). Based on the lists of pathways and networks potentially enriched in metastatic PTC (Figure 2A). RSPO4 showed to be a promising candidate. RSPO4 is a gene encoding a secretory protein that is a member of the R-spondin family of proteins [34]. Within this family, they share a common domain organization consisting of a signal peptide, cysteine-rich/furin-like domain, thrombospondin domain and a C-terminal basic region. There is no published study showing RSPO4 upregulation in aggressive thyroid cancer. The status and function of this protein is unknown in PTC. Though it has been found to cause anonychia congenital, a condition characterized by the absence of fingernails and toenails [34]. Public databases, literature and stringDB was explored to find genes associated with RSPO4 in order to connect with the information from our results ( Figure 3A). In addition to RSPO4, other family members (specifically RSPO1 and RSPO3, not RSPO2) were differentially expressed in metastatic PTC, but only RSPO4 showed the highest statistical significance that offered an unprecedented opportunity for further investigation ( Figure 3B). RSPO proteins are known to bind to leucine-rich repeat-containing G-protein-coupled receptor (LGRs). In our study, it was confirmed that LGR4, LGR5, and LGR6 similarly had high expression ( Figure 3B). is unknown in PTC. Though it has been found to cause anonychia congenital, a condition characterized by the absence of fingernails and toenails [34]. Public databases, literature and stringDB was explored to find genes associated with RSPO4 in order to connect with the information from our results ( Figure 3A). In addition to RSPO4, other family members (specifically RSPO1 and RSPO3, not RSPO2) were differentially expressed in metastatic PTC, but only RSPO4 showed the highest statistical significance that offered an unprecedented opportunity for further investigation ( Figure 3B). RSPO proteins are known to bind to leucine-rich repeat-containing G-protein-coupled receptor (LGRs). In our study, it was confirmed that LGR4, LGR5, and LGR6 similarly had high expression ( Figure 3B).

RSPO4 Can Serve as a Robust Biomarker of Metastatic PTC
To validate the differentially expressed genes that can discriminate between metastatic and non-metastatic PTC from our study, we used a public cancer database (The Cancer Genome Atlas-TCGA, Figure 4) to confirm our findings in an independent data set. The TCGA contains 496 PTC samples of which 240 are primary non-metastatic tumors (M0), 6 metastatic PTC (M1), and 165 cases were undefined regarding the metastatic status (MX). The survival analysis showed that high BRAF (Figure 4) expression was correlated with decreased survival rates both in metastatic and non-metastatic. (Figure 4). Though in the TCGA dataset RSPO4 was also highly expressed in the metastatic samples (p < 0.005).

RSPO4 Can Serve as a Robust Biomarker of Metastatic PTC
To validate the differentially expressed genes that can discriminate between metastatic and non-metastatic PTC from our study, we used a public cancer database (The Cancer Genome Atlas-TCGA, Figure 4) to confirm our findings in an independent data set. The TCGA contains 496 PTC samples of which 240 are primary non-metastatic tumors (M0), 6 metastatic PTC (M1), and 165 cases were undefined regarding the metastatic status (MX). The survival analysis showed that high BRAF (Figure 4) expression was correlated with decreased survival rates both in metastatic and non-metastatic. (Figure 4). Though in the TCGA dataset RSPO4 was also highly expressed in the metastatic samples (p < 0.005).

RSPO4 Overexpression Occurs Selectively in BRAF V600E PTC
Since both RSPO4 and BRAF gene was highly expressed in our metastatic cohort (Figure 5A), we used Integrative Genome Viewer (IGV-https://software.broadinstitute.org/software/igv/ accessed on 1 January, 2022) to screened for common mutations described in thyroid cancer to identify those impacted on our target genes identified to discriminate between non-metastatic and metastatic PTC cases. Among our 20 samples, 10 harbored an A to T substitution at position 140,453,136 on chromosome 7, that is harboring a BRAF V600E classical mutation ( Figure 5B). Interestingly, BRAF mutated samples showed overexpression of RSPO4 (FPKM > 2 for BRAF V600E while FPKM < 1 for BRAF wt). Overall, the relative fold-change between BRAF V600E and BRAF wt was about 10 and statistically highly significant (p = 6.27 × 10 −7 ). This result was confirmed by qPCR and IHC in 170 PTC samples ( Figure 5C). In addition, the independent TCGA platform comparing BRAF V600E to BRAF wt papillary PTC found similar gene expression patterns ( Figure 5). Our study shows that BRAF V600E mutation is associated with RSPO4 at transcript and protein level.

RSPO4 Overexpression Occurs Selectively in BRAF V600E PTC
Since both RSPO4 and BRAF gene was highly expressed in our metastatic cohort ( Figure 5A), we used Integrative Genome Viewer (IGV-https://software.broadinstitute. org/software/igv/, accessed on 1 January, 2022) to screened for common mutations described in thyroid cancer to identify those impacted on our target genes identified to discriminate between non-metastatic and metastatic PTC cases. Among our 20 samples, 10 harbored an A to T substitution at position 140,453,136 on chromosome 7, that is harboring a BRAF V600E classical mutation ( Figure 5B). Interestingly, BRAF mutated samples showed overexpression of RSPO4 (FPKM > 2 for BRAF V600E while FPKM < 1 for BRAF wt). Overall, the relative fold-change between BRAF V600E and BRAF wt was about 10 and statistically highly significant (p = 6.27 × 10 −7 ). This result was confirmed by qPCR and IHC in 170 PTC samples ( Figure 5C). In addition, the independent TCGA platform comparing BRAF V600E to BRAF wt papillary PTC found similar gene expression patterns ( Figure 5). Our study shows that BRAF V600E mutation is associated with RSPO4 at transcript and protein level.

BRAF Inactivation and RSPO4 Overexpression Influence PTC Cell Proliferation, Motility and Invasion
To evaluate the impact of BRAF-RSPO4 in metastatic and non-metastatic conditions, both genes (BRAF and RSPO4) were deleted in different regions of the gene by using CRISPR technology ( Figure 6A,B) as we previously described [35]. RSPO4 protein and mRNA expression was decreased after BRAF knockout when compared to wild-type (Figure 6A-D) in our panel of metastatic and non-metastatic human thyroid cancer cell lines (THJ-16T, 8505C and TPC-1). The CRISPR-RSPO4 and CRISPR-BRAF knockout carried out discrete changes in cellular morphology, however the knockout condition significantly impacted on the phenotype of metastatic and non-metastatic thyroid cancer cell lines by significantly decrease the proliferation, invasion, and migration status ( Figure 6E-H).

BRAF Inactivation and RSPO4 Overexpression Influence PTC Cell Proliferation, Motility and Invasion
To evaluate the impact of BRAF-RSPO4 in metastatic and non-metastatic conditions, both genes (BRAF and RSPO4) were deleted in different regions of the gene by using CRISPR technology ( Figure 6A,B) as we previously described [35]. RSPO4 protein and mRNA expression was decreased after BRAF knockout when compared to wild-type ( Figure 6A-D) in our panel of metastatic and non-metastatic human thyroid cancer cell lines (THJ-16T, 8505C and TPC-1). The CRISPR-RSPO4 and CRISPR-BRAF knockout carried out discrete changes in cellular morphology, however the knockout condition significantly impacted on the phenotype of metastatic and non-metastatic thyroid cancer cell lines by significantly decrease the proliferation, invasion, and migration status ( Figure 6E-H).

RSPO4 Inhibition in PTC Cells Prevents the Activation of Proteins Involved in Focal Adhesion Signaling
Metastatic and non-metastatic thyroid cancer cell lines carrying CRISPR-RSPO4 knockout, and their controls were used to understand how the gene modulation could influence the phenotypic changes related to cell movement and invasion. As expected, inhibition of RSPO4 induced a downregulation of Wnt/catenins and an upregulation of E-cadherins. Interestingly, RSPO4 downregulation resulted in a potent inhibition of activated forms of mutiple phosphorylated proteins involved in focal adhesion signaling ( Figure 7A-C). These include Y397 (p-PTK2), the autophosphorylation site of focal adhesion kinase (PTK2/FAK) critical for PTK2 activation and interaction with downstream partners such as 118Y-p-Paxillin, an adapter proteins and partner of PTK2 in focal adhesion signaling and regulation of focal adhesion turnover during cancer cell locomotion; and 416Y-p-Src, which regulate Src multiple downstream signaling pathways, including Src interaction with PTK2 during in cell migration and invasion.

RSPO4 Inhibition in PTC Cells Prevents the Activation of Proteins Involved in Focal Adhesion Signaling
Metastatic and non-metastatic thyroid cancer cell lines carrying CRISPR-RSPO4 knockout, and their controls were used to understand how the gene modulation could influence the phenotypic changes related to cell movement and invasion. As expected, inhibition of RSPO4 induced a downregulation of Wnt/catenins and an upregulation of Ecadherins. Interestingly, RSPO4 downregulation resulted in a potent inhibition of activated forms of mutiple phosphorylated proteins involved in focal adhesion signaling (Figure 7A-C). These include Y397 (p-PTK2), the autophosphorylation site of focal adhesion kinase (PTK2/FAK) critical for PTK2 activation and interaction with downstream partners such as 118Y-p-Paxillin, an adapter proteins and partner of PTK2 in focal adhesion signaling and regulation of focal adhesion turnover during cancer cell locomotion; and 416Y-p-Src, which regulate Src multiple downstream signaling pathways, including Src interaction with PTK2 during in cell migration and invasion.

Discussion
Genomic testing of thyroid cancer nodules is becoming an important component in decision-making in clinical care, reducing the need for diagnostic surgery and improving accuracy of cancer risk assessment. Despite advances in molecular testing and the discovery of new promising therapeutics, effective treatments for advanced radioactive iodine (RAI)refractory metastatic differentiated thyroid cancer (DTC) are still lacking. This study focused on the identification of novel biomarkers for aggressive PTC using a cohort of PTC cases with over than 10 years fellow-up and characterized by an unusual invasiveness, including extensive lymph node metastases and higher incidence of locoregional recurrence. 2034 differentially expressed genes (DEGs) were identified by high-throughput RNA sequencing (RNA-Seq) to be associated with PTC metastatic profile. The pathway enrichment analysis was used to gain mechanistic insight into the gene panel generated from the genomescale (omics) experiment. Considering the top 50 genes specifically deregulated in BRAF mutated PTC, we identified markers involved in multiple signaling pathways, including Wnt, focal adhesion, and cell-matrix interaction. Deregulation of members of these pathways have been identified to contribute to cancer metastasis in other tumor types. [36][37][38] and some are tightly regulated by the tumor surrounding microenvironment [39][40][41].
Our study reveals that PTC with mutated BRAF-(V600E) (but not BRAF-wild-type) presented a consistent overexpression of RSPO4, a member of the R-Spondin family of secretory proteins involved in the regulation of Wnt signaling. This family is comprised of RSPO1, RSPO2, RSPO3, and RSPO4; these members share overall 60% sequence homology and similar domain organization [34]. Despite the high structural homology between these members, the expression patterns and phenotypes in knockout mice have demonstrated striking differences [42]. In our cohort, from the four family members, only RSPO4 was significantly upregulated. This is the first report identifying RSPO4 gene and protein overexpression in metastatic PTC. Earlier studies demonstrated that the stimulation of RSPO1-4 greatly potentiates the activity of Wnt ligands in both the canonical (B-catenin-dependent) and noncanonical (B-catenin-independent) pathways with major effects on the actin cytoskeleton, leading to loss of focal adhesions [43][44][45][46]. Furthermore, our study identified RSPO4 overexpression to be associated with BRAF V600E mutation, supporting a cooperation of these events to promote influence cell-matrix interaction, and cell motility and invasion. Interestingly, downregulation of RSPO4 in cell lines induced a strong inhibition of activated forms of the focal adhesion-associated proteins including Fak, Paxillin, and Src which are major signaling molecules involved in the regulation cancer cell stiffness and cell locomotion/migration [47][48][49][50]. However, future studies are necessary to establish molecular mechanisms by which RSPO4 regulate these multiple intracellular targets and to exploit their potential therapeutic utility to develop inhibitors tailored for patients with advanced PTC that can also be exploited in combination with standard radiotherapy/chemotherapy to overcome resistance/recurrence.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cells12010139/s1, Table S1: Genes differentially expressed between metastatic (ME) and non-metastatic (NM) PTC patients.  Informed Consent Statement: Written informed consent has been obtained from the patients.

Data Availability Statement:
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.