An integrated multiomics-based risk model for predicting the prognosis and immune microenvironment status in patients with papillary thyroid carcinoma
Original Article

An integrated multiomics-based risk model for predicting the prognosis and immune microenvironment status in patients with papillary thyroid carcinoma

Nizhen Xu1 ORCID logo, Zheng Shi2, Mingjie Zheng3, Deguang Zhang1, Qiqi He4, Chu Zhu5, Zhenlei Zhang6,7, Gaofei He1, Junjie Chu1, Jinxi Jiang1, Xiaoxiao Lu1, Xiujun Cai8,9,10,11

1Department of Thyroid and Head & Neck Surgery, Institute of Micro-Invasive Surgery of Zhejiang University, Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou, China; 2School of Medicine, Shaoxing University, Shaoxing, China; 3Department of Cardiology, The First Affiliated Hospital, School of Medicine, Zhejiang University, Hangzhou, China; 4Department of Nursing, Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou, China; 5Department of Radiation Oncology, Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou, China; 6Department of Orthopaedic Surgery, Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou, China; 7Zhejiang Key Laboratory of Mechanism Research and Precision Repair of Orthopaedic Trauma and Aging Diseases, Zhejiang University School of Medicine, Hangzhou, China; 8Key Laboratory of Laparoscopic Technology of Zhejiang Province, Department of General Surgery, Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou, China; 9Zhejiang Minimal Invasive Diagnosis and Treatment Technology Research Center of Severe Hepatobiliary Disease, Zhejiang Research and Development Engineering Laboratory of Minimally Invasive Technology and Equipment, Hangzhou, China; 10Zhejiang University Cancer Center, Hangzhou, China; 11Liangzhu Laboratory, Zhejiang University Medical Center, Hangzhou, China

Contributions: (I) Conception and design: X Cai, N Xu; (II) Administrative support: X Cai; (III) Provision of study materials or patients: X Cai, N Xu, M Zheng, D Zhang; (IV) Collection and assembly of data: N Xu, Z Shi, Q He, C Zhu, Z Zhang, G He; (V) Data analysis and interpretation: N Xu, Z Shi, M Zheng, J Chu, J Jiang, X Lu; (VI) Manuscript writing: All authors; (VII) Final approval of manuscript: All authors.

Correspondence to: Xiujun Cai, MD. Key Laboratory of Laparoscopic Technology of Zhejiang Province, Department of General Surgery, Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, No. 3 East Qingchun Road, Shangcheng District, Hangzhou 310016, China; Zhejiang Minimal Invasive Diagnosis and Treatment Technology Research Center of Severe Hepatobiliary Disease, Zhejiang Research and Development Engineering Laboratory of Minimally Invasive Technology and Equipment, Hangzhou, China; Zhejiang University Cancer Center, Hangzhou, China; Liangzhu Laboratory, Zhejiang University Medical Center, Hangzhou, China. Email: srrsh_cxj@zju.edu.cn.

Background: Papillary thyroid carcinoma (PTC) is the most common type of primary endocrine malignancy. The tumor immune microenvironment (TIME) and genetic alterations play crucial roles in the progression of PTC. With the advance in research, there has been a heightened focus on re-evaluating molecular targeted therapies and identifying novel targets through molecular biology-based approaches. This study aimed to identify robust prognostic biomarkers and to construct a reliable risk model for patients with PTC by integrating multiomics data.

Methods: RNA-sequencing (RNA-seq) data from six Gene Expression Omnibus (GEO) datasets and genome-wide association study (GWAS) data were integrated to identify differentially expressed genes (DEGs) and facilitate Mendelian randomization (MR) analysis for causal inference. The CIBERSORT algorithm was employed to evaluate immune cell infiltration. A prognostic risk model was constructed via least absolute shrinkage and selection operator (LASSO) Cox regression and validated in The Cancer Genome Atlas (TCGA) and GEO cohorts. Functional experiments, including Cell Counting Kit-8 (CCK-8), transwell, and wound-healing assays, were conducted to investigate the role of the key gene ALOX15B in PTC cells.

Results: We identified seven PTC-associated genes (ALOX15B, TIAM1, TMC6, GPX3, RAP1GAP, JUN, and PAPSS2) through expression quantitative trait loci and MR analysis. A robust three-gene risk model (ALOX15B, RAP1GAP, and JUN) was established. Patients in the high-risk group exhibited significantly poorer progression-free survival (PFS), which was confirmed to be an independent prognostic factor by multivariate analysis [hazard ratio =1.355, 95% confidence interval (CI): 1.052–1.746; P=0.02]. The high-risk group was characterized by an immunosuppressive TIME—including decreased CD8+ T-cell and increased regulatory T-cell abundance—higher tumor mutational burden, and a higher frequency of BRAF mutations. Conversely, the low-risk group showed a higher prevalence of NRAS mutations. In addition, functional assays in vitro revealed that ALOX15B markedly enhanced the proliferative, migratory, and invasive capacities of PTC cells.

