Biomarker Identification and Risk Prediction Model Development for Differentiated Thyroid Carcinoma Lung Metastasis Based on Primary Lesion Proteomics

Thyroid cancer stands as the most prevalent malignant tumor within the endocrine system, with differentiated thyroid carcinoma (DTC), encompassing papillary thyroid carcinoma and follicular thyroid carcinoma, constituting more than 90% of cases (1). The annual incidence of DTC is on the rise (2), with approximately 15% of patients experiencing distant metastasis (3). Notably, the lungs represent the primary site of metastasis, profoundly impacting patent prognosis and survival (4). Clinical management of these cases poses considerable challenges, underscoring the need to develop effective methods for risk stratification and metastasis prediction.

Currently, the bulk of studies on risk factors for lung metastasis in DTC (LMDTC) rely on the analysis of patient’s clinical characteristics. For instance, Sugino and colleagues (5) conducted a retrospective study involving 171 patients with DTC under the age of 19, revealing gender, lymph node metastasis, extrathyroidal invasion, and the number of lymph node metastases as independent risk factors for lung metastasis occurrence. Although there have been reports on biomarker studies for LMDTC in recent years, most derive conclusions from mining databases such as Gene Expression Omnibus and The Cancer Genome Atlas. A few studies have employed genomic and transcriptomic analysis techniques to study small sample sizes. For example, some studies have identified low miR-199a-3p expression and miR-424-5p overexpression as linked to distant DTC metastasis (6, 7). MiR-222 has been found to promote papillary thyroid carcinoma invasion and metastasis by targeting PPP2R2A (8). Wang and colleagues (9) developed a combined prediction model integrating gene expression profiles and clinical pathologic features, identifying DLX5, COX6B2, and LYPD1 as key genes for predicting distant DTC metastasis. However, the tissue sample sizes in these studies did not exceed 20 cases.

Proteins, serving as the agents of life processes, reflect the ongoing events within organisms, representing the outcome of gene expression. Therefore, investigating proteins is of paramount importance. Proteomics, which centers on the study of the proteome, employs sensitive and specific proteomic techniques to explore proteins linked to tumor occurrence and development (10–12). Although proteomic approaches are widely used in tumor research, reports on proteomic studies related to DTC metastasis remain scarce.

To identify potential biomarkers for predicting LMDTC risk at the protein level, we conducted 4D label–free proteomic analysis and dual validation using parallel reaction monitoring (PRM) and IHC on formalin-fixed, paraffin-embedded (FFPE) samples obtained from primary lesions of patients with LMDTC and non-lung metastatic DTC (NMDTC). Our analysis revealed eight hub proteins as potential biomarkers for predicting LMDTC risk. Furthermore, we augmented the sample size and developed single and combined LMDTC risk prediction models based on clinical features and hub protein expression. Internal validation using the bootstrap method with 2,000 resamples demonstrated superior performance of the combined prediction model over the single prediction model. Subsequently, we employed the combined prediction model to construct a nomogram. This study aims to aid in the risk stratification of lung metastasis in patients with DTC and provide a pivotal foundation for developing clinical personalized treatment plans and follow-up strategies for high-risk patients.

Primary lesion FFPE samples were collected from 181 patients with DTC sourced from the Yunnan Cancer Hospital biobank (Supplementary Table S1). Initially, a subset of patients was randomly selected for proteomic analysis and small sample PRM and IHC validation to identify potential biomarkers for predicting LMDTC, comprising 44 cases of LMDTC and 25 cases of NMDTC (Supplementary Table S2). Subsequently, the sample size was expanded to construct LMDTC risk prediction models, encompassing 81 cases of LMDTC and 100 cases of NMDTC. The clinical diagnostic criteria for LMDTC were adhered to (13, 14), including (i) abnormal findings on iodine 131 whole-body imaging, (ii) elevated serum thyroglobulin (TG) levels with supporting evidence from chest CT or PET/CT, or (iii) confirmation of lung lesions as metastatic DTC via pathology after needle biopsy or surgical resection. Patient clinical characteristics such as gender, age, preoperative laboratory tests, and postoperative pathology were also recorded. Written informed consent was obtained from all patients before data collection, and the study adhered to ethical guidelines outlined in the Declaration of Helsinki, approved by the Ethics Committee of Yunnan Cancer Hospital (ID: KYCS2023-094).

