Insights into serum metabolic biomarkers for early detection of incident diabetic kidney disease in Chinese patients with type 2 diabetes by random forest

Baseline characteristics of study participants

An overview of this study design was shown in Figure 1. For the 164 diabetic patients in the discovery and validation cohort, their median diabetic duration was 8 years, median eGFR was 76 (43–104) mL/min/1.73 m², and median UACR was 80 (ranges from 10 to 842) mg/g Cr. Approximately 67% of them had a history of ≥ one diabetic microvascular or macrovascular complication (Supplementary Table 1). The median eGFR and UACR of 30 healthy participants were 99 (95–112) mL/min/1.73 m² and 2.7 (2.3–3.9) mg/g Cr., respectively. The baseline characteristics of the two cohorts were compared, the patients group showed lower eGFR and an increase of UACR, RRI, systolic blood pressure (SBP), serum creatinine, urea, uric acid, cystatin C and urinary β2-microglobulin (β2-MG) concentrations along with DKD severity (Supplementary Table 2). We noticed that Stage 1a patients showed an enlargement of kidney size in compare with healthy participants by renal ultrasound images and testing body surface-area (BSA) related renal volumes (Supplementary Figure 1); however, for patients at Stage 2–4, their BSA related renal volumes were gradually decreased with DKD progressed (Supplementary Figure 1), which is consistent with previous findings that kidneys tended to be smaller in the most advanced stages of CKD [9]. Our results suggested that patients were likely to have abnormal kidney hypertrophy and enlargement at Stage 1a, followed by gradual renal atrophy and volume depletion at later stages.

Figure 1. The pipeline of this study. Three independent cohorts were recruited to perform the metabolite biomarker in discovery, validation and follow-up groups, respectively.

Characterization of metabolites in study participants

A total of 7480 compounds were detected and we found that the trend of MS-detected sCr among all the participants showed remarkable consistency with clinically measured sCr (Supplementary Figure 2), suggesting that UPLC-Orbitrap-MS is accurate and efficient for high-throughout metabolites detection. Step-wise filtering was performed based on two criteria: removing metabolites with unstable signals and only retaining metabolites with significant different concentration levels between healthy control and patients, 80 candidates (72 metabolites and 8 ratios) were screened out for a next-step analysis (Supplementary Tables 3 and 4, Supplementary Figure 3). The distribution plot of preprocessed data was shown as Supplementary Figure 4. The fold changes of these metabolites among different stages were calculated and shown in Supplementary Tables 5 and 6.

The 72 metabolites are classified to 7 categories: sulfate metabolites, amino acids, organic acids, acylcarnitine, purine derivatives, steroids and monosaccharides. Metabolomic network based on the 72 metabolites in the discovery set from Stage 0–4 was shown in Supplementary Figure 4. Comparing with healthy group (Stage 0), merely 3 down- and 6 up- regulated metabolites were found in Stage 1a; nevertheless, it increased to 7 down- and 40 up- regulated metabolites in Stage 4. The SAMs were enriched in Tryptophan metabolism (hsa00380) and Phenylalanine metabolism (hsa00360) pathways, indicating that amino acid metabolism disruption is a dominate signature of DKD (Figure 2). Among these significant altered metabolites (SAMs), 1,5-anhydro-D-glucitol (1,5-AG) was remarkably reduced in stages 1–4 compared to healthy group (fold change = −26.5 to −2.60, Supplementary Table 5). As demonstrated by previous studies that 1,5-AG is a potential biomarker for monitoring the progression of diabetes [10, 11], we therefore separately tested the correlation of two clinical glycemic markers – fast blood glucose (FBG) and hemoglobin A1C (HbA1c) with 1,5-AG in our cohort. It showed that 1,5-AG has strong negative correlation with HbA1c and FBG in stage 1a-3 patients (r ranges were −0.95 to −0.64 and −0.87 to −0.42, respectively); however, abnormal correlation was observed between 1,5-AG and FBG in stage 4 discovery sets, with r = 0.25 (Supplementary Figure 5A, 5B). Correlation between 1,5-AG and HbA1c was stronger among stage 1a-3 patients than stage 0–4 (Supplementary Figure 5A, 5B), suggesting that 1,5-AG may serve better as a potential glycemic marker in stage 1a-3 DKD patients than late stage.