Conclusions: A novel three-gene signature was developed and was demonstrated to be an independent prognostic indicator for patients with PTC. This model effectively reflects the intrinsic characteristics of the TIME and genomic instability and can inform risk stratification and therapeutic targeting. Due to its oncogenic activity, ALOX15B may represent a viable therapeutic target in PTC.

Keywords: Thyroid cancer; biomarkers; tumor immune microenvironment (TIME); Mendelian randomization (MR)


Submitted May 02, 2026. Accepted for publication May 28, 2026. Published online Jun 26, 2026.

doi: 10.21037/gs-2026-0267


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.

Figure 1 Identification of candidate genes associated with PTC through integrated transcriptome and MR analysis. (A,B) PCA of six PTC transcriptome datasets (GSE129562, GSE29265, GSE3678, GSE50901, GSE6004, and GSE97001) before (A) and after (B) batch effect correction, showing successful normalization of the data. (C) Volcano plot of DEGs from the integrated GEO datasets. Significantly upregulated (red) and downregulated (blue) genes are highlighted. (D) Heatmap of the top 50 upregulated and 50 downregulated DEGs in 37 paired PTC tissues and matched normal tissues. (E,F) Venn diagram of the intersection between MR-derived PTC-associated genes (n=189) and DEGs (n=450), yielding seven candidate genes. (G) MR forest plot displaying the effects of the seven candidate genes on PTC risk. (H) Chromosomal locations of the seven candidate genes. CI, confidence interval; DEG, differentially expressed gene; FC, fold change; GEO, Gene Expression Omnibus; MR, Mendelian randomization; NSNP, number of single-nucleotide polymorphisms; OR, odds ratio; PC, principal component; PCA, principal component analysis; PTC, papillary thyroid carcinoma.

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

Figure 2 Immune microenvironment remodeling and derivation of a risk model in PTC. (A) Relative abundance of immune cells in PTC and normal tissues (CIBERSORT). (B) Box plot showing the composition of 22 immune cell types in PTC and normal tissues. (C) Correlation matrix depicting the relationships among the 22 immune cell types (red: positive correlations; blue: negative correlations). (D,E) Validation of the seven candidate genes expression in the independent (D) GEO (GSE33630) and (E) TCGA datasets. Box plots show the expression levels in PTC and normal tissues. (F) Correlation heatmap analyzing the interactions among the seven candidate genes. (G) The tuning parameters (log λ) of the seven candidate genes were selected to cross-verify the error curve. According to the minimal criterion and 1-SE criterion, perpendicular imaginary lines were drawn at the optimal value. (H) The LASSO coefficient profile of the seven candidate genes. Data are presented as the mean ± SD. *, P<0.05; ***, P<0.001. 1-SE, one standard error rule; GEO, Gene Expression Omnibus; LASSO, least absolute shrinkage and selection operator; PTC, papillary thyroid carcinoma; SD, standard deviation; TCGA, The Cancer Genome Atlas.

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.

Figure 3 The risk model as an independent prognostic indicator for PTC in the TCGA. (A) Distribution of the risk model-based risk score. (B) Kaplan-Meier survival curves comparing the PFS between the low- and high-risk groups. (C) Time-dependent ROC curves analyzing the predictive accuracy of the risk model for 1-, 3-, and 5-year PFS. (D,E) Forest plots from (D) univariate and (E) multivariate Cox regression analyses of the risk score and clinical parameters for PFS. (F) Time-dependent C-index comparing the prognostic performance of the risk score against other clinical variables. (G) The nomogram integrating the risk model and clinical features for predicting 1-, 3-, and 5-year PFS. (H) The calibration curves for the nomogram the comparing predicted and observed PFS at 1, 3, and 5 years. ***, P<0.001. AUC, area under the curve; C-index, concordance index; CI, confidence interval; PFS, progression-free survival; PTC, papillary thyroid carcinoma; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

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.

Figure 4 Characterization of the tumor immune microenvironment based on the risk model in the TCGA. (A) TME difference in the low- and high-risk groups. (B) Relative abundance of immune cells across risk groups. (C) Box plot showing the composition of 22 immune cell types in the low- and high-risk groups. (D) TMB comparison between the risk groups. (E,F) Mutational landscape of the top 15 driver genes in the high-risk (E) and low-risk (F) groups. (G) Kaplan-Meier analysis of PFS between the low- and high-TMB groups. (H) PFS of patients classified according to TMB status and risk model from the Kaplan-Meier curve analysis. *, P<0.05; ***, P<0.001. H, high; L, low; PFS, progression-free survival; TCGA, The Cancer Genome Atlas; TMB, tumor mutation burden; TME, tumor microenvironment.

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.

