Matched analysis of circulating selenium with the breast cancer selenotranscriptome: a multicentre prospective study

Introduction Low serum selenium and altered tumour RNA expression of certain selenoproteins are associated with a poor breast cancer prognosis. Selenoprotein expression stringently depends on selenium availability, hence circulating selenium may interact with tumour selenoprotein expression. However, there is no matched analysis to date. Methods This study included 1453 patients with newly diagnosed breast cancer from the multicentric prospective Sweden Cancerome Analysis Network – Breast study. Total serum selenium, selenoprotein P and glutathione peroxidase 3 were analysed at time of diagnosis. Bulk RNA-sequencing was conducted in matched tumour tissues. Fully adjusted Cox regression models with an interaction term were employed to detect dose-dependent interactions of circulating selenium with the associations of tumour selenoprotein mRNA expression and mortality. Results 237 deaths were recorded within ~ 9 years follow-up. All three serum selenium biomarkers correlated positively (p < 0.001). All selenoproteins except for GPX6 were expressed in tumour tissues. Single cell RNA-sequencing revealed a heterogeneous expression pattern in the tumour microenvironment. Circulating selenium correlated positively with tumour SELENOW and SELENON expression (p < 0.001). In fully adjusted models, the associations of DIO1, DIO3 and SELENOM with mortality were dose-dependently modified by serum selenium (p < 0.001, p = 0.020, p = 0.038, respectively). With increasing selenium, DIO1 and SELENOM associated with lower, whereas DIO3 expression associated with higher mortality. Association of DIO1 with lower mortality was only apparent in patients with high selenium [above median (70.36 µg/L)], and the HR (95%CI) for one-unit increase in log(FPKM + 1) was 0.70 (0.50–0.98). Conclusions This first unbiased analysis of serum selenium with the breast cancer selenotranscriptome identified an effect-modification of selenium on the associations of DIO1, SELENOM, and DIO3 with prognosis. Selenium substitution in patients with DIO1-expressing tumours merits consideration to improve survival. Supplementary Information The online version contains supplementary material available at 10.1186/s12967-023-04502-y.


Introduction
Breast cancer remains a significant global health challenge, with an estimated 2.3 million new cases and 685,000 deaths annually [1].Discovery of novel prognostic factors may improve prognosis by identifying high risk women early and by personalizing or intensifying therapy regimens [2].
Recently, there has been growing scientific interest in the role of the essential trace element selenium (Se) in risk, progression and prognosis of breast cancer [3].Several large epidemiological studies have reported an independent association of low dietary intake or marginal serum levels of Se with a distinctly poor prognosis [4][5][6][7][8].Se affects various physiological pathways by being translationally incorporated into the products of 25 human selenoprotein genes, that mainly act in antioxidative defense systems, quality and structure control, and by activating or inactivating thyroid hormones [9][10][11].The biosynthesis of selenoproteins depends on their RNA expression level and on Se availability, and hereby correlates with Se intake and serum Se concentrations [12].Accordingly, marginal Se intake correlates to low serum concentrations of the Se transporter selenoprotein P (SELENOP), resulting in suboptimal systemic Se supply [13].This deficit is e.g.translated in kidney into relatively low biosynthesis and secretion, and hence lower activity of the plasma glutathione peroxidase 3 (GPx3).Accordingly, low serum Se and SELENOP concentrations as well as low serum GPx3 activity have been associated with higher risk of mortality and recurrence after breast cancer diagnosis [4].At the same time, tumour tissue gene expression levels of several selenoproteins, such as iodothyronine deiodinases (DIO), glutathione peroxidases (GPx) or thioredoxin-reductases (TXNRD) have also been reported as prognostic factors for breast cancer in large-scale genomic profiling studies [14][15][16].
Despite the evidence that both circulating Se and tumour gene expression of selenoproteins are associated with breast cancer prognosis, there is a lack of data from a matched analysis of serum Se biomarkers and tumour selenoprotein expression.Se availability constitutes a key factor for RNA stability of certain selenoprotein transcripts, and for translation, i.e., Se controls the rate of biosynthesis of the different selenoproteins from a given RNA expression level [17,18].Therefore, the association between selenoprotein transcript expression and breast cancer prognosis may be modified by Se availability, and a sufficiently high Se status may be required for translating differences in RNA expression levels into the corresponding gene products and disease-modifying selenoprotein activities.
To test this hypothesis, we simultaneously analysed three complementary biomarkers of Se status (total serum Se, SELENOP, GPx3) and conducted bulk RNAsequencing of tumours of 1453 patients with a new diagnosis of primary invasive breast cancer.The patients were followed for ~ 9 years.The main aim was to assess whether the association between tumour RNA expression of selenoprotein genes with breast cancer prognosis is modified by circulating Se and selenoprotein levels.