Figure 2. KEGG pathway analysis of all SAMs. (A) An overview view of pathway analysis; (B) Table of the matched pathway with p-values from pathway enrichment analysis and pathway impact values from the pathway topology analysis using MetaboAnalyst 4.0 and KEGG database (Hits ≥ 2); (C, D) Simplified pathways of tryptophan metabolism and phenylalanine metabolism with the change trends of metabolites and their ratios at Stage 4 compared with the normal group.

Given that DKD is one of the consequences induced by diabetes, we hypothesized that decreased levels of 1,5-AG may be relevant with DKD development. However, our results showed that 1,5-AG exhibited non-significant correlation with neither eGFR nor UACR (Supplementary Figure 5C), indicating that diabetic progression has limited contribution to DKD development (scatter plot of 1,5-AG levels among healthy controls and different stages of patients was shown in Supplementary Figure 5D).

Identification of candidate DKD biomarkers (CDBs)

Receiver Operating Characteristic (ROC) curve analysis and Spearman’s coefficient coexpression analysis were used to evaluate the power of each metabolite as well as the combinations of every 2–9 compounds in DKD staging. The metabolites that closely associated (|r| ≥ 0.6) with eGFR in all participants were shown in Supplementary Figure 6. Strict rank coefficient cut-off values of 0.8 for Stage 0–4 (all participants) and Stage 1–4 (all patients), and 0.7 for Stage 1 and 2 (early-stage patients) were applied to identify biomarkers that closely correlated with eGFR progressive, four metabolites were screened out consist of SAdo ((M-H)^- = 382.1005 at 2.89 min), pseudouridne ((M-H)^- = 243.0622 at 0.93 min), ADT ((M-H)^- = 367.1497 at 2.21 min) and L,L-TMAP ((M+H)⁺ = 229.1546 at 1.06 min) (The regression plots among the four metabolites, eGFR and UACR see Figure 3). Since few studies have investigated the basic features of SAdo, the demonstration of its peak identification was plotted and calibrated (Supplementary Figures 7 and 8). To alleviate the bias induced by sex, age, SBP and UACR, partial correlation analysis (PCA) was performed among the four metabolites, serum cystatin C, MS-detected serum creatinine (MS-sCr) and log(eGFR). The PCA showed that log(MS-sCr), log(pseudouridine) and log(L,L-TMAP) were strikingly correlated with log(eGFR) (|r| > 0.9) in Stage 0–4 (Table 1). In addition, we evaluated the association between the four metabolites and kidney function related factors such as UACR, urinary β2- microglobulin, renal resistive index and the decrease of total BSA related renal volume, and found they were closely related as well (|r| > 0.5, see Table 1). The correlation among interested metabolites, total BSA-related renal volume and renal resistive index was calculated as well (Supplementary Table 7). In conclusion, our results suggested that the four metabolites are possibly involved in DKD progression and have potential to be utilized as candidate DKD biomarkers (CDBs).

Figure 3. Linear regression analysis among CDBs, sCr, eGFR (A) and UACR (B) in all stages and T2D patients after log10 transformation, respectively. Four metabolites showed similar strong predictive power with MS-detected sCr as their R2 of the equations were all above 0.85. R2 of UACR prediction model were close to 0.5.

Evaluation of CDB’s capacity in staging DKD by random forest (RF)