FFPE samples were lysed using an appropriate amount of SDT lysate, followed by protein quantification using the bicinchoninic acid method. Samples underwent SDS-PAGE and FASP enzymatic hydrolysis, with desalting of peptide fragments accomplished using a C18 cartridge. Following lyophilization of peptide fragments, they were reconstituted with 40 μL of 0.1% formic acid solution and quantified (OD280). The prepared samples then underwent 4D label–free proteomics analysis and PRM verification.

Samples were analyzed using a timsTOF Pro mass spectrometer (Bruker, Bremen, Germany) equipped with a CaptiveSpray ion source. Database searches were conducted using MaxQuant software (V1.6.17.0; RRID: SCR_014485), employing the LFQ algorithm for label-free quantitative analysis. For qualitative matching, MaxQuant utilized the FDR principle with a set FDR (accuracy of peptide and protein level identification) of 1%, requiring the identified protein to contain at least one unique peptide.

Differential analysis of quantitative results involved screening data with at least half of the repeated experimental data in the sample group that were non-null values [FC_{LMDTC/NMDTC,P} = Mean (R_LMDTC,P)/Mean (R_NMDTC,P) (R: relative protein amount, P: protein)]. To amplify protein differences, fold change (FC) had to be log base 2 transformed. Proteins meeting the screening criteria, i.e., expression difference >1.2-fold and t test <0.05, were classified as differentially expressed proteins (DEPs).

Protein bioinformatics analysis included Gene Ontology (GO; RRID:SCR_002811) and Kyoto Encyclopedia of Genes and Genomes (KEGG; RRID:SCR_012773) pathway annotation and enrichment analyses, along with localization prediction analyses of DEPs using WoLF PSORT (wolfpsort.hgc.jp; RRID:SCR_002472). Hierarchical clustering heatmaps were generated using Matplotlib software (MatPlotLib; RRID:SCR_008624) to classify both the samples and protein expression dimensions simultaneously. ID was used to establish protein–protein interaction (PPI) networks with direct or indirect interactions with the target proteins in the STRING (string-db.org; RRID:SCR_005223) database. AnyChart software (V8.11.0.1934) was used to generate the PPI network analysis results. Weighted gene co-expression network analysis (WGCNA) was conducted using the R package WGCNA to screen for hub proteins, which are then used to construct gene co-expression networks among different samples based on the interconnectivity of genes within modules and the association between modules and phenotypes.

PRM validation was performed for candidate hub proteins identified from proteomic analysis using a nanoliter flow NanoElute system for separation, followed by mass spectrometry analysis employing the PRM-PASEF mode of the timsTOF Pro mass spectrometer. Data analysis was conducted using SpectroDive 11.0. Supplementary Table S3 details the parameter settings of the peptides.

TNMplot (TNMplot.com) is a survival analysis tool rooted in tumor–node–metastasis staging, utilizing data sourced from public databases such as The Cancer Genome Atlas and Gene Expression Omnibus. Currently, tumor–node–metastasis staging and corresponding survival data for more than 30 different cancer types are available. We employed the TNMplot database to align and analyze the selected DEPs. Subsequently, receiver operating characteristic (ROC) analysis of DEPs was conducted using MedCalc Statistics software (V20.1.0.0; RRID:SCR_015044). The predictive performance of DEPs in determining the risk of LMDTC was evaluated via the AUC.