SCAN-B study design
The Sweden Cancerome Analysis Network Breast Study (SCAN-B) is a multicentre population based prospective study that intends to identify novel biomarkers for early identification of patients with a poor prognosis based on serum and tumour tissue (genomics) as a matrix.In brief, the study, which is still ongoing since 2010 enrols patients with a newly diagnosed or suspected primary breast cancer in the catchment area of Southern Sweden at multiple cancer centres.The study was approved by the Regional Ethical Review Board in Lund, Sweden (Registration numbers 2009/658, 2010/383, 2012/58, 2013/459, and 2015/277) and registered under the ClinicalTrials.govID NCT02306096.All enrolled participants have given written informed consent for inclusion in the study and analyses/procedures, which were conducted as an integrative part of clinical routine, as described before [19].

Follow up and covariate assessment
All clinical data was collected by standardized procedures by clinicians and referred to the Swedish National Quality Registry for Breast Cancer (NKBC).Vital status was obtained from the Swedish National Population Registry, which maintains records for all Swedish citizens.Covariates included patients' characteristics, tumour characteristics, and information on treatment procedures, as described before [5,19].These variables comprised age, sex, and information on menopausal status of the patients was reported.Histopathological information comprised tumour size, the side of involved breast, histopathological subtype, histological grade (Nottingham Histological Grade, NHG), Ki67 expression, HER2 expression, ER expression, PGR expression, and number of involved lymph nodes.Information on therapy comprised the surgical procedure regarding the breast, regarding the axillary lymph nodes, application of endocrine or immune therapy, radiation or chemotherapy.Mode of diagnosis, i.e. either clinical or by screening was also reported.The almost full completeness (99.9%) and over 90% validity of NKBC with regard to information on vital status, and other have been externally validated before [20].

Serum selenium biomarkers
Laboratory analyses and results for quantification of Se biomarkers in this study have been described before [4,5].In brief, total Se was measured using total reflection X-ray fluorescence (TXRF) spectroscopy, SELENOP was measured using a validated commercial ELISA (sele-nOtest, selenOmed GmbH, Berlin, Germany), and activity of glutathione peroxidase 3 was measured using an established coupled enzyme reaction.All analyses were conducted by scientists and technicians blinded to any clinical information, with samples arranged in a randomized order with regard to order of enrolment in the study, as reported before [4,5].

