The experimental samples were collected from the First Affiliated Hospital of Nanjing Medical University from June 2015 to June 2021. Blood samples were collected from 6:00 am to 8:30 am after overnight fasting and kept under 4°C before stored at –80°C within 6 hours after plasma isolation (14). Subjects included in this study were free of metabolic abnormalities such as hypoproteinemia, weight loss, and negative nitrogen balance. This prospective study was approved by the Ethics Committee of the First Affiliated Hospital of Nanjing Medical University (No. 2016-SRFA-149), and informed consent was obtained from all subjects.

1, 2-¹³C2-Myristic acid, methyl myristate, methoxamine hydrochloride (purity 98%) and pyridine (≥99.8% GC) were purchased from Sigma-Aldrich (St. Louis, MO, USA). N-methyl-trimethylsilyltrifluoroacetamide (MSTFA) and 1% trimethylchlorosilane (TMCS) were provided by Pierce Chemical (Rockford, IL, USA). Methanol and n-heptane were HPLC grade and obtained from Merck (Darmstadt, Germany).

Plasma samples were processed, extracted, and derived in accordance with our previously developed methods (15). 50.0 μL plasma was added into 200.0 μL methanol (containing 1, 2-¹³C₂-Myristic acid, 5.0 µg/mL). The specimens were vigorously extracted for 5.0 min and centrifuged at 20000×g for 10.0 min at 4°C. A 100.0 μL aliquot of the resulting supernatant was transferred to a GC vial and evaporated to dryness in a Speed-Vac concentrator (Thermo Fisher Scientific, Savant™ SC250EXP, Holbrook, USA). The dried plasma samples were then methoxymated, where 30.0 μL of methoxyamine pyridine solution (10.0 mg/mL) was added to the residue and incubated for 16 h at room temperature. Then the samples were trimethylsilylated for another 1.0 h by adding 30.0 μL of MSTFA with 1% TMCS as a catalyst. At last, 30.0 μL n-heptane with methyl myristate (15.0 µg/mL) was added to each GC vial as an external standard to monitor the stability of GC/MS.

A 0.5 μL sample aliquot was injected into gas chromatography coupled to a mass spectrometer (Shimadzu GCMS-QP2010 Ultra, Kyoto, Japan) in split mode (split ratio 8:1). It was equipped with an Rtx-5MS capillary column (0.25 mm × 30 m × 0.25µm, Restek, PA, USA). The injector temperature was set at 250°C. Helium was used as the carrier gas at a 1.5 mL/min flow rate. The column temperature was initially kept at 80°C for 3.0 min, then raised to 300°C at a rate of 20°C/min, and held for 5.0 min. The mass spectrometer ion source temperature and interface temperature were both 220°C, and ions were generated by a 70-eV electron beam at a current of 3.2 mA. The mass spectra were acquired over the mass range of 50-700 m/z in a full scan mode, with each run of 19.0 min. Following the abovementioned procedure, the quality control (QC) samples were prepared for the plasma pool. All samples were randomly selected for GC-MS analysis to diminish the systematic variations.

The metabolites were by comparing the mass spectrum and retention indexes for the analyte with the corresponding values from the literature and various libraries [e.g., Mainlib and Public in the National Institute of Standards and Technology (NIST) library 2.0 (2008) and Wiley 9 (Wiley-VCH Verlag GmbH & Co KGaA, Weinheim, Germany)]. Some standard compounds were also utilized to identify the metabolites.

After normalization against the internal standard, all the semiquantitative data from GC/MS were log₁₀-transformed. The transformed data were imported into SIMCA-P 14.1 software (Germany, Sartorius, Goettingen) and pre-processed for multivariate statistical analysis using par scaling. Principal component analysis (PCA) and partial least square to latent structure discriminant analysis (PLS-DA) models were built and plotted to show the clustering or separation of samples from different groups. PLS-DA models were constructed and plotted to show the clustering or separation of samples from different groups. The goodness of fit for the PLS-DA models was evaluated using three quantitative parameters: R²X, R²Y and Q². R²X and R²Y are the explained variations, and Q² is the predicted variation, with a higher level of R²Y and Q²Y indicating the model’s better fit and predictive performance. To avoid the classification obtained by supervised learning methods being chance and to test whether the model reproduces well and whether the data in the model are over-fitted, the validity of the built model was examined by a 7-fold cross-check and replacement test (200 times, cross-validation). The intercept of the R² and Q² regression lines to the axes was used to measure overfitting, and the model was valid when the intercept of Q² was negative.