Tissue sections underwent deparaffinization and hydration, followed by antigen retrieval and incubation with 3% hydrogen peroxide in darkness. After blocking, primary Abs, including SUCLG1 (Proteintech Cat# 14923-1-AP; RRID: AB_2197177), SUCLG2 (Proteintech Cat# 14240-1-AP; RRID: AB_2197187), DLAT (Proteintech Cat# 13426-1-AP; RRID: AB_2091774), IDH3B (Abcam Cat# ab135444; RRID: AB_2819360), ACSF2 (Proteintech Cat# 16140-1-AP; RRID: AB_2878222), ACO2 (Abcam Cat# ab129069; RRID: AB_11144142), CYCS (Proteintech Cat# 10993-1-AP; RRID: AB_2090467), and VDAC2 (Proteintech Cat# 11663-1-AP; RRID: AB_2304144), were added and incubated overnight at 4°C. Following incubation with goat anti-rabbit IgG (Fitzgerald Industries International Cat# 41R-1079; RRID: AB_10814969) at room temperature for 1 hour, sections were stained with DAB chromogen and hematoxylin and then dehydrated and mounted.

Patients with DTC meeting the inclusion criteria were enrolled, and their clinical characteristics were collected. Semi-quantitative analysis of eight hub proteins using IHC was performed (Supplementary Table S4). Data analysis utilized R4.3.1 software, with each experiment repeated at least thrice. Categorical data were presented as percentages (%) and compared using the χ² test, whereas continuous data were expressed as mean ± SD and compared using t test and one-way ANOVA. A P-value <0.05 was considered statistically significant.

Univariate and multivariate analyses were employed to explore the correlation between clinical characteristics, the eight hub proteins, and LMDTC. Variables with a P-value <0.05 in univariate analyses were included in a multivariate stepwise backward logistic regression to fit three risk prediction models for LMDTC based on clinical features, hub protein expression, and a combination of both. Subsequently, the bootstrap method with 2,000 resamples was utilized for internal validation of the three models. Model performance was primarily evaluated using a confusion matrix, ROC analysis, calibration curve, DCA curve, as well as NRI and IDI, with the best-performing model used to construct a nomogram.

The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (ProteomeXchange; RRID:SCR_004055; https://proteomecentral.proteomexchange.org) via the iProX partner repository with the dataset identifier PXD048443.

Proteomic analysis between the LMDTC and NMDTC groups revealed 418 significantly different proteins.

Mass spectrometry analysis confirmed 43,947 peptide sequences, with our standardization processing pipeline identifying 5,281 proteins (Supplementary Fig. S1). Principal component analysis demonstrated sample overlap across different groups, indicating a high degree of similarity between these samples (Supplementary Fig. S2). Pearson’s correlation analysis further supported this high degree of correlation between the two groups (Fig. 1A). Differential expression analysis yielded 418 DEPs, comprising 186 upregulated and 232 downregulated proteins (Fig. 1B). Particularly, STING1, HMGN5, TESC, PBXIP1, MXRA7, and VWA1 were significantly upregulated, whereas SEC11C, RNF213, HCLS1, HINT2, RAC2, and SOD2 were significantly downregulated in the LMDTC group compared with the NMDTC group (Fig. 1C).

GO and KEGG enrichment analyses revealed that most of the DEPs were located in the mitochondrial matrix, which were mainly involved in various metabolic pathways, the tricarboxylic acid (TCA) cycle, oxidative phosphorylation, and mitochondrial ATP synthesis.

Unsupervised clustering analysis indicated significant differences in protein expression between the LMDTC and NMDTC groups (Supplementary Fig. S3). We also performed GO and KEGG enrichment analyses of all DEPs. GO enrichment analysis indicated involvement in the AA catabolic process, the TCA cycle, mitochondrial ATP synthesis, and fatty acid oxidation. The GO molecular function enrichment results revealed enrichment in the extracellular matrix, SUMO activating enzyme, propionyl-CoA carboxylase, and several enzymes in ATP synthesis. In terms of GO cellular components, the majority of DEPs were localized in the mitochondrial matrix, lysosomal lumen, mitochondrial nucleoid, and other cellular structures related to energy metabolism (Fig. 1D). KEGG pathway enrichment results further validated the role of these DEPs in multiple AA metabolism processes, the TCA cycle, oxidative phosphorylation, and fatty acid metabolism (Fig. 1E). The top 50 DEPs with the largest difference in upregulation and downregulation were selected as target proteins for direct PPI network analysis. Notably, ATP5J, NDUFA6, COX7A2, and CYCS proteins were identified as the hub nodes in the network (Fig. 1F).

Thirty candidate hub proteins were preliminarily screened out through a comprehensive analysis of important factors such as FC, pathway enrichment, correlation degree of protein interaction, and researchers’ attention.

We performed initial screening for candidate hub proteins through five methods. (i) According to the FC of differential expression and combined with relevant literature, we selected the top three upregulated (STING1, HMGN5, and PBXIP1) and top four downregulated proteins (RAC2, HINT2, HCLS1, and RNF213; Fig. 1B and C). (ii) First, we performed KEGG pathway enrichment analysis for 418 DEPs, and the top 20 enriched pathways were used to construct a pathway–gene network and ranked based on the degree (Supplementary Table S5). There were 25 DEPs in the network and were downregulated in the LMDTC group. Then, we constructed a PPI network using the 418 DEPs and filtered proteins with degree higher than 40 in the network (Supplementary Table S6). There were eight shared proteins in the set of 25 DEPs in the top 20 KEGG pathway and the set of proteins with degree >40 (Fig. 1G). (iii) We used WGCNA to construct a gene co-expression network of proteomics data (Supplementary Fig. S4). The results showed that proteins in the blue module significantly correlated with tumor metastasis (Fig. 1H), and 116 hub proteins in the blue module were screened according to MM > 0.8 and GS > 0.2 (Fig. 1I). The Venn diagram of the 116 hub proteins obtained by WGCNA, 25 DEPs in the top 20 KEGG pathways, and proteins with degree >40 in the PPI network showed there were five shared proteins (Fig. 1J). Based on ROC analysis, VDAC1, VDAC2, SDHA, and CYCS could better distinguish metastasis and non-metastasis (AUC > 0.70, P < 0.05; Fig. 1K). (iv) According to the GO and KEGG enrichment analyses, DEPs were significantly enriched in the TCA cycle and were mostly distributed in mitochondria. We speculate that LMDTC is correlated to the TCA cycle, thus we selected 13 proteins enriched in the TCA cycle (see “KEGG Analysis Table of Differential Proteins” in the Supplementary Material). (v) In addition, our group has also carried out research on ferroptosis and tumor metastasis, and we are also interested in whether ferroptosis is involved in LMDTC. We showed the overlap of the 418 DEPs and ferroptosis-related proteins (acquired from the FerrDb V2 database, http://www.zhounan.org/ferrdb/current/), and there were nine shared DEPs (Fig. 1L). Using the above five methods, we try to take into account the important factors for screening hub proteins such as differential change fold, pathway enrichment, protein interaction degree, and researchers’ attentions, so as to screen hub proteins more comprehensively and objectively. Finally, a total of 30 candidate hub proteins were used for the following analysis.

Eight hub proteins, namely, SUCLG1/2, DLAT, IDH3B, ACSF2, ACO2, CYCS, and VDAC2, were identified using the dual validation of PRM and IHC, highlighting their potential as biomarkers for predicting LMDTC.

To further verify the expression of candidate hub proteins, we collected FFPE samples of primary lesions from eight patients with LMDTC and five patients with NMDTC for PRM experiments. The experimental results revealed that 17 of the 30 candidate hub proteins were consistent with the proteomic analysis and were statistically significant (Supplementary Table S7). These 17 proteins were downregulated in the LMDTC (Fig. 2A and B). Similarly, GO enrichment analysis revealed that these 17 hub proteins were mainly involved in pyruvate metabolism and the TCA cycle and correlated to cell respiration and energy metabolism (Fig. 2C). KEGG enrichment analysis further revealed the correlation of these hub proteins to the TCA cycle, cell respiration, and material metabolism (Fig. 2D). Additionally, we used the STRING database to analyze the interactions of these 17 hub proteins. The PPI network demonstrated the interaction of the 15 proteins with each other, except for TIMM9 and HINT2 (Fig. 2E).

The screened 15 hub proteins were analyzed using the TNMplot database (TNMplot.com), revealing eight hub proteins that were consistent with the expression trend in the database (gradually decreasing with the occurrence and development of tumors) and statistically significant: SUCLG1/2, DLAT, IDH3B, ACSF2, ACO2, CYCS, and VDAC2 (Fig. 2F). We performed ROC analysis of these eight hub proteins, revealing high AUC values (>0.90) between the LMDTC group and the NMDTC group (Fig. 2G). Among them, IDH3B and ACSF2 had the highest AUC of 0.975, and the AUC for the remaining six proteins ranged from 0.925 to 0.952. Therefore, the eight proteins were considered to possess a high predictive power for the risk of LMDTC.

We collected 10 pairs of primary lesion tissues from LMDTC and NMDTC and verified the expression of the eight hub proteins using IHC. The expression of the eight hub proteins was significantly lower in the LMDTC group than that in the NMDTC group (Fig. 2H and I), which was consistent with proteomic analysis and PRM verification. Furthermore, proteomic analysis and dual-verification by PRM/IHC highlighted their potential as biomarkers for predicting LMDTC.

The performance of the combined prediction model was better than the single prediction model, suggesting that larger primary tumor size, higher TG levels, invasion of surrounding tissues, N1b staging, and lower expression of SUCLG1, DLAT, CYCS, and ACSF2 could increase the risk of LMDTC.

We collected clinical characteristics and semi-quantitative IHC results of the eight hub proteins from an expanded sample of 81 patients with LMDTC and 100 patients with NMDTC. The baseline data of all patients are presented in Table 1. Protein expressions were categorized as negative (−), weakly positive (+), positive (++), and strongly positive (+++; Supplementary Fig. S5). Logistic regression analyses on the risk factors of lung metastasis in patients with DTC were performed (Table 2). Univariable logistic regression analysis revealed that multifocal, tumor size, TG, adjacent tissue infiltration, T stage, N stage, SUCLG1/2, DLAT, IDH3B, ACSF2, ACO2, CYCS, and VDAC2 were significantly associated with LMDTC (P < 0.05). Multivariable logistic regression analysis further revealed that tumor size, TG, SUCLG1, and DLAT were independent risk factors for the occurrence of lung metastasis in patients with DTC. Specifically, the larger the tumor diameter and the higher the TG level, the lower the expression of SUCLG1 and DLAT and the higher the risk of lung metastasis in patients with DTC. Subsequently, we incorporated the six clinically significant features identified in the univariate analysis, the eight hub proteins, and the combination of the two (14 variables in total) into a stepwise backward logistic regression to fit three LMDTC risk prediction models.

The combined prediction model demonstrated the highest true positive rate and true negative rate (Fig. 3A), with an accuracy of 0.958, a sensitivity of 0.920, a specificity of 0.952, and an F1 score of 0.939, all of which were superior to the single prediction model (based on clinical characteristics analysis or the expression of hub proteins, Supplementary Fig. S6A and S6B). These findings suggest that the combined prediction model can better differentiate whether patients with DTC will experience lung metastasis. The AUC values of the three models were 0.986, 0.937, and 0.953, respectively (Fig. 3B; Supplementary Fig. S6C–S6E), underscoring the better ability of the combined model to differentiate whether patients with DTC will experience lung metastasis compared with the single prediction model. Subsequent calibration curve analysis demonstrated good consistency between the predicted metastasis results of the three models and the actual occurrence of lung metastasis (Fig. 3C; Supplementary Fig. S6F–S6G), with the combined prediction model exhibiting lower Brier scores (Fig. 3D). Decision curve analysis indicated the highest clinical benefit with the combined prediction model (Fig. 3E). Additionally, combined NRI and IDI (Supplementary Table S8) suggested superior predictive ability of the combined prediction model over the single prediction model. A nomogram based on the combined prediction model (Fig. 3F) highlighted factors influencing the risk of lung metastasis in patients with DTC, such as large primary tumor size, high TG levels, surrounding tissue invasion, N1b staging and lower SUCLG1, DLAT, CYCS, and ACSF2 expression.

The prognosis of DTC is often good, but the overall survival rate of patients is greatly reduced when distant metastasis occurs. Approximately 15% of DTC cases develop distant metastasis. Due to the high hemodynamics, the most common organ for distant metastasis of DTC is the lung (15), accounting for more than 70% of distant metastatic DTC cases (3), underscoring its diagnostic and therapeutic challenges in clinical work. Evaluating the risk of metastasis and improving the prognosis through the identification of high-efficiency biomarkers and the establishment of predictive models is vital in the domain of LMDTC. This study identified 418 DEPs, with significant differences between LMDTC and NMDTC groups, using the proteomic analysis of DTC primary lesions. GO and KEGG enrichment analyses revealed that the majority of DEPs were located in the mitochondrial matrix, which was primarily involved in various metabolic pathways, the TCA cycle, oxidative phosphorylation, and mitochondrial ATP synthesis. A comprehensive analysis of important factors such as FC, pathway enrichment, correlation degree of protein interaction, and researchers’ attention revealed 30 candidate hub proteins via a preliminary screening. Subsequently, eight hub proteins, namely SUCLG1/2, DLAT, IDH3B, ACSF2, ACO2, CYCS, and VDAC2, were identified using the dual validation of PRM and IHC. Moreover, these findings highlight the potential of these hub proteins as biomarkers for predicting LMDTC (Supplementary Table S9).

Furthermore, a majority of the potential biomarkers were enzymes vital to pyruvate metabolism and the TCA cycle pathway, which is the key metabolic pathway of glucose, AAs, and fatty acids in human cells and the central pathway of oxidative phosphorylation, ensuring that the requirements of cellular bioenergy, substance synthesis, and redox balance are met (16). Abnormal proliferation is a hallmark of tumor cells, wherein cells increase the production of biological macromolecules and energy to meet the higher metabolic demands (17, 18). Tumor cells, even in the presence of sufficient oxygen for mitochondrial oxidative phosphorylation, tend to metabolize glucose to lactic acid and utilize products such as asparagine, glutamine, and fatty acids to meet their increased metabolic demands rather than induce the more energy-efficient TCA cycle pathway (19–21). Thus, we speculate that if the TCA cycle is blocked, tumor cells can meet the higher metabolic demands through the glycolysis pathway. Studies report that when tumor cells survive under glycolysis conditions, copper-induced cell death is also reduced (22), thereby increasing tumor proliferation and invasive ability. Recently, several studies demonstrated that multiple components in the TCA cycle play crucial roles in regulating energy uptake and decomposition. For instance, the abnormal regulation of isocitric dehydrogenase (IDH), α-oxoglutarate dehydrogenase complex, and fumarate hydratase promotes the occurrence and development of tumors. Therefore, the role of pyruvate metabolism and the TCA cycle pathway in tumor metabolism and disease progression should be explored in this context (20, 23–25).

Acetyl coenzyme A (A-CoA) is an essential substrate for the TCA cycle. It catalyzes the oxidative decarboxylation of pyruvate through the pyruvate dehydrogenase complex (PDHC) and generates NADH, which enters the electron transfer chain (ETC) as a biohydrogen carrier and electron donor and eventually combines with oxygen to form water and drive ATP production (26, 27). DLAT, an important part of the PDHC, is involved in the generation of A-CoA and can link glycolysis with the TCA cycle (28). If the expression of DLAT is reduced, A-CoA will be reduced, resulting in the inhibition of the TCA cycle. Recent studies have shown that DLAT is related to copper-induced cell death. Cu⁺ can directly bind lipoylated DLAT and undergo polymerization, which produces toxic stress and leads to cell death. This new cell death mode is termed cuproptosis. Inhibition of the TCA cycle leads to the accumulation of A-CoA, which reduces DLAT activity and inhibits cuproptosis through competitive feedback (22, 29). Thus, the inhibition of the pyruvate metabolics–TCA cycle pathway is hypothesized to inhibit cuproptosis through the feedback inhibition of DLAT expression, ultimately promoting tumor progression.

In addition to PDHC catalyzing the oxidative decarboxylation of pyruvate to generate NADH, other reactions in the TCA cycle can also generate NADH, including IDH catalyzing the oxidative decarboxylation of isocitric acid to generate α-ketoglutaric acid (α-KG) and NADH, α-oxoglutarate dehydrogenase complex catalyzing the oxidation of α-KG decarboxylation to generate succinyl-coenzyme A and NADH, and L-malate dehydrogenase catalyzing the oxidative decarboxylation of malic acid to generate oxaloacetic acid and NADH. NADH generated above enters the ETC through respiratory chain complex I (30). Additionally, succinate dehydrogenase in the TCA cycle catalyzes the oxidative decarboxylation of succinate to generate fumaric acid and reduce flavine adenine dinucleotide, wherein the latter enters the ETC through respiratory chain complex II (31). These reactions connect the TCA cycle with the ETC. Additionally, succinate-CoA ligase subunits α and β (SUCLG1/2) hydrolyze succinyl-coenzyme A to succinic acid, and ACO2 catalyzes the isomerization of citrate to isocitric acid, indirectly affecting the production of NADH and flavine adenine dinucleotide. If the expression of IDH3B, SUCLG1/2, and ACO2 is downregulated, it reduces the TCA cycle intermediates and inhibits mitochondrial oxidative phosphorylation. Moreover, the reduction or interruption in ETC activity leads to the impairment of cellular oxygen utilization, resulting in hypoxia.

Tumor cells exhibit various adaptive mechanisms in hypoxic niches (32), with hypoxia-inducible factor-1α (HIF1α) being a key factor in the tumor hypoxia response (33). In hypoxic environments, HIF1α accumulates and migrates to the nucleus, in which it binds to the hypoxia response element located in the 5′ regulatory regions of its downstream target genes, thereby activating the transcriptional regulation of multiple enzymes (34–36). These enzymes promote glycolysis and lactic acid production while inhibiting PDHC activity and glucose transfer from the TCA cycle to the glycolytic pathway (25). Certain enzymatic reactions in the TCA cycle affect the accumulation of HIF. For example, mutant IDH converts α-KG to 2-hydroxyglutarate (37) and the accumulation of 2-hydroxyglutarate inhibits prolyl hydroxylase domain proteins and increases HIF1α production (38). The stable expression of HIF1α can upregulate the expression of target genes such as VEGF and promote tumor angiogenesis, invasion, and metastasis (39, 40). In addition to affecting the function of HIF, the accumulation of fumaric acid promotes tumor proliferation by reducing intracellular reactive oxygen species (ROS) levels via supporting the stabilization of NF-E2–related factor 2, a key regulator of the cellular antioxidant defense system (41).

The products of the pyruvate metabolism and TCA cycle pathway can generate ROS through the ETC (27, 42). ROS induces apoptosis mediated by CYCS (43) and ferroptosis with the participation of iron (44). Thioester formation of CoA mediated by ACSF2 catalyzes the initial reaction of fatty acid metabolism and is vital for mitochondrial lipid peroxidation and ROS generation in ferroptosis (44). Decreased ACSF2 expression suppresses ferroptosis by lowering ROS levels, thereby enhancing tumor invasion and metastasis. VDAC2 acts as a channel for Fe²⁺ to enter the outer mitochondrial membrane to induce ferroptosis (44–46); moreover, it promotes the binding of Bax to mitochondria to induce apoptosis (47). Similarly, the downregulation of VDAC2 promotes the invasion and migration of tumor cells by inhibiting ferroptosis and apoptosis.

It can be concluded that the mechanism of action of the eight hub proteins in tumor progression is not isolated and unrelated but functions as a network that can be regulated by each other (Fig. 4). The regulation of metastasis of these proteins in renal cell carcinoma, colorectal cancer, liver cancer, and ovarian cancer has been reported (26, 48–53). As a biomarker of LMDTC, whether it also acts through the above mechanism needs to be further studied.

Furthermore, we developed and internally validated three LMDTC risk prediction models based on clinical features and the expression of hub proteins, with the combined prediction model demonstrating superior performance compared with the single prediction model. Therefore, we selected the combined prediction model to construct a nomogram, which indicated that larger primary tumor size, higher TG levels, surrounding tissue invasion, N1b staging, and lower SUCLG1, DLAT, CYCS, and ACSF2 expression increase the risk of lung metastasis in patients with DTC.

This study has certain limitations. The lower rate of LMDTC resulted in a relatively small sample size, which could impact the accuracy of the predictive model. Future research should involve randomized, large sample, and multi-center studies to further optimize the predictive model. Furthermore, systematic basic research utilizing the eight hub proteins is necessary to elucidate their functions and mechanisms of action. Overall, this study aids in stratifying the risk of lung metastasis in patients with DTC and offers a solid foundation for developing clinical personalized treatment plans and follow-up strategies for high-risk patients.

T. Chen reports grants from the National Natural Science Foundation of China, the Project of Regional Science Fund, and Yunnan Fundamental Research Projects during the conduct of the study; as well as a patent for a method and system for constructing a risk prediction model for lung metastasis of differentiated thyroid cancer pending. No disclosures were reported by the other authors.

X. Peng: Data curation, software, formal analysis, investigation, visualization, writing–original draft, writing–review and editing. H. Zhao: Data curation, software, formal analysis, methodology. L. Ye: Resources, investigation. F. Hou: Data curation, software, formal analysis. Z. Yi: Data curation, investigation. Y. Ren: Resources. L. Lu: Data curation, investigation. F. Chen: Investigation. J. Lv: Investigation. Y. Wang: Data curation, software, investigation. H. Cai: Data curation, investigation. X. Zheng: Data curation, investigation. Q. Yang: Resources. T. Chen: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

This work was supported by the National Natural Science Foundation of China (Grant No. 81960426, T. Chen), Yunnan Province Applied and Basic Research Foundation (Grant No. 202401AT070008, T. Chen), and Ten Thousand Talent Plans for Young Top-notch Talents of Yunnan Province (T. Chen).