To evaluate whether CDBs can be applied for staging DKD, RF was employed to assess the classification power of the four CDBs as RF is a powerful supervised classification technique for decision making via building large numbers of decision tree models and merging all predictions from these trees to get an accurate and unprogressed prediction [12]. It exhibited that CDBs can specifically differentiate Stage 1a patients from Stage 1b-4, Stage 1a from Stages 1b-2, Stage 1b from Stage 2, and Stages 1a-1b from Stages 2-4 (average AUC > 0.700, Table 2). Among them, pseudouridine and SAdo achieve better performance than MS-sCr in all the staging process (Supplementary Table 8). Multiple combinations of the 4 metabolites and MS-sCr were generated and used for assessing their ability for DKD stratification. Any one of the four candidate DKD biomarkers combined with MS-sCr can gain higher AUC than MS-sCr alone (Table 2). Among these models, the No. 7 model (MS-sCr + pseudouridine + L,L-TMAP) ranks the best in phasing Stage 1a from the rest with average AUC > 0.9 (Table 2).

Comparison of CDBs between male and female patients

Increasing evidence suggested that sexual difference is a significant factor related with DKD progression, which leads to a complex personalized approach for DKD diagnosis and treatment in clinical practice [13, 14]. To investigate the association between the four CDBs and sexual difference, we compared their levels between male and female patients. For SAdo, pseudouridine and ADT at early and later stages, no significant differences was observed (Figure 4). Multiple linear regression analysis which included sex as a covariate showed that CDBs show insignificant sex dependence with eGFR, suggesting that CDBs can be utilized in both male and female patients (Supplementary Table 9).

Figure 4. The MS-detected CDBs. (A–D) and sCr (E) were evaluated and compared between male and female participants at the early and later stages’ patients. Unlike with sCr, three of the four CDBs include SAdo, pseudouridine and ADT displayed non-significant differences between male and female patients at early or later stages. F, women; M, men. p value was calculated by Student’s t-test and Mann-Whitney U according to the data normality. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001, respectively. Horizontal and error bars in the scatter plots represent mean ± SEM.

Predict eGFR using CDBs signatures

We hypothesized that the use of a combinations of multiple biomarkers may be more sensitive and specific than sCr in evaluating the kidney function of diabetic patients. To test the potential of CDBs in predicting DKD, non-parametric methods include univariate linear regression (ULR) and multivariate linear regression (MvLR) were applied to calculate the association among eGFR, UACR, creatinine and four CDBs. The ULR analysis using Stage 0–4 data found a high linear relationship between each CDB and log(eGFR) (training R² = 0.87–0.95, root mean square errors (RMSEs) = 0.08–0.13; predictive R² = 0.91–0.95) which was very close to MS-detected sCr (training R² = 0.95, RMSE = 0.11; predictive R² = 0.95) (Figure 5), suggested that CDBs are good covariates to be applied for eGFR prediction. A stepwise MvLR analysis was performed to test the effects of four CDBs and covariates (sex, age, SBP and UACR) in calculating CDB-predicted eGFR (BeGFR) using data from Stage 0–4 and Stage 0–2, respectively. Among all the individuals (Stage 0–4), MS-sCr, pseudouridine, L,L-TMAP and sex are the most significant variables, and were therefore considered as confounding covariants for BeGFR estimation. The predictive outcome R² was optimized from 0.971 (log(MS-detected sCr) and sex as covariates) to 0.987 (Table 3).

Figure 5. Univariate linear regression plots of BeGFR against MDRD eGFR using the four CDBs. (A–D) and MS-detected sCr (E) for all participants at Stages 0–4 after log₁₀ transformation. Univariate linear regression analysis of each selected metabolites with all participants’ log(MDRD eGFR) resulted in a high linear relationship (training R² = 0.85–0.94, root mean square errors (RMSEs) = 0.08–0.13; predictive R² = 0.91–0.95), which was similar with that of MS-detected sCr (training R² = 0.95, RMSE = 0.11; predictive R² = 0.95). ^***p < 0.001 β, unstandardized coefficient of linear regression. ^*R² was calculated based on the log(BeGFR) against log (eGFR) using the equation of the model and data of the discovery cohort. ^†R² was measured based on that using the equation of the model of the discovery cohort and data of the validation cohort.

