STRA6 regulates tumor immune microenvironment and is a prognostic marker in BRAF-mutant papillary thyroid carcinoma

Background BRAF mutation is one of the most common genetic alterations contributing to the initiation and progression of papillary thyroid carcinoma (PTC). However, the prognostic value of BRAF mutation for PTC is limited. Novel markers are needed to identify BRAF-mutant patients with poor prognosis. Methods Transcriptional expression data were downloaded from the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets. Pathway enrichment was performed by Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis and gene set enrichment analysis (GSEA). Protein-protein interaction networks were predicted by the GeneMANIA. The correlation between STRA6 expression and immune infiltration was analyzed by tumor immune estimation resource (TIMER) and tumor-immune system interaction database (TISIDB). Immunohistochemistry was used to detect the STRA6 protein expression level of PTC. Infiltration of regulatory T cells (Tregs) and CD8+ T cells in tumor samples were analyzed by fluorescent multiplex immunohistochemistry. Results In BRAF-mutant PTC, STRA6 was extremely upregulated and predicted unfavorable survival, which was an independent risk factor for increased mortality risk. Bioinformatic analyses indicated that STRA6 might activate the MAPK pathway synergistically with BRAFV600E. The expression of STRA6 was associated with immune infiltrates and T cell exhaustion. Fluorescent multiplex immunohistochemistry showed that STRA6 increased Tregs abundance and decreased CD8+ T cells infiltration in PTC. Moreover, STRA6 promoted epithelial-mesenchymal transition via increased cancer-associated fibroblasts infiltration. Conclusions Our study demonstrates STRA6 may serve as a prognostic marker for BRAF-mutated PTC, which may drive thyroid carcinogenesis via activation of oncogenic pathway and regulation of tumor immunosuppressive microenvironment.


Introduction
Thyroid carcinoma (TC) is the most common endocrine malignancy with an increasing incidence worldwide (1). Most of TCs arise from thyroid follicular epithelial cells, among which 80% are papillary thyroid carcinoma (PTC) (2). Initiation and progression of TC involve multiple genetic alterations, of which the most common oncogenic mutation is BRAF V600E (3,4). Although BRAF V600E mutation is associated with poorer recurrence-free probability (5), it was suggested that not all of the PTCs with BRAF V600E mutation belonged to high-risk recurrence in American Thyroid Association Management Guidelines (6) or even were associated with aggressive clinicopathological outcomes (7,8). Thus, novel markers are needed to identify BRAF-mutant PTC patients with poor prognosis.
STRA6, functions as a Vitamin A transporter and cytokine receptor (9), has been reported in multiple cancers. STRA6 activated the Wnt/b-catenin signaling in gastric cancer and served as a prognostic biomarker (10,11). Upregulation of STRA6 in colon cancer was correlated with poor oncologic prognosis and transduced JAK2-STAT signaling, which contributes to colon tumorigenesis (12, 13). Moreover, STRA6 polymorphisms were correlated with clinical characteristics and might be potential markers in locally-advanced and metastatic non-small cell lung cancer (14). However, the prognostic value of STRA6 in PTC remains unknown.
Tumor microenvironment, composed of tumor cells, immune cells, cytokines, etc., plays an essential role in cancer metastasis, progression and recurrence. The escape of tumor cells often occurs in an immunosuppressive microenvironment (15). It was reported that regulatory T cells (Tregs) were enriched in primary thyroid tumors (16) and had higher infiltration in lymph node metastases (17). Besides, the density of tumor-associated macrophages was increased in advanced thyroid cancer, which correlated with decreased survival rate of patients (18). Recent study revealed that STRA6 was associated with infiltration of antigen-presenting cells in stomach adenocarcinoma, which could serve as a potential mRNA vaccine candidate (19). Nevertheless, the engagement of STRA6 in tumor immune microenvironment of TC remains elusive.
In this study, we demonstrate STRA6 is upregulated in BRAFmutant PTC. High STRA6 expression is associated with unfavorable prognosis for individuals with BRAF V600E . Immune infiltration and pathway enrichment were investigated in samples with high STRA6 expression. STRA6 may contribute to the progression of BRAF-mutant PTC via regulation of immunosuppressive tumor microenvironment and activation of oncogenic pathway.

