Longitudinal multi-omics of host-microbe dynamics in prediabetes

Abstract

Type 2 diabetes mellitus (T2D) is a growing health problem, but little is known about its early disease stages, its effects on biological processes or the transition to clinical T2D. To understand the earliest stages of T2D better, we obtained samples from 106 healthy individuals and individuals with prediabetes over approximately four years and performed deep profiling of transcriptomes, metabolomes, cytokines, and proteomes, as well as changes in the microbiome. This rich longitudinal data set revealed many insights: first, healthy profiles are distinct among individuals while displaying diverse patterns of intra- and/or inter-personal variability. Second, extensive host and microbial changes occur during respiratory viral infections and immunization, and immunization triggers potentially protective responses that are distinct from responses to respiratory viral infections. Moreover, during respiratory viral infections, insulin-resistant participants respond differently than insulin-sensitive participants. Third, global co-association analyses among the thousands of profiled molecules reveal specific host-microbe interactions that differ between insulin-resistant and insulin-sensitive individuals. Last, we identified early personal molecular signatures in one individual that preceded the onset of T2D, including the inflammation markers interleukin-1 receptor agonist (IL-1RA) and high-sensitivity C-reactive protein (CRP) paired with xenobiotic-induced immune signalling. Our study reveals insights into pathways and responses that differ between glucose-dysregulated and healthy individuals during health and disease and provides an open-access data resource to enable further research into healthy, prediabetic and T2D states.

Conflict of interest statement

M.S. is a cofounder of Personalis, Qbio, SensOmics, January, Filtricine and Akna and advisor for Genapsys. M.S. is also an inventor on provisional patent number 62/757,629: ‘Methods for evaluation and treatment of glycemic dysregulation and applications thereof. ’ S.M.S.-F.R., K.C., W.Z. and T. Mishra are also listed as inventors. All other authors declare no competing interests.

Figures

Fig. 1. Summary of study design, cohort details and data.

a, Samples were collected from 106 participants for nearly four years, with quarterly healthy visits and dense sampling during stress events. b, Summary of visits in each category and categorization of participants as insulin-sensitive or insulin-resistant. c, Sample sources, types of multi-omic assay and number for each data type profiled per visit. d, Conceptual illustration of the data structure for downstream analysis. Supplementary Table 3 lists the number of visits and assays profiled for each participant.

Fig. 2. Variances observed among health visits.

We characterized 624 healthy baselines. a, Top, intra-class correlation (ICC) levels in analytes from each ’ome (gut microbes, median 32.65; gut microbial genes, median 19.19, P = 7.261 × 10−15 by Wilcoxon-rank test, two-sided). Bottom, separation pattern by multidimensional scaling (MDS) among participants defined by the top 30 most personally distinct analytes. Fourteen participants are shown in different colours, with healthy visits presented as dots bound by contours denoting personal variable space. b, Expression of analytes significantly associated with SSPG in healthy baselines. Red, increased expression; blue, decreased expression. Pink, clinical laboratory tests; red, cytokines; grey, metabolites; orange, proteins; dark green, gut microbes at the genus level per row of the heatmap. HDL, high-density lipoprotein; TGL, triglycerides; NEUTAB, neutrophil absolute count; uncl., unclassified; WBC, white blood cell count; EGFR, estimated glomerular filtration rate; GROA, also as CXCL1, growth-regulated alpha protein. Undefined abbreviations are protein symbols.

Fig. 3. Dynamics of differential responses during RVI.

a, Top, examples of the most significantly enriched integrated canonical pathways that showed differential changes in response to RVI according to the AUC test (two-sided). All P values were corrected using the Benjamini–Hochberg (B–H) method for multiple hypothesis correction. The Ingenuity Pathway Analysis (IPA) Z-score measures the match between expected relationship direction and observed gene expression. The ratio presents the percentage of molecules enriched from the total number of molecules within a pathway. Bottom, a simplified example of enriched neuroinflammation signalling showing molecules that were significantly upregulated (red) or downregulated (green) in response to RVI. The full pathway is shown in Extended Data Fig. 4d. b, Two time-series clusters (left and middle) that show up- and downregulation of significant ’omics molecules during RVI (AUC test, two-sided) with their top enriched pathways (black), clinical laboratory tests (red), gut microbial communities (green), and nasal microbial communities (blue). Trend lines are colour-encoded with red shades denoting high membership values of genes belonging to the time-series cluster. Right, a general time-series pattern was observed regarding viral load in relation to both nasal and gut microbial changes. The y axis denotes the percentage of maximal abundance. ALKP, alkaline phosphatase; INSU, insulin; INSF, insulin, fasting; LDLHDL, low-density to high-density lipoprotein ratio; ALCRU, aluminum/creatinine ratio random, urine; BASO, basophil in percentage; RDW, red blood cell distribution width; EOS, eosinophil in percentage; CO2, carbon dioxide. Virus 1 and Virus 2 denote generically for different viral species that changed in individuals. c, Heatmap showing different responses to RVI in the insulin-sensitive (77 RVI and 62 healthy categorized time points) and insulin-resistant groups (33 RVI and 17 healthy categorized time points) at the integrated pathway level based on the AUC test (two-sided) over the entire course of RVI and paired t-test (two-sided) for stage-wise comparison. Left, B–H corrected pathway P value enrichment heatmap. Right, activation Z-score heatmap. Heatmap bottom annotation: AUC, EE, EL and RE represent analysis over the entire course of RVI (AUC), at EE, EL and RE stages as compared to healthy baseline.