Renal function of early-stage DKD (Stage 1 and 2) is reversible and manageable [15, 16]; however, most DKD patients are asymptomatic and indolent [17, 18]. Considered that early detection is of great vital for lifetime benefits for DKD patients, we specifically tested the predictive potency of CDBs and eGFR in early-stage participants. Surprisingly, pseudouridine and L,L-TMAP can enhance the predictive power of MS-detected sCr at early stages’ patients and healthy participants (Stage 0–2), predictive R² in the validation datasets was significantly improved from 0.70 (log(MS-detected creatinine) and sex as covariates) to 0.82 (Table 3), demonstrating that the two CDBs are potential biomarkers for the early detection of DKD. The best model for BeGFR estimation at the early stage (training R² = 0.7733, RMSEs = 0.0513) is: log(BeGFR) = −0.675 log(MS-detected sCr) −0.467 log(pseudourdine) + 0.101 (if male) + 1.559. Our results indicated that the combination of multiple biomarkers achieves better performance than standard sCr.

Follow-up study and prognostic assessment

Due to the limitations for purchasing the commercial standards of L-L-TMAP, only three CDBs include ADT, SAdo and pseudouridine and some clinical indexes were measured in the follow-up cohort. The association were assessed between 7 variates (sex, age, eGFR, sCr, ADT, SAdo and pseudouridine) as well as their combinations with DKD progression. The first-time collected serum samples from 106 subjects of Stage 2 and 3 DKD patients were examined to determine the concentration of the metabolites. In Stage 2 unprogressed and Stage 2 progressed groups, for each single variate, ADT ranks the top prognostic power (average AUC = 0.9184) than sCr alone (average AUC = 0.9133). Surprisingly, the combinations of “UACR + ADT + sCr” and “UACR + ADT + age + sex” were extremely associated with the future progression of DKD (AUC ranges from 0.9592 to 1) (Figure 6A, Supplementary Table 10). The same method was applied in Stage 3 unprogressed and Stage 3 progressed patients. Three single variates, including UACR, pseudouridine and ADT, gained the top 3 strongest association with DKD progression (average AUC values = 0.8889, 0.8302 and 0.8117, respectively). For the Stage 3 patients, either of the four metabolites (ADT, sCr, SAdo and pseudouridine) combined with UACR can optimize average AUC value ≥ 0.9012; among the combinations, “UACR + ADT + pseudouridine + SAdo + Sex” and “UACR + pseudouridine + sCr” achieve the best performance with an average AUC = 0.929 (Figure 6A, Supplementary Table 11). Ignoring the initial DKD phasing status, “UACR + ADT + SAdo” is the best combination for DKD prognostic assessment (average AUC = 0.9502).

Figure 6. Evaluation of the prognostic performance of CDBs in follow-up cohort. (A) The distribution of AUC values using single and combinations of variate(s) in follow-up patients. With stratified random sampling and random forest, AUC of distinction between and progressed patients were calculated 100 times with single and multiple variables. Results of AUC average and standard deviation indicated that ADT_SAdo_UACR (AUC average: 0.9502; CI: 0.9062–0.9805) manifested the best prediction, followed with UACR_ADT_SAdo_sCr (AUC average: 0.9482; CI:0.9248-0.9805) and ADT_UACR (AUC average: 0.9443; CI: 0.9141–0.9727). (B) Comparisons of the three CDB levels and UACR between “progressed” and “unprogressed” groups in all the follow-up individuals. (C) In stage 2 patients, levels of three CDBs and UCAR were remarkably different in “progressed” vs. “unprogressed”, ^**p < 0.01 and ^***p < 0.001, respectively. (D) Risk sores of three CDBs and UCAR in DKD progression by logistic regression analysis.