Selenoprotein gene expression in patient's tumours
The detailed protocols for RNA-sequencing have been described before, i.e. either the protocol as shown before [19] or using the Illumina stranded TruSeq mRNA procedure.Protocols were established using the Illumina Neo-Prep system or the KingFisher system.For the purpose of this study, gene expression values in Fragments Per Kilobase per Million reads (FPKM) was generated.An established analysis pipeline was used to extract FPKM values by alignment and estimation of gene expression data.The pipeline has been described in detail before [21], and involves the tools picard tools, trimmomatic, bowtie, hisat2, stringtie, dbSNP56 and GENCODE.Genes were annotated based on gene and transcript definitions contained in GENCODE Release 27.After adding an offset of 1, FPKM data were log-transformed for the analyses in this study.Single cell RNA sequencing data was accessed from the Single Cell Portal of the Broad Institute (https:// singl ecell.broad insti tute.org/ single_ cell), and included 100,064 single cells from 26 primary breast cancer tumours [22].Access and visualization was conducted on March 21st 2023.

Statistical analyses
Descriptive statistics were presented as median along with interquartile range (IQR) for continuous variables, and as frequencies along with percent for categorical data.Correlation matrices were generated to depict the relationship between RNA expression of selenoproteins and serum Se biomarkers.Non-parametric Spearman's rank correlation test was applied to compute Spearman's R and p values for the correlations, and cut-off for p-values were adjusted in correlation matrices taking into account the number of genes tested (0.05:23).
Linear multivariable Cox proportional hazards models were employed to calculate hazard ratios along with 95% confidence intervals (CI) for one increment in log(FPKM + 1) in gene expression for each gene.Models were adjusted for established clinical predictors of breast cancer prognosis, i.e. age at diagnosis (years), menopausal status (pre-menopausal, post-menopausal, uncertain), tumour size (mm), NHG (I, II, III), lymph node involvement (at least 4, 1 to 3, submicrometastasis (< 0.2 mm) or no involvement), HER2 expression (positive, negative), ER expression (positive, negative), PGR expression (positive, negative), histological type (ductal, lobular, ductal + lobular/other, other), laterality (right or left breast).Regression analyses were conducted in the entire cohort, and separately in the low and high Se group, based on each different biomarker, while median level of the cohort served as unbiased cut-off.In order to detect potential effect modification by Se biomarkers, an interaction term with the continuous Se biomarker variable was added.An interaction was considered statistically significant in case of p interaction < 0.05.Significant interactions between continuous variables were visualized using contour plots, and by visualizing the association based on tertiles of Se biomarker.Missing variables made up only a small portion (0.4%) of all values contained in variables in the adjusted models (Additional file 1: Fig. S1), therefore Cox regression models were computed using complete cases.
All analyses were conducted using the R language (The R Foundation for Statistical Computing, version 4.3.0) on the RStudio environment (RStudio, PBC, version 2022.2.3.492).

Results
Based on availability of tissue and RNA-sequencing data, as well as serum sampling and Se status assessment, a total of 1453 patients with complete RNA-sequencing and data on serum Se biomarkers were included in the final analyses.A detailed description of the study flow chart is included in Additional file 1: Fig. S2.Follow-up time comprised 9,701 years in total, corresponding to a mean follow-up time of 6.68 years, and 237 deaths were recorded in this time frame.

Baseline patient and tumour characteristics
Baseline patient characteristics and tumour characteristics as well as applied therapy regimens according to vital status during the study are presented in Table 1.Patients were divided based on whether they died during the follow-up period or survived.Patients that died over the course of the follow-up were older at time of diagnosis, more frequently post-menopausal, had larger tumours, more lymph nodes involved, lower serum Se and SELE-NOP concentrations and a lower serum GPx3 activity.The association of serum Se biomarkers with prognosis have been assessed in this cohort previously, displaying dose-dependent associations of low serum Se with a poor prognosis, independent of various confounders [5].Table 2 depicts the mode of clinical diagnosis, and the treatment regimens applied in the study cohort, comprising both adjuvant and surgical methods.

Serum selenium and tumour selenoprotein expression
Figure 1 A depicts the study design.Selenoprotein mRNA expression levels within the tumours are displayed in Fig. 1B.As GPX6 was not expressed in the tumour and SELENOV displayed only a low expression, both genes were excluded from further analyses.Analyses included the deiodinase family involved in thyroid hormone regulation (DIO1-3), the glutathione peroxidases involved in antioxidative defence (GPX1-4), thioredoxin reductases involved in cellular redox regulation (TXNRD1-3), selenoproteins located within the endoplasmic reticulum (e.g.SELENOM, SELENOF etc.), as well as the other selenoproteins with specific functions [9]. Figure 1 C displays the correlation between tumour selenoprotein expression and serum Se and selenoprotein levels.Within the group of selenoprotein genes, the highest correlation was between MSRB1 and SEPHS2 (R = 0.57), while SELENOI and DIO3 displayed the most prominent negative correlation (R = -0.34).All three serum biomarkers correlated positively, with serum Se and SELENOP displaying the   highest correlation (R = 0.59).Serum Se biomarkers were mostly not correlated with selenoprotein gene expression levels, except for a weak correlation between serum Se and SELENOP with SELENOW (R = 0.18 and R = 0.13, respectively) (Fig. 1D), and between serum Se and SELE-NON (R = 0.082).Figure 1E displays the selenoprotein gene expression in different cells in breast cancer.

Selenoprotein gene expression and survival based on selenium biomarkers
Figure 2 displays the association of selenoprotein mRNA expression of each gene with survival, in the whole cohort and in patients with low or high serum Se levels.There were significant interactions between serum Se with DIO1, DIO3, and SELENOM, p < 0.001, p = 0.020, and p = 0.038, respectively.Association of DIO1 with lower mortality was only apparent in patients with high Se [above median (70.36 µg/L)], HR (95%CI) for one-unit increase in log(FPKM + 1) was 0.70 (0.50-0.98).
The complex interaction between serum Se and DIO1, DIO3 and SELENOM is depicted in Fig. 3. Figure 3 A displays lower hazard ratios for a simultaneous increase in serum Se and DIO1 expression.This relationship is emphasized in Fig. 3B, which shows a decreased hazard ratio with increasing DIO1 levels, however only in patients with a relatively high Se level, i.e. residing in the 2nd or 3rd tertile.On the other hand, Fig. 3C displays an inverse interaction of Se with DIO3 levels, where increasing serum Se and simultaneous increase in DIO3 are associated with an elevated hazard ratio.Accordingly, DIO3 is associated with higher mortality in patients with high Se only, i.e., solely in the 3rd tertile (Fig. 3D).The interaction of serum Se with SELENOM was similar to DIO1, where SELENOM associated with a lower mortality rate with increasing Se levels (Fig. 3E and F).
In further sensitivity analyses, interactions of serum Se with DIO1, DIO3 and SELENOM were tested after further adjusting for treatment methods used.Endocrine therapy, immune-, chemo-and radiotherapy as well as surgical procedures regarding the breast and axilla were added to the fully adjusted models one by one, and all combined, and nearly no changes were observed in p values for interaction (Additional file 1: Table S1).
An analysis of The Cancer Genome Atlas Program (TCGA) data displayed no overall associations of DIO1, DIO3, SELENOM, SELENOW, and SELENON with survival, when not incorporating serum Se (Additional file 1: Fig. S5), highlighting the need for consideration of both serum Se and tumour selenoprotein expression in order to ascertain effects of selenoprotein mRNA expression on prognosis.

Discussion
In this large multicentric prospective study, the first matched analysis of circulating serum Se levels with the breast cancer selenotranscriptome was performed.As expected, serum Se levels were mostly unrelated to selenoprotein mRNA expression levels.Serum Se dosedependently interacted with the association of DIO1, DIO3, and SELENOM and survival.With increasing serum Se, DIO1 and SELENOM associated with lower, and DIO3 expression associated with higher mortality.These opposing Se-dependent effects of DIO1 and DIO3 imply a mechanism of action involving alteration of local thyroid hormone status, in agreement with prior data and experiments with model systems [23,24].Selenium substitution particularly in patients with DIO1 expressing tumours may improve survival, which should be considered for testing as an adjuvant therapy in randomized controlled trials.
Aberrations in the genome of breast cancer lead to alterations in expression of various genes, and hence to alterations in protein levels.Some of these proteins, such as HER2, ER, PGR are involved in key-regulatory mechanisms of breast cancer progression [25,26].Therefore, assessment of RNA-expression as indirect measure of actual protein levels constitutes an increasingly popular tool to predict prognosis [27][28][29].For most proteins, RNA-expression has been shown to be a reliable proxy for actual protein levels [27].Expression of selenoproteins however, is subject to more complex regulation involving a key limiting Se-dependent step in the translational process [30].The incorporation of the characteristic Sec residue(s) via recoding of the UGA stop codon to a sense codon is an important regulator of translation during the biosynthesis of functional selenoproteins, which is mainly regulated by both transcript abundance and dietary intake of Se [31,32].Thus, RNA-sequencing of selenoproteins may reliably reflect true protein expression only in patients where Se intake is sufficiently high, ensuring saturated levels of circulating Se needed for optimal supply of tissues and high intracellular Se concentrations.In statistical terms, this hypothesis is described as a testable interaction between serum Se and gene expression levels in relation to survival.Fig. 2 A Analysis scheme.B Cox regression models in the whole cohort and in low and high selenium subsets, divided according to median selenium concentration of the cohort, i.e. 70.36 µg/L.All models were adjusted for age, tumour size, histological grade, lymph node involvement, expression of HER2/ER/PGR-Receptor, laterality of the tumour, and histological type.P for interaction was tested by adding an interaction term between serum selenium and the gene of interest, marked by purple asterisk In line with this hypothesis, we observed strong interaction effects, particularly for DIO1 and DIO3.Deiodinases are the most important regulators of thyroid hormone activity, which is essential for cellular proliferation and differentiation, and hence implicated in cancer progression and cancer mortality [14,33,34].Conflictingly, both increased and decreased circulating thyroid hormone levels have been linked to breast cancer survival, which may be attributed to the complex nature of thyroid hormone transport, metabolism and action [35][36][37][38].Local regulators, such as thyroid hormone transporters, receptors and deiodinases play a critical role in thyroid hormone metabolism and action.Recent studies have demonstrated that DIO3 is a prognostic factor in breast cancer [14], and that e.g.low thyroid hormone receptor alpha 2 expression is associated with higher breast cancer mortality [35].In addition to promoting proliferation, thyroid hormone action also mediates cellular differentiation in physiological processes [39].Active thyroid hormones have been found to induce differentiation into a more benign phenotype in hepatic cancer cells [40].In murine models of basal cell carcinoma, Dio3 knockdown with concomitant increase in local T3 led to a five-fold decrease in tumour growth [41].Our findings are in agreement with these studies, and suggest a Se-mediated potentiation of the favourable effects of DIO1 and the unfavourable effects of DIO3 on breast cancer survival.These findings indicate potential involvement of local thyroid hormone action in breast cancer progression, as DIO3 is the primary thyroid hormone inactivating enzyme, and DIO1 plays a crucial role in the deiodination of T4 to T3 [33].Further epidemiological and mechanistic studies are necessary to investigate this hypothesis in more detail.
The association of SELENOM with mortality was also modified by circulating Se levels, potentiating the favourable association with survival.SELENOM encodes selenoprotein M, which is located in the endoplasmic reticulum and involved in protein folding [42].Although little is known about the specific roles of selenoprotein M, its functional homolog selenoprotein F (SELENOF) has been shown to be involved in cancer progression [43,44].In line with our findings of favourable effects of SELENOM on survival, SELENOF was recently described as a tumour suppressor in breast cancer, and enhancing its expression reduced tumour growth both in vivo and in murine breast cancer models [45].Similar to SELE-NOF, and in line with our findings higher expression of SELENOM has been shown to be a protective prognostic factor in other cancer types, such as gastric cancer and cholangiocarcinoma [46,47].
Interactions with circulating Se were observed for three of the 23 tested selenoprotein genes only.One likely explanation for this finding is based on the organ-and selenoprotein-specific hierarchical regulation of selenoprotein expression, which serves to provide regular functioning of essential tissues in a Se deficient state.Hence, some of the selenoproteins are preferentially sustained in case of Se deficiency, while others display a more responsive decrease in activity and reduced expression when the supply is low [48].Although this hierarchy is further regulated specific to different tissues, DIO1 has been shown to be one of the most responsive selenoproteins in the liver, displaying 95% decrease in activity in rats with severe Se deficiency [49].This is in line with our findings, as the interaction between DIO1 and circulating Se levels was the most prominent.Another explanation is the heterogenous distribution of selenoprotein gene expression within cells residing in the tumour and its microenvironment, as outlined in Fig. 1E.While DIO1 as an instance is mostly abundant in cancer cells, DIO2 shows nearly no expression in malignant cells, but rather in cancer associated fibroblasts, which may explain the lack of interaction for this member of the deiodinase family in our study.
Previously in the cohort used in this study, serum Se, serum SELENOP concentrations, activity of the serum GPx3 as well as novel autoantibodies targeting SELENOP were shown to be independent predictors of survival after breast cancer diagnosis, outperforming established clinical prognostic markers [4,5].Other large breast cancer studies are in line with these findings [6][7][8], and the effects seem to be also consistent for some other cancer entities [50][51][52][53][54].In this study, we aimed to further exploring the potential mechanisms of action underlining the consistent associations, and to examine whether patients with tumours displaying certain selenoprotein (See figure on next page.)Fig. 3 A Contour plot of the interaction between DIO1 expression and serum selenium concentrations on mortality.B Cox regression models depicting the association of DIO1 expression with mortality according to 10th, 50th and 90th quantiles of circulating selenium.C Contour plot of the interaction between DIO3 expression and serum selenium concentrations on mortality.D Cox regression models depicting the association of DIO3 expression with mortality according to 10th, 50th and 90th quantiles of circulating selenium.E Contour plot of the interaction between SELENOM expression and serum selenium concentrations on mortality.F Cox regression models depicting the association of SELENOM expression with mortality according to 10th, 50th and 90th quantiles of circulating selenium.All models were adjusted for age, tumour size, histological grade, lymph node involvement, expression of HER2/ER/PGR-Receptor, laterality of the tumour, and histological type gene expression profiles are particularly likely to benefit from a higher Se status.Our results indicate a potential mechanism of action in local thyroid hormone action due to significant and inverse interactions with DIO1 and the thyroid hormone inactivating DIO3.Clinically, our results indicate that Se-deficient patients with DIO1 expressing tumours may distinctly benefit from Se substitution, whereas patients with high DIO3 may not.
To the best of our knowledge, this study represents the first attempt to investigate the relationship between serum Se levels and RNA levels of selenoproteins in tumour tissues simultaneously in relation to survival, providing a more comprehensive understanding of the complex relationship between Se, selenoproteins, and cancer progression.Another noteworthy strength includes the large sample size as well as the large number of covariates assessed by physicians and pathologists.This allowed for conducting complex interaction analyses that require large sample sizes.The nearly complete database enabled studying effects independent of established clinical prognostic factors.The primary outcome, all-cause mortality derives from the Swedish National Registry, and the covariates were extracted from NKBC, which exhibited over 99.9% completeness and over 90% validity in an independent validation study conducted at time of participant recruitment of this study [20].SCAN-B has been fully integrated into clinical routine, without interfering with clinical decision making, which increases generalizability of the results.Importantly, Se status was measured by multiple biomarkers, all linearly correlating, which ensures a correct quantification of the main exposure.All Se measurements were made in a doubleblinded fashion, reducing risk of bias and increasing internal validity of the results.
Despite the strengths, a significant limitation is the observational nature, which precludes the ability to infer causality.Although adjustment was done for most important potential confounders, there may be residual confounding.A set of patients were excluded due to missing variables, although missingness was very low (Additional file 1: Fig. S1).The study's population sample from Sweden may limit its generalizability to other populations and ethnicities other than European.Further studies in other populations are necessary to confirm the findings.

Conclusion
Our unbiased analysis of circulating Se levels and the breast cancer selenotranscriptome revealed that Se modifies (potentiates) the associations of type 1 deiodinase expression with a favourable prognosis.On a mechanistical aspect, the study contributes to a growing evidence of an importance of thyroid hormones on cancer progression.Specifically, the favourable effects of DIO1 and opposing effects of DIO3 together argue for a beneficial effect of in-tissue local thyroid hormone action.
Clinically, the results provide a translational value, as serum Se measurement and (histo)pathological assessment of DIO1 expression in tumours of breast cancer patients undergoing surgery can be integrated into clinical routine to stratify patients according to potential benefit from Se substitution.As our data provides observational evidence, however, this potential improved survival remains to be tested in well-designed RCTs.Nevertheless, this approach can readily be used in clinical routine in order to provide (surrogate) prognostic value.

Fig. 1 A
Fig. 1 A Study scheme.B Gene expression of selenoproteins in tumour samples of 1453 patients.C Spearman's correlation matrix of tumour gene expression of selenoprotein genes and circulating selenium biomarker concentrations.D Spearman's correlation of serum biomarkers with each other and serum biomarkers with SELENOW expression in the tumour.E Single cell RNA-expression of selenoprotein genes in different cells in breast cancer

Table 1
Baseline patients characteristics according to vital status

Table 2
Therapy regimens according to vital status