Purpose:

The rising global high incidence of differentiated thyroid carcinoma (DTC) has led to a significant increase in patients presenting with lung metastasis of DTC (LMDTC). This population poses a significant challenge in clinical practice, necessitating the urgent development of effective risk stratification methods and predictive tools for lung metastasis.

Experimental Design:

Through proteomic analysis of large samples of primary lesion and dual validation employing parallel reaction monitoring and IHC, we identified eight hub proteins as potential biomarkers. By expanding the sample size and conducting statistical analysis on clinical features and hub protein expression, we constructed three risk prediction models.

Results:

This study identified eight hub proteins—SUCLG1/2, DLAT, IDH3B, ACSF2, ACO2, CYCS, and VDAC2—as potential biomarkers for predicting LMDTC risk. We developed and internally validated three risk prediction models incorporating both clinical characteristics and hub protein expression. Our findings demonstrated that the combined prediction model exhibited optimal predictive performance, with the highest discrimination (AUC: 0.986) and calibration (Brier score: 0.043). Application of the combined prediction model within a specific risk threshold (0–0.97) yielded maximal clinical benefit. Finally, we constructed a nomogram based on the combined prediction model.

Conclusions:

As a large sample size study in LMDTC research, the identification of biomarkers through primary lesion proteomics and the development of risk prediction models integrating clinical features and hub protein biomarkers offer valuable insights for predicting LMDTC and establishing personalized treatment strategies.

Translational Relevance

The increasing incidence of differentiated thyroid carcinoma (DTC) has led to a rise in the number of patients with lung metastatic DTC (LMDTC), posing a significant challenge in clinical management. Consequently, there is a critical need for effective methods to stratify risk and predict metastatic potential in this patient population. Existing literature on predicting LMDTC is limited, with studies often limited by small tissue sample sizes. This study addresses this gap by employing a large sample of DTC primary lesions and utilizing 4D label–free proteomic analysis. Through dual validation of results via parallel reaction monitoring and IHC, the study identified eight hub proteins as potential biomarkers for predicting LMDTC. This approach has significant advantages in terms of sample size and technology advancement. Moreover, by expanding the sample size and conducting statistical analysis of clinical features and hub proteins, this study is the first to establish a combined prediction model for LMDTC risk. Further optimization of this model holds promise for enhancing risk stratification in patients with DTC, enabling the formulation of individualized clinical treatment plans and follow-up strategies for high-risk patients.

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

Sample collection

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

Sample preparation

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.

