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 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).

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/).

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.

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.

Fluorescent multiplex immunohistochemistry (mIHC) was applied to characterize immune infiltration in PTC. After incubation with primary antibody, the slides were incubated with an HRP-conjugated secondary antibody and fluorophore-conjugated Tyramide Signal Amplification reagent (Cat. 10003100100, Panovue). Multiplex panel design was as follows: FoxP3 (diluted 1:100, ab215206, abcam), Opal 570 TSA (diluted 1:200); CD8 (diluted 1:100, 70306S, CST), Opal 690 TSA (diluted 1:200). Finally, the slides were stained with DAPI (4’, 6-diamidino-2-phenylindole) and mounted with anti-fated reagent.

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.

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.

We investigated the relationship between STRA6 expression and clinicopathological features in the TCGA cohort. High STRA6 expression was associated with BRAF mutation (P = 0.001) and distant metastasis (P = 0.019) ( Table 1 ). Then we divided PTC patients into two groups according to the genotype of BRAF. Intriguingly, high STRA6 expression was significantly associated with aggressive clinicopathological outcomes in BRAF-mutant PTC, but not in BRAF WT groups. For individuals with BRAF mutation, rates of distant metastasis were 2/130 (1.5%) vs. 3/26 (11.5%) (P = 0.033) in low-STRA6 group vs. high-STRA6 group, respectively. The recurrence rates were significantly higher in high-STRA6 PTC than low-STRA6 PTC [11/40 (27.5%) vs. 20/212 (9.4%), P = 0.003]. As for mortality, it was robustly higher in PTC with high STRA6 expression than those with low STRA6 expression [5/40 (12.5%) vs. 4/212 (1.9%), P = 0.004] ( Table 2 ). Kaplan-Meier analysis showed higher STRA6 expression significantly predicted poorer recurrence free survival (RFS) (P = 0.0007) and overall survival (OS) (P = 0.0006) ( Figures 2A, B ).

Furthermore, we explored the hazard ratios (HRs) of STRA6-related 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 ). STRA6 overexpression was a risk factor for mortality with an HR being 7.12 (1.91-26.56) (P = 0.003), which became more significant with an HR being 19.02 (3.48-103.94) (P = 0.001) after multivariate adjustment ( Figure 2D ). Therefore, individuals with high STRA6 expression show increased mortality risk in BRAF-mutated PTC.

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.

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 confirm the role of STRA6 in Tregs and CD8+ T cells infiltration, 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.

(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 TNFα/NF-κB 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.