In general clinical practice, eGFR was employed as a popular reporter for grading diseases as its levels reflect the status of renal function decline; while UACR was mostly used as the predictive biomarker for disease’s progression [7]. Hence, we only compared the changes of UACR and CDBs between progressed and unprogressed groups in all follow-up individuals, it showed that the levels of ADT, pseudouridine and AUCR were significant different between the two categories while SAdo showed slightly but non-significant differences (Figure 6B). Among stage 2 incidents, concentrations of ADT, pseudouridine, Sado and AUCR were remarkably different in progressed patients comparing with the unprogressed (Figure 6C). Further logistic regression analysis indicated that ADT and pseudouridine are risk factors for DKD development (Figure 6D). With every single increase of standard deviation (SD) of ADT, the risk of DKD progression is enhanced for 2.151 folds; for pseudouridine, the risk scores are 1.741 folds (Figure 6D). To better evaluate the associations between the levels of risk factors and the DKD progressing, we used the duration of time for progressing to later stage in follow-up individuals for survival curve analysis. It showed that patients with higher levels of either pseudouridine or ADT had significant less survival probabilities (p < 0.05), which is similar to UACR (Supplementary Figure 9). Taken together, our cross-sectional study indicated that abnormal metabolism is involved in DKD progression and our follow-up study validated the predictive power of CDBs in DKD development.

Study design and participants

The 194 serum samples were collected spanning five DKD stages. As these two groups of samples were subjected to untargeted metabolites analysis at different times, we analyzed them separately to remove batch effect and separately referred them as the discovery and validation sets (Figure 1). Participants were required to cease taking unnecessary medications and fasted for 8 hours before serum collection. We examined all participants’ clinical parameters, reviewed and recorded their medical history, medical complications and dietary habits (see Supplementary Table 1). Based on the ADA and KDIGO criteria [33, 38] for diagnosis of diabetes, participants were firstly classified into healthy group (stage 0) and diabetic group. The healthy and diabetic groups are age and sex matched (Supplementary Table 2). We followed the MDRD formula: “eGFR (mL/min/1.73 m²) = 186 × (serum creatinine)^−1.154 × (age in years)^−0.203 × 0.742 (if female) × 1.210 (if African American)” to calculate the eGFR values [39]. According to the eGFR levels, diabetic individuals were further stratified to five types, which are stage 1a-b and stage 2 to 4 (criteria see Supplementary Materials and Methods); the healthy group participants were regarded as stage 0. The renal resistive index (RRI) was calculated as (peak systolic velocity - end diastolic velocity)/peak systolic velocity derived from the kidney doppler ultrasonography. Additional methods for measuring eGFR values were shown in Supplementary Table 12 [22]. All the DKD-related clinical parameters were measured in The Fourth Clinical Medical College of Guangzhou University of Chinese Medicine (Shenzhen Traditional Chinese Medicine Hospital) followed by the standard procedures.

Measurements of serum untargeted metabolites

Serum samples and an equal volume of quality control (QC) samples (Supplementary Table 3) were deproteinated with cold methanol that contains internal standards. Ultra-Performance Liquid Chromatography-Orbitrap-Mass Spectrometry (UPLC-Orbitrap-MS) analysis was conducted on a Waters ACQUITY UPLC system coupled to a Thermo Scientific Orbitrap Fusion Lumos Tribrid mass spectrometer for mass spectrometry (MS) analysis. For detailed steps of UPLC-Orbitrap-MS, please see the Supplementary Materials and Methods.

Untargeted metabolites analysis