4D label–free quantitative proteomic analysis and identification of differentially expressed proteins

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 [FCLMDTC/NMDTC,P = Mean (RLMDTC,P)/Mean (RNMDTC,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).

Bioinformatics analysis

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 based on LC-MS/MS

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 and receiver operating characteristic analysis

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.

IHC

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.

Construction of the LMDTC risk prediction model

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 χ2 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.

Data availability

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

Figure 1.

Thirty candidate hub proteins were identified between the LMDTC and NMDTC groups based on proteomic analysis. A, Heatmap illustrating Pearson’s correlation analysis in all samples. The color bar indicates Pearson’s correlation coefficients. B, Volcano plot displaying significant DEPs between LMDTC and NMDTC groups. C, The top six proteins with the most pronounced upregulation and downregulation in expression. Significance levels: *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. D, Top 10 enriched GO terms for GO biological processes (blue), GO molecular functions (yellow), and GO cellular components (red). The point size reflects the number of proteins. E, Top 20 enriched KEGG pathways. F, Direct PPI network of the top 50 DEPs. The point size represents the degree of the node. G, Venn diagram showing proteins in the top 20 enriched KEGG pathways (25 DEPs) and proteins with a degree >40 in the PPI network. H, Module–trait relationships. The color scale (red–blue) represents the strength of the correlation between the module and the trait. I, Scatter plot of the module membership and gene significance in the blue module. J, Venn diagram of 116 hub proteins in the blue module, 25 DEPs in the top 20 enriched pathways, and proteins with a degree >40 in the PPI network. K, ROC curves and AUC of five DEPs is shown in Fig. 1J. L, Venn diagram of 418 DEPs and ferroptosis-related genes (red indicates suppressor genes, and green indicates driver genes).

Figure 1.

Thirty candidate hub proteins were identified between the LMDTC and NMDTC groups based on proteomic analysis. A, Heatmap illustrating Pearson’s correlation analysis in all samples. The color bar indicates Pearson’s correlation coefficients. B, Volcano plot displaying significant DEPs between LMDTC and NMDTC groups. C, The top six proteins with the most pronounced upregulation and downregulation in expression. Significance levels: *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001. D, Top 10 enriched GO terms for GO biological processes (blue), GO molecular functions (yellow), and GO cellular components (red). The point size reflects the number of proteins. E, Top 20 enriched KEGG pathways. F, Direct PPI network of the top 50 DEPs. The point size represents the degree of the node. G, Venn diagram showing proteins in the top 20 enriched KEGG pathways (25 DEPs) and proteins with a degree >40 in the PPI network. H, Module–trait relationships. The color scale (red–blue) represents the strength of the correlation between the module and the trait. I, Scatter plot of the module membership and gene significance in the blue module. J, Venn diagram of 116 hub proteins in the blue module, 25 DEPs in the top 20 enriched pathways, and proteins with a degree >40 in the PPI network. K, ROC curves and AUC of five DEPs is shown in Fig. 1J. L, Venn diagram of 418 DEPs and ferroptosis-related genes (red indicates suppressor genes, and green indicates driver genes).

Close modal

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

Figure 2.

Eight hub proteins were identified through the dual validation of PRM and IHC, showing potential as biomarkers for predicting LMDTC. A, PRM validation of 30 candidate hub proteins, with 17 downregulated hub proteins statistically significant (green point). B, Boxplot of the top nine hub proteins with the highest differential expression between the LMDTC group (M) and the NMDTC group (NM). Significance levels: *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001; ****, P-value < 0.0001. C, Top 10 enriched GO terms for GO biological processes (blue), GO molecular functions (yellow), and GO cellular components (red). The point size indicates the number of proteins. D, Top 20 enriched KEGG pathways. E, PPI network of the 17 hub proteins. The 15 proteins with correlations to each other were used for the following analysis. F, Comparison of the TNMplot database, indicating eight hub proteins with expression trends consistent with it. G, ROC analysis of the eight hub proteins showed high predictive efficiency (AUC > 0.90, P < 0.05). H, IHC staining of the eight hub proteins. Scale bar, 0.100 mm (magnification, 200×). I, Bar plot of semi-quantitative IHC analysis. Significance levels: *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001.

Figure 2.

Eight hub proteins were identified through the dual validation of PRM and IHC, showing potential as biomarkers for predicting LMDTC. A, PRM validation of 30 candidate hub proteins, with 17 downregulated hub proteins statistically significant (green point). B, Boxplot of the top nine hub proteins with the highest differential expression between the LMDTC group (M) and the NMDTC group (NM). Significance levels: *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001; ****, P-value < 0.0001. C, Top 10 enriched GO terms for GO biological processes (blue), GO molecular functions (yellow), and GO cellular components (red). The point size indicates the number of proteins. D, Top 20 enriched KEGG pathways. E, PPI network of the 17 hub proteins. The 15 proteins with correlations to each other were used for the following analysis. F, Comparison of the TNMplot database, indicating eight hub proteins with expression trends consistent with it. G, ROC analysis of the eight hub proteins showed high predictive efficiency (AUC > 0.90, P < 0.05). H, IHC staining of the eight hub proteins. Scale bar, 0.100 mm (magnification, 200×). I, Bar plot of semi-quantitative IHC analysis. Significance levels: *, P-value < 0.05; **, P-value < 0.01; ***, P-value < 0.001.

Close modal

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.

Table 1.

Baseline clinical characteristics of patients and the expression of hub proteins.

NameLevelsDistant metastasis (N = 81)No distant metastasis (N = 100)P
Sex Female 51 (63%) 75 (75%) 0.112 
 Male 30 (37%) 25 (25%)  
Age Mean ± SD 43.4 ± 15.8 43.2 ± 10.5 0.941 
BMI Mean ± SD 23.0 ± 4.1 23.7 ± 3.1 0.185 
Benign_thyroid_lesions No 30 (37%) 39 (39%) 0.907 
 Yes 51 (63%) 61 (61%)  
Multifocal No 13 (16%) 39 (39%) 0.001 
 Yes 68 (84%) 61 (61%)  
AFP Mean ± SD 2.9 ± 1.7 3.1 ± 1.5 0.271 
WBC Mean ± SD 6.5 ± 1.7 6.2 ± 1.5 0.301 
RBC Mean ± SD 5.0 ± 0.5 4.9 ± 0.5 0.686 
Platelet Mean ± SD 254.6 ± 87.5 246.1 ± 71.0 0.481 
Hemoglobin Mean ± SD 147.7 ± 24.6 145.1 ± 16.1 0.428 
Lymphocyte Mean ± SD 2.1 ± 0.6 2.1 ± 0.6 0.971 
Tumor_size Mean ± SD 2.4 ± 1.5 0.7 ± 0.7 <0.001 
Neutrophils Mean ± SD 3.8 ± 1.5 3.7 ± 1.2 0.502 
Eosinophils Mean ± SD 0.2 ± 0.2 0.1 ± 0.1 0.152 
Basophil Mean ± SD 0.1 ± 0.6 0.0 ± 0.1 0.184 
CD3 Mean ± SD 68.2 ± 12.7 69.2 ± 9.9 0.547 
CD4 Mean ± SD 36.1 ± 9.9 36.9 ± 8.9 0.596 
CD8 Mean ± SD 23.9 ± 10.4 22.5 ± 6.8 0.328 
CD4_to_CD8_cell_ratio Mean ± SD 2.0 ± 1.3 1.9 ± 1.0 0.455 
NK Mean ± SD 13.9 ± 8.9 14.6 ± 7.9 0.61 
Blood_sugar Mean ± SD 5.1 ± 1.2 5.0 ± 0.8 0.624 
LDH Mean ± SD 169.7 ± 25.9 162.6 ± 28.7 0.088 
Ca Mean ± SD 2.4 ± 0.3 2.3 ± 0.3 0.551 
Mean ± SD 4.1 ± 0.6 4.1 ± 0.5 0.696 
Mean ± SD 1.2 ± 0.2 1.2 ± 0.2 0.164 
TG <10 11 (13.6%) 53 (53%) <0.001 
 ≥10 and <30 15 (18.5%) 32 (32%)  
 ≥30 55 (67.9%) 15 (15%)  
TSH Mean ± SD 3.6 ± 3.2 3.3 ± 2.2 0.371 
TGAb Mean ± SD 253.8 ± 833.4 85.0 ± 189.0 0.078 
TPOAb Mean ± SD 50.9 ± 116.3 69.0 ± 115.9 0.299 
Infiltrated_the_adjacent_tissue No 37 (45.7%) 91 (91%) <0.001 
 Yes 44 (54.3%) 9 (9%)  
T_stage T1/2 30 (37%) 89 (89%) <0.001 
 T3/4 51 (63%) 11 (11%)  
N_stage N0 13 (16%) 74 (74%) <0.001 
 N1 68 (84%) 26 (26%)  
SUCLG1 Medium positive 42 (51.9%) 17 (17%) <0.001 
 Strong positive 6 (7.4%) 71 (71%)  
 Weakly positive 33 (40.7%) 12 (12%)  
DLAT Medium positive 20 (24.7%) 37 (37%) <0.001 
 Negative/weakly positive 50 (61.7%) 15 (15%)  
 Strong positive 11 (13.6%) 48 (48%)  
IDH3B Medium positive 28 (34.6%) 36 (36%) <0.001 
 Strong positive 13 (16%) 44 (44%)  
 Weakly positive 40 (49.4%) 20 (20%)  
ACSF2 Medium positive 20 (24.7%) 12 (12%) 0.011 
 Negative/weakly positive 19 (23.5%) 15 (15%)  
 Strong positive 42 (51.9%) 73 (73%)  
SUCLG2 Medium positive 24 (29.6%) 33 (33%) 0.028 
 Negative/weakly positive 45 (55.6%) 38 (38%)  
 Strong positive 12 (14.8%) 29 (29%)  
ACO2 Medium positive 32 (39.5%) 36 (36%) 0.044 
 Negative/weakly positive 24 (29.6%) 17 (17%)  
 Strong positive 25 (30.9%) 47 (47%)  
CYCS Medium positive 22 (27.2%) 38 (38%) 0.008 
 Strong positive 19 (23.5%) 35 (35%)  
 Weakly positive 40 (49.4%) 27 (27%)  
VDAC2 Medium positive 34 (42%) 24 (24%) <0.001 
 Negative/weakly positive 21 (25.9%) 16 (16%)  
 Strong positive 26 (32.1%) 60 (60%)  
NameLevelsDistant metastasis (N = 81)No distant metastasis (N = 100)P
Sex Female 51 (63%) 75 (75%) 0.112 
 Male 30 (37%) 25 (25%)  
Age Mean ± SD 43.4 ± 15.8 43.2 ± 10.5 0.941 
BMI Mean ± SD 23.0 ± 4.1 23.7 ± 3.1 0.185 
Benign_thyroid_lesions No 30 (37%) 39 (39%) 0.907 
 Yes 51 (63%) 61 (61%)  
Multifocal No 13 (16%) 39 (39%) 0.001 
 Yes 68 (84%) 61 (61%)  
AFP Mean ± SD 2.9 ± 1.7 3.1 ± 1.5 0.271 
WBC Mean ± SD 6.5 ± 1.7 6.2 ± 1.5 0.301 
RBC Mean ± SD 5.0 ± 0.5 4.9 ± 0.5 0.686 
Platelet Mean ± SD 254.6 ± 87.5 246.1 ± 71.0 0.481 
Hemoglobin Mean ± SD 147.7 ± 24.6 145.1 ± 16.1 0.428 
Lymphocyte Mean ± SD 2.1 ± 0.6 2.1 ± 0.6 0.971 
Tumor_size Mean ± SD 2.4 ± 1.5 0.7 ± 0.7 <0.001 
Neutrophils Mean ± SD 3.8 ± 1.5 3.7 ± 1.2 0.502 
Eosinophils Mean ± SD 0.2 ± 0.2 0.1 ± 0.1 0.152 
Basophil Mean ± SD 0.1 ± 0.6 0.0 ± 0.1 0.184 
CD3 Mean ± SD 68.2 ± 12.7 69.2 ± 9.9 0.547 
CD4 Mean ± SD 36.1 ± 9.9 36.9 ± 8.9 0.596 
CD8 Mean ± SD 23.9 ± 10.4 22.5 ± 6.8 0.328 
CD4_to_CD8_cell_ratio Mean ± SD 2.0 ± 1.3 1.9 ± 1.0 0.455 
NK Mean ± SD 13.9 ± 8.9 14.6 ± 7.9 0.61 
Blood_sugar Mean ± SD 5.1 ± 1.2 5.0 ± 0.8 0.624 
LDH Mean ± SD 169.7 ± 25.9 162.6 ± 28.7 0.088 
Ca Mean ± SD 2.4 ± 0.3 2.3 ± 0.3 0.551 
Mean ± SD 4.1 ± 0.6 4.1 ± 0.5 0.696 
Mean ± SD 1.2 ± 0.2 1.2 ± 0.2 0.164 
TG <10 11 (13.6%) 53 (53%) <0.001 
 ≥10 and <30 15 (18.5%) 32 (32%)  
 ≥30 55 (67.9%) 15 (15%)  
TSH Mean ± SD 3.6 ± 3.2 3.3 ± 2.2 0.371 
TGAb Mean ± SD 253.8 ± 833.4 85.0 ± 189.0 0.078 
TPOAb Mean ± SD 50.9 ± 116.3 69.0 ± 115.9 0.299 
Infiltrated_the_adjacent_tissue No 37 (45.7%) 91 (91%) <0.001 
 Yes 44 (54.3%) 9 (9%)  
T_stage T1/2 30 (37%) 89 (89%) <0.001 
 T3/4 51 (63%) 11 (11%)  
N_stage N0 13 (16%) 74 (74%) <0.001 
 N1 68 (84%) 26 (26%)  
SUCLG1 Medium positive 42 (51.9%) 17 (17%) <0.001 
 Strong positive 6 (7.4%) 71 (71%)  
 Weakly positive 33 (40.7%) 12 (12%)  
DLAT Medium positive 20 (24.7%) 37 (37%) <0.001 
 Negative/weakly positive 50 (61.7%) 15 (15%)  
 Strong positive 11 (13.6%) 48 (48%)  
IDH3B Medium positive 28 (34.6%) 36 (36%) <0.001 
 Strong positive 13 (16%) 44 (44%)  
 Weakly positive 40 (49.4%) 20 (20%)  
ACSF2 Medium positive 20 (24.7%) 12 (12%) 0.011 
 Negative/weakly positive 19 (23.5%) 15 (15%)  
 Strong positive 42 (51.9%) 73 (73%)  
SUCLG2 Medium positive 24 (29.6%) 33 (33%) 0.028 
 Negative/weakly positive 45 (55.6%) 38 (38%)  
 Strong positive 12 (14.8%) 29 (29%)  
ACO2 Medium positive 32 (39.5%) 36 (36%) 0.044 
 Negative/weakly positive 24 (29.6%) 17 (17%)  
 Strong positive 25 (30.9%) 47 (47%)  
CYCS Medium positive 22 (27.2%) 38 (38%) 0.008 
 Strong positive 19 (23.5%) 35 (35%)  
 Weakly positive 40 (49.4%) 27 (27%)  
VDAC2 Medium positive 34 (42%) 24 (24%) <0.001 
 Negative/weakly positive 21 (25.9%) 16 (16%)  
 Strong positive 26 (32.1%) 60 (60%)  
Table 2.

Univariable/multivariable logistic regression analysis of risk factors for lung metastasis in patients with DTC.

Dependent: GroupLevelsDistant metastasisNo distant metastasisOR (univariable)OR (multivariable)
(N = 81)(N = 100)
Sex Female 51 (63%) 75 (75%)   
 Male 30 (37%) 25 (25%) 0.57 (0.30–1.07, P = 0.082)  
Age Mean ± SD 43.4 ± 15.8 43.2 ± 10.5 1.00 (0.98–1.02, P = 0.938)  
BMI Mean ± SD 23.0 ± 4.1 23.7 ± 3.1 1.06 (0.97–1.15, P = 0.173)  
Benign_thyroid_lesions No 30 (37%) 39 (39%)   
 Yes 51 (63%) 61 (61%) 0.92 (0.50–1.68, P = 0.787)  
Multifocal No 13 (16%) 39 (39%)   
 Yes 68 (84%) 61 (61%) 0.30 (0.15–0.61, P = 0.001) 1.28 (0.10–15.69, P = 0.845) 
AFP Mean ± SD 2.9 ± 1.7 3.1 ± 1.5 1.11 (0.92–1.35, P = 0.272)  
WBC Mean ± SD 6.5 ± 1.7 6.2 ± 1.5 0.91 (0.75–1.09, P = 0.300)  
RBC Mean ± SD 5.0 ± 0.5 4.9 ± 0.5 0.88 (0.47–1.63, P = 0.684)  
Platelet Mean ± SD 254.6 ± 87.5 246.1 ± 71.0 1.00 (0.99–1.00, P = 0.470)  
Hemoglobin Mean ± SD 147.7 ± 24.6 145.1 ± 16.1 0.99 (0.98–1.01, P = 0.410)  
Lymphocyte Mean ± SD 2.1 ± 0.6 2.1 ± 0.6 0.99 (0.61–1.60, P = 0.970)  
Tumor_size Mean ± SD 2.4 ± 1.5 0.7 ± 0.7 0.22 (0.13–0.36, P < 0.001) 0.20 (0.06–0.63, P = 0.006) 
Neutrophils Mean ± SD 3.8 ± 1.5 3.7 ± 1.2 0.93 (0.74–1.15, P = 0.491)  
Eosinophils Mean ± SD 0.2 ± 0.2 0.1 ± 0.1 0.25 (0.04–1.65, P = 0.150)  
Basophil Mean ± SD 0.1 ± 0.6 0.0 ± 0.1 0.33 (0.04–2.73, P = 0.301)  
CD3 Mean ± SD 68.2 ± 12.7 69.2 ± 9.9 1.01 (0.98–1.04, P = 0.534)  
CD4 Mean ± SD 36.1 ± 9.9 36.9 ± 8.9 1.01 (0.98–1.04, P = 0.594)  
CD8 Mean ± SD 23.9 ± 10.4 22.5 ± 6.8 0.98 (0.95–1.02, P = 0.306)  
CD4_to_CD8_cell_ratio Mean ± SD 2.0 ± 1.3 1.9 ± 1.0 0.90 (0.69–1.18, P = 0.443)  
NK Mean ± SD 13.9 ± 8.9 14.6 ± 7.9 1.01 (0.97–1.05, P = 0.608)  
Blood_sugar Mean ± SD 5.1 ± 1.2 5.0 ± 0.8 0.92 (0.68–1.25, P = 0.610)  
LDH Mean ± SD 169.7 ± 25.9 162.6 ± 28.7 0.99 (0.98–1.00, P = 0.090)  
Ca Mean ± SD 2.4 ± 0.3 2.3 ± 0.3 0.72 (0.24–2.16, P = 0.555)  
Mean ± SD 4.1 ± 0.6 4.1 ± 0.5 0.90 (0.52–1.55, P = 0.695)  
Mean ± SD 1.2 ± 0.2 1.2 ± 0.2 0.31 (0.06–1.55, P = 0.154)  
TG <10 11 (13.6%) 53 (53%)   
 ≥10 and <30 15 (18.5%) 32 (32%) 0.44 (0.18–1.08, P = 0.074) 0.04 (0.00–0.82, P = 0.036) 
 ≥30 55 (67.9%) 15 (15%) 0.06 (0.02–0.13, P < 0.001) 0.02 (0.00–0.35, P = 0.007) 
TSH Mean ± SD 3.6 ± 3.2 3.3 ± 2.2 0.95 (0.85–1.06, P = 0.358)  
TGAb Mean ± SD 253.8 ± 833.4 85.0 ± 189.0 1.00 (1.00–1.00, P = 0.102)  
TPOAb Mean ± SD 50.9 ± 116.3 69.0 ± 115.9 1.00 (1.00–1.00, P = 0.299)  
Infiltrated_the_adjacent_tissue No 37 (45.7%) 91 (91%)   
 Yes 44 (54.3%) 9 (9%) 0.08 (0.04–0.19, P < 0.001) 0.06 (0.00–4574334.86, P = 0.760) 
T_stage T1/2 30 (37%) 89 (89%)   
 T3/4 51 (63%) 11 (11%) 0.07 (0.03–0.16, P < 0.001) 4.32 (0.00–329669706.67, P = 0.875) 
N_stage N0 13 (16%) 74 (74%)   
 N1 68 (84%) 26 (26%) 0.07 (0.03–0.14, P < 0.001) 0.21 (0.02–2.52, P = 0.218) 
SUCLG1 Medium positive 42 (51.9%) 17 (17%)   
 Strong positive 6 (7.4%) 71 (71%) 29.24 (10.69–79.94, P < 0.001) 572.72 (20.66–15876.37, P < 0.001) 
 Weakly positive 33 (40.7%) 12 (12%) 0.90 (0.38–2.14, P = 0.809) 6.72 (0.89–50.97, P = 0.065) 
DLAT Medium positive 20 (24.7%) 37 (37%)   
 Negative/weakly positive 50 (61.7%) 15 (15%) 0.16 (0.07–0.36, P < 0.001) 0.27 (0.03–2.50, P = 0.247) 
 Strong positive 11 (13.6%) 48 (48%) 2.36 (1.01–5.53, P = 0.048) 19.63 (1.51–254.95, P = 0.023) 
IDH3B Medium positive 28 (34.6%) 36 (36%)   
 Strong positive 13 (16%) 44 (44%) 2.63 (1.19–5.81, P = 0.017) 1.90 (0.13–26.99, P = 0.634) 
 Weakly positive 40 (49.4%) 20 (20%) 0.39 (0.19–0.81, P = 0.011) 0.41 (0.06–3.09, P = 0.389) 
ACSF2 Medium positive 20 (24.7%) 12 (12%)   
 Negative/weakly positive 19 (23.5%) 15 (15%) 1.32 (0.49–3.52, P = 0.585) 0.76 (0.05–10.90, P = 0.843) 
 Strong positive 42 (51.9%) 73 (73%) 2.90 (1.29–6.51, P = 0.010) 10.10 (0.89–114.82, P = 0.062) 
SUCLG2 Medium positive 24 (29.6%) 33 (33%)   
 Negative/weakly positive 45 (55.6%) 38 (38%) 0.61 (0.31–1.21, P = 0.160)  
 Strong positive 12 (14.8%) 29 (29%) 1.76 (0.75–4.13, P = 0.196)  
ACO2 Medium positive 32 (39.5%) 36 (36%)   
 Negative/weakly positive 24 (29.6%) 17 (17%) 0.63 (0.29–1.38, P = 0.247)  
 Strong positive 25 (30.9%) 47 (47%) 1.67 (0.85–3.30, P = 0.139)  
CYCS Medium positive 22 (27.2%) 38 (38%)   
 Strong positive 19 (23.5%) 35 (35%) 1.07 (0.50–2.30, P = 0.869) 0.18 (0.01–2.47, P = 0.200) 
 Weakly positive 40 (49.4%) 27 (27%) 0.39 (0.19–0.80, P = 0.010) 0.21 (0.03–1.36, P = 0.102) 
VDAC2 Medium positive 34 (42%) 24 (24%)   
 Negative/weakly positive 21 (25.9%) 16 (16%) 1.08 (0.47–2.49, P = 0.858) 1.08 (0.09–13.43, P = 0.953) 
 Strong positive 26 (32.1%) 60 (60%) 3.27 (1.63–6.56, P < 0.001) 3.60 (0.42–30.91, P = 0.242) 
Dependent: GroupLevelsDistant metastasisNo distant metastasisOR (univariable)OR (multivariable)
(N = 81)(N = 100)
Sex Female 51 (63%) 75 (75%)   
 Male 30 (37%) 25 (25%) 0.57 (0.30–1.07, P = 0.082)  
Age Mean ± SD 43.4 ± 15.8 43.2 ± 10.5 1.00 (0.98–1.02, P = 0.938)  
BMI Mean ± SD 23.0 ± 4.1 23.7 ± 3.1 1.06 (0.97–1.15, P = 0.173)  
Benign_thyroid_lesions No 30 (37%) 39 (39%)   
 Yes 51 (63%) 61 (61%) 0.92 (0.50–1.68, P = 0.787)  
Multifocal No 13 (16%) 39 (39%)   
 Yes 68 (84%) 61 (61%) 0.30 (0.15–0.61, P = 0.001) 1.28 (0.10–15.69, P = 0.845) 
AFP Mean ± SD 2.9 ± 1.7 3.1 ± 1.5 1.11 (0.92–1.35, P = 0.272)  
WBC Mean ± SD 6.5 ± 1.7 6.2 ± 1.5 0.91 (0.75–1.09, P = 0.300)  
RBC Mean ± SD 5.0 ± 0.5 4.9 ± 0.5 0.88 (0.47–1.63, P = 0.684)  
Platelet Mean ± SD 254.6 ± 87.5 246.1 ± 71.0 1.00 (0.99–1.00, P = 0.470)  
Hemoglobin Mean ± SD 147.7 ± 24.6 145.1 ± 16.1 0.99 (0.98–1.01, P = 0.410)  
Lymphocyte Mean ± SD 2.1 ± 0.6 2.1 ± 0.6 0.99 (0.61–1.60, P = 0.970)  
Tumor_size Mean ± SD 2.4 ± 1.5 0.7 ± 0.7 0.22 (0.13–0.36, P < 0.001) 0.20 (0.06–0.63, P = 0.006) 
Neutrophils Mean ± SD 3.8 ± 1.5 3.7 ± 1.2 0.93 (0.74–1.15, P = 0.491)  
Eosinophils Mean ± SD 0.2 ± 0.2 0.1 ± 0.1 0.25 (0.04–1.65, P = 0.150)  
Basophil Mean ± SD 0.1 ± 0.6 0.0 ± 0.1 0.33 (0.04–2.73, P = 0.301)  
CD3 Mean ± SD 68.2 ± 12.7 69.2 ± 9.9 1.01 (0.98–1.04, P = 0.534)  
CD4 Mean ± SD 36.1 ± 9.9 36.9 ± 8.9 1.01 (0.98–1.04, P = 0.594)  
CD8 Mean ± SD 23.9 ± 10.4 22.5 ± 6.8 0.98 (0.95–1.02, P = 0.306)  
CD4_to_CD8_cell_ratio Mean ± SD 2.0 ± 1.3 1.9 ± 1.0 0.90 (0.69–1.18, P = 0.443)  
NK Mean ± SD 13.9 ± 8.9 14.6 ± 7.9 1.01 (0.97–1.05, P = 0.608)  
Blood_sugar Mean ± SD 5.1 ± 1.2 5.0 ± 0.8 0.92 (0.68–1.25, P = 0.610)  
LDH Mean ± SD 169.7 ± 25.9 162.6 ± 28.7 0.99 (0.98–1.00, P = 0.090)  
Ca Mean ± SD 2.4 ± 0.3 2.3 ± 0.3 0.72 (0.24–2.16, P = 0.555)  
Mean ± SD 4.1 ± 0.6 4.1 ± 0.5 0.90 (0.52–1.55, P = 0.695)  
Mean ± SD 1.2 ± 0.2 1.2 ± 0.2 0.31 (0.06–1.55, P = 0.154)  
TG <10 11 (13.6%) 53 (53%)   
 ≥10 and <30 15 (18.5%) 32 (32%) 0.44 (0.18–1.08, P = 0.074) 0.04 (0.00–0.82, P = 0.036) 
 ≥30 55 (67.9%) 15 (15%) 0.06 (0.02–0.13, P < 0.001) 0.02 (0.00–0.35, P = 0.007) 
TSH Mean ± SD 3.6 ± 3.2 3.3 ± 2.2 0.95 (0.85–1.06, P = 0.358)  
TGAb Mean ± SD 253.8 ± 833.4 85.0 ± 189.0 1.00 (1.00–1.00, P = 0.102)  
TPOAb Mean ± SD 50.9 ± 116.3 69.0 ± 115.9 1.00 (1.00–1.00, P = 0.299)  
Infiltrated_the_adjacent_tissue No 37 (45.7%) 91 (91%)   
 Yes 44 (54.3%) 9 (9%) 0.08 (0.04–0.19, P < 0.001) 0.06 (0.00–4574334.86, P = 0.760) 
T_stage T1/2 30 (37%) 89 (89%)   
 T3/4 51 (63%) 11 (11%) 0.07 (0.03–0.16, P < 0.001) 4.32 (0.00–329669706.67, P = 0.875) 
N_stage N0 13 (16%) 74 (74%)   
 N1 68 (84%) 26 (26%) 0.07 (0.03–0.14, P < 0.001) 0.21 (0.02–2.52, P = 0.218) 
SUCLG1 Medium positive 42 (51.9%) 17 (17%)   
 Strong positive 6 (7.4%) 71 (71%) 29.24 (10.69–79.94, P < 0.001) 572.72 (20.66–15876.37, P < 0.001) 
 Weakly positive 33 (40.7%) 12 (12%) 0.90 (0.38–2.14, P = 0.809) 6.72 (0.89–50.97, P = 0.065) 
DLAT Medium positive 20 (24.7%) 37 (37%)   
 Negative/weakly positive 50 (61.7%) 15 (15%) 0.16 (0.07–0.36, P < 0.001) 0.27 (0.03–2.50, P = 0.247) 
 Strong positive 11 (13.6%) 48 (48%) 2.36 (1.01–5.53, P = 0.048) 19.63 (1.51–254.95, P = 0.023) 
IDH3B Medium positive 28 (34.6%) 36 (36%)   
 Strong positive 13 (16%) 44 (44%) 2.63 (1.19–5.81, P = 0.017) 1.90 (0.13–26.99, P = 0.634) 
 Weakly positive 40 (49.4%) 20 (20%) 0.39 (0.19–0.81, P = 0.011) 0.41 (0.06–3.09, P = 0.389) 
ACSF2 Medium positive 20 (24.7%) 12 (12%)   
 Negative/weakly positive 19 (23.5%) 15 (15%) 1.32 (0.49–3.52, P = 0.585) 0.76 (0.05–10.90, P = 0.843) 
 Strong positive 42 (51.9%) 73 (73%) 2.90 (1.29–6.51, P = 0.010) 10.10 (0.89–114.82, P = 0.062) 
SUCLG2 Medium positive 24 (29.6%) 33 (33%)   
 Negative/weakly positive 45 (55.6%) 38 (38%) 0.61 (0.31–1.21, P = 0.160)  
 Strong positive 12 (14.8%) 29 (29%) 1.76 (0.75–4.13, P = 0.196)  
ACO2 Medium positive 32 (39.5%) 36 (36%)   
 Negative/weakly positive 24 (29.6%) 17 (17%) 0.63 (0.29–1.38, P = 0.247)  
 Strong positive 25 (30.9%) 47 (47%) 1.67 (0.85–3.30, P = 0.139)  
CYCS Medium positive 22 (27.2%) 38 (38%)   
 Strong positive 19 (23.5%) 35 (35%) 1.07 (0.50–2.30, P = 0.869) 0.18 (0.01–2.47, P = 0.200) 
 Weakly positive 40 (49.4%) 27 (27%) 0.39 (0.19–0.80, P = 0.010) 0.21 (0.03–1.36, P = 0.102) 
VDAC2 Medium positive 34 (42%) 24 (24%)   
 Negative/weakly positive 21 (25.9%) 16 (16%) 1.08 (0.47–2.49, P = 0.858) 1.08 (0.09–13.43, P = 0.953) 
 Strong positive 26 (32.1%) 60 (60%) 3.27 (1.63–6.56, P < 0.001) 3.60 (0.42–30.91, P = 0.242) 

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.

Figure 3.

Performance comparison of three risk prediction models for LMDTC. A, Confusion matrix of the combined prediction model based on clinical characteristics and hub proteins. B, Comparison of the ROC curves of the three models. C, Calibration curve of the combined prediction model. D, Comparison of the calibration curves of the three models. E, Comparison of the DCA curves of the three models. F, Nomogram of the combined prediction model for lung metastasis risk in DTC based on clinical characteristics and the expression of hub proteins. All analyses were performed using the bootstrap method, with resampling = 2,000 for internal validation.

Figure 3.

Performance comparison of three risk prediction models for LMDTC. A, Confusion matrix of the combined prediction model based on clinical characteristics and hub proteins. B, Comparison of the ROC curves of the three models. C, Calibration curve of the combined prediction model. D, Comparison of the calibration curves of the three models. E, Comparison of the DCA curves of the three models. F, Nomogram of the combined prediction model for lung metastasis risk in DTC based on clinical characteristics and the expression of hub proteins. All analyses were performed using the bootstrap method, with resampling = 2,000 for internal validation.

Close modal

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 (1921). 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, 2325).

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 (3436). 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 Fe2+ to enter the outer mitochondrial membrane to induce ferroptosis (4446); 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, 4853). As a biomarker of LMDTC, whether it also acts through the above mechanism needs to be further studied.