Fig. 4. Temporal multi-omics responses to immunization.

a, Heatmap showing different responses during RVI (n = 156 categorized time points) and immunization (n = 117 categorized time points, AUC test, two-sided, over the entire course of events and paired t-test, two-sided, for stage-wise comparison), same positions/labels as in Fig. 3. b, The ROC curve for classifiers designed to discriminate RVI time points from healthy baselines at different discrimination thresholds in testing data based on an LR model. 95% confidence intervals are shown as coloured ribbons. AUC scores of ROC curves are listed on the right for classification accuracy.

Fig. 5. Correlational networks capture multi-omics association structures that differ between insulin-resistant and insulin-sensitive groups.

a, Gut microbial associations are different in insulin-sensitive from insulin-resistant participants (insulin-sensitive n = 190 healthy visits, insulin-resistant n = 184 healthy visits, as not all visits had stool sampling). For each gut microorganism genus, significant associations (by SparCC) to other genera are displayed with their correlation coefficients as dots and the total number of significant associations (q < 0.05) as the bar in the middle. b, Examples of microbial-cytokines correlations (by CLR+rmcorr) that are significant in insulin-sensitive but not insulin-resistant participants. Longitudinal measurements and the correlation trend line are coloured per individual and q values are indicated at the top right of each comparison.

Fig. 6. Personal features that precede the onset of T2D revealed by longitudinal tracking.

a, Left, scatter plot showing the number of outliers among thousands of molecules for each participant; right, the percentage of outliers contributed by each ome. b, Top, collection timeline for participant ZNED4XZ in our study relative to onset of T2D, with green arrows pointing to eight healthy baselines (H), purple for one stress visit (S), red for four antibiotics visits (Ax) and blue for one immunization visit (Im). Bottom, heatmap showing levels of selected immune cytokines together with A1C (%) and fasting glucose (mg dl−1) across time. Red denotes increased expression; blue represents decreased expression, with shades corresponding to levels.

Extended Data Fig. 1. Summary of visits sampled and technical normalization across the study.

a, Visit chart along the time axis illustrated for each subject. b, Number of healthy baselines across the study per subject. c, Number of event visits (top, RVI; bottom, immunization) since the onset of events. Visits with defined days from onset of events were counted. d, Comparison of before and after normalization of cytokine data in order to correct batch effects during the assay. The after-normalization data show overall reduced coefficient of variation (CV). Red line, 25% CV. An example is shown for correcting five samples from two batches. e, Comparison of before (left) and after (right) LOESS normalization of metabolite data using l-lysine as an example. l-Lysine signal drift with time was independently corrected by fitting a LOESS curve to the l-lysine signal measured in QCs. QCs were injected every ten biological samples and consisted of an equimolar mixture of 150 random samples from the study. Red, l-lysine signals from QCs; green, biological samples in HILIC ESI(+) MS mode.

Extended Data Fig. 2. Variance among healthy baselines (n = 624) observed for each measurement.

a, Examples of measurements with their normalized values and corresponding IQR either outlined by boxes or quantified by numbers below. b, Scatter plots show diverse association patterns between IQR and expression or abundance level for measurements across all healthy visits in different ’omes. c, Scatter plots of gut and nasal mean abundance versus corresponding ICC for each measurement across all healthy baselines. Each measurement is presented as one grey dot; darkness of grey reflects the degree of overlap of nearby dots. d, MDS plots for each ’ome. MDS1 and MDS2 were shown as the x and y axis, respectively. Only 14 subjects with at least ten healthy baselines were used for plotting to achieve clear visualization of the clustering patterns. Subjects were given the same colour code across plots. Individual separation scores (Ind_score) calculate the average distance between any two individuals across all analytes in each ’ome.

