Identication of Oncosuppressing Effects of Seven Selenoproteins in Thyroid Cancer: Implications for Distinct Roles of Selenoproteins in Tumorigenesis

Background: Low selenium levels are associated with increased incidence and advanced stage of thyroid cancers. In response to changes in selenium levels, a hierarchy of selenoprotein biosynthesis allows tissue-specic ne-tuning of the 25 selenoproteins. To determine the role of individual selenoproteins on thyroid carcinogenesis, we carried out a multiomic data mining study. Methods: The expressions of individual selenoproteins and their correlation with prognosis in thyroid cancers were analyzed using Oncomine, GEPIA and Kaplan–Meier plotter platforms. Co-expression analysis using cBioportal database was carried out to identify genes that are correlated with selenoproteins. GO and KEGG enrichment were further performed for genes correlated with selenoproteins that are clinically signicant. Results: DIO1, GPX3, SELENOP, SELENOM SELENOS, SELENOO and SELENOV were signicantly downregulated and with a poor prognosis in thyroid cancers. 35 biological processes including GO: 0045926 (negative regulation of growth) and 23 biological processes including GO:0001525 (angiogenesis) and GO:0007155 (cell adhesion) were enriched in DIO1-positively and DIO1-negatively correlated genes, respectively. 41 biological processes including GO: 0045926 (negative regulation of growth) and 23 biological processes including GO:0007165 (signal transduction) and GO: 0000165 (MAPK cascade) were enriched in GPX3-positively and GPX3-negatively correlated genes, respectively. The antitumor effects of SELENOM and SELENOS might be attributed to their protection against endoplasmic reticulum (ER) stress. SELENOO was revealed to be correlated with ER stress, mitochondria translation and telomere maintenance. Biological processes of SELENOV-correlated genes were enriched in oxidation-reduction process and ER calcium ion homeostasis. Moreover, cell adhesion and angiogenesis were also shown to be negatively regulated by SELENOV, providing an anti-metastatic effect similar as DIO1. Conclusion: This study conrms the distinct roles of the 25 selenoproteins of

While involved in various biological reactions, a notable feature of selenoproteins is their participation in cellular oxidoreductase reactions. Since oxidative stress may cause mutagenesis and genomic variability in cells, selenoproteins have long been thought to protect against cancer progression by its antioxidant role [5] . Indeed, interfering with selenoprotein biosynthesis in mice by altering Sec tRNA coding gene, Trsp, could accelerate tumor development and progression in prostate cancer [6] , colon cancer [7] and breast cancer [8] . However, large clinical trials have generated controversial results concerning selenium supplementation for cancer prevention [9][10][11] . This inconsistent result could be explained by the organspeci c and context-speci c role of selenoproteins contributed by a hierarchical regulation of selenoprotein biosynthesis [12] . However, the potential mechanisms of individual selenoproteins in tissuespeci c cancer development and progression remains obscure.
To determine the role of individual selenoproteins on pathogenesis of thyroid cancer, we carried out a multiomic data mining study and identi ed DIO1, GPX3, SELENOP, SELENOM, SELENOS, SELENOO and SELENOV were signi cantly downregulated and with a poor prognosis in thyroid cancers. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment were further performed for genes that are correlated with these selenoproteins to explore the potential mechanisms of their oncosuppressing role in thyroid cancers. This study provides useful information for revealing the currently unknown functions of selenoproteins and their roles in thyroid carcinogenesis.

Oncomine Analysis
Data regarding mRNA expression of 25 selenoproteins in thyroid cancers was retrieved from Oncomine (https://www.oncomine.org/resource/login.html) [13] . Primary lters for "thyroid gland carcinoma" and "cancer vs. normal analysis" were used. 12 analyses with 1147 specimens including 574 normal thyroid tissues and 573 thyroid cancer tissues were available. The database was queried for each selenoproteins individually with the lter settings as p-value < 0.01, fold change > 2, and gene ranking in the top 10%.

GEPIA Dataset
GEPIA dataset (http://gepia.cancer-pku.cn/) provides access to a large collection of data from The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) projects [14] . In this study, differential analysis of selenoprotein mRNA expression with a cut-off p-value of 0.01 in 512 thyroid cancer tissues and 337 normal thyroid tissues was validated. The relationships between expression of each selenoproteins with pathological stages of thyroid cancers were analyzed and a cut-off p-value of 0.05 was used.

The Kaplan-Meier Plotter Analysis
To analyze the impact of selenoproteins on survival, 502 samples from patients with thyroid cancers who underwent an average follow-up period of 240 months were retrieved from the Kaplan-Meier plotter (http://kmplot.com). They were divided into two groups according to median levels of selenoprotein expression. Kaplan-Meier survival plots were constructed to compare the two patient groups and calculate the hazard ratio (HR) with 95% con dence intervals (CIs), and the log rank p-values.
cBioPortal Analysis Genes that are correlated with each selenoproteins were retrieved in a papillary thyroid cancer (PTC) dataset (496 samples, TCGA Cell 2014) from cBioPortal database (https://www.cbioportal.org/) using the co-expression tool [15] . These genes harbor a correlation with Spearman's higher than 0.4 or lower than − 0.4.
Correlation between different selenoproteins in thyroid cancers were also analyzed using the cocorrelation tool on cBioportal database.

Gene Ontology and Pathway Enrichment Analysis
Functions and pathways enrichment of genes that are correlated with each selenoproteins were analyzed using DAVID (https://david.ncifcrf.gov/) [16] and KEGG (https://genome.jp/kegg/) [17] online databases with p-value < 0.05 as the cut-off criterion.

Selenoprotein mRNA Expression Analysis in Thyroid Cancers
To explore the expression pattern of selenoproteins in thyroid cancers, a meta-analysis of transcriptomic pro ling data using Oncomine database was performed (Fig. 1). We observed DIO2, GPX3, SELENOP, SELENOM, SELENOS and SELENOO were underexpressed and GPX1 was overexpressed in thyroid cancer tissues. Since samples retrieved from distinct databases are different, we used another online platform, GEPIA, to con rm the results above. As shown in Fig. 2, a signi cantly low expression of DIO1, DIO3, GPX3, SELENOO, and SELENOV and a considerably high expression of GPX1 were observed in thyroid cancer tissues. The expressions of DIO2, SELENOP, SELENOM, and SELENOS were also reduced in thyroid cancer tissues, although no statistical signi cance was observed.

The Prognostic Analysis of Selenoproteins in Thyroid Cancers
To investigate the extent to which selenoprotein expression is associated with prognosis, we detected the expression of each selenoprotein individually in different stages of thyroid cancer using GEPIA platform. As shown in Fig. 3, a negative correlation with pathological stages was observed in DIO1, DIO2, DIO3, GPX3, TXNRD2, SELENOK, SELENOP, SELENOS, SELENOV, SEPHS2 and MSRB1, while a positive correlation was seen in GPX4.
Since the overall 5-year survival rate of thyroid cancer patients is generally well with 98% [18] , we thus compare the recurrence-free survival (RFS) of patients with high and low expression of each selenoprotein by Kaplan-Meier plotter analysis. As shown in Table 1 and Figure S1, low expression of DIO1, GPX2, GPX3, TXNRD2, SELENOF, SELENOK, SELENOM, SELENOO, SELENOP, SELENOS, SELENOV, SELENOW and, and high expression of SELENOI and SEPHS2 were associated with worse RFS. In order to identify selenoproteins that are clinically signi cant in thyroid cancers, we integrated the above selenoprotein expression data and prognostic data and found DIO1, GPX3, SELENOP, SELENOM SELENOS, SELENOO and SELENOV were downregulated and with a poor prognosis in thyroid cancers. These results indicated an anti-tumor role of these selenoproteins in thyroid cancers, and we thus performed further studies to predict their biological functions and related pathways.
Gene Ontology and Pathway Enrichment of Genes Correlated with DIO1, GPX3, SELENOP, SELENOM SELENOS, SELENOO and SELENOV in Thyroid Cancers To identify genes that are correlated with DIO1, GPX3, SELENOP, SELENOM SELENOS, SELENOO and SELENOV, we performed co-expression analysis (RNA Seq V2 RSEM) using cBioportal database. The correlated genes of DIO1, GPX3, SELENOP, SELENOM, SELENOS, SELENOO and SELENOV with a Spearman's correlation higher than 0.4 or lower than − 0.4 were listed in Table S1-S7, respectively. The top 200 genes on the correlation ranks were selected and GO analysis were performed using DAVID. As for DIO1, 35 biological processes including GO: 0045926 (negative regulation of growth) were enriched in the top 200 DIO1-positive-correlated genes and 23 biological processes including GO:0001525 (angiogenesis) and GO:0007155 (cell adhesion) were enriched in the top 200 DIO1-negatively-correlated genes by GO analysis (Table S8). The top 5 processes on each list were shown in Fig. 4A. In addition, KEGG analysis identi ed 7 pathways related to DIO1 in thyroid cancers, including hsa04978: Mineral absorption, hsa04918: Thyroid hormone synthesis, hsa04922: Glucagon signaling pathway, hsa04977: Vitamin digestion and absorption, hsa04020: Calcium signaling pathway, hsa00980: Metabolism of xenobiotics by cytochrome P450, and hsa05205: Proteoglycans in cancer (Table 2). With regards to SELENOP, GO and KEGG enrichment failed to identify any functions and pathways except for the biological process GO: 0018105 (peptidyl-serine phosphorylation), which is involved in selenoprotein biosynthesis (data not shown).

Correlations of Selenoproteins with Each Other
Next, we calculated the Spearman's rank correlations of selenoproteins with each other by analyzing their mRNA expression (RNA Seq V2 RSEM) using cBioPortal online database for thyroid cancers. The results indicated a signi cant and positive correlations (r ≥ 0.6) in the following selenoproteins: DIO1 with SELENOV; GPX1 with GPX4, SELENOW, and SELENOH; TXNRD2 with SELENOO; SELENOK with SELENOS (Fig. 4G).

Discussion
Low serum and tissue levels of selenium are associated with increased incidence and advanced stage of thyroid cancers [2] . In response to changes in selenium supply, a hierarchy of selenoprotein biosynthesis allows tissue-speci c ne-tuning of the 25 selenoproteins [12] . However, the role of individual selenoproteins on thyroid tumorigenesis and metastasis remains unclear. Here, we found an underexpression of DIO1, GPX3, SELENOP, SELENOM SELENOS, SELENOO and SELENOV in thyroid cancers and their downregulation was associated with a poor prognosis. The connection of a decreased expression to poor tumor prognosis strongly indicates a potential onco-suppressing role of these selenoproteins in thyroid cancers.
DIOs play an important role in the regulation of thyroid hormone metabolism and intracellular action. DIO2 convert thyroxine (T4) to the active hormone, triiodothyronine (T3) through outer ring deiodination; while DIO3 inactivate THs through inner ring deiodinase. DIO1 involves in both inner and outer ring deiodinations and is viewed as a scavenger enzyme in iodine recycling from the inactive iodothyronines [19] . Since thyroid hormones regulate cellular growth, development and differentiation, it is not surprising that expression of these deiodinases altered in various cancers, such as breast cancers, renal cancers, liver cancers, lung cancers, brain tumors, and thyroid cancers [20][21] . Mechanistic studies to determine how exactly these deiodinases modify tumor growth, however, are still lacking. In the current study, a tumor suppressing role of DIO1 was established and GO analysis showed DIO1-negative-correlated genes were mainly enriched in biological processes including angiogenesis and cell adhesion. The antimetastatic effect of DIO1 was in fact comparable to that observed in renal cancer cells, in which ectopic expression of DIO1 inhibit proliferation and migration through regulating the expression of integrins and their ligands, including collagens, bronectin, laminin, and vitronectin [22] . In addition, biological process GO: 0045926 (negative regulation of growth), which involves 5 members of metallothionein (MT) family, was identi ed to be positively related with DIO1. MTs are a family of low molecular and cysteine-rich proteins regulating heavy metal homeostasis and oxidation-reduction reactions [23] . In thyroid cells, defects in MT1G and MT1M expression has been linked to malignant transformation [24][25][26] . Yet, the speci c mechanisms how DIO1 and MTs are related with each other remain unclear. Zinc plays a key role in thyroid hormone metabolism specially by regulating deiodinases enzyme activities [27] and exerts anticancer role via many cellular processes such as antioxidant reactions, transcriptional regulations, DNA and RNA repair [28] . Since MTs are able to bind both Zinc and Selenium, it is possible that defects in DIO1 and MTs expression are attributed to changes of these physiology heavy metals in the context of thyroid cancers.
GPX1-4 and GPX6 are selenium-containing GPX family members that catalyze hydroperoxide reduction using GSH as reductant. Among the GPX family members, GPX3 is the only one that can be secreted into the plasma, and thus plays a critical role in reduction of extracellular ROS stress [29] . Downregulation of GPX3 by hypermethylation has been well established in thyroid carcinogenesis, while expression of other GPXs in thyroid cancer remain controversial [30] . Consistently, we found a decreased GPX3 expression connected with a poor prognosis of thyroid cancer, whereas the connection was not observed in other GPX family members. In thyroid, GPX3 is detected in the follicular lumen to protect the thyroid follicular cells from possible damage by excess H2O2, an oxidative agent required for thyroid hormone synthesis [31] . GO analysis revealed signal transduction was enriched in GPX3-negative-correlated genes, which may be explained by the accumulation of hydroperoxides, particularly H2O2, by GPX3 downregulation during transformation of thyroid cancer cells [32] . In addition, MAPK signaling, a major and highly prevalent oncogenic pathway in thyroid cancers, was identi ed to be negatively correlated with GPX3 [33] . GPX3 has been reported to inhibit cell growth and metastasis in prostate cancer through downregulation of c-Met, which is a tyrosine kinase receptor initiating MAPK activation [34][35] . Thus, suppression of MAPK cascade by GPX3 might be mediated by c-Met, in which further studies are needed to clarify the precise mechanism. Moreover, the positive correlation of MTs and GPX3 might be explained by their responses to changes in selenium levels as discussed above.
SELENOM and SELENOS are ER-resident selenoproteins regulating protein folding and degradation, consistent with the results of our GO analysis [36] . In cancer cells, enhanced protein synthesis triggers ER stress and the sustained ER stress promotes tumorigenesis and metastasis [37] . Accordingly, the antitumor effects of SELENOM and SELENOS may be attributed to their protection against ER stress in thyroid cells. In addition, GO analysis revealed that SELENOM and SELENOS regulate mitochondrial translation responsible for the maintenance of energetic balance. Currently, there is no evidence implicating their expression in mitochondria. Since a direct interaction of ER and mitochondria has been reported [38] , it is possible that SELENOM and SELENOS are somehow involved in the interconnectedness of the ER and mitochondria.
To date, functions of SELENOO and SELENOV are largely unknown. Our data showed SELENOO is involved in ER stress, mitochondria translation and telomere maintenance, similar as SELENOM and SELENOS. The tumor suppressing role of SELENOO, thus, probably connects with its potential role in regulation of ER functions. SELENOV was recently reported to protect against reactive oxygen and nitrogen species-mediated ER stress and oxidative injury [39] , which is likely one of the mechanisms for its anti-cancer role. Consistently, our data showed SELENOV is related with biological processes including oxidation-reduction process and ER calcium ion homeostasis. Moreover, cell adhesion and angiogenesis were also shown to be negatively regulated by SELENOV, providing an anti-metastatic effect similar as DIO1. In accordance with our correlation analysis, functions of SELENOV and DIO1 might be signi cantly correlated with each other. However, the precise mechanisms and implications in cancer progression still remain to be clari ed.
Whereas well known for its role in selenium transport, SELENOP has antioxidant capacity exerting its "double-edged sword" effects in the process of cell transformation [40] . After produced in the liver, SELENOP is secreted into plasma and travel to targeted organs participating in synthesis of other selenoproteins [41] . Here, we showed SELENOP was downregulated in thyroid cancer tissues and GO analysis failed to identify any other biological processes except for selenoprotein biosynthesis. Thus, a low expression of SELENOP in thyroid cancer tissues mainly indicates a reduced local supply of selenium. However, selenium supplementation for thyroid cancer prevention might not be effective and even could be harmful, since expressions of selenoproteins are not all reduced and some even increased in the context of thyroid cancer. Thus, targeted regulation of individual selenoproteins, such as DIO1, GPX3, SELENOP, SELENOM SELENOS, SELENOO and SELENOV, may serve as a potential therapeutic option for thyroid cancer.

Competing of Interests
The authors declare that there is no con ict of interest regarding the publication of this paper. Author's contributions ZY and CP analyzed and interpreted the data regarding the expression of selenoproteins. DXY did prognosis analysis of selenoproteins in thyroid cancer. ZY retrieved the genes that are correlated with the selenoproteins. LHJ and WY performed GO and KEGG analysis for selenoprotein-related genes. FJ and LS was a major contributor in writing the manuscript. All authors read and approved the nal manuscript. Meta-analysis of gene expression pro ling for the 25 selenoproteins in thyroid cancers using Oncomine database. With a cut-off p-value <0.01 and fold change > 2, underexpressed selenoproteins were shown in blue and overexpressed selenoproteins were shown in red.

Figure 2
Selenoprotein mRNA expression analysis in thyroid cancers using GEPIA platform. With a cut-off p-value of 0.01, underexpressed selenoproteins were shown in blue and overexpressed selenoproteins were shown in red (T: thyroid cancer tissues; N: non-malignant thyroid tissues; * p ≤0.01).

Figure 3
Prognostic analysis of selenoproteins in thyroid cancers using GEPIA platform. With a cut-off p-value of 0.05, selenoproteins negatively correlated with pathological stages were shown in blue, selenoproteins positively correlated with pathological stages were shown in red.