UPLC-Orbitrap-MS data from the discovery and validation cohorts were analyzed separately. Data were firstly processed by Progenesis QI 2.3 software (Nonlinear Dynamics, Waters, Milford, MA, USA) for peak detection and alignment, and then subjected to Matlab (MathWorks, Natick, MA, USA) for exclusion of unreliable features with missing rates >40% and missing value imputation [40, 41]. We performed baseline correction via cubic spline interpolation to align the baseline levels of data obtained at different times [42] (Supplementary Figures 10–84). Unstable signals with a coefficient of variation (CV%) >30% across the QC samples were filtered out. Compounds were identified upon matching their mass to charge ratio (m/z) and mass fragmentation patterns against available reference standards and Human Metabolome Database (hmdb.ca) (Supplementary Table 13) [43]. To gain a unique view of DKD, only metabolites that repeat the same mass fragmentation pattern, retention time and show the same trend of significant statistical differences in both discovery and validation sets were kept for further investigation.

Evaluation of CDBs’ performance on DKD disease stage classification using random forest

Random forest (RF) algorithm is applicable for evaluating the performance of metabolites on differentiating disease stages [44]. We used the RF package scikit-learn [45] that was implemented by Python (version 3.8) to evaluate the classification power of candidate DKD biomarkers (n_estimators = 10–100). Area under curve (AUC) using 1–9 metabolite models was determined by RF, respectively. To avoid overfitting, samples in each cohort were randomly and evenly divided into training and testing sets for the establishment and performance evaluation of the models. This step was repeated 100 times to obtain the mean AUC values using the testing set of the two cohorts.

Significant altered metabolites analysis

One-way ANOVA followed by Fisher’s LSD post-hoc test was used to identify significant altered metabolites (SAMs). Metabolites with p-value < 0.05 and false discovery rate (FDR) < 0.1 among any two of the stages were regarded as SAMs.

Metabonomic networks analysis

Networks of metabolites (Pubmed ID listed in Supplementary Table 13) were generated through MetaMapp [46] with default parameters. CytoScape [47] was used to visualize the networks. Their metabolic pathway output was generated on the basis of their KEGG reaction pairs while their chemical and structural relationships were constructed by their Tanimoto similarity.

Linear regression of CDBs for log(eGRF) estimation

Using CDBs and covariates such as sex, age, SBP and UACR, univariate linear regression (ULR) and multivariate linear regression (MvLR) were utilized to determine log(eGFR). Discovery group was used as the training set. The unstandardized regression coefficients (β) of the training set’s model were applied to generate equations for log(eGFR) prediction in testing set -- the validation cohort. Variables that contributed to the model with p < 0.05 were selected by stepwise linear regression analysis for the best model in favor of a simpler model.

Follow-up study and targeted metabolites prognostic assessment

To identify the progression of DKD, we recruited an extra cohort of 106 patients at stage 2 and 3 for an average 3.9 years of follow-up study. Patients’ serums were collected since the starting point and the sCr, Bun and eGFR were measured every 3 months. The serum samples collected at the first time were subjected to targeted metabolomic analysis for the quantification of ADT, SAdo and pseudouridine by an ultra-high performance liquid chromatography (Shimadzu, Kyoto, Japan) coupled with the AB SCIEX Q-Trap 5500 triple quadrupole mass spectrometer (AB SCIEX, Toronto, Canada). During the follow-up period, patients who remained at their original stage were regarded as the “unprogressed” group; those who progressed to later stages accompanying with a 25% drop in eGFR [48] were “progressed” group. Random forest algorithm [20] were used to evaluate covariates that associated with DKD progression. The prognostic power of biomarkers during DKD stage progression was assessed via performing 100 iterations for each variate. Logistic regression analysis was used to evaluate the risk scores of interested metabolites by SPSS. The average AUC values were calculated and compared among all the variates. Detailed methods of targeted metabolite analysis with UPLC-QQQ-MS/MS and random forest analysis can be found in the Supplementary Materials and Methods (Supplementary Tables 14–16).

Data availability

The raw data were submitted to National Genomics Data Center database under the bioproject PRJCA013833 (https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA013833). All the data analyzed during the current study are available from the corresponding author on reasonable request.