To assess the differences among the groups, analysis of variance (ANOVA) and multiple comparisons (LSD) were used. The independent-sample t-test and Mann-Whitney U test were used to analyze normally and non-normally distributed data. OR calculations and ROC curve analysis were performed using SPSS 26.0 (SPSS Inc., Chicago, IL, USA), bar graphs were produced using GraphPad Prism 8.0, and heatmap and pathway analysis were performed using the online software MetaboAnalyst (https://www.metaboanalyst.ca/).

All individuals were randomly divided into train and validation datasets in a ratio of seven to three. Hyperparameters of seven machine learning models [including XGBoost, AdaBoost, Random Forest (RF), Gaussian Naive Bayes (GNB), Multilayer Perceptron (MLP), Support Vector Machine (SVM), and k-Nearest Neighbors (KNN)] were optimized using 10-fold cross-validation. Ten groups were created randomly from the training set. In each iteration of the 10-fold cross-validation method, nine groups were randomly selected for training and the remaining groups were used as test sets. Thus, the test sets for each group were selected sequentially, which ensured that the evaluation results did not overlap. Then, to minimize errors due to unreasonable test set selection, the results of the 10 evaluations were averaged.

The ROC curve analysis, calibration curve, decision curve analysis (DCA), accuracy, F1-score, sensitivity and specificity were used to assess the model’s performance. Model discrimination was assessed with ROC analysis, and the accuracy of its predictions was assessed with AUC. The calibration curve showed the calibration and the extent to which the model’s predictions deviated from actual events. Clinical utility and net benefit were assessed with DCA, which allows estimating the net benefit by calculating the difference between the true positive rate and the false positive rate, weighted by the odds ratio of the selected risk threshold probabilities.

All analyses were performed with R software (version 4.0) and Python version 3.7.

RNA-sequencing expression (level 3) profiles and corresponding clinical information for lung adenocarcinoma were downloaded from the TCGA dataset (https://portal.gdc.com). The current-release (V8) GTEx datasets were obtained from the GTEx data portal website (https://www.gtexportal.org/home/datasets). Statistical analyses were performed using R software v4.0.3 (R Foundation for Statistical Computing, Vienna, Austria). P-value <0.05 was considered statistically significant. All analyses were performed using the online website HOME for Researchers (https://www.home-for-researchers.com/static/index.html#/).

A549 and BEAS-2B cell lines were purchased from the Type Culture Center, Chinese Academy of Sciences (Shanghai, China). These cell lines were grown in RPMI-1640 medium supplemented with 10% (v/v) fetal bovine serum and 100 U/mL penicillin and streptomycin at 37°C and 5% CO2. GC/MS and metabolomics analysis methods were the same as previously reported (16).

A total of 1083 subjects were included in this study, including 670 lung adenocarcinoma (LUAD), 135 benign lung nodules (BLN) and 278 healthy controls (HC). Benign lung nodules mainly include pulmonary hamartomas, hemangiomas, and inflammatory pseudotumors. LUAD and BLN participants were newly diagnosed and had not undergone anti-cancer treatment, including radiation therapy, chemotherapy, surgical intervention, or medication administration. All patients included in this study were diagnosed by pathological examination. The distribution of subjects is shown in Figure 1A and Supplementary Table 1 .

GC/MS analysis of the plasma samples aligned the metabolites in typical chromatograms ( Supplementary Figure 1 ). Deconvolution of the GC/MS chromatograms produced 103 independent peaks from the plasma samples, 54 of which were authentically identified as metabolites ( Supplementary Table 2 ). Quantitative data were acquired for each metabolite in the plasma samples of the HC, BLN and LUAD cases.

The PCA and PLS-DA score plots ( Figures 1B, C ) demonstrated good clustering of pooled QC samples, indicating reproducibility of the assay and consistent instrument performance throughout the experiment. The supervised PLS-DA model showed that the samples in the HC and lung nodule groups (including LUAD and BLN) were distributed in different quadrants ( Figure 1C ), indicating significant metabolic differences between the two groups. Similarly, the overlap between LUAD and BLN indicates similar metabolic models. R²X, R²Y and Q² of the PLS-DA model were 0.628, 0.552 and 0.527, respectively, which suggested that the PLS-DA model had good adaptability and predictability. The permutation plot ( Supplementary Figure 2 ) demonstrated that the PLS-DA models were valid: the Q² regression line had a negative intercept, and all of the permuted R² values to the left of the intercept were lower than the original point to the right. These results suggest that both variability and similarity in metabolic patterns are present in the three groups.

Based on statistical analysis ( Table 1 ), 43 and 42 discriminant metabolites were found to differentiate LUAD and BLN patients from healthy controls, respectively. Similarly, LUAD cases primarily showed different metabolomic patterns from BLN patients. According to the statistical analysis, 26 distinct metabolites were identified between LUAD and BLN.

Of the 42 metabolites differentiating BLN from HC, the levels of glycolate, creatinine and aminomalonic acid were higher in BLN, while 4-hydroxybutanoic acid was lower, and all the above metabolites deviated further in LUAD ( Figure 1D ). These findings indicate that the above metabolites are involved in the development of lung nodules (from averagely minimal damage to lung cancer). Although glycine, phosphate and fructose-6-phosphate showed significant differences between the LUAD and BLN groups, they had no significant difference between the BLN and HC groups ( Figure 1E ). It is, therefore, suggested that these markers are associated with LUAD. In addition, BLN cases showed deviations in lysine, myo-inositol, and arachidonic acid levels, whereas the LUAD and HC groups did not show any significant differences ( Figure 1F ), suggesting that they are markers associated with BLN.

Metabolic pathway analysis showed that patients with lung nodules were more affected by the following metabolic pathways: aminoacyl-tRNA biosynthesis, alanine, aspartate and glutamate metabolism, glyoxylate and dicarboxylate metabolism, citrate cycle, arginine biosynthesis ( Figures 1G, H ). The differential metabolic pathways between LUAD patients and BLN were: aminoacyl-tRNA biosynthesis, alanine, aspartate and glutamate metabolism, glyoxylate and dicarboxylate metabolism, glutathione metabolism, arginine biosynthesis ( Figure 1I ).

The ROC analysis showed monopalmitin, succinate, glutamine, alpha-tocopherol, malate, asparagine and pyroglutamate performed well for HC and BLN differentiation (AUC ≥ 0.80) ( Supplementary Table 3 ). Meanwhile, succinate and creatinine performed well for the differentiation of HC and LUAD (AUC ≥ 0.80) ( Supplementary Table 4 ). However, each metabolite performed poorly in differentiating LUAD and BLN (AUC < 0.65) ( Supplementary Table 5 ).

ORs values were calculated to assess the role of the above metabolites as risk factors for predicting BLN/LUAD occurrence ( Table 2 ). Glycolate, oxalate, glyceric acid, β-alanine, aminomalonic acid, creatinine, cystine and arachidonic acid were found to be risk factors for BLN (ORs>1.0, P<0.05). In the meantime, glycolate, oxalate, glycine, glyceric acid, aminomalonic acid and creatinine were risk factors for LUAD (ORs>1.3, P<0.05). These results suggest that high levels of glycolate, oxalate, glyceric acid, aminomalonic acid and creatinine were independent risk factors of lung nodules. Whereas glycine may be a LUAD-specific risk factor, similarly, β-alanine, cystine and arachidonic acid may be BLN-specific risk factors.

Analyzing the differential metabolites between LUAD and BLN ( Table 2 ), we found that alanine, phosphate, proline, glycine, serine, threonine, aminomalonic acid, malate, pyroglutamate, creatinine, pyrophosphoric acid, asparagine, glutamine, ornithine, lysine, palmitic acid, myo-inositol, monopalmitin 18 substances may be risk factors for the development of BLN to LUAD (ORs>1.1, P<0.05).

The occurrence and progression of lung adenocarcinoma are progressive processes over time. Identifying the metabolic changes during lung adenocarcinoma progression is highly important for diagnosing, treating, and prognosis. The pTNM staging is the most commonly used method of tumor staging, which mainly consists of three parts: T-primary tumor size, N-lymph node metastasis and M-distant metastasis. We set to stage 0, stage I and stage II as early-stage, stage III and stage IV as advanced-stage, considering the progression of LUAD and pTNM staging.

Statistical analysis showed that the metabolism in lung adenocarcinoma changes with tumorigenesis and progression ( Table 3 ). Thirty-nine differential metabolites were found between early-stage LUAD and HC groups, and thirteen between advanced-stage LUAD and early-stage LUAD ( Figure 2A ). Notably, plasma asparagine, myo-inositol, ornithine, pyrophosphoric acid, threonine, and glutamine levels gradually decreased with increasing pTNM staging, while oxalate was elevated ( Figure 2B ). OR values suggested that oxalate may be a risk factor for progressive exacerbation of LUAD disease ( Supplementary Table 6 ). These findings suggest that the above metabolites may be implicated in the development and progression of LUAD (from early to the advanced stage).

The correlation between metabolites and tumor size was analyzed with Spearman. Ten metabolites (alanine, valine, phosphate, leucine, asparagine, glutamine, ornithine, lysine, palmitic acid, and linoleic acid) were negatively correlated with tumor size ( Supplementary Figure 3 ). With the appearance of lymphatic metastases, plasma levels of oxalate, glyceric acid, nonanoic acid, and arachidonic acid increased in LUAD, and phosphate, proline, glycine, serine, methionine, creatinine, pyrophosphoric acid, asparagine, glutamine, ornithine, lysine, uric acid, and myo-inositol were decreased ( Figure 2C ; Supplementary Table 7 ). In addition, high plasma levels of oxalate and glyceric acid may be risk factors for lymphatic metastasis in lung adenocarcinoma ( Supplementary Table 7 ). Significant differences were observed in the glycolate, oxalate, glycine, and arachidonic acid levels between LUAD with distant metastasis and non-distant metastases ( Supplementary Table 8 ). High levels of oxalate and arachidonic acid may be prognostic biomarkers for distant metastasis in lung adenocarcinoma ( Figure 2D ; Supplementary Table 8 ).

In summarizing the changes in metabolites during lung adenocarcinoma progression, it is clear that patients’ metabolic profiles are changing. As the pTNM stages advanced and tumor metastasis, plasma oxalate levels increased, indicating that oxalate is closely associated with lung adenocarcinoma progression. Notably, OR values suggested that plasma oxalate may be a risk marker of the progression in LUAD ( Figure 2E ).

To improve the differential diagnostic performance between HC, BLN and LUAD, we developed AI-based prediction models using seven machine learning algorithms (including XGBoost, RF, AdaBoost, MLP, SVM, KNN and GNB). Glycolate, oxalate, glyceric acid, aminomalonic acid and creatinine were used to build predictive models for lung nodules and healthy controls. The results are shown in Figure 3 ; Supplementary Table 9 . Compared with others, the RF model performs best in both the training and validation set (accuracy =1.00 and 0.85; F1-score =1.00 and 0.87; AUC=1.00 and 0.89, respectively). The calibration curves showed high consistency between the validation cohorts’ predicted and observed survival probability. The DCA results showed that the RF model had an excellent net benefit across the whole range of threshold probabilities. Meanwhile, the RF models also showed excellent differentiation effects for distinguishing BLN or LUAD from HC, respectively ( Supplementary Figures 4 , 5 ). These results suggested that the RF model has the best diagnostic accuracy and applicability for distinguishing lung nodules from healthy controls.

4-hydroxybutanoic acid, monopalmitin, myo-inositol, β-alanine, oxalate, alanine, fructose-6-phosphate, glycolate, phosphate, and aminomalonic acid were screened out to construct prediction models for LUAD and BLN ( Supplementary Figure 6 ). The results are shown in Figure 4 ; Supplementary Table 10 . Compared with others, the RF model performs best in both the training and validation set (accuracy =1.00 and 0.73; F1-score =1.00 and 0.79; AUC=1.00 and 0.76, respectively). The calibration curves showed high consistency between the validation cohorts’ predicted and observed survival probability. The DCA results showed that the RF model had an excellent net benefit across the whole range of threshold probabilities. These results suggested that the RF model can distinguish LUAD from BLN with high performance and accuracy.

A sample set was obtained from the Affiliated Hospital of Nanjing University Medical School in Nanjing, China for the purpose of external validation. This sample set comprised of 33 healthy controls, 23 patients with benign lung nodules, and 77 patients diagnosed with LUAD. The RF model previously developed was assessed on this distinct sample set, with the outcomes presented in Figure 5 . The model exhibited high accuracy in predicting healthy controls and benign lung nodules, achieving an AUC of 0.99 ( Figure 5A ). The calibration curve and decision curve analysis further validated the model’s precision and substantial net benefit ( Figures 5B, C ). In patients with benign lung nodules and lung adenocarcinoma, the model demonstrated an AUC of 0.866, indicating its relevance to these specific subgroups ( Figures 5D, F ). A noteworthy finding of the study was the correlation between elevated plasma oxalate levels and advanced pTNM stage and tumor metastasis in LUAD patients, as evidenced by an independent external validation set, underscoring the significant association between oxalate and LUAD progression ( Figures 5G–I ).

The primary sources of oxalate in humans are exogenous dietary intake (20-50%) and endogenous synthesis (50-80%) (17). The subjects collected in this study were mainly from Jiangsu Province, China, and there was slight dietary variation between subjects. Although ascorbic acid has been shown to have an essential effect on oxalate production, the process proceeds mainly through a non-enzymatic reaction, and the mechanism of action is unknown (18). In addition, we have ruled out the possibility of additional ascorbic acid supplementation in LUAD patients during treatment.

Glyoxalate is the primary precursor of endogenous oxalate synthesis (19). TCGA and GTEx database results showed significant differences in the expression of enzymes and transporters related to oxalate metabolism (LDHA, GRHPR, AGT, DAO, HAO2 and SLC26A1) in tumor tissues of lung adenocarcinoma patients ( Figures 6A, B ). Univariate and multivariate Cox regression analyses showed ( Figure 6C ) that LDHA, a key enzyme in the biosynthesis of oxalate, was significantly associated with the overall survival (OS) of LUAD patients and was a separately available prognostic indicator (HR = 1.66, 95% CI = 1.34 - 2.07, P < 0.001). A risk score was then calculated for each patient based on LDHA expression levels and risk factors, and patients were classified into low-risk and high-risk ( Figure 6D ). The heatmap revealed that high-risk patients tended to express LDHA genes at high levels, and low-risk patients tended to express LDHA genes at low levels ( Figure 6D ). Survival curves showed that patients with low-risk scores significantly had longer survival times than those with high-risk scores ( Figure 6E ). ROC analysis showed that LDHA could predict prognosis ( Figure 6F ).

The expression of LDHA in tumor tissues of lung adenocarcinoma increases with the advancement of the pTNM stage and lymphatic metastasis, suggesting that LDHA may be associated with abnormal immune function ( Figures 7A–C ). The European prospective investigation into cancer and nutrition (EPIC) array was used to analyze the correlation between immune cells and LDHA in the TCGA dataset. LDHA was negatively correlated with the expression of B cells, T cell CD8+, and endothelial cells ( Figures 7D–F ) and positively correlated with the expression of NK cells and uncharacterized cells ( Figures 7G, H ). In conclusion, LDHA is positively associated with the progression and prognosis of lung adenocarcinoma. The performance of LDHA further demonstrates the strong potential of oxalate as a risk factor of tumor progression in lung adenocarcinoma.

In order to assess the LDHA expression levels in LUAD, the protein expression of the LDHA gene was analyzed through immunohistochemistry, utilizing data from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). The results indicated a significantly higher expression of LDHA in lung adenocarcinoma tissue compared to normal lung tissue ( Figure 7I ).

The relative content of oxalate and the expression of LDHA were further examined and analyzed in A549 and BEAS-2B cells. The results, depicted in Figures 7J, K , indicated a significantly higher oxalate content in A549 cells compared to BEAS-2B cells, accompanied by elevated LDHA expression. These findings suggest that increased LDHA and oxalate levels may represent a specific alteration in LUAD.