Figure 4.

Schematic diagram depicting the potential mechanisms by which hub proteins may influence tumor metastasis. (The eight hub proteins are represented by purple ellipses.)

Figure 4.

Schematic diagram depicting the potential mechanisms by which hub proteins may influence tumor metastasis. (The eight hub proteins are represented by purple ellipses.)

Close modal

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

Note: Supplementary data for this article are available at Clinical Cancer Research Online (http://clincancerres.aacrjournals.org/).

1.
Xu
S
,
Huang
H
,
Huang
Y
,
Qian
J
,
Wang
X
,
Xu
Z
, et al
.
Comparison of lobectomy vs total thyroidectomy for intermediate-risk papillary thyroid carcinoma with lymph node metastasis
.
JAMA Surg
2023
;
158
:
73
9
.
2.
Lim
H
,
Devesa
SS
,
Sosa
JA
,
Check
D
,
Kitahara
CM
.
Trends in thyroid cancer incidence and mortality in the United States, 1974–2013
.
JAMA
2017
;
317
:
1338
48
.
3.
Hirsch
D
,
Levy
S
,
Tsvetov
G
,
Gorshtein
A
,
Slutzky-Shraga
I
,
Akirov
A
, et al
.
Long-term outcomes and prognostic factors in patients with differentiated thyroid cancer and distant metastases
.
Endocr Pract
2017
;
23
:
1193
200
.
4.
Tahara
M
,
Kiyota
N
,
Hoff
AO
,
Badiu
C
,
Owonikoko
TK
,
Dutcus
CE
, et al
.
Impact of lung metastases on overall survival in the phase 3 SELECT study of lenvatinib in patients with radioiodine-refractory differentiated thyroid cancer
.
Eur J Cancer
2021
;
147
:
51
7
.
5.
Sugino
K
,
Nagahama
M
,
Kitagawa
W
,
Ohkuwa
K
,
Uruno
T
,
Matsuzu
K
, et al
.
Distant metastasis in pediatric and adolescent differentiated thyroid cancer: clinical outcomes and risk factor analyses
.
J Clin Endocrinol Metab
2020
;
105
:
dgaa545
.
6.
Liu
C
,
Xing
M
,
Wang
L
,
Zhang
K
.
miR-199a-3p downregulation in thyroid tissues is associated with invasion and metastasis of papillary thyroid carcinoma
.
Br J Biomed Sci
2017
;
74
:
90
4
.
7.
Liu
X
,
Fu
Y
,
Zhang
G
,
Zhang
D
,
Liang
N
,
Li
F
, et al
.
miR-424-5p promotes anoikis resistance and lung metastasis by inactivating hippo signaling in thyroid cancer
.
Mol Ther Oncolytics
2019
;
15
:
248
60
.
8.
Huang
Y
,
Yu
S
,
Cao
S
,
Yin
Y
,
Hong
S
,
Guan
H
, et al
.
MicroRNA-222 promotes invasion and metastasis of papillary thyroid cancer through targeting protein phosphatase 2 regulatory subunit B alpha expression
.
Thyroid
2018
;
28
:
1162
73
.
9.
Wang
W
,
Shen
C
,
Zhao
Y
,
Sun
B
,
Bai
N
,
Li
X
.
Identification and validation of potential novel biomarkers to predict distant metastasis in differentiated thyroid cancer
.
Ann Transl Med
2021
;
9
:
1053
.
10.
Bose
U
,
Wijffels
G
,
Howitt
CA
,
Colgrave
ML
.
Proteomics: tools of the trade
.
Adv Exp Med Biol
2019
;
1073
:
1
22
.
11.
Karczewski
KJ
,
Snyder
MP
.
Integrative omics for health and disease
.
Nat Rev Genet
2018
;
19
:
299
310
.
12.
Peng
L
,
Cantor
DI
,
Huang
C
,
Wang
K
,
Baker
MS
,
Nice
EC
.
Tissue and plasma proteomics for early stage cancer detection
.
Mol Omics
2018
;
14
:
405
23
.
13.
Guidelines Working Committee of Chinese Society of Clinical Oncology
.
Guidelines of Chinese society of clinical oncology (CSCO) differentiated thyroid cancer
.
J Cancer Control Treat
2021
;
34
:
1164
201
.
14.
Haugen
BR
,
Alexander
EK
,
Bible
KC
,
Doherty
GM
,
Mandel
SJ
,
Nikiforov
YE
, et al
.
2015 American thyroid association management guidelines for adult patients with thyroid nodules and differentiated thyroid cancer: the American thyroid association guidelines task force on thyroid nodules and differentiated thyroid cancer
.
Thyroid
2016
;
26
:
1
133
.
15.
Liu
X
,
Fu
Y
,
Zhang
G
,
Zhang
D
,
Liang
N
,
Li
F
, et al
.
miR-424-5p promotes anoikis resistance and lung metastasis by inactivating hippo signaling in thyroid cancer
.
Mol Ther Oncolytics
2019
;
15
:
248
60
.
16.
Anderson
NM
,
Mucka
P
,
Kern
JG
,
Feng
H
.
The emerging role and targetability of the TCA cycle in cancer metabolism
.
Protein Cell
2018
;
9
:
216
37
.
17.
Thompson
CB
.
Metabolic enzymes as oncogenes or tumor suppressors
.
N Engl J Med
2009
;
360
:
813
15
.
18.
Shorthouse
D
,
Bradley
J
,
Critchlow
SE
,
Bendtsen
C
,
Hall
BA
.
Heterogeneity of the cancer cell line metabolic landscape
.
Mol Syst Biol
2022
;
18
:
e11006
.
19.
Law
RC
,
Nurwono
G
,
Park
JO
.
A parallel glycolysis provides a selective advantage through rapid growth acceleration
.
Nat Chem Biol
2024
;
20
:
314
22
.
20.
Pavlova
NN
,
Thompson
CB
.
The emerging hallmarks of cancer metabolism
.
Cell Metab
2016
;
23
:
27
47
.
21.
Jin
HR
,
Wang
J
,
Wang
ZJ
,
Xi
MJ
,
Xia
BH
,
Deng
K
, et al
.
Lipid metabolic reprogramming in tumor microenvironment: from mechanisms to therapeutics
.
J Hematol Oncol
2023
;
16
:
103
.
22.
Tsvetkov
P
,
Coy
S
,
Petrova
B
,
Dreishpoon
M
,
Verma
A
,
Abdusamad
M
, et al
.
Copper induces cell death by targeting lipoylated TCA cycle proteins
.
Science
2022
;
375
:
1254
61
.
23.
Eng
C
,
Kiuru
M
,
Fernandez
MJ
,
Aaltonen
LA
.
A role for mitochondrial enzymes in inherited neoplasia and beyond
.
Nat Rev Cancer
2003
;
3
:
193
202
.
24.
Yan
H
,
Parsons
DW
,
Jin
G
,
McLendon
R
,
Rasheed
BA
,
Yuan
W
, et al
.
IDH1 and IDH2 mutations in gliomas
.
N Engl J Med
2009
;
360
:
765
73
.
25.
Shen
T
,
Wang
H
,
Tang
B
,
Zhu
G
,
Wang
X
.
The impact of RNA binding proteins and the associated long non-coding RNAs in the TCA cycle on cancer pathogenesis
.
RNA Biol
2023
;
20
:
223
34
.
26.
Arnold
PK
,
Finley
LWS
.
Regulation and function of the mammalian tricarboxylic acid cycle
.
J Biol Chem
2023
;
299
:
102838
.
27.
Ryan
DG
,
O’Neill
LAJ
.
Krebs cycle reborn in macrophage immunometabolism
.
Annu Rev Immunol
2020
;
38
:
289
313
.
28.
Arnold
PK
,
Jackson
BT
,
Paras
KI
,
Brunner
JS
,
Hart
ML
,
Newsom
OJ
, et al
.
A non-canonical tricarboxylic acid cycle underlies cellular identity
.
Nature
2022
;
603
:
477
81
.
29.
Xie
J
,
Yang
Y
,
Gao
Y
,
He
J
.
Cuproptosis: mechanisms and links with cancers
.
Mol Cancer
2023
;
22
:
46
.
30.
Kory
N
,
Uit de Bos
J
,
van der Rijt
S
,
Jankovic
N
,
Güra
M
,
Arp
N
, et al
.
MCART1/SLC25A51 is required for mitochondrial NAD transport
.
Sci Adv
2020
;
6
:
eabe5310
.
31.
Martínez-Reyes
I
,
Chandel
NS
.
Mitochondrial TCA cycle metabolites control physiology and disease
.
Nat Commun
2020
;
11
:
102
.
32.
Shen
T
,
Xia
W
,
Min
S
,
Yang
Z
,
Cheng
L
,
Wang
W
, et al
.
A pair of long intergenic non-coding RNA LINC00887 variants act antagonistically to control Carbonic Anhydrase IX transcription upon hypoxia in tongue squamous carcinoma progression
.
BMC Biol
2021
;
19
:
192
.
33.
Zhou
L
,
Jiang
J
,
Huang
Z
,
Jin
P
,
Peng
L
,
Luo
M
, et al
.
Hypoxia-induced lncRNA STEAP3-AS1 activates Wnt/β-catenin signaling to promote colorectal cancer progression by preventing m6A-mediated degradation of STEAP3 mRNA
.
Mol Cancer
2022
;
21
:
168
.
34.
Kim
JW
,
Tchernyshyov
I
,
Semenza
GL
,
Dang
CV
.
HIF-1-mediated expression of pyruvate dehydrogenase kinase: a metabolic switch required for cellular adaptation to hypoxia
.
Cell Metab
2006
;
3
:
177
85
.
35.
Wu
Q
,
You
L
,
Nepovimova
E
,
Heger
Z
,
Wu
W
,
Kuca
K
, et al
.
Hypoxia-inducible factors: master regulators of hypoxic tumor immune escape
.
J Hematol Oncol
2022
;
15
:
77
.
36.
Lv
R
,
Liu
X
,
Zhang
Y
,
Dong
N
,
Wang
X
,
He
Y
, et al
.
Pathophysiological mechanisms and therapeutic approaches in obstructive sleep apnea syndrome
.
Signal Transduct Target Ther
2023
;
8
:
218
.
37.
Zhao
J
,
Yao
K
,
Yu
H
,
Zhang
L
,
Xu
Y
,
Chen
L
, et al
.
Metabolic remodelling during early mouse embryo development
.
Nat Metab
2021
;
3
:
1372
84
.
38.
Dang
L
,
White
DW
,
Gross
S
,
Bennett
BD
,
Bittinger
MA
,
Driggers
EM
, et al
.
Cancer-associated IDH1 mutations produce 2-hydroxyglutarate
.
Nature
2009
;
462
:
739
44
.
39.
Martin
F
,
Linden
T
,
Katschinski
DM
,
Oehme
F
,
Flamme
I
,
Mukhopadhyay
CK
, et al
.
Copper-dependent activation of hypoxia-inducible factor (HIF)-1: implications for ceruloplasmin regulation
.
Blood
2005
;
105
:
4613
19
.
40.
Wu
Y
,
Terekhanova
NV
,
Caravan
W
,
Naser Al Deen
N
,
Lal
P
,
Chen
S
, et al
.
Epigenetic and transcriptomic characterization reveals progression markers and essential pathways in clear cell renal cell carcinoma
.
Nat Commun
2023
;
14
:
1681
.
41.
DeNicola
GM
,
Karreth
FA
,
Humpton
TJ
,
Gopinathan
A
,
Wei
C
,
Frese
K
, et al
.
Oncogene-induced Nrf2 transcription promotes ROS detoxification and tumorigenesis
.
Nature
2011
;
475
:
106
9
.
42.
Mills
EL
,
Kelly
B
,
Logan
A
,
Costa
ASH
,
Varma
M
,
Bryant
CE
, et al
.
Succinate dehydrogenase supports metabolic repurposing of mitochondria to drive inflammatory macrophages
.
Cell
2016
;
167
:
457
70.e13
.
43.
Liu
X
,
Kim
CN
,
Yang
J
,
Jemmerson
R
,
Wang
X
.
Induction of apoptotic program in cell-free extracts:requirement for dATP and cytochrome c
.
Cell
1996
;
86
:
147
57
.
44.
Battaglia
AM
,
Chirillo
R
,
Aversa
I
,
Sacco
A
,
Costanzo
F
,
Biamonte
F
.
Ferroptosis and cancer: mitochondria meet the “iron maiden” cell death
.
Cells
2020
;
9
:
1505
.
45.
Yagoda
N
,
von Rechenberg
M
,
Zaganjor
E
,
Bauer
AJ
,
Yang
WS
,
Fridman
DJ
, et al
.
RAS-RAF-MEK-dependent oxidative cell death involving voltage-dependent anion channels
.
Nature
2007
;
447
:
864
8
.
46.
Yang
Y
,
Luo
M
,
Zhang
K
,
Zhang
J
,
Gao
T
,
Connell
DO
, et al
.
Nedd4 ubiquitylates VDAC2/3 to suppress erastin-induced ferroptosis in melanoma
.
Nat Commun
2020
;
11
:
433
.
47.
Chin
HS
,
Li
MX
,
Tan
IKL
,
Ninnis
RL
,
Reljic
B
,
Scicluna
K
, et al
.
VDAC2 enables BAX to mediate apoptosis and limit tumor development
.
Nat Commun
2018
;
9
:
4976
.
48.
Yang
W
,
Wang
Y
,
Huang
Y
,
Yu
J
,
Wang
T
,
Li
C
, et al
.
4-Octyl itaconate inhibits aerobic glycolysis by targeting GAPDH to promote cuproptosis in colorectal cancer
.
Biomed Pharmacother
2023
;
159
:
114301
.
49.
Chen
Y
,
Gu
D
,
Wen
Y
,
Yang
S
,
Duan
X
,
Lai
Y
, et al
.
Identifying the novel key genes in renal cell carcinoma by bioinformatics analysis and cell experiments
.
Cancer Cell Int
2020
;
20
:
331
.
50.
Zhang
Z
,
Zhu
H
,
Li
Q
,
Gao
W
,
Zang
D
,
Su
W
, et al
.
Gene expression profiling of tricarboxylic acid cycle and one carbon metabolism related genes for prognostic risk signature of colon carcinoma
.
Front Genet
2021
;
12
:
647152
.
51.
Li
N
,
Li
H
,
Wang
Y
,
Cao
L
,
Zhan
X
.
Quantitative proteomics revealed energy metabolism pathway alterations in human epithelial ovarian carcinoma and their regulation by the antiparasite drug ivermectin: data interpretation in the context of 3P medicine
.
EPMA J
2020
;
11
:
661
94
.
52.
You
X
,
Tian
J
,
Zhang
H
,
Guo
Y
,
Yang
J
,
Zhu
C
, et al
.
Loss of mitochondrial aconitase promotes colorectal cancer progression via SCD1-mediated lipid remodeling
.
Mol Metab
2021
;
48
:
101203
.
53.
Zhao
Z
,
Liu
M
,
Xu
Z
,
Cai
Y
,
Peng
B
,
Liang
Q
, et al
.
Identification of ACSF gene family as therapeutic targets and immune-associated biomarkers in hepatocellular carcinoma
.
Aging (Albany NY)
2022
;
14
:
7926
40
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.

Supplementary data