Transcriptional data
Transcriptional expression data and clinicopathological information of the Cancer Genome Atlas (TCGA) cohort were downloaded from UCSC Xena (https://xena.ucsc.edu/). GSE199023 was downloaded from Gene Expression Omnibus (GEO) datasets (https://www.ncbi.nlm.nih.gov/). The RNA sequencing data of the FAH-SYSU cohort was obtained from our previous study (20).

Tumor immune estimation
Tumor immune estimation resource (TIMER) 2.0 (http://timer. cistrome.org/) is a web server for comprehensive analysis of immune infiltrates with diverse algorithms. It was applied to analyze STRA6 expression level in pan-cancer based on the TCGA dataset. The correlation of STRA6 expression between BRAF mutation status was further explored by Timer 2.0. It was also used to analyze the correlation of immune infiltrates with STRA6 expression and BRAF mutation status. The relation between immunomodulators and STRA6 expression was evaluated via tumor-immune system interaction database (TISIDB) (http://cis.hku.hk/TISIDB/index.php).

Differential expression analysis and functional annotations
NetworkAnalyst (https://www.networkanalyst.ca/) is a platform for comprehensive gene expression profiling analysis. We divided the BRAF-mutant PTC from the TCGA cohort into two groups according to STRA6 expression. After uploading a gene expression table, differently expressed genes (DEGs) were analyzed by the platform and followed with Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Additionally, principal component analysis (PCA) between two groups was performed on Clustvis website (https://biit. cs.ut.ee/clustvis/). Protein-protein interaction (PPI) was predicted by GeneMANIA (http://genemania.org/).

Gene set enrichment analysis
GSEA software was downloaded from the website (http://www. gsea-msigdb.org/gsea/index.jsp). We performed analysis in 50 patients with high STRA6 expression and 50 patients with low STRA6 expression of the TCGA cohort. Hallmark gene sets were used to explore the mechanism of STRA6 in thyroid carcinoma progression.

Immunohistochemistry
The study was approved by the Institutional Research Ethics Committee of The First Affiliated Hospital of Sun Yat-sen University (no. [2021] 109). 30 pairs of PTC tissues were collected from the First Affiliated Hospital of Sun Yat-sen University with informed consent from patients. Paraffin embedded tissues were deparaffinized with xylene for two times and dehydrated with a series of ethanol. Antigen retrieval was performed by boiling the slides in EDTA buffer (C1034, Solarbio). After cooling, the slides were blocked with 20% goat serum at room temperature for 30 min and incubated with STRA6 primary antibody (diluted 1:200, H00064220-D01P, Novus) at 4°C overnight. Followed by rinse with PBS for 3 times, the slices were incubated with secondary antibody at room temperature for 30 min. Detection of protein expression was performed with a 3,3'-diaminobenzidine (DAB) peroxidase substrate kit (K5007; Dako) under a microscope. Subsequently, the slides were counterstained with Harris hematoxylin, followed by dehydration in graded alcohol and incubation in xylene. The immunohistochemistry (IHC) staining score was calculated with the following formula: IHC score = extent score × intensity score. The staining extent was scored as 0 (1%-4% positive cells), 1(5%-25% positive cells), 2 (26%-50% positive cells), 3 (51%-75% positive cells), 4 (≥75% positive cells). The staining intensity was scored as 0-3, corresponding to no staining, weak staining, intermediate staining and strong staining, respectively.

Statistical analysis
Experimental data was analyzed with GraphPad Prism 8.0 and presented as mean ± standard deviation (SD). Comparison between two groups was analyzed with Mann-Whitney U test or Student's t test. The correlation between STRA6 expression and clinicopathological features in the TCGA cohort was analyzed with a chi-squared test. Kaplan-Meier analysis with a log-rank test was used to evaluate the prognostic value of STRA6 in BRAF-mutant PTC. Univariate Cox regression analysis was performed to investigate the correlation between variables and survival. After multivariate adjustment for age, gender, tumor size, multifocality, extrathyroidal extension and lymph node metastasis, multivariate Cox regression analysis further accessed whether high STRA6 expression was an independent risk factor for recurrence and mortality. A P value of <0.05 was considered statistically significant.

STRA6 is upregulated in BRAF-mutated PTC.
To unveil the expression pattern of STRA6 across different cancers, we analyzed the mRNA expression level of STRA6 in the TCGA cohort. STRA6 was upregulated in most cancers including thyroid carcinoma ( Figure 1A). Notably, the expression of STRA6 was positively correlated with BRAF mutation in TC and melanoma ( Figure S1A). Compared with BRAF wild-type (WT) PTC, STRA6 was significantly upregulated in BRAF-mutated samples ( Figure 1B). Based on a public dataset GSE199023, zebrafish thyroid tumors transfected with BRAF V600E oncogenes showed higher expression of STRA6 as compared to those transfected with CCDC6-RET fusion ( Figure 1C). We further analyzed the mRNA expression in our cohort based on distinct molecular subtypes. STRA6 was remarkably upregulated in BRAF subtype compared with other PTC subtypes and benign nodules ( Figure 1D). Overexpression of STRA6 protein level in BRAF V600E PTC was also confirmed by IHC ( Figure 1E). Taken together, STRA6 is upregulated in thyroid cancer and exhibits extremely higher expression in BRAF-mutant PTC.
High STRA6 expression is an independent risk factor for mortality in BRAF-mutated PTC Furthermore, we explored the hazard ratios (HRs) of STRA6related risk (Table 3). In the entire TCGA cohort, univariate Cox regression analysis showed that upregulation of STRA6 was an independent risk factor for recurrence and mortality, with HRs being 2.98 (1.62-5.48) (P < 0.001) and 4.67 (1.69-12.87) (P = 0.003), respectively. After multivariate adjustment for age, gender, tumor size, multifocality, extrathyroidal extension and lymph node metastasis, HRs of recurrence and mortality remained significantly (P = 0.013 and P = 0.017). Then we divided patients into two groups according to BRAF mutation status. In BRAF wild-type patients, neither recurrence rates nor mortality rates was associated with STRA6 expression. In BRAF-mutant patients, STRA6-related recurrence risk was 3.32 (1.59-6.96) (P =0.001) and became 1.97 (0.81-4.81) (P = 0.135) after multivariate adjustment ( Figure 2C).   Figure 2D). Therefore, individuals with high STRA6 expression show increased mortality risk in BRAF-mutated PTC.

STRA6 plays oncogenic and immunoregulatory roles in BRAF-mutated PTC progression
To explore the mechanism of STRA6 contributing to the poor prognosis in BRAF-mutant PTC, we divided patients into two groups according to STRA6 expression. Principal component analysis (PCA) showed separation between high-STRA6 PTC and low-STRA6 PTC ( Figure 3A). Differential expression genes (DEGs) between two groups were shown in Figure 3B. Of which, most genes were involved in malignant tumor progression, such as OBP2B, TNNT1, KLK6, KRT17 etc. Moreover, chemokine signaling pathway and cytokine-cytokine receptor interaction were enriched with the upregulated genes in high-STRA6 groups, indicating STRA6 may engage in the regulation of the immune interaction ( Figure 3C). We also found that STRA6 expression was associated with immunomodulators in the TCGA cohort ( Figures S1B, C), which indicated that STRA6 may function as an immunoregulator in TC progression. Additionally, the MAPK signaling pathway, which has emerged as one of the frequently activated pathways in thyroid tumorigenesis, was also remarkably enriched with the upregulated genes in high-STRA6 PTC. Downregulation of genes enriched in thyroid hormone synthesis further suggested the oncogenic role of STRA6 in TC progression ( Figure 3C). PPI networks of STRA6 were constructed by GeneMANIA. It was predicted that 20 proteins were binding to STRA6. Of note, STRA6 interacted with several proteins modulated by the MAPK signaling pathway, including DIO3, BMP4 and DUSP6 ( Figure 3D). Therefore, our results demonstrate that STRA6 may contribute to BRAF-mutated PTC progression by regulating oncogenic pathways and immune microenvironment.

STRA6 is correlated with immune infiltration in thyroid carcinoma
Given immune infiltration takes part in tumor progression, we performed Gene Set Enrichment Analysis (GSEA) with transcriptome data from the TCGA cohort. The results showed that samples with high STRA6 expression enriched in the IL2-STAT5 signaling pathway ( Figure 4A). STRA6 expression was positively correlated with the mRNA level of IL2, STAT5a and STAT5b ( Figure 4B). IL2 is a critical immunomodulatory cytokine, which increases Tregs and induces CD8+ T cell exhaustion via activation of STAT5. Thus, we evaluated the association between STRA6 expression and immune cell infiltration in PTC samples. After purity adjustment, upregulation of STRA6 was correlated with abundance of Tregs. However, STRA6 expression was negatively correlated with the infiltration of CD8+ T cells ( Figure 4C). We also assessed the expression of exhausted T cells markers, including LAG3, PDCD1G2, CTLA4, HAVCR2, PRDM1 and GZMB in the TCGA cohort. STRA6 was positively correlated with most T cell exhaustion markers ( Figure 4D). Pan-caner analysis also revealed the role of STRA6 in inducing T cell exhaustion ( Figure 4E). Tregs were significantly abundant in high-STRA6 PTC, while CD8+ T cells had a lower infiltration level ( Figure 4F). The infiltration of other immune cells was also correlated with STRA6 expression, including NK cells, neutrophils, MDSCs, macrophages, myeloid dendritic cells, monocytes and B cells ( Figure S2A).
To confirmtheroleofSTRA6inTregsandCD8+Tcellsinfiltration, we performed fluorescent multiplex immunohistochemistry (mIHC) in PTC with different STRA6 expression. Samples with high STRA6 expression were mainly infiltrated with FoxP3+ Tregs while CD8+ T cells were mainly observed in samples with low STRA6 expression (Figures 5A,  B). Adjacent thyroid tissues were less infiltrated with immune cells ( Figure 5C). Upregulation of STRA6 was associated with lower CD8 + cells and higher Tregs infiltrates ( Figure 5D), which was consistent with the infiltration level of immune cells in PTC with BRAF mutation (Figures S3A, B). Therefore, STRA6 may limit antitumor immunity and induce immunosuppressive tumor microenvironment in TC. Frontiers in Endocrinology frontiersin.org STRA6 promotes epithelial-mesenchymal transition via increased cancer-associated fibroblasts infiltration in thyroid carcinoma (CAF) Cancer-associated fibroblast is one of the important components in tumor microenvironment. CAFs infiltration was positively associated with the expression of STRA6 ( Figure 6A). We also analyzed the expression level of CAFs markers in PTC. COLA1, COLA2 and COL3A1 were all upregulated in high-STRA6 samples ( Figure 6B). Compared with low-STRA6 group, CAFs were abundant in high-STRA6 group ( Figure 6C). Higher infiltration level of CAFs was also observed in BRAF-mutant PTC ( Figure S3C). CAFs play important roles in cancer development via regulation of metastasis. The results of GSEA showed that the TNFa/NF-kB signaling pathway and Epithelial-Mesenchymal Transition (EMT) hallmark gene set was significantly enriched in high-STRA6 PTC (Figures 6D, E). Consistently, the mRNA expression of STRA6 was positively correlated with the expression of EMT markers, including Vimentin (r = 0.1446), MMP7 (r = 0.5345), MMP9 (r = 0.3594), MMP14 (r = 0.4984) ( Figure 6F). Taken together, STRA6 may promote TC metastasis via CAF-mediated EMT process.

Discussion
As the most common genetic alteration in PTC, the prevalence of BRAF mutation differs in races. Compared with western populations (3), the incidence of BRAF mutation was much higher in Chinese population (4,21), which was detected in about 75% PTC patients. However, BRAF mutation may not be a single and independent prognostic marker for PTC. Studies revealed that BRAF mutation may not totally associate with clinicopathological characteristics or survival (22,23). In our study, we found STRA6 is associated with BRAF genotype. Upregulation of STRA6 predicts poor RFS and OS in BRAFmutant PTC, which is an independent risk factor for mortality. Previous study showed that BRAF V600E mutation was significantly correlated with PTC-related mortality in an unadjusted analysis and the association was no longer significant after adjusting for multiple clinical factors (24). Thus, our study indicated that STRA6 may act as a prognostic marker, which may provide individualized diagnosis and therapeutic strategy for BRAFmutant PTC. We further explored the mechanism of STRA6 in BRAF-mutant PTC progression. Interestingly, when we divided BRAF-mutant PTC into two groups according to STRA6 expression, the MAPK signaling pathway is significantly enriched with the upregulated genes. STRA6 also interacts with several proteins regulated by the MAPK pathway, including BMP4, DIO3, and DUSP6. It was reported that BMP4 can inhibit heat-induced apoptosis by enhancing the activation of ERK (25). DIO3 is increased in PTC through activation of the MAPK pathway (26,27). DUSP6 is a member of the MAPK phosphatase family and plays pro-tumorigenic role in TC (28,29). Moreover, BRAF mutation drives aberrant activation of the MAPK pathway (30) and promotes TC tumorigenesis. Dual activation of the MAPK pathway by STRA6 and BRAF mutation, may contribute to the poor survival of BRAF-mutant PTC with high STRA6 expression.
Pathway analyses showed that STRA6 involves in the regulation of the immune interaction and activation of the IL2-STAT5 pathway. Our study also revealed that STRA6 expression is correlated with the infiltration of multiple immune cells. Notably, upregulation of STRA6 results in higher Tregs abundance and lower CD8+ T cells infiltration. Additionally, STRA6 is correlated with T cell exhaustion. In support of our results, IL-2 is mainly produced by activated CD4+ T cells and induces the expansion of Tregs (31). Persistent activation of the IL2-STAT5 pathway triggers CD8+ T cells into exhausted status (32). Decreased infiltration of CD8+ T cells and a low CD8/regulatory T cell ratio in tumors are critical for unfavorable prognosis of cancers (33,34). Our results also showed that high STRA6 expression is associated with increased CAFs and promotes EMT process, which further explains the worse prognosis in BRAF-mutant PTC. CAFs are one of the most abundant stromal components in the tumor microenvironment. By secreting TGF-b, CCL2, CCL5, IL-6 etc., CAFs recruit immunosuppressive cells into the tumor stroma, which contributes to cancer metastasis and progression (35). The secretion of various cytokines promotes aggressive phenotypes in tumor cells via the activation of EMT process (36). We also confirmed BRAF V600E mutation was associated with suppressive immune cell infiltration (37). In hence, immune suppression mediated by STRA6 and BRAF contributes to PTC progression synergistically.
Taken together, STRA6 is overexpressed in TC and associated with BRAF mutation. Upregulation of STRA6 predicts unfavorable prognosis for BRAF-mutant PTC. High STRA6 expression is correlated with decreased CD8+ T cells as well as abundant Tregs and CAFs, which induces immunosuppressive microenvironment in TC. Moreover, STRA6 may play an oncogenic role via the MAPK pathway and EMT process. In conclusion, STRA6 might be a prognostic marker for BRAF-mutant PTC and individualized therapeutic strategies for the patients are needed.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement
The studies involving human participants were reviewed and approved by the Institutional Research Ethics Committee of The First Affiliated Hospital of Sun Yat-sen University. The patients/ participants provided their written informed consent to participate in this study.