Extended Data Fig. 3. Factors that contributed to the variance among healthy baselines.

a, b, Examples of measurements correlated with the day factor in our study. Rmcorr correlation is used for the association calculations. a, The association between the level of ALKP and the day factor was significant both in all 106 subjects and in 27 subjects sampled for more than 900 days. b, The association between the level of pregnenolone sulfate and the day factor was significant only in 27 subjects sampled for more than 900 days but not when all 106 subjects were included. Individuals were coloured differently for comparison. c, The expression pattern of measurements that differed significantly between insulin-resistant (n = 35, red) and insulin-sensitive (n = 31, green) subjects in healthy baselines (Wilcoxon-rank test, two-sided). On the left: pink, clinical laboratory tests; grey, metabolites; dark green, gut microorganisms.

Extended Data Fig. 4. Differential responses to RVIs.

a, Comparison of AUC test (two-sided) performance with LR analysis with the time covariate. Left, distribution plot of P values resulting from the AUC test. Right, distribution plot of P values resulting from the LR analysis. The pie chart shows the number of differentially expressed genes (q < 0.1) for both methods and the overlapping number. b, Top, comparison of results from paired t-test (two-sided) with DESeq2 method (Wald test) for stage pairwise analyses (for example, EE versus personal healthy baseline). Seventy-one per cent of transcripts identified by DESeq2 method overlap with paired t-test results. Bottom, P value distribution for paired t-test (left) and DESeq2 (right). Both methods show a homogenous distribution of P values with enrichment for significant ones. Supplementary Table 39 lists molecules found significant by paired t-test or DESeq2 method. c, Illustration of IL10 pathway enriched in response to RVI. Multi-omic molecules that changed significantly during the course of RVI by AUC test are highlighted in red (upregulated) or green (downregulated). d, An example of enriched neuroinflammation signalling showing molecules significantly upregulated (red) or downregulated (green) in response to RVI. Simplified pathway graph is shown in Fig. 3a (bottom). e, Illustration of insulin signalling (left) and cardiac hypertrophy (right) pathways, which were enriched in response to RVI. Molecules that changed significantly during the course of RVI by AUC test are highlighted in dark pink (upregulated) or green (downregulated).

Extended Data Fig. 5. Trend analysis of differential changes in response to RVI.

n = 156 RVI categorized time points and n = 89 healthy time points based on AUC test, two-sided. a, Four time-course patterns clustered by differential ’omic molecules that changed in response to RVI. b, The elbow method showed the optimal number of clusters (K = 4). c, Example of co-occurring changes in nasal and stool bacteria during upper respiratory rhinovirus infection. In subject ZOZOW1T, a high load of human Rhinovirus was found in RVI EE stage, accompanied by an increase in Moraxella (80%) in nasal samples. The abundances of Rhinovirus and Moraxella in nasal samples decreased in RVI EL and continued to decrease during the RE stage. Rhinovirus could not be detected at the next healthy time point and Moraxella showed a relatively low abundance (5%). Bacteroides, the dominant genus in the stool samples, decreased by 40% from EE to EL and returned to the level of the previous healthy time point at the RE stage. Rhinovirus was identified by nasal transcriptome and quantified by mapping the transcriptome reads to the virome references database. Moraxella and Bacteroides in nasal and stool samples were identified and quantified by 16S rRNA gene sequencing. d, Acute phase response pathway representation upon RVI in the insulin-sensitive group (left; Z-score 3.9 at EE stage) and the insulin-resistant group (right; Z-score 1.2 at EE stage). Upregulated molecules are highlighted in red, downregulated ones in green.

Extended Data Fig. 6. Trend analysis of differential changes in response to RVI and immunization.

Data include 117 immunization (VAC) and 59 healthy categorized time points based on AUC test, two-sided. a, Four time-course patterns clustered by differential ’omic molecules that changed in response to immunization. b, The elbow method shows the optimal number of clusters (K = 4). c, Heatmap of immune molecules that changed significantly in response to RVI and immunization. The heatmap uses log2(baseline normalized read counts) to measure changes in expression in every category (EE, EL, RE) compared to the personal baseline (–H and +H). Red, upregulation; blue, downregulation. d, NF-κB pathway upregulation after RVI (Z-score = +1.96). e, NF-κB pathway downregulation upon immunization (Z-score = –4.9). Red, upregulated; green, downregulated.

Extended Data Fig. 7. Trend analysis of differential changes in response to respiratory viral infections (RVI, n = 156) and immunization (VAC, n = 117).

a, Heatmap of cytokine expression over the course of RVI and VAC. The heatmap uses log2(baseline normalized read counts) to measure changes in expression in every category (EE, EL, RE) compared to personal baseline (–H and +H). Red, upregulation; blue, downregulation. b–d, Heatmaps of nasal KO genes (b), nasal microbial taxa (c) and gut microbial taxa (d) over the course of RVI and VAC. The heatmaps use the log2(baseline normalized read counts) to measure changes in values in every category (EE, EL, RE) compared to the personal baseline (–H and +H).

Extended Data Fig. 8. Receiver operating characteristic (ROC) curves for classifiers designed to separate RVI time points from healthy baselines at different discrimination thresholds.

a, b, Performance in the training (top) and test sets (bottom) for LR (a) and SVM (b) models. Plots show the true positive rate (y axis) versus the false positive rate (x axis) for each ’ome and all ’omes combined. AUC scores for ROC curves are listed on the right for classification accuracy. For binary classification between healthy and RVI data, combined multi-omes had the highest prediction performance, followed by the metabolome, compared to others in the test cohort. c, ROC curves and AUC scores based on all-pairs testing using the LR model for classifying RVI events. The plot shows the true positive rate (y axis) versus the false positive rate (x axis) for each pairwise combination. Bottom right, summary of all pairwise combinations of multi-omes and their respective AUC scores in percentages.

Extended Data Fig. 9. ROC curves for classifiers designed to separate immunization time points from healthy baselines at different discrimination thresholds.

a, b, Performance in the training (top) and test sets (bottom) for LR (a) and SVM (b) models. Plots show the true positive rate (y axis) versus the false positive rate (x axis) for each ’ome and all ’omes combined. AUC scores of ROC curves are listed on the right for classification accuracy. For binary classification between healthy and immunization, combined multi-omes had the highest prediction performance, followed by metabolomes, compared to others in the test cohort. c, ROC curves and AUC scores based on all-pairs testing using the LR model for classifying immunization events. The plot shows the true positive rate (y axis) versus the false positive rate (x axis) for each pairwise combination. Bottom right, summary of all pairwise combinations of multi-omes and their respective AUC scores in percentages.

Extended Data Fig. 10. Associations among multi-omic molecules.

Insulin-resistant n = 215, insulin-sensitive n = 238, all n = 624 including baselines from individuals with unknown insulin sensitivity, Rmcorr correlation. a, Comparison of associations between A1C and fasting glucose level (GLU) for between-individual associations (at the cohort level, top) and within-individual associations (at individual level, bottom) in insulin-resistant and insulin-sensitive groups. Correlation coefficients and adjusted P values (B–H corrected) are shown in each comparison. For comparison, data from individuals are coloured differently, with a single colour for all time points from the same individual. At the cohort level, the median fasting plasma glucose per participant correlated positively with the median haemoglobin A1C level for both insulin-sensitive and insulin-resistant groups. However, at the individual level, there was no significant correlation. This is expected, as fasting glucose is very variable and can be affected by diet or lifestyle within hours, whereas A1C approximates average glucose levels over about 3 months. b, Associations with HSCRP are shown either at the cohort level (left, between-individual correlation) or at individual level (right, within-individual correlation) for comparison between insulin-resistant and insulin-sensitive individuals. Only significant associations (q < 0.05) in insulin-resistant or insulin-sensitive subjects at individual level are used as examples. c, Significant (by SparCC) association patterns between gut microorganisms at the genus level in insulin-sensitive (left) and insulin-resistant individuals (right, insulin-sensitive n = 190, insulin-resistant n = 184, as not all visits included stool sampling). Red, positive correlation; blue, inverse correlation. d, An example of gut microorganisms: genus Butyricimonas and its significant associations with host metabolites (by CLR+rmcorr) that are specific to insulin-sensitive or insulin-resistant groups. Black dots denote the correlation coefficients of the genus related to corresponding metabolites on the left, with their 95% confidence intervals shown as either green ribbons for insulin-sensitive-specific or red ribbons for insulin-resistant-specific associations. The P values of those associations are shown on the right.

Extended Data Fig. 11. Personal features reveal molecular markers before the onset of diseases.

a, Left, scatter plot showing the number of outliers among thousands of molecules for each participant; right: the percentage of outliers contributed by each ’ome. The top part of the plot is also shown in Fig. 6a. b) Z-scores of outlier molecules in participant ZJTKAE3. c, Metabolite–cytokine profiles of the cohort. All healthy baselines (H1–H8) of participant ZNED4XZ are highlighted in red. d, Correlational network among molecules that were highly associated with IL1RA and enriched in xenobiotic-induced immune response pathways. Metabolites, orange ovals; transcripts, blue rectangles; proteins, yellow squares. Positive associations, green lines; negative associations, red lines.