An integrated multiomics-based risk model for predicting the prognosis and immune microenvironment status in patients with papillary thyroid carcinoma
Highlight box
Key findings
• A three-gene signature, consisting of ALOX15B, RAP1GAP, and JUN, was developed and confirmed to be an independent predictor of progression-free survival in patients with papillary thyroid carcinoma (PTC). High-risk patients were found to have an immunosuppressive tumor immune microenvironment, higher tumor mutational burden, and a higher frequency of BRAF mutations. ALOX15B was discovered to promote PTC cell proliferation, migration, and invasion in vitro.
What is known and what is new?
• Although the majority of patients with PTC have a good prognosis, a portion experience recurrence or disease progression, and this heterogeneity cannot be fully accounted for by current clinicopathological assessment.
• This study employed multiomics analysis, Mendelian randomization, external validation, and functional experiments to establish a robust three-gene prognostic model. The findings suggested an association of poor prognosis with immune suppression and genomic alterations, and ALOX15B was identified as a potential functional driver in PTC.
What is the implication, and what should change now?
• This signature may improve molecular risk stratification in patients with PTC and support further evaluation of ALOX15B as a biomarker and potential therapeutic target. Prospective clinical validation is needed before the novel signature and biomarker can be applied routinely in clinic.
Introduction
Papillary thyroid carcinoma (PTC) represents the most prevalent type of thyroid cancer, accounting for approximately 80% of all thyroid malignancies (1-3). It primarily affects individuals in their third to fifth decades of life, with a higher incidence among women than among men. Although the majority of patients with PTC experience favorable outcomes following standard treatment, a subset undergo disease recurrence, progression, or even metastasis, leading to a significantly worse prognosis (4-7). The current clinicopathological staging systems, such as the American Joint Committee on Cancer (AJCC) tumor-node-metastasis (TNM) classification, remain insufficient for precisely predicting the heterogeneous clinical outcomes of patients with PTC. Consequently, there is a pressing need to identify molecular biomarkers capable of enhancing risk prediction and further clarifying the biological processes underlying PTC progression.
Advances in molecular biology and public human genetic databases have provided new opportunities for identifying disease-associated genes and candidate therapeutic targets in PTC (8,9). Genome-wide association studies (GWASs) have identified multiple single-nucleotide polymorphisms (SNPs) associated with thyroid carcinoma (10,11). However, many GWAS-identified variants are located in noncoding or intergenic regions, which limits their direct interpretation and their ability to pinpoint functional driver genes (12).
Mendelian randomization (MR), especially when integrated with expression quantitative trait loci (eQTL) data, can help prioritize genes with potential causal relevance to disease risk (13-19). By combining genetic association data with gene expression profiles, this approach may provide additional evidence for linking genetic variation, gene expression, and disease susceptibility. Therefore, integrated analysis of GWAS, eQTL, and transcriptomic data represents a useful strategy for identifying candidate genes involved in PTC pathogenesis.
The development and prognosis of PTC are influenced not only by genetic alterations but also by the dynamic tumor immune microenvironment (TIME) (20-24). Infiltrating immune cells can either suppress or promote tumor progression, and their composition is increasingly recognized as an important correlate of clinical outcome (25-27). However, the immune infiltration landscape of PTC and its association with molecular risk features remain incompletely understood. High-throughput sequencing and bioinformatic methodologies provide opportunities to explore these associations and to construct molecular models for prognosis prediction.
In this study, we employed an integrated approach combining genomics and transcriptomics to identify robust genetic determinants of PTC pathogenesis and prognosis. We first performed differential expression analysis on multiple Gene Expression Omnibus (GEO) datasets and applied MR analysis using GWAS data to infer causal relationships between gene expression and PTC risk. Furthermore, we characterized the immune infiltration patterns in PTC using the CIBERSORT algorithm and developed a prognostic gene signature using least absolute shrinkage and selection operator (LASSO) Cox regression. Finally, in vitro functional experiments were conducted to validate the role of ALOX15B in regulating proliferation, migration, and invasion in PTC cells. The aim of this work was to establish a comprehensive molecular framework for understanding PTC progression and to generate hypotheses for future prognostic and mechanistic studies. We present this article in accordance with the TRIPOD reporting checklist (available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0267/rc).
Methods
Data source
Gene expression profiles and corresponding clinical information related to PTC in humans were collected via microarray dataset analysis. All materials were sourced from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which provides downloadable gene expression data and platform probe annotations. To ensure the quality of the dataset, filtering criteria were applied. This included a requirement for a minimum of eight samples, with at least four PTC and four matched nontumor controls. Additionally, all samples were required to be free from chemical treatment or genetic manipulation, and only datasets with publicly accessible raw data or array gene expression profiles in the GEO were considered. The RNA-sequencing (RNA-seq) transcriptome data, clinical information, and mutation data of patients with PTC were collected from the Genomic Data Commons Data Portal within The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/).
Differential expression analysis
Initial processing of the six individual datasets (GSE129562, GSE29265, GSE3678, GSE50901, GSE6004, and GSE97001) was completed in R version 4.4.0 (The R Foundation for Statistical Computing, Vienna, Austria). Following the integration of datasets, a batch correction was conducted on a total of 74 samples comprising 37 PTC and 37 matched nontumor control samples. Differential expression screening was implemented through the limma package in R with classical Bayesian data analysis, with the criteria of P<0.05 and |log2 fold change (FC)| >1 being applied to screen for significant differentially expressed genes (DEGs). Result visualization included volcano plots and hierarchical clustering heatmaps created with the R pheatmap package. Normalization and standardization procedures were conducted via expression matrices and probe annotation files from the GEO. To address batch effects and improve visualization, principal component analysis (PCA) was conducted via the prcomp function in R. This analysis aided in the identification of key genes distinguishing PTC samples from normal controls, providing further validation of the findings.
eQTL analysis of exposure data
To detect genetic variations associated with gene expression, an analysis of eQTL was carried out through use of both transcriptome and genotype information from a variety of study populations. The summary eQTL data used in this investigation were obtained from the GWAS Catalog (https://gwas.mrcieu.ac.uk/). Highly correlated SNPs (P<5e−08) were identified as instrumental variables via the TwoSampleMR software package in R. The criteria for linkage disequilibrium were set at r2<0.001 and clumping distance of 10,000 kb. An F-test value >10 was used as a filter to exclude SNPs that show weak associations or failed to adequately explain phenotypic variation.
Determination of outcome data
Outcome data for this study were derived from the MRC Integrative Epidemiology Unit Open GWAS project repository (https://gwas.mrcieu.ac.uk). The specific GWAS ID used for the analysis was ieu-a-1082. It is worth noting that all GWAS statistical summaries cited in this study are freely downloadable to the public. All original analyses received ethical approval, and all procedures were conducted in accordance with relevant ethical standards.
MR analysis
To evaluate the association between specific genes and PTC, we implemented MR analysis with the “TwoSampleMR” R package, applying the inverse variance weighting (IVW) approach as the primary analytical method. In addition, we performed multiple sensitivity analyses using MR-Egger, simple mode, weighted median, and weighted mode methodologies. The three-criteria method was used to identify disease-related genes: first, the IVW method was used to select genes with P<0.05; second, genes were further prioritized based on the consistency of effect estimates [odds ratios (ORs)] across the five MR methods; finally, pleiotropic genes with P<0.05 were excluded.
Identification of the differential genes associated with PTC
An overlap analysis of these disease-related genes and DEGs between the PTC samples and control samples was conducted. Following this, the intersected genes were individually subjected to MR analysis to determine their causal relationships with PTC. Heterogeneity tests, pleiotropy tests, and leave-one-out sensitivity analyses were performed throughout the study to assess the consistency and dependability of the outcomes. In order to visually portray and support these findings, scatter plots, forest plots, and funnel plots were generated and examined.
GEO and TCGA database verification
The dataset GSE33630 and TCGA-PTC were analyzed with R software to validate differences in PTC-related genes between PTC samples and control samples. The results were then compared with our MR analysis. The preprocessing methods used were the same as those described above.
Prognostic model
Machine learning algorithms were applied to discern the feature genes among the differential genes related to PTC. The glmnet R package was used to implement the LASSO method. Subsequently, receiver operating characteristic (ROC) curves for the candidate genes were generated in the TCGA-PTC dataset via the survivalROC R package. A risk score was derived from the predictive gene signature, and patients were stratified into high- and low-risk groups based on the median value. Kaplan-Meier analysis was employed to produce progression-free survival (PFS) curves, and differences between risk groups were assessed via the log-rank test. To evaluate whether the gene signature independently predicted prognosis, univariate and multivariate Cox regression analyses were performed. A nomogram incorporating the risk score, age, gender, and TNM stage was constructed, and its predictive accuracy for 1-, 3-, and 5-year PFS was assessed.
Immune infiltration analysis of the model
The relative abundance of 22 immune cell subsets within TCGA-PTC specimens was evaluated via the CIBERSORT algorithm. Subsequently, statistical comparisons between sample groups were conducted via the Wilcoxon test (P<0.05). The distribution of immune cell scores was visualized through box plots.
Clinical samples
In this study, thyroid cancer and matched adjacent paracarcinoma tissue samples were obtained from surgical specimen archives at Sir Run Run Shaw Hospital. Notably, none of the patients received chemotherapy or radiotherapy prior to undergoing the radical operation. To ensure accurate diagnoses, all tissue samples included in the study were independently assessed by at least two senior pathologists in compliance with the World Health Organization’s criteria. Following surgery, the samples were promptly preserved in liquid nitrogen to maintain their integrity. Prior to participating in the study, each patient provided written informed consent, ensuring that they were fully aware of the research procedures and implications. This study was approved by the Ethics Committee of Sir Run Run Shaw Hospital (approval No. 20260302-065). This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.
Cell culture
The cell lines used in the study, including HEK-293T, human thyroid follicular epithelial Nthy-ori3-1, and various human PTC cell lines (TPC-1, B-CPAP, K1, and KTC-1) were procured from American Type Culture Collection (ATCC; Manassas, VI, USA). Cells were maintained in Dulbecco’s Modified Eagle Medium (DMEM) containing 10% fetal bovine serum (Gibco, Thermo Fisher Scientific, Waltham, MA, USA), 100 U/mL of penicillin, and 100 µg/mL of streptomycin (Invitrogen, Thermo Fisher Scientific). Cultures were incubated at 37 ℃ in a humidified atmosphere with 5% CO2, and regular mycoplasma testing was conducted to verify the absence of contamination.
Lentivirus infection
The vector pMSCV/ALOX15B, which is designed to overexpress human ALOX15B, was developed by subcloning the coding sequence of human ALOX15B, which had been amplified through polymerase chain reaction (PCR), into the pMSCV vector sourced from Clontech Laboratories (Takara Bio, Kusatsu, Japan). This strategic construction allows for a significant increase in the expression levels of ALOX15B within the host cells. To achieve the silencing of the native ALOX15B expression, we proceeded to clone three distinct short hairpin RNA (shRNA) oligonucleotides. These oligonucleotides include shRNA-1 (AUCCAGACCAAUGUCAUUAAU), shRNA-2 (UGGAAGACUUGGAGCUCAAUA), and shRNA-3 (CUCUGAAAUCAUCGGUAUCUA), which were incorporated into the vector pSuper-retro-puro to create pSuper-retro-ALOX15B-shRNA(s). The medium that included the lentiviral particles was collected, and ultracentrifugation was conducted to concentrate the lentivirus. Infection with lentivirus was carried out for 48 hours and was followed by selection via 2 µg/mL of puromycin (Gibco). The cells that remained after antibiotic selection were combined to create stable mass transfectants.
RNA extraction and reverse transcription quantitative PCR (RT-qPCR)
Total RNA was extracted from cells and tissues via TRIzol reagent (Invitrogen) according to the manufacturer’s protocol. For complementary DNA (cDNA) synthesis, we utilized the Mir-X miR First-Strand Synthesis Kit (Takara Bio), and quantitative PCR analysis (qPCR) was conducted via SYBR Premix Ex Taq II (Takara Bio). Actin served as an internal standard control. Relative RNA expression data were analyzed through the comparative Ct method (2−ΔΔCt method). We designed all primers and confirmed their specificity using the National Center for Biotechnology Information (NCBI) primer blast tools. All primer sequences we procured from Tsingke Biological (Beijing, China).
Western blotting analysis and antibodies
Cells were lysed on ice for 30 minutes in radioimmunoprecipitation assay buffer (cat. No. BL509A; Biosharp, Hefei, China) supplemented with a protease inhibitor cocktail (cat. No. BL6049A; Biosharp). Equivalent quantities of total protein were placed on sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) gels for electrophoresis at 120 V for 1.5 hours and subsequently transferred onto 0.2 µm polyvinylidene fluoride (PVDF) membranes (MilliporeSigma, Burlington, MA, USA) at 4 ℃ for 1.5 hours. The membranes were obstructed with a solution of Tris-buffered saline with Tween 20 (TBST) and 5% skimmed milk powder for 1 hour at room temperature and were then exposed to primary antibodies overnight at 4 ℃. Following thorough rinsing with TBST, the membranes were treated with horseradish peroxidase-conjugated secondary antibody for 1.5 hours at room temperature, which was accompanied by three subsequent rinses with TBST. Finally, the protein bands were visualized with an enhanced chemiluminescence kit (cat. No. P0018AS; Beyotime Biotechnology, Shanghai, China) with an Amersham Imager 600 (GE HealthCare, Chicago, IL, USA).
Cell Counting Kit-8 (CCK-8) proliferation assay
Transfected cells were seeded into 96-well plates at consistent densities. At days 0, 1, 2, 3, 4, and 5, each well received 10 µL of CCK-8 reagent. Following a 2-hour incubation at 37 ℃, absorbance at 450 nm was measured under a spectrophotometer (Thermo Fisher Scientific, Vantaa, Finland).
Wound-healing assay
After transfection, cells were evenly distributed in six-well plates and cultured at 37 ℃ and 5% CO2 to form a monolayer of cells. When the cells reached 50–60% confluence, the center of each well was scraped with a 10 µL pipette tip. This was followed by three washes with phosphate-buffered saline to remove any detached or dead cells, and fresh serum-free medium was added to maintain cell growth. Cell migration was observed at 0 and 24 hours after injury, and images were obtained in the same field of view with an inverted microscope. Wound closure was quantified by measuring the distance between wound edges and normalizing to the 0-hour control. Experiments were performed in at least three independent replicates.
Transwell migration assays
Transwell plates (8 µm; Labselect, Labgic, Hefei, China) were used for migration experiments. For migration assays, 2×105 cells were seeded in 200 µL of serum-free medium in the upper chamber, while 600 µL of complete medium was placed in the lower chamber. After 24 hours of incubation, cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. Images were acquired with a Primovert inverted light microscope (Zeiss, Oberkochen). Migrated cells were quantified by counting at least three randomly selected fields per well.
Immunohistochemistry
Following antigen retrieval with DAKO Target Retrieval Solution (Agilent Technologies, Santa Clara, CA, USA) for 30 minutes, the sections were incubated with specific antibodies via the avidin-biotin complex system (Vector Laboratory, Newark, CA, USA) and visualized with 3,3'-diaminobenzidine serving as the substrate. Nuclei were counterstained with hematoxylin.
Statistical analyses
All statistical analyses were performed with R software version 4.4.0 and SPSS version 21.0 (IBM Corp., Armonk, NY, USA) software. Data were obtained from a minimum of three independent experiments and are presented as the mean ± standard deviation (SD). Data were analyzed with the unpaired Student’s t-test unless otherwise stated. Relationships between categorical variables were evaluated with the Chi-squared test. The Kaplan-Meier method and log-rank test were used to analyze the prognosis of patients with PTC. P<0.05 was considered statistically significant.
Results
Identification of the differential genes associated with PTC
To examine the potential involvement of messenger RNAs (mRNAs) in the pathogenesis of PTC, we performed an integrated analysis of RNA-seq data from 37 matched pairs of PTC and adjacent normal tissues, obtained from six bioinformatics datasets (GSE129562, GSE29265, GSE3678, GSE50901, GSE6004, and GSE97001). An overview of the included datasets is provided in Table S1. We used R version 4.4.0 to correct and merge the expression values of each gene across the respective datasets and employed PCA to eliminate batch effects. As shown in Figure 1A, the batch effects were evident in the six PTC gene datasets. After correction, all samples demonstrated acceptable uniformity, as shown by the PCA visualized in Figure 1B. We performed differential expression analysis on 37 pairs of PTC and adjacent tissues after batch correction. Using log2 FC and adjusted P values as thresholds, we identified 228 upregulated and 222 downregulated DEGs. Detailed information for these significant DEGs is provided in supplementary material 1 (available at https://cdn.amegroups.cn/static/public/gs-2026-0267-1.xls). The volcano plot of the merged GEO dataset is shown in Figure 1C. The top 50 upregulated and the top 50 downregulated DEGs are shown in Figure 1D.
SNP data were downloaded from the GWAS Catalog website, and based on P values, r2 values, clumping distance >10,000 kb, and an F-test value >10, we selected a total of 26,152 SNPs as instrumental variables (supplementary material 2 available at https://cdn.amegroups.cn/static/public/gs-2026-0267-2.xls). We then performed MR analysis, filtering genes with P<0.05 using the IVW method, further screening genes based on the consistency of OR values and excluding pleiotropic genes with P<0.05. Ultimately, 189 genes related to PTC were selected (supplementary material 3 available at https://cdn.amegroups.cn/static/public/gs-2026-0267-3.xls). We intersected these 189 PTC-related genes with DEGs, finally obtaining three PTC-related upregulated genes (TIAM1, ALOX15B, and TMC6) and four PTC-related downregulated genes (GPX3, RAP1GAP, JUN, and PAPSS2), as shown in Figure 1E,1F. To assess the causal effect of each gene on the disease, we performed MR-Egger, simple mode, weighted median, and weighted mode for further validation on these seven PTC-related genes. The results showed that all three up-regulated PTC-related genes were significantly and positively correlated with the outcomes of PTC [ALOX15B: OR =1.647, 95% confidence interval (CI): 1.120–2.420, P=0.01; TIAM1: OR =1.270, 95% CI: 1.001–1.611, P=0.049; TMC6: OR =1.250, 95% CI: 1.021–1.530, P=0.03] and all four downregulated PTC-related genes were significantly and negatively correlated with the outcomes of PTC (GPX3: OR =0.848, 95% CI: 0.721–0.998, P=0.047; JUN: OR =0.795, 95% CI: 0.653–0.967, P=0.02; PAPSS2: OR =0.779, 95% CI: 0.608–1.000, P=0.049; RAP1GAP: OR =0.895, 95% CI: 0.810–0.989, P=0.03) (Figure 1G). Statistical assessments detected no significant heterogeneity or pleiotropy among the PTC-related genes (P>0.05), indicating these factors would not substantially influence the findings. The subsequent leave-one-out sensitivity analysis revealed that the effect sizes of the included instrumental variables were nearly identical to the overall effect size, affirming the robustness of the analysis. Additional information about each gene, from scatter plots, forest plots, funnel plots, and leave-one-out sensitivity analyses, can be found in supplementary material 4 (available at https://cdn.amegroups.cn/static/public/gs-2026-0267-4.pdf). Chromosomal locations of the seven candidate genes are graphically displayed in Figure 1H.
Estimation of immune cell infiltration in PTC
We used the CIBERSORT algorithm to evaluate the relative proportions of 22 types of immune cells in normal thyroid and PTC tissue. We found that specific immune cell subtypes were significantly different between normal thyroid and PTC group samples (Figure 2A). Specifically, the proportion of γδ T cells was significantly decreased in PTC, and the proportions of resting dendritic cells and resting mast cells were significantly increased in PTC (Figure 2B). We further performed correlation analysis on these immune cells, and the results showed that the highest positive correlation was found between γδ T cells and M1 macrophages [correlation (cor) =0.6], and the highest negative correlation was found between naïve B cells and memory B cells (cor =−0.44) (Figure 2C).
Validation of the differential analysis and construction of a risk model
To further validate the findings from the MR analysis, we selected datasets from TCGA and GEO (GSE33630). The expression of TIAM1, TMC6, and ALOX15B were significantly higher in PTC than in normal thyroid tissues (P<0.001). In addition, the expression of PAPSS2, RAP1GAP, JUN and GPX3 were all significantly lower in PTC than in normal thyroid tissues (P<0.001) (Figure 2D,2E). The expression levels of three upregulated genes and four downregulated genes in the TCGA and GEO datasets (GSE33630) were consistent with the findings from our MR analysis, thereby reinforcing the robustness of the MR results. We also analyzed the correlation among the seven PTC-related genes to gain a deeper understanding of their interactions. The strongest positive correlation was between TIAM1 and TMC6 (cor =0.85), while the strongest negative correlation was observed between RAP1GAP and both TIAM1 and TMC6 (cor =−0.78) (Figure 2F). LASSO Cox regression analysis is a widely used multiple linear regression method that improves the predictive accuracy and interpretability of statistical models. It also allows for simultaneous variable selection and regularization, which is particularly suitable for high-dimensional data with poor correlation and salient predictive values. By optimizing feature selection and avoiding overfitting, the LASSO method effectively identifies the most available predictive markers to support the prediction of prognostic indicators for clinical results. A dashed vertical line represented the first rank value of log λ, which had the minimum segment likelihood deviation. Therefore, we used a TCGA dataset and identified three genes that were used to create a risk model for evaluating the prognosis of individuals with PTC (Figure 2G,2H). The risk model was constructed by combining the levels of expression from three specific genes, with their respective regression coefficients being input into the LASSO algorithm. Specifically, the risk score was calculated as follows: risk score = (0.175 × ALOX15B expression) + (−0.083 × RAP1GAP expression) + (−0.021 × JUN expression).
The risk model as an independent prognostic indicator
To assess the independent predictive capability of the risk model for the outcomes of patients with PTC, we generated a risk score for each patient according to the model. Subsequently, we categorized the patients into high- and low-risk groups, employing the median risk score as the cutoff point. The distribution of risk grades between the low- and high-risk groups can be clearly seen in Figure 3A. According to the Kaplan-Meier analysis, there was a significant difference in PFS between the two groups. Specifically, patients classified in the high-risk group exhibited a shorter PFS in comparison to those in the low-risk group (Figure 3B). The ROC curve served as the metric to evaluate the predictive performance of the risk model. The area under the curve (AUC) values for assessing the model’s accuracy in predicting the 1-, 3-, and 5-year PFS were 0.760, 0.700, and 0.708, respectively (Figure 3C). These findings demonstrate that the risk model possesses a potential capability to accurately forecast the prognosis of patients with PTC. We carried out both univariate and multivariate Cox regression analyses. In the univariate Cox regression analysis, the hazard ratio (HR) and the 95% CI for risk scores were 1.415 and 1.115–1.796 (P=0.004), respectively (Figure 3D). In the subsequent multivariate Cox regression analysis, the HR and 95% CI for risk scores were 1.355 and 1.052–1.746, respectively (P=0.02) (Figure 3E). These findings indicate that the risk model operates independently of clinicopathologic parameters, including gender, age, and TNM stage. In order to better evaluate the uniqueness of risk scores in forecasting the results of patients with PTC, we additionally examined the conformity index of the risk score. The findings indicated that the concordance index of the risk score exceeded that of alternative clinical variables over time, suggesting that the risk grade could more effectively predict the prognosis of PTC (Figure 3F). This suggests that utilizing the risk score may be more beneficial in accurately assessing and predicting outcomes for patients with PTC. A nomogram was created to forecast PFS over 1, 3, and 5 years through the combination of risk grade and clinical risk features. When compared with clinical variables, the prognostic model’s risk grade exhibited superior predictive performance within the nomogram (Figure 3G). Examination of correlation diagrams indicated excellent alignment between the anticipated and actual rates of PFS over 1-, 3-, and 5-year periods (Figure 3H). Furthermore, patients with PTC classified with high-risk scores demonstrated a significantly shorter PFS compared to those with low-risk scores across various clinical subgroups. This pattern was observed consistently in the age group, the gender group, the N grade group, and the TNM stage group (Figure S1). Consequently, these findings validated the ability of the risk model to independently predict the prognosis of patients with PTC.
Estimation of the TIME
The TIME in PTC samples was analyzed according to the risk model. We found that the immune score and Estimation of Stromal and Immune Cells in Malignant Tumor Tissues (ESTIMATE) score were significantly different between the low- and high-risk groups (Figure 4A). We next investigated the relative abundance of 22 types of immune cells in the two groups. Unsurprisingly, we discovered that the proportion of CD8 T cells was significantly lower in high-risk group, while the proportions of regulatory T cells (Tregs) were significantly higher (Figure 4B,4C). We next calculated the tumor mutation burden (TMB) scores utilizing the somatic mutation data derived from TCGA. The results revealed that the TMB in the high-risk group was significantly higher than that in the low-risk group, indicating a correlation between the risk model and the TMB (Figure 4D). We further annotated the mutations according to the variant effect predictor. We found that the 15 driver genes with the highest alteration frequency were significantly different between the low- and high-risk groups. Specifically, the high-risk group had a significantly higher mutation frequency of BRAF (73%) than did the low-risk group (47%), while the low-risk group had significantly higher mutation frequency of NRAS (11%) than did the high-risk group (5%) (Figure 4E,4F). Patients with PTC were then categorized into high-TMB (H-TMB) and low-TMB (L-TMB) groups based on the median TMB score. According to the Kaplan-Meier analysis, there was a significant difference in PFS between the two groups. Specifically, patients classified in the H-TMB group exhibited a shorter PFS in comparison to those in the L-TMB group (Figure 4G). Patients with H-TMB and L-TMB in the high-risk groups (defined as H-TMB + high risk and L-TMB + high risk, respectively) showed a poorer PFS compared to patients diagnosed with H-TMB and L-TMB in the low-risk groups (H-TMB + low risk and L-TMB + low risk, respectively). Notably, patients classified as L-TMB in the high-risk group (L-TMB + high risk) had even poorer PFS outcomes compared to patients identified as H-TMB in the low-risk group (H-TMB + low risk) (Figure 4H). These findings suggest that the risk score is associated with both mutation burden and computationally inferred immune features, but further experimental validation is required to determine their biological and clinical implications.
Interference of ALOX15B expression modulated the proliferation, migration, and invasion of PTC cells
To investigate the functional role of ALOX15B, RT-qPCR analysis of 12 paired PTC and adjacent normal tissue specimens revealed a significant upregulation of ALOX15B expression in PTC tumors compared to matched normal tissues (Figure 5A). Furthermore, immunohistochemical analysis confirmed the expression profile of ALOX15B in PTC specimens (Figure 5B). Among 60 samples from patients with PTC, 48 (80%) exhibited positive ALOX15B immunostaining, while the remaining 12 cases (20%) showed negative expression. Next, we examined ALOX15B in conventional PTC cell lines and found marked upregulation in KTC-1 and B-CPAP as compared to normal thyroid follicular epithelial cells (Figure 5C). We established stable KTC-1 and B-CPAP cell lines expressing shRNA targeting ALOX15B, and the expression level of ALOX15B was confirmed by Western blotting and RT-qPCR (Figure 5D,5E). CCK-8 assays revealed that ALOX15B silencing significantly suppressed the proliferation of PTC cells (Figure 5F). Furthermore, transwell and wound-healing assays demonstrated that ALOX15B knockdown markedly inhibited the migratory and invasive capacities of PTC cells (Figure 5G,5H). To experimentally validate these observations, we established PTC cell lines with stable ALOX15B overexpression, and transcriptional upregulation was confirmed by Western blotting (Figure 5I). As anticipated, in contrast to the effects of ALOX15B silencing, overexpression of ALOX15B markedly enhanced the proliferation, migration, and invasion of PTC cells (Figure 5J-5L). Taken together, these findings suggest that the in vitro manipulation of ALOX15B expression modulates both the proliferative and metastatic capacity of PTC cells.
Discussion
PTC is generally considered to be an indolent malignant tumor of the endocrine system (28-30). Although most patients with PTC have a favorable prognosis, significant heterogeneity between individuals exists. Recurrence severely impacts patients’ quality of life, and the 5-year postoperative recurrence rate is as high as 30%, posing a substantial burden on healthcare systems and society (31,32). Therefore, elucidating the mechanisms underlying PTC progression and identifying reliable biomarkers for predicting and monitoring disease behavior are urgently needed to improve clinical management.
In this study, we integrated multiomics data and bioinformatic approaches to systematically investigate the molecular mechanisms and immune microenvironment features driving PTC development. By harmonizing six GEO datasets and applying rigorous bioinformatic criteria, we identified 450 DEGs. Among these, seven were further indicated by MR analysis to be causative factors in PTC pathogenesis. Specifically, ALOX15B, TIAM1, and TMC6 were upregulated and may likely function as oncogenes, whereas GPX3, RAP1GAP, JUN, and PAPSS2 were downregulated and may act as tumor suppressors. The consistent expression patterns of these genes in TCGA and GSE33630 datasets reinforce the reliability of our findings.
We also employed the CIBERSORT algorithm to examine the differences in immune cell infiltration patterns between normal thyroid tissue and the tumor microenvironment (TME) of PTC. The relative proportion of γδ T cells was decreased in PTC tissues, whereas resting dendritic cells and resting mast cells were increased. Because γδ T cells can contribute to antitumor immune surveillance and cytokine production, their reduced abundance may suggest weakened immune surveillance in the PTC microenvironment (33-37). Increased resting dendritic cells may reflect an immature or less activated antigen-presenting state, and mast cells have been associated with angiogenesis and tissue remodeling in several tumor contexts (38-40). Nevertheless, these immune findings were derived from bulk transcriptomic deconvolution and should be interpreted as computational associations rather than direct experimental evidence of immune suppression or immune escape. Further validation using multiplex immunohistochemistry, flow cytometry, or single-cell transcriptomics is required to confirm these immune alterations and clarify their functional relevance in PTC.
The three-gene risk model constructed via LASSO Cox regression showed independent prognostic value in the TCGA-PTC cohort. In multivariate Cox regression analysis, the risk score remained significantly associated with PFS after adjustment for clinicopathological variables, and the nomogram integrating the risk score with clinical parameters showed good agreement between predicted and observed outcomes. Time-dependent ROC analysis revealed AUCs of 0.760, 0.700, and 0.708 for 1-, 3-, and 5-year PFS, respectively. Therefore, although this three-gene signature may provide additional molecular information for risk stratification in PTC, it is not sufficient to serve as a standalone clinical decision-making tool at this stage. In addition, although the expression patterns of the candidate genes were validated in an independent GEO dataset, the prognostic model was mainly developed and evaluated using the TCGA-PTC cohort. An independent external cohort with complete survival or PFS information was not available in the present study. Therefore, the prognostic performance of this three-gene signature requires further validation in larger, multicenter, prospective cohorts before it can be translated into routine clinical risk stratification. Furthermore, the high-risk group was characterized by an elevated TMB and a more immunosuppressive microenvironment, evidenced by the reduced infiltration of CD8⁺ T cells and an increase in the abundance of Tregs. These features suggest that the aggressive phenotype associated with high-risk scores may be driven not only by genomic instability but also by immune evasion mechanisms. Notably, BRAF mutations were significantly more prevalent in the high-risk group, whereas NRAS mutations were more common in the low-risk group. This finding aligns with established evidence that BRAF-mutant PTC tends to exhibit more aggressive clinical behavior (41-44), while RAS-driven tumors are generally associated with indolent disease progression (45,46), further supporting the biological relevance of the risk model.
Notably, ALOX15B, which encodes a key enzyme involved in lipid metabolism, has been implicated in various malignancies but remains poorly characterized in PTC (47). Mechanistically, ALOX15B may influence PTC cell behavior through lipid mediator-dependent signaling. As a lipoxygenase involved in polyunsaturated fatty acid metabolism, ALOX15B can participate in the generation of bioactive lipid metabolites, including 15-HETE. Previous studies in other tumor types have linked lipoxygenase-derived lipid mediators to PI3K/Akt signaling, STAT3 activation, extracellular matrix remodeling, and matrix metalloproteinase regulation (48-50). These pathways may partly explain the enhanced proliferative, migratory, and invasive phenotypes observed after ALOX15B overexpression in PTC cells. ALOX15B-derived lipid mediators have also been implicated in inflammatory and immune-regulatory processes in other disease contexts (51). However, these mechanisms remain speculative in the context of PTC and require direct validation through pathway inhibition, transcriptomic profiling, lipidomic analysis, and rescue experiments.
However, several limitations of this study should be acknowledged. First, although the expression patterns of the candidate genes were validated in TCGA and GSE33630 datasets, the prognostic model was mainly developed and evaluated in the TCGA-PTC cohort. An independent external cohort with complete survival or PFS information was not available, and prospective multicenter validation is required before clinical application. Second, immune infiltration was estimated using CIBERSORT based on bulk transcriptomic data; therefore, conclusions regarding immune suppression or immune escape should be interpreted cautiously and validated using direct experimental approaches. Third, although our gain- and loss-of-function experiments support a functional role of ALOX15B in PTC cells, the downstream signaling pathways were not directly examined. Further mechanistic studies are needed to determine whether ALOX15B regulates PTC progression through lipid mediator production, PI3K/Akt or STAT3 signaling, epithelial-mesenchymal transition, or remodeling of the TME. Finally, in vivo validation in animal models will be necessary to confirm the functional role of ALOX15B and the translational value of the proposed prognostic model.
Conclusions
In conclusion, our study identified PTC-associated genes through integrated transcriptomic and MR analyses and established a three-gene signature with moderate prognostic value. The risk score was associated with distinct genomic alterations and computationally inferred immune features. Functional experiments further suggested that ALOX15B promotes PTC cell proliferation, migration, and invasion in vitro. These findings provide potential biomarkers and hypotheses for future mechanistic and translational studies, but further external validation and mechanistic investigation are required.
Acknowledgments
None.
Footnote
Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0267/rc
Data Sharing Statement: Available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0267/dss
Peer Review File: Available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0267/prf
Funding: This work was supported by
Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://gs.amegroups.com/article/view/10.21037/gs-2026-0267/coif). The authors have no conflicts of interest to declare.
Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved
Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.
References
- Zheng G, Chen S, Ma W, et al. Spatial and Single-Cell Transcriptomics Unraveled Spatial Evolution of Papillary Thyroid Cancer. Adv Sci (Weinh) 2025;12:e2404491. [Crossref] [PubMed]
- Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin 2020;70:7-30. [Crossref] [PubMed]
- Castellanos LE, Zafereo ME, Sturgis EM, et al. Pediatric Papillary Thyroid Carcinoma: Outcomes After Surgery Without Adjuvant Radioactive Iodine. J Clin Endocrinol Metab 2025;110:e208-17. [Crossref] [PubMed]
- Donaldson LB, Yan F, Morgan PF, et al. Hobnail variant of papillary thyroid carcinoma: a systematic review and meta-analysis. Endocrine 2021;72:27-39. [Crossref] [PubMed]
- Boucai L, Zafereo M, Cabanillas ME. Thyroid Cancer: A Review. JAMA 2024;331:425-35. [Crossref] [PubMed]
- Ruan X, Huang Y, Zeng Y, et al. The SOX12-YBX1-LDHA signaling axis drives metastasis in papillary thyroid carcinoma. Cell Death Dis 2025;16:474. [Crossref] [PubMed]
- Bischoff LA, Ganly I, Fugazzola L, et al. Molecular Alterations and Comprehensive Clinical Management of Oncocytic Thyroid Carcinoma: A Review and Multidisciplinary 2023 Update. JAMA Otolaryngol Head Neck Surg 2024;150:265-72. [Crossref] [PubMed]
- Su WM, Gu XJ, Dou M, et al. Systematic druggable genome-wide Mendelian randomisation identifies therapeutic targets for Alzheimer's disease. J Neurol Neurosurg Psychiatry 2023;94:954-61. [Crossref] [PubMed]
- King EA, Davis JW, Degner JF. Are drug targets with genetic support twice as likely to be approved? Revised estimates of the impact of genetic support for drug mechanisms on the probability of drug approval. PLoS Genet 2019;15:e1008489.
- Zhou W, Brumpton B, Kabil O, et al. GWAS of thyroid stimulating hormone highlights pleiotropic effects and inverse association with thyroid cancer. Nat Commun 2020;11:3981. [Crossref] [PubMed]
- Hwangbo Y, Park YJ. Genome-Wide Association Studies of Autoimmune Thyroid Diseases, Thyroid Function, and Thyroid Cancer. Endocrinol Metab (Seoul) 2018;33:175-84. [Crossref] [PubMed]
- Zhang Z, Fang T, Chen L, et al. Multi-omics Mendelian randomization integrating GWAS, eQTL, and mQTL data identified genes associated with breast cancer. Am J Cancer Res 2024;14:1433-45. [Crossref] [PubMed]
- Weith M, Beyer A. The next step in Mendelian randomization. Elife 2023;12:e86416. [Crossref] [PubMed]
- Spiga F, Gibson M, Dawson S, et al. Tools for assessing quality and risk of bias in Mendelian randomization studies: a systematic review. Int J Epidemiol 2023;52:227-49. [Crossref] [PubMed]
- Elhage KG, Kranyak A, Jin JQ, et al. Mendelian Randomization Studies in Atopic Dermatitis: A Systematic Review. J Invest Dermatol 2024;144:1022-37. [Crossref] [PubMed]
- Li Y, Sundquist K, Zhang N, et al. Mitochondrial related genome-wide Mendelian randomization identifies putatively causal genes for multiple cancer types. EBioMedicine 2023;88:104432. [Crossref] [PubMed]
- Zhang Y, Wang M, Li Z, et al. An overview of detecting gene-trait associations by integrating GWAS summary statistics and eQTLs. Sci China Life Sci 2024;67:1133-54. [Crossref] [PubMed]
- Jia F, Wei Z, Kong X, et al. Causal Associations Between Lifestyle Habits and Risk of Benign Prostatic Hyperplasia: A Two-Sample Mendelian Randomization Study. J Gerontol A Biol Sci Med Sci 2024;79:glad187. [Crossref] [PubMed]
- Lukas E, van de Weijer M, Bergstedt J, et al. Causal inference in the field of arrhythmia: An introduction to mendelian randomization. Heart Rhythm 2025;22:203-16. [Crossref] [PubMed]
- Li X, Jian J, Zhang A, et al. The role of immune cells and immune related genes in the tumor microenvironment of papillary thyroid cancer and their significance for immunotherapy. Sci Rep 2024;14:18125. [Crossref] [PubMed]
- Li C, Wang P, Dong Z, et al. Single-cell transcriptomics analysis reveals that the tumor-infiltrating B cells determine the indolent fate of papillary thyroid carcinoma. J Exp Clin Cancer Res 2025;44:91. [Crossref] [PubMed]
- Liao Z, Cheng Y, Zhang H, et al. A novel prognostic signature and immune microenvironment characteristics associated with disulfidptosis in papillary thyroid carcinoma based on single-cell RNA sequencing. Front Cell Dev Biol 2023;11:1308352. [Crossref] [PubMed]
- Liddy DD, Zhang Z, Shirlekar K, et al. Impact of age on genomic alterations and the tumor immune microenvironment in papillary thyroid cancer. Endocr Relat Cancer 2024;31:e230341. [Crossref] [PubMed]
- Zheng X, Sun R, Wei T. Immune microenvironment in papillary thyroid carcinoma: roles of immune cells and checkpoints in disease progression and therapeutic implications. Front Immunol 2024;15:1438235. [Crossref] [PubMed]
- Ferkel SAM, Holman EA, Sojwal RS, et al. Tumor-Infiltrating Immune Cells in Colorectal Cancer. Neoplasia 2025;59:101091. [Crossref] [PubMed]
- Xu T, Zhang H, Yang BB, et al. Tumor-infiltrating immune cells state-implications for various breast cancer subtypes. Front Immunol 2025;16:1550003. [Crossref] [PubMed]
- Klobuch S, Seijkens TTP, Schumacher TN, et al. Tumour-infiltrating lymphocyte therapy for patients with advanced-stage melanoma. Nat Rev Clin Oncol 2024;21:173-84. [Crossref] [PubMed]
- Hou ZK, Zhao J, Zhang M, et al. Preoperative Identification of Papillary Thyroid Carcinoma Subtypes and Lymph Node Metastasis via Deep Learning-Assisted Surface-Enhanced Raman Spectroscopy. ACS Nano 2025;19:21807-19. [Crossref] [PubMed]
- Wang L, Zhang L, Ma R, et al. Semaglutide Reprograms Macrophages via the GLP-1R/PPARG/ACSL1 Pathway to Suppress Papillary Thyroid Carcinoma Growth. J Clin Endocrinol Metab 2025;110:2777-89. [Crossref] [PubMed]
- Wang JR, Zafereo ME, Cabanillas ME, et al. The Association Between Thyroid Differentiation Score and Survival Outcomes in Papillary Thyroid Carcinoma. J Clin Endocrinol Metab 2025;110:356-63. [Crossref] [PubMed]
- Papini P, Rossi L, Matrone A, et al. Prophylactic central neck dissection in clinically node-negative papillary thyroid carcinoma: 10-year impact on surgical and oncologic outcomes. Surgery 2025;181:109258. [Crossref] [PubMed]
- Kim BC, Pak SJ, Kwon D, et al. Clinically Significant Central Lymph Node Metastasis is Not Common in Patients with Clinically N0 Papillary Thyroid Carcinoma on Both Ultrasonography and CT. Thyroid 2025;35:415-23. [Crossref] [PubMed]
- Muñoz-Ruiz M, Llorian M, D'Antuono R, et al. IFN-γ-dependent interactions between tissue-intrinsic γδ T cells and tissue-infiltrating CD8 T cells limit allergic contact dermatitis. J Allergy Clin Immunol 2023;152:1520-40. [Crossref] [PubMed]
- Fagundes BO, de-Sousa TR, Victor JR. Gamma-delta (γδ) T cell-derived cytokines (IL-4, IL-17, IFN-γ and IL-10) and their possible implications for atopic dermatitis development. Int J Dermatol 2023;62:443-8. [Crossref] [PubMed]
- Ibusuki A, Kawai K, Nitahara-Takeuchi A, et al. TCR signaling and cellular metabolism regulate the capacity of murine epidermal γδ T cells to rapidly produce IL-13 but not IFN-γ. Front Immunol 2024;15:1361139. [Crossref] [PubMed]
- Li J, Cao Y, Liu Y, et al. Multiomics profiling reveals the benefits of gamma-delta (γδ) T lymphocytes for improving the tumor microenvironment, immunotherapy efficacy and prognosis in cervical cancer. J Immunother Cancer 2024;12:e008355. [Crossref] [PubMed]
- Kang I, Kim Y, Lee HK. γδ T cells as a potential therapeutic agent for glioblastoma. Front Immunol 2023;14:1273986. [Crossref] [PubMed]
- Bhandarkar V, Dinter T, Spranger S. Architects of immunity: How dendritic cells shape CD8(+) T cell fate in cancer. Sci Immunol 2025;10:eadf4726. [Crossref] [PubMed]
- Rodrigues PF, Wu S, Trsan T, et al. Rorγt-positive dendritic cells are required for the induction of peripheral regulatory T cells in response to oral antigens. Cell 2025;188:2720-2737.e22. [Crossref] [PubMed]
- Song X, Jiao J, Qin J, et al. Single-cell transcriptomics reveals the heterogeneity and function of mast cells in human ccRCC. Front Immunol 2024;15:1494025. [Crossref] [PubMed]
- Wen SS, Wu YJ, Wang JY, et al. BRAF(V600E)/p-ERK/p-DRP1(Ser616) Promotes Tumor Progression and Reprogramming of Glucose Metabolism in Papillary Thyroid Cancer. Thyroid 2024;34:1246-59. [Crossref] [PubMed]
- Xie H, Wei B, Shen H, et al. BRAF mutation in papillary thyroid carcinoma (PTC) and its association with clinicopathological features and systemic inflammation response index (SIRI). Am J Transl Res 2018;10:2726-36.
- Li P, Liu Y, Wei T, et al. Effect and Interactions of BRAF on Lymph Node Metastasis in Papillary Thyroid Carcinoma With Hashimoto Thyroiditis. J Clin Endocrinol Metab 2024;109:944-54. [Crossref] [PubMed]
- Perampalam S, Wu K, Gild M, et al. The Association between Lymphocytic Thyroiditis and Papillary Thyroid Cancer Harboring Mutant BRAF: A Systematic Review and Meta-Analysis. Thyroid 2024;34:1082-93. [Crossref] [PubMed]
- Kim YH, Yoon SJ, Kim M, et al. Integrative Multi-omics Analysis Reveals Different Metabolic Phenotypes Based on Molecular Characteristics in Thyroid Cancer. Clin Cancer Res 2024;30:883-94. [Crossref] [PubMed]
- Moposita Molina FD, Agnes G, Zandoná MR, et al. Prevalence and diagnostic reliability of BRAF, RAS mutations, and RET/PTC rearrangements in a Latin American public health service population with thyroid nodular disease. PLoS One 2025;20:e0329407. [Crossref] [PubMed]
- Fang W, Liu J, Zhang F, et al. A novel cholesterol metabolism-related ferroptosis pathway in hepatocellular carcinoma. Discov Oncol 2024;15:7. [Crossref] [PubMed]
- Soumya SJ, Binu S, Helen A, et al. 15(S)-HETE-induced angiogenesis in adipose tissue is mediated through activation of PI3K/Akt/mTOR signaling pathway. Biochem Cell Biol 2013;91:498-505. [Crossref] [PubMed]
- Yang L, Ma C, Zhang L, et al. 15-Lipoxygenase-2/15(S)-hydroxyeicosatetraenoic acid regulates cell proliferation and metastasis via the STAT3 pathway in lung adenocarcinoma. Prostaglandins Other Lipid Mediat 2018;138:31-40. [Crossref] [PubMed]
- Souza FDC, Ferreira MT, Colquhoun A. Influence of Lipoxygenase Inhibition on Glioblastoma Cell Biology. Int J Mol Sci 2020;21:8395. [Crossref] [PubMed]
- Snodgrass RG, Zezina E, Namgaladze D, et al. A Novel Function for 15-Lipoxygenases in Cholesterol Homeostasis and CCL17 Production in Human Macrophages. Front Immunol 2018;9:1906. [Crossref] [PubMed]