Figure 5 Functional characterization of ALOX15B in PTC progression. (A) RT-qPCR validation of upregulated ALOX15B in PTC tissues compared with adjacent normal tissues (n=12 pairs). (B) Representative IHC images of ALOX15B expression, showing positive staining in PTC tissue and negative staining in normal thyroid tissue. (C) ALOX15B expression in normal thyroid follicular epithelial cells (Nthy-ori 3-1) and PTC cell lines (K1, KTC-1, B-CPAP, and TPC-1). (D,E) Validation of ALOX15B knockdown efficiency in KTC-1 and B-CPAP cells by Western blotting (D) and RT-qPCR (E). (F) Cell proliferation assessed by CCK-8 assay in PTC cells after ALOX15B knockdown. (G) Wound-healing assay showing the migration ability of PTC cells after ALOX15B knockdown. Representative images at 0 and 24 h are shown. The relative migration distance was quantified. (H) Representative images (left) and quantification (right) of cell migration assessed by transwell assays following ALOX15B knockdown (0.1% crystal violet stained). (I) Validation of ALOX15B overexpression efficiency in PTC cells by Western blotting. (J) Cell proliferation assessed by CCK-8 assay in PTC cells with ALOX15B overexpression. (K) Wound-healing assay showing the migration ability of PTC cells overexpressing ALOX15B. Representative images at 0 and 24 h are shown. The relative migration distance was quantified. (L) Representative images (above) and quantification (below) of cell migration assessed by transwell assays after ALOX15B overexpression (0.1% crystal violet stained). Data are presented as the mean ± SD. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001. CCK-8, Cell Counting Kit-8; H&E, hematoxylin and eosin; IHC, immunohistochemistry; OD, optical density; PTC, papillary thyroid carcinoma; RT-qPCR, reverse-transcription quantitative polymerase chain reaction.

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 the Natural Science Foundation of Zhejiang Province (No. LQN25H160022), the National Natural Science Foundation of China (No. 82500383), and the Zhejiang Province Medicine and Health Technology Plan Project (No. 2025KY885).

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. This study was approved by the Ethics Committee of Sir Run Run Shaw Hospital (Approval No. 20260302-065), with written informed consent obtained from each patient. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

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

  1. 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]
  2. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin 2020;70:7-30. [Crossref] [PubMed]
  3. 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]
  4. 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]
  5. Boucai L, Zafereo M, Cabanillas ME. Thyroid Cancer: A Review. JAMA 2024;331:425-35. [Crossref] [PubMed]
  6. 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]
  7. 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]
  8. 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]
  9. 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.
  10. 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]
  11. 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]
  12. 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]
  13. Weith M, Beyer A. The next step in Mendelian randomization. Elife 2023;12:e86416. [Crossref] [PubMed]
  14. 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]
  15. 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]
  16. 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]
  17. 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]
  18. 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]
  19. 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]
  20. 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]
  21. 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]
  22. 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]
  23. 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]
  24. 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]
  25. Ferkel SAM, Holman EA, Sojwal RS, et al. Tumor-Infiltrating Immune Cells in Colorectal Cancer. Neoplasia 2025;59:101091. [Crossref] [PubMed]
  26. 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]
  27. 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]
  28. 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]
  29. 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]
  30. 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]
  31. 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]
  32. 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]
  33. 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]
  34. 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]
  35. 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]
  36. 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]
  37. Kang I, Kim Y, Lee HK. γδ T cells as a potential therapeutic agent for glioblastoma. Front Immunol 2023;14:1273986. [Crossref] [PubMed]
  38. 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]
  39. 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]
  40. 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]
  41. 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]
  42. 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.
  43. 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]
  44. 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]
  45. 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]
  46. 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]
  47. 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]
  48. 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]
  49. 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]
  50. Souza FDC, Ferreira MT, Colquhoun A. Influence of Lipoxygenase Inhibition on Glioblastoma Cell Biology. Int J Mol Sci 2020;21:8395. [Crossref] [PubMed]
  51. 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]
Cite this article as: Xu N, Shi Z, Zheng M, Zhang D, He Q, Zhu C, Zhang Z, He G, Chu J, Jiang J, Lu X, Cai X. An integrated multiomics-based risk model for predicting the prognosis and immune microenvironment status in patients with papillary thyroid carcinoma. Gland Surg 2026;15(6):173. doi: 10.21037/gs-2026-0267

Download Citation