Propensity Score-Matched Analysis to Identify Pathways Associated with Loss of Sodium Iodide Symporter in Papillary Thyroid Cancer

Sodium iodide symporter (NIS) expression in thyroid follicular cells plays an important role in normal physiology and radioactive iodine therapy for thyroid cancer. Loss of NIS expression is often seen in thyroid cancers and may lead to radioiodine refractoriness. To explore novel mechanisms of NIS repression beyond oncogenic drivers, clinical and RNA-seq data from the thyroid cancer dataset of The Cancer Genome Atlas were analyzed. Propensity score matching was used to control for various genetic background factors. We found that tumoral NIS expression was negatively correlated with tumor size. Additionally, low NIS expression was the only factor associated with recurrence-free survival in a Cox multivariate regression analysis. After matching for clinicopathologic profiles and driver mutations, the principal component analysis revealed distinct gene expressions between the high and low NIS groups. Gene set enrichment analysis suggested the downregulation of hedgehog signaling, immune networks, and cell adhesions. Positively enriched pathways included DNA replication, nucleotide excision repair, MYC, and Wnt/β-catenin pathways. In summary, we identified several potential targets which could be exploited to rescue the loss of NIS expression and develop redifferentiation strategies to facilitate radioactive iodine therapy for thyroid cancer.


Introduction
Sodium iodide symporter (NIS), encoded by the SLC5A5 gene on chromosome 19p13.11, is a transmembrane glycoprotein that mediates active transport of iodide from the bloodstream into thyroid follicular cells [1]. NIS-mediated radioactive iodine (RAI) uptake is the cornerstone of the diagnosis and treatment of thyroid cancer with RAI. For patients with differentiated thyroid cancer, RAI therapy is associated with a lower risk of all-cause and cancer-specific mortality [2]. When thyroid cancer becomes RAI-refractory, the 10-year survival rate of patients with thyroid cancer drops to less than 20%, and the life expectancy decreases to 3 to 5 years [3].
Impaired RAI uptake or accumulation accounts for the refractoriness of RAI therapy. This parallels the dedifferentiation process during tumor progression, with a decrease in, or loss of, NIS expression and/or targeting to the plasma membrane [4]. Several oncogenic pathways are associated with the loss of NIS expression [5]. The most well-known molecular link to NIS loss in thyroid cancer is the BRAF V600E mutation-the most common genetic alteration in papillary thyroid cancer [6]. A selective BRAF inhibitor, dabrafenib, has been shown to stimulate RAI uptake in clinical trials [7]. However, novel therapeutic targets are still urgently needed to tailor NIS redifferentiation and further improve patient outcomes.
In the present study, we used clinical and molecular data from the dataset of The Cancer Genome Atlas for Thyroid Cancer (TCGA-THCA) to explore potential pathways associated with the loss of NIS expression. To eliminate the confounding effects of different genetic backgrounds, a propensity score-matched approach was adopted to provide a balanced and innovative comparison.

Data Acquisition
Clinicopathologic information, Mutation Annotation Format files, and normalized level 3 RNA-seq data were downloaded from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/; accessed on 11 June 2021). The American Joint Committee on Cancer (AJCC) stage, American Thyroid Association (ATA) recurrence risk group, recurrence events, and patients' vital status were obtained and updated. Gene expression levels were acquired at RNA-seq by expectation maximization (RSEM) values [8]. The thyroid differentiation score (TDS) and BRAF-RAS score (BRS) were referenced from the initial TCGA-THCA global analysis [9]. A lower or negative BRS represents more mitogenactivated protein kinase (MAPK) output by oncogenic BRAF V600E. Additionally, the BRAF-RAS class was based on the gene expression profile, which could be classified as BRAF-like (resembling those with the BRAF V600E mutation) or RAS-like (resembling those with RAS mutations).
A formal power analysis was not conducted. Clinicopathologic parameters were compared between groups using the chi-square test, Fisher's exact test, or the Cochran-Armitage trend test for categorical variables and the Mann-Whitney U test for continuous variables. Correlations between parameters were investigated using Spearman's rank correlation coefficients. Survival curves were generated using the Kaplan-Meier method and compared using the log-rank statistic [10]. Variables with p < 0.1 in univariate analysis were subjected to multivariate Cox regression modeling.

Propensity Score Matching
Propensity score matching was performed to adjust for differences in baseline covariates and driver mutations that might have a major impact on gene expression profiling. A 1:1 matching scheme was utilized with the nearest-neighbor matching method without a replacement. The propensity score was assessed using logistic regression, which included clinicopathologic profiles (age, sex, female, tumor subtype, tumor size, multifocality, extrathyroidal extension, lymph node metastasis, and stage) and driver mutations (BRAF, RAS, and TERT promoter mutations) [11]. Balance in covariate distribution between groups after the matching process was examined by post-matching analysis.

Principal Component Analysis
Principal component analysis was performed to visually explore the distribution pattern of different groups according to their gene expression profiles [12]. The percent explained variance was examined using eigenvalue decomposition. The first two components providing an adequate representation of the data were used to produce the biplot.

Pathway Analysis
To prevent the problem of arbitrary cutoff thresholds in over-representation analysis, gene set enrichment analysis (GSEA) was performed [13]. The GSEA strategy requires considering entire genes in the enrichment analysis without selecting significant genes. The enrichment score was calculated from the rank order of all genes in the annotation category. The Molecular Signatures Database v7.5.1 was retrieved from the GSEA website (https://www.gsea-msigdb.org/; accessed on 7 February 2022). The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Pathway Interaction Database (PID) databases were examined to identify potential pathways associated with NIS loss.

Results
A total of 500 patients with papillary thyroid cancer were included in the analysis before matching. As expected, the NIS expression level was closely correlated with TDS (Spearman's rho = 0.542, p < 0.001). The age at diagnosis was not associated with the NIS expression level (p = 0.728), but there was a negative correlation between the tumor size and NIS expression (Spearman's rho = −0.153, p = 0.001, Figure 1a). No difference in overall survival was observed between the groups (log-rank 0.172). However, patients in the low NIS group had a significantly shorter recurrence survival time (p < 0.001). Extrathyroidal extension, lymph node metastasis, AJCC st TERT promoter mutation, and NIS expression status were associated with recurrence survival in univariate analysis ( Table 2). In multivariate analysis, low NIS expression the only factor associated with recurrence-free survival, with a hazard ratio of 3.136 0.014). Following matching for tumor size, AJCC stage, and driver mutations, 121 pati in each group were selected for transcriptome analyses. The groups after matc showed good balance in clinicopathologic features (Supplementary Table S1). Princ component analysis of transcriptome profiling showed distinct gene expression betw Data are expressed as the number (percentage) or median (interquartile range). After a median split, patients were assigned to either the high or low NIS groups. Patients in the low NIS group had larger tumor size, more advanced extrathyroidal extension, and a higher risk of recurrence (Table 1). A higher prevalence of the BRAF V600E mutation was seen in the low NIS group, consistent with a positive correlation between BRS and NIS expression (Spearman's rho = 0.303, p < 0.001, Figure 1b). Of interest, TERT promoter mutations were more frequent in the low NIS group than in the high NIS group.
No difference in overall survival was observed between the groups (log-rank p = 0.172). However, patients in the low NIS group had a significantly shorter recurrence-free survival time (p < 0.001). Extrathyroidal extension, lymph node metastasis, AJCC stage, TERT promoter mutation, and NIS expression status were associated with recurrence-free survival in univariate analysis (Table 2). In multivariate analysis, low NIS expression was the only factor associated with recurrence-free survival, with a hazard ratio of 3.136 (p = 0.014).
Following matching for tumor size, AJCC stage, and driver mutations, 121 patients in each group were selected for transcriptome analyses. The groups after matching showed good balance in clinicopathologic features (Supplementary Table S1). Principal component analysis of transcriptome profiling showed distinct gene expression between the high and low NIS groups, although BRAF-like and RAS-like classes remained separable ( Figure 2).  Pathway analysis was performed by GSEA to compare gene expression between the two groups. KEGG pathways enriched in the low NIS group included downregulation of hedgehog signaling, immune networks, and cell adhesions, while nucleotide excision repair and DNA replication were upregulated (Table 3). Consistently, negatively enriched PID pathways were the ones primarily relevant to T-cell immunity (Table 4). In contrast, MYC, LKB1, and Wnt/β-catenin pathways were positively enriched in the low NIS group. To further corroborate these findings, we analyzed the association between the expression of NIS and key effectors in the entire TCGA-THCA patient cohort. As shown in Figure 3, NIS expression was positively correlated with GLI1 (Spearman's rho = 0.386, p < 0.001) and was negatively correlated with CTNNB1 (Spearman's rho = −0.138, p = 0.002).

Discussion
Refractoriness to RAI therapy represents an unmet therapeutic need in thyroid cancer. Early identification of patients who are refractory to RAI therapy not only saves patients from ineffective radiation exposure but also expedites the implementation of alternative treatment plans. An updated meta-analysis reported that extrathyroidal extension, BRAF V600E mutation, TERT promoter mutations, and high-risk histological subtypes were associated with RAI refractoriness [14]. In this study, we substantiated that the loss of NIS expression was associated with high-risk clinicopathologic features of papillary thyroid cancer. Consistent with a previous report [15], we found that NIS expression was inversely correlated with tumor size. This finding suggests that the NIS loss occurs during the process of tumor progression, and molecular pathways that are in the relatively late stages probably play some role in NIS downregulation. The recognition of molecular mechanisms that are involved in the loss of NIS expression and RAI refractoriness may shed some light on the puzzles in our understanding of thyroid dedifferentiation processes.
Age is an important prognostic factor in thyroid cancer, and thyroid cancer is the only cancer type that incorporates age in the tumor-node-metastasis stage. The thyroid tissue of children has smaller follicles and higher expression of NIS, pendrin, and dual oxidases [16]. Additionally, younger patients were more likely to have RAI-avid distant metastases [17]. It was therefore reasoned that physiological and pathological changes associated with aging can influence NIS expression and RAI responsiveness. However, we found that tumor NIS expression was unrelated to patient age. Our results are in agreement with a previous study from Portuguese researchers [18]. The negative impact of aging on prognosis likely results from metabolic alterations and immune dysregulation [19].
We previously highlighted that the low expression of thyroid differentiation genes was associated with disease recurrence [20]. In the present study, low NIS expression was identified as a predictor of recurrence-free survival, independent of oncogenic drivers. MAPK and PI3K/AKT pathways are well-characterized mechanisms leading to NIS repression [5]. BRAF mutations are particularly enriched in RAI-refractory metastases [21]. Without adequately matching, the results of comparisons would be substantially confounded by these classic pathways. The strength of our approach lies in balancing the genetic discrepancy between the high and low NIS groups.
In the low NIS group, GSEA indicated the downregulation of hedgehog signaling. Although experimental evidence suggests that the hedgehog pathway promotes the selfrenewal of thyroid cancer stem cells [22], hedgehog pathway activity is generally reduced in papillary thyroid cancer [23]. We recently reported that ethacridine, a TAZ activator, induced differentiation, despite TAZ being conventionally considered to preserve the stemness properties of progenitor cells [24]. The blockage of hedgehog signaling by cyclopamine could result in markedly decreased type 3 deiodinase expression [25]. Whether a cause-andeffect relationship exists between hedgehog signaling and NIS downregulation requires further investigation.
The Wnt/β-catenin pathway is also essential for stem cell maintenance. A previous study showed that β-catenin overexpression decreased cytomembrane localization of NIS and RAI uptake ability in FTC-133 cells [26]. However, in rat thyroid follicular cell lines, β-catenin as downstream of thyrotropin and insulin-like growth factor 1 increased the expression of Pax8 and NIS [27]. It should be noted that differential β-catenin expression between benign and malignant thyroid neoplasms was seen in cytoplasmic or nuclear localization (indicating accumulation of the Wnt canonical pathway) but not membranous localization (non-Wnt signaling) [28]. Furthermore, cytoplasmic/nuclear β-catenin expression was strongly associated with galectin-3 expression. Galectin-3 is a carbohydrate-binding protein commonly overexpressed in thyroid cancer and can serve as a diagnostic biomarker of thyroid cancer. Recently, we reported that galectin-3 inhibitors could suppress anoikis resistance and invasive capacity in thyroid cancer cells [29]. Activation of the Wnt/β-catenin pathway along with NIS downregulation is a fascinating topic for future work.
A non-pump function of intracellular non-membranous NIS may have a tumorigenic role [30]. In the current study, we observed upregulation of the DNA replication pathway and MYC pathway in the low NIS group. Network analysis revealed that MYC is one of the essential proteins that regulate NIS expression in anaplastic thyroid cancer [31]. Interestingly, a BRD4 inhibitor, JQ1, was able to suppress MYC and enhance NIS expression [32]. MYC-stimulated cell proliferation may potentiate genotoxic stress induced by oncogenic drivers. NIS repression can also result from DNA damage involving ATM-mediated mechanisms [33]. We noted that dermatopontin is among the top-ranked differentially expressed genes among both the high and low NIS groups. Dermatopontin impeded the proliferation of thyroid cancer cells through MYC repression [34]. It will be intriguing to determine whether MYC transcriptionally regulates NIS expression.
Taken together, through a propensity score-matched analysis, our transcriptome analysis identified several novel pathways that could serve as potential targets in future studies to reverse the loss of NIS expression in thyroid cancer. More research is warranted to determine the translatability of these findings for developing redifferentiation strategies for the management of patients with RAI-refractory thyroid cancer.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/cimb44040101/s1, Table S1: Clinicopathologic features of the TCGA-THCA patient cohort, grouped by the expression of sodium iodide symporter (NIS) after propensity score matching.