Metabolomics analysis of type 2 diabetes remission identifies 12 metabolites with predictive capacity: a CORDIOPREV clinical trial study | BMC Medicine

Aim and objectives

The aim is to characterize the metabolic profile of non-RE and RE patients in serum samples at the baseline to identify biomarkers of interest to assist in the diagnosis, monitoring, and treatment of the disease. These biomarkers can be particularly useful to predict who will recover from T2DM following a dietary intervention.

The specific objectives of our study are to analyse the metabolomics differences before the dietary intervention and identify biomarkers of interest to assist in the prediction of T2DM recovery.

Study design and participants

This study was developed within the framework of the CORDIOPREV (CORonary Diet Intervention with Olive oil and cardiovascular PREVention) study, registered at Clinicaltrials.gov (number NCT00924937). This study is an ongoing controlled, single-blind, and randomized trial, with 1002 CHD patients. The trial protocol and subsequent revisions were approved by the Reina Sofia University Hospital Ethics Committee, following the Helsinki Declaration and good clinical practices. All patients signed a written informed consent to participate in the study.

Patients’ recruitment took place between November 2009 and February 2012, mostly at Reina Sofia University Hospital, Córdoba, Andalusia, Spain, with contributions from other hospitals in Córdoba and Jaen areas, in Andalusia, Spain. Complete details of the study methods, rationale, inclusion criteria, cardiovascular risk factors, and baseline characteristics are found elsewhere [13]. In brief, eligible participants, in the age range of 20 to 75 years, had established CHD with no clinical events in the previous 6 months. They all had at least a 5-year life expectancy and no other concurrent major diseases and were willing to participate in a long-term monitoring study [13].

In this work, 183 patients, from the CORDIOPREV study (https://www.cordioprev.es/index.php/es/) diagnosed with diabetes, underwent a dietary interventional study where participants were offered either an LF or MED diet for 5 years (Fig. 1). In our study, blood samples were taken during fasting (time 0) and 120 min after a glucose boost.

Oral glucose tolerance test

The patients underwent an OGTT at the baseline and, once a year, every year during the dietary intervention. Before the test, patients had fasted (from food/drugs) for 12 h and were asked to refrain from smoking and alcohol intake during the preceding 7 days. They were also asked to avoid strenuous physical activity, a day before the test. At 8:00 A.M., patients were admitted to the laboratory to perform the oral glucose tolerance test (OGTT) (75 g flavoured glucose load, Trutol 75; Custom Laboratories, Baltimore, MD, USA). Blood samples were taken at times corresponding to 0, 30, 60, 90, and 120 min to determine the glucose and insulin concentrations [14]. The insulin sensitivity index (ISI) was calculated from the OGTT using the following formula: ISI = 10.000 ÷ √([fasting plasma glucose X fasting plasma insulin] × [mean glucose in OGTT × mean insulin in OGTT]) [14]. HOMA-IR was calculated as described by Song et al. [15]. Beta-cell function was calculated using the disposition index (DI) as follows: DI = ISI × [AUC30 min insulin/AUC30 min glucose], where AUC30 min is the area under the curve between baseline and that at 30 min of the OGTT for insulin (pmol/l) and glucose (mmol/l) measurements, calculated by the trapezoidal method [16]. The indices used to determine tissue-specific insulin resistance (IR) were the hepatic insulin resistance index (HIRI, fasting plasma insulin × fasting plasma glucose) and the muscle insulin sensitivity index (MISI, (dG/dt)/mean of plasma insulin) [17]. Insulinogenic Index (IGI) was calculated by measuring plasma insulin at 30 min − fasting plasma insulin (mU/L)/(plasma glucose at 30 min − fasting plasma glucose(mg/dL) [18].

The adipose tissue (AT) insulin resistance index (Adipo-IR) was determined according to the formula: Adipo-IR = fasting plasma NEFA (mM) × fasting plasma insulin (pmol/L), which has been found as a suitable and useful method in clinical practice to estimate AT insulin sensitivity [19].

Randomization and masking

The process of randomization has been reported elsewhere [13]. Briefly, this is based on the following variables: sex (male, female), age (under and over 60 years old), and previous myocardial infarction (yes, no). Eight different groups were created to represent all the possible combinations of the above factors. Therefore, eight different blocks were created to assign the diets (bloc randomization). Dietitians were the only members of the intervention team to be aware of the dietary group of each participant.

Dietary assessment

The participants were randomized to consume two diets: the Med diet or an LF diet [13]. The LF diet consists of < 30% total fat (< 10% saturated fat, 12–14% MUFA fat, and 6–8% PUFA fat), 15% protein, and a minimum of 55% carbohydrates. The Med diet consists of a minimum of 35% of calories as fat (22% MUFA fat, 6% PUFA fat, and < 10% saturated fat), 15% proteins, and a maximum of 50% carbohydrates [20]. Neither energy restriction, nor physical activity was specifically encouraged. In both diets, the cholesterol content was adjusted to < 300 mg/d.

The Mediterranean and low-fat diets were designed to provide a wide variety of foods, including vegetables, fruit, cereals, potatoes, legumes, dairy products, meat, and fish. The participants in both intervention groups received the same intensive dietary counselling. The nutritionists administered personalized individual advice every 6 months. In addition, quarterly group education sessions were held with up to 20 participants per session; separate group sessions were performed every 3 months, and dietary counselling by phone was done every 2 months [20]. At the beginning of the study, and every 6 months afterwards, each patient had a face-to-face interview with a nutritionist to complete a 137-item semi-quantitative food frequency questionnaire (validated in Spain [21]). The dietary evaluation was calculated by the 14-item Med Diet Adherence Screener, which was used for measuring adherence to the Med diet [22]. Moreover, a 9-item dietary adherence screener was used to measure adherence to the LF diet guidelines. A more detailed report on dietary adherence has been published recently by our research group [20].

Diabetes remission criteria

Remission required the following: (i) the absence of glucose-lowering treatment and was defined by levels of HbA1c < 6·5%, (ii) a fasting plasma glucose < 126 mg/dl, and (iii) a 2-h plasma glucose in the 75 g OGTT < 200 mg/dl maintained for at least 2 years. This agrees with the American Diabetes Association (ADA) diagnosis criteria [23].

Sample preparation

Plasma samples (100 μL) were immersed in bath ice and treated with 300 μL of 3:1 (v/v) methanol–acetonitrile (MeOH–ACN). The treated samples were vortexed for 2 min and subsequently cooled at − 20 °C for 3 min. Centrifugation was carried out for 15 min at 4 °C and 13,800 × g in a thermostatic centrifuge Thermo Sorvall Legend Micro 21 R from Thermo (Thermo Fisher Scientific, Bremen, Germany), and the supernatant phase was isolated. This phase was dried by evaporation and reconstituted with 60 μL of 3:1 (v/v) MeOH–ACN. All samples were processed in a 1200 Series LC system (Agilent Technologies, Waldbronn, Germany) coupled to an Agilent 6530 high-resolution QTOF mass spectrometer equipped with a dual electrospray ionization source.

LC–QTOF MS/MS analysis

A Poroshell 120 EC-C18 column (50 mm × 2.1 mm i.d., 2.7 μm particle size, from Agilent), kept at 25 °C, was used to carry out the chromatographic division. The mobile phases consisted of (A) 0.1% formic acid in deionized water and (B) 0.1% formic acid in acetonitrile. The protocol used for the elution consisted of 0–2 min, 5% B; 2–11 min and the percentage of mobile phase B was modified from 0 to 100%. The final percentage was held for 6 min. Five minutes post-run was included to equilibrate the column. The flow rate was maintained at 0.4 mL/min. The injected sample volume was 5.0 μL and the injector needle was washed 10 times with 70% methanol between injections. Therefore, the needle seat was flushed for 15 s at a flow rate of 4 mL/min, with 70% methanol, to avoid cross-contamination between samples. The autosampler was kept at 4 °C to increase sample stability. The settings of the electrospray ionization source, which was operated in negative and positive ionization modes, were as follows: capillary voltage ± 3.5 kV, Q1 voltage 130 V, N2 pressure in the nebulizer 35 psi; N2 flow rate and temperature as drying gas 10 L min–1 and 325 °C, respectively. MS/MS data were acquired in both polarities, using the centroid mode at a rate of 2.5 spectra s–1 in extended dynamic range mode (2 GHz). Accurate mass spectra in the MS scan were acquired in the m/z range 40–1100 and the MS/MS mode in the m/z range 30–1100. The instrument gave a typical resolution of 18,000 full width at half maximum (FWHM) at m/z 118.0862 and 35,000 FWHM at m/z 922.0098. The instrument was calibrated and tuned as recommended by the manufacturer. To assure the desired mass resolution, continuous internal calibration was performed during analyses by using the signals at m/z 121.0509 (protonated purine) and m/z 922.0098 [protonated hexakis(1H,1H,3H-tetrafluoropropoxy) phosphazine or HP-921] in the positive ion mode, while in the negative ion mode, ions with m/z 119.0362 (proton abstracted purine) and m/z 966.0007 (formate adduct of HP-921) were used. The collision energy was set at 20 V for the whole run. The analytical samples were injected in auto MS/MS acquisition mode to obtain fragmentation information from a maximum of two precursors selected per cycle with an exclusion window of 0.1 min after 2 consecutive selections of the same precursor.

Data processing

The MassHunter Workstation software (version B7.00 Qualitative Analysis, Agilent Technologies, Santa Clara, CA, USA) was used to process all the data obtained by LC–QTOF in data-dependent acquisition MS/MS mode. Treatment of raw data files started with the extraction of potential molecular features (MFs) with the suited algorithm included in the software. For this aim, the extraction algorithm considered all ions exceeding 500 counts for both polarities with a maximum charge state of 2 for the obtained chromatograms. The count cut-off value was established considering the chromatographic background noise. Additionally, only MFs defined by two or more ions were considered, with a tolerance for the isotopic distribution of 0.0025 m/z for peak spacing tolerance, plus 7.0 ppm in mass accuracy. Only the following potential ions and adducts were considered in positive (H + , Na + , K + , NH4 +) and negative ionization (H − , HCOO + , Cl +) modes. Furthermore, a potential neutral loss by dehydration was also included to identify features corresponding to the same potential metabolite.

Identification of metabolites was supported on MS and MS/MS information by using METLIN MS and MS/MS databases (http://metlin.scripps.edu), the Human Metabolome Database (HMDB, 3.6 version), and the LIPID MAPS website ((http://www.lipidmaps.org); in all cases, the MFs obtained in the previous step were used. A database with all identified metabolites was used to perform a targeted compound extraction analysis using a tolerance window of 0.8 min and 5 ppm mass accuracy. This step was performed with Profinder Analysis (version B8.00, Agilent Technologies, Santa Clara, CA, USA). A table with the peak area of all identified compounds in the different samples injected was obtained as a result.

Statistical analysis

Metabolites showing in at least 80% of the samples were selected for further analysis. To allow predictive modelling, imputation was carried out when needed, substituting missing values by half the smallest value of the appropriate metabolite.

LC–MS data (polar and apolar) was imported into Matlab (R2015a, Mathworks UK) and analysed using the statistic toolbox and algorithms from Korrigan Toolbox version 0.1 (Korrigan Sciences Ltd, UK). In Matlab, matrices were log 10 normalized. The biostatistical pipeline for the multivariate statistical analysis considered a preliminary unsupervised principal component analysis (PCA), followed by a supervised pairwise O-PLS DA [24, 25], which identifies the specific modulations driven by the appropriate predictor (i.e. individuals who returned from T2DM versus those who did not).

To assess the predictive power of the O-PLS DA models, R² (the explained variance) was calculated. This parameter evaluates the model maximizing variance given by the endogenous variables. The Q², or goodness of prediction, assesses the predictive relevance of the model and is based on a matrix partition technique that ignores part of the data (in our case a seventh part each time), estimates the model parameters, and predicts the omitted parts using the estimates obtained previously. Q² greater than 0 means the model has predictive value. In addition, the overfitting of the model (the difference between R²Y and Q²Y) was also considered, and only models with less than 50% overfitting were further considered. Model parameters and associated metabolites were reported and used for a Cox proportional hazard model, GLM, and ROC calculations in R (version 4.0.5 (2021–03-31, https://www.r-project.org) using the packages ‘caret’ and ‘pROC’. Unadjusted Cox proportional-hazard models calculated the hazard ratio (HR) of every metabolite previously identified in the O-PLS DA model within a 95% confidence interval (CI). This unadjusted Cox allowed the identification of the betas for every metabolite. This information was used to calculate the patient’s likelihood of recovering from diabetes by running a Cox analysis adjusted for sex, age, body mass index (BMI), HDL, triglycerides, and intensity of statin therapy based on the tertiles resulting from the multiplication of the betas previously obtained by the abundance of each metabolite for every patient. Finally, generalized linear models were calculated for (i) all the clinical variables (i.e. sex, age, BMI, HDL, triglycerides, and intensity statin therapy), (ii) all the metabolites, (iii) the glycated haemoglobin, (iv) the clinical variables and the 12 metabolites, and (v) the clinical variables and the glycated haemoglobin. ROC analyses were carried out for these three models, and AUC, sensitivity, specificity, accuracy, and threshold were estimated for the models. Finally, DeLong analysis was used to compare whether the AUCs of these models were or not significantly different between them.