Introduction

Diabetes mellitus is considered to be a refractory disease causing a significant socio-economic burden. Type 2 diabetes mellitus (T2DM) is the dominant type of diabetes mellitus characterized by metabolic disorders, insulin resistance, and deficiency of insulin secretion (1, 2). The pathogenesis of this chronic disease is complex and genetic and environmental factors, such as sugar and lipid intake, gut microbiota, many metabolites, and even air pollutants, contribute to its increase in prevalence (3, 4). Accumulating evidence has linked gut microbiota with T2DM development by variety ways. The gut microbiota can impact the integrity of the intestinal epithelial barrier, mediate insulin resistance, as well as regulate the function of mitochondria (5, 6, 7). They can regulate local or systemic immunity and inflammation, which also contributes to the development of T2DM (8). Moreover, various gut microbial metabolites, such as short-chain fatty acids, bile acid, and tryptophan-derived metabolites, have been reported to be closely related to the pathogenesis of T2DM (5, 6, 9, 10). However, the interactions between gut microbiota and its host with T2DM have not yet been fully characterized. In pathological conditions, the dysregulation of the host can lead to changes in gut microbiota composition. In turn, the microbiota plays a regulatory role to participate in the development of T2DM. Despite the complicate interactions between host and microbiota in the context of T2DM, some studies suggest antidiabetic interventions targeting the gut microbiota such as fecal microbiota transplanting (FMT) can be applied as a clinical treatment of T2DM (11, 12, 13).

Meanwhile, dysfunction of lipid metabolism contributes to T2DM development by inducing lipotoxicity in humans and animal models (14, 15). The long-chain fatty acids (LCFAs) are the principal lipid components naturally occurring in animal fats and vegetable oil, as well as the main metabolites of fat. LCFAs such as palmitic acid (PA, C16:0), palmitoleic acid (C16:1N7), and oleic acid (C18:1N9), are reported to have strong association with T2DM. In the last few decades, there is increasing evidence that frequent consumption of LCFAs contributes to metabolic diseases such as obesity and T2DM due to the high PA content (16, 17). PA is a saturated fatty acid and its increase in serum is a significant contributing factor in T2DM development (18, 19, 20). The suggested mechanisms by which PA mediates T2DM include increasing diacylglycerol and ceramide synthesis (21), mitochondrial and endoplasmic reticulum stress (22, 23), and activation of pro-inflammatory pathways (24). Nevertheless, whether gut microbiota involving in PA mediated T2DM and the interactions between gut microbiota and LCFA metabolites in T2DM development are still unclear.

Spontaneous development of T2DM in non-human primates (e.g. macaques) is highly similar to human T2DM, such as insulin resistance in early stages and later abnormal glucose tolerance and T2DM development follows the same pathological changes of pancreatic islets and complications. In fact, T2DM macaques and humans have amyloidosis of pancreatic islets but it is not observed in mouse models (25, 26), suggesting that T2DM macaques are the best animal model that can simulate and reproduce human T2DM and its complications (27). However, previous studies indicated that naturally occurring spontaneous T2DM macaques in captive populations were rare even if individuals were given a high lipid and high sugar diet (28, 29). After nasally fed cynomolgus macaques with a high-fat dietary emulsion for 12 months, the macaques did not experience significant increases in fasting blood glucose and glycosylated hemoglobin (29). The controversial effects of high fat on the T2DM development are worthy of further investigations to better understand this complex disease.

Here, this study identified spontaneous development of T2DM in individuals (hereafter spontaneous T2DM macaques) from a large group of rhesus macaques. These spontaneous T2DM macaques have never been treated with anti-diabetic drugs and therefore provide valuable models for pathogenesis investigation of T2DM. Based on the macaque model, we used multi-omics techniques to address the interactions between gut microbiota, host gene expression, and fecal metabolites and the development of T2DM. Our results demonstrated that gut microbiota and LCFA metabolites played important roles in the pathogenesis of T2DM. We validated the increased content of plasma PA and activation of inflammation in the T2DM macaques. In addition, we successfully induced prediabetes in mice by transplanting fecal microbiota from T2DM macaques into mice in conjunction with a diet high in PA. We also revealed the specific structure of gut microbiota that promoted T2DM development by regulating the absorption of excess PA, which provides key evidence of the causative effect of gut microbiota on the pathogenesis of T2DM.

Results

Metagenome and 16S sequencing demonstrate alterations on gut microbiota in spontaneous T2DM macaques

We identified eight spontaneous T2DM macaques out of 1698 individuals from long-term glucose monitoring in a captive population (Table S1). The fasting plasma glucose (FPG), fasting plasma insulin (FPI), and HOMA-IR levels in the T2DM macaques were significantly higher than in the control group (p<0.01), suggesting insulin resistance in T2DM macaques. However, glycosylated hemoglobin A1c (HbA1c), triglycerides (TG), total cholesterol (TC), high-density lipoprotein cholesterol (HDL), low-density lipoprotein cholesterol (LDL), and BMI did not significantly differ from the controls (p>0.05, Table 1). Comparison of gut microbiota between the two groups based on 16S rRNA amplicon found that both Shannon and Simpson index decreased in T2DM macaques but it was not significantly different from controls (p>0.05; Figures 1A and B). PCoA indicated that the microbiota composition of spontaneous T2DM macaques was different from the control group (Figure 1C). As one of the most dominant families in both groups, Lachnospiraceae showed significantly higher abundance in the T2DM group (14.48%) than in the control group (6.66%). While the abundance of the Lactobacillaceae family was significantly greater in the control group (25.95%) compared to spontaneous T2DM macaques (20.88%) (Figure 1D). A total of 21 microbes were identified as differentially microbes between the T2DM group and control group, where ten microbes including Erysipelotrichaceae and five members in Lachnospiraceae family (Ruminococcus gnavus, Lachnospira, Coprococcus sp. Dorea longicatena, and Roseburia) were significantly greater in the T2DM group (Figure 1E).

The changes in gut microbiota in spontaneous T2DM macaques.

(A) Alpha diversity estimates (Shannon index) between T2DM and control groups (ns, not significant, two-tailed t-test, n=8).

(B) Alpha diversity estimates (Simpson index) between T2DM and control groups (ns, not significant, two-tailed t-test, n=8).

(C) Principal Coordinate Analysis (PCoA) (n=8).

(D) Differential analysis of gut microbial composition in T2DM and control groups (n=8).

(E) LEfSe analysis between T2DM and control groups (n=8).

(F) Differential analysis of gut microbial function in T2DM and control groups (n=5). The pathways with red color were associated with T2DM and inflammation. Error bar is mean with ± standard deviation (s.d.).

(G) Differential analysis of gut microbial CAZy enzyme in T2DM and control groups (n=5). CBMs: carbohydrate-binding module (p=0.022, two-tailed t-test); GTs: Glycosyl Transferases (p=0.013, two-tailed t-test); PLs: Polysaccharide Lyases (p=0.017, two-tailed t-test); AA: Auxiliary activity enzymes (ns, not significant, two-tailed t-test); GH: Glycoside hydrolases (ns, not significant, two-tailed t-test); CE: Carbohydrate esterases (p=0.039, two-tailed t-test). For all boxplots: centre lines, upper and lower bounds show median values, 25th and 75th quantiles; upper and lower whiskers show the largest and smallest non-outlier values. In c, ellipses represent the 95% confidence intervals.

Physiological and biochemical parameters of Control and T2DM group

The metagenome results showed that the T2DM group had a higher abundance of Erysipelorictchaeae, Eubacterium rectale, Lachnospiraceae, Negativicutes, Blautia, and Coriobacteriia than the control group (Figure S1A). Functional enrichment demonstrated a total of 74 KEGG pathways with significant differences between the two groups. These pathways were mainly associated with T2DM and inflammation, including type II diabetes mellitus, glucagon signaling pathway, starch and sucrose metabolism, lipid metabolism, and toll-like receptor signaling pathways (Figure 1F). Among the six CAZy enzyme families, the Carbohydrate-binding module (CBMs), Glycosyl Transferases (GTs), Carbohydrate esterases (CEs), and Polysaccharide Lyases (PLs) families were significantly upregulated in the T2DM group, indicating significant changes in carbohydrate metabolism (Figure 1G). The results of metagenome and 16S sequencing demonstrated significant alterations in the composition and function of gut microbiota in spontaneous T2DM macaques, with a greater abundance of microbes associated with T2DM and fewer beneficial microbes.

Fecal metabolome and blood transcriptome reveals dysfunction of fatty acid β oxidation and tryptophan metabolism in T2DM macaques

The UHPLC-MS based metabolome analysis on fecal samples of T2DM and control macaques identified 1564 metabolites belonging to various types of secondary metabolites, with lipids and lipid-like molecules being most abundant (31.1%) (Figure S1B). We found 64 significantly differential metabolites between the two groups using a combined multidimensional statistical analysis (OPLS-DA) and univariate statistical analysis (T-test) (VIP>1, p<0.05) (Figures 2A and B; Table S2). Among them, muscon, indole-3-acetaldehyde, and serotonin were significantly lower in the T2DM group (Figure 2C) and are associated with anti-inflammatory activity (30, 31, 32). Notably, 64 differential metabolites were significantly enriched in one pathway: tryptophan metabolism (Figure 2D), and we identified two significantly different metabolites, indole-3-acetaldehyde and serotonin in this pathway. Both metabolites are microbiota-derived tryptophan metabolites and the ligands for AhR (31, 32). The lower concentration of AhR ligands may lead to the development of inflammation and metabolic syndromes in humans (33). In addition, we found the contents of many acylcarnitine metabolites were significantly higher in the T2DM group (VIP>1, p<0.1), including l-propionylcarnitine, hexanoyl-l-carnitine, (r)-butyrylcarnitine, and isovaleryl-l-carnitine (Table S2). Among them, l-propionylcarnitine, a kind of acylcarnitines, was the most upregulated metabolite in the T2DM macaques (FC: 16.19) (Figure 2C). Acylcarnitines are products of incomplete oxidation of LCFAs, which can activate pro-inflammatory signaling pathways and ultimately inhibit insulin activity (34). Our results indicated incomplete LCFAs β oxidation in spontaneous T2DM macaques, while similar characteristics were also observed in insulin resistant and T2DM humans (35, 36).

The alterations of fecal metabolites and gene expression in spontaneous T2DM macaques.

(A) Orthogonal partial least squares discriminant analysis (OPLS-DA) score plots based on the metabolic profiles.

(B) Volcano plots of metabolomics (p<0.05, two-tailed t-test).

(C) Fecal metabolites with significant differences between T2DM and control groups (VIP>1, p<0.05, two-tailed t-test).

(D) Enrichment analysis of the differentially abundant pathways between T2DM and control groups (p<0.05, two-tailed t-test).

(E) Non-metric multidimensional scaling (NMDS) analysis between T2DM and control groups (p=0.019, two-tailed t-test).

(F) Volcano plots of DEGs (log fold change≥1, p<0.05, two-tailed t-test).

(G) The GO and KEGG pathway enrichment analyses (p<0.05, two-tailed t-test).

(H) Weighted Gene Co-Expression Network Analysis (WGCNA).

(I) Venn analysis between hub genes and DEGs. In A and E, ellipses represent the 95% confidence intervals. Data shown are from 8 individuals per group.

Blood transcriptome analysis was consistent with metabolome results, indicating dysfunction of fatty acid β oxidation and inflammation in T2DM macaques. Gene expression in T2DM macaques exhibited significant differences from the controls (Figure 2E), and a total of 161 differentially expressed genes (DEGs) (26 upregulated and 135 downregulated) were identified in T2DM macaques at a FDR level of 0.05 (Figure 2F). Enrichment analysis of the DEGs was linked to diabetes, fatty acid metabolism, and inflammation, such as diabetic cardiomyopathy, glucose homeostasis, fatty acid metabolic process, and chemokine production (Figure 2G). We also identified 26 differential enrichment pathways between the two groups by aggregate fold change (AFC), and most were associated with insulin resistance and inflammation, including insulin resistance, PI3K-Akt signaling pathway, bacterial invasion of epithelial cells, and NOD-like receptor signaling pathway (Figure S2A). As shown in the insulin resistance pathway, expression of IL6 and IRS1 were upregulated, while INSR was downregulated in T2DM macaques (Figure S2B). The WGCNA analysis identified four modules that were related to T2DM. Darkgreen module and brown module were significantly positively correlated with T2DM, while green module and darkred module were negatively correlated with T2DM (Figure 2H). Genes in the darkgreen module and brown module related to lipolysis and inflammation were significantly upregulated in T2DM macaques, and genes in the green module and darkred module related to fatty acid metabolism and insulin secretion were significantly downregulated (Figure S2C). With the cut-off (|kME>0.8|), a total of 546 genes in the four modules were identified as hub genes (Table S3). Of these 546, 102 hub genes were also DEGs (Figure 2I). Several of these genes have been reported as correlated with T2DM in humans, including IGF2BP2, LEPR, RAP1A, SESTRIN 3, and ITLN1 (37, 38, 39, 40, 41). In addition, ACSM3, HADHB, and EFHB are involved in fatty acid β oxidation and inflammation (42, 43, 44).

Validation of LCFAs accumulation and inflammation in T2DM macaques

Collectively, metabolome and transcriptome results indicated dysfunction in fatty acid β oxidation and tryptophan metabolism in T2DM macaques, which may lead to LCFA accumulation and inflammation. To support this conclusion, we performed targeted medium-and long-chain fatty acid mass spectrometry of plasma and examined serum inflammatory cytokines in the macaques. A total of 34 fatty acids were detected, and among the five types of fatty acids, the concentration of saturated fatty acid (SFA) was significantly greater in the T2DM group (p<0.05, Figure 3A), while other types of fatty acids were not significantly different between the two groups (p>0.05, Figures 3B-E). In particular, concentrations of PA, palmitoleic acid, and oleic acid were significantly higher in the T2DM group than control group (p<0.05 and VIP>1, Figures 3F and G). Increased content of PA, palmitoleic acid, and oleic acid was also found in human T2DM (45, 46). PA was an important metabolite mediating insulin resistance through three main mechanisms, being increased diacylglycerol and ceramide synthesis, mitochondrial and endoplasmic reticulum stress, and activation of pro-inflammatory pathways through membrane receptors (21, 22, 23, 24). Analysis on the serum inflammatory cytokines found that IL-1β was significantly higher in the T2DM group (p<0.05, Figure 3H), but TNF-α and IL-6 levels showed no significant difference between the two groups (p>0.05, Figures 3I and J). IL-1β is a major player in a variety of autoinflammatory diseases and a key promoter of T2DM systemic and tissue inflammation (47). Moreover, blood routine examination showed an increase of WBC number, NEU percentage and NEU number and a decrease of LYM percentage and LYM number, also indicating the inflammation in the T2DM macaques (p<0.05, Table 2).

LCFAs accumulation and inflammation occurred in spontaneous T2DM macaques.

(A-E) The contents of SFA (A, p=0.038), MUFA (B), PUFA (C), N3 (D), and N6 (E) in plasma (ns, not significant, two-tailed t-test).

(F) The univariate analysis by two-tailed t-test, error bar is mean with ± s.d.

(G) The multidimensional analysis by VIP value (VIP>1).

(H-J) The contents of serum inflammatory cytokines, including IL-1β (H, p=0.032), TNF-α (I) and IL-6 (J) (ns, not significant, two-tailed t-test).

(K) Correlation analysis between DEGs, differential metabolites, and differential microbes using Spearman rank correlation (|r|>0.5, adj p<0.05). For all boxplots: centre lines, upper and lower bounds show median values, 25th and 75th quantiles; upper and lower whiskers show the largest and smallest non-outlier values. Data shown are from 7 individuals per group.

blood routine examination of Control and T2DM group

To investigate effect of gut microbiota on the LCFAs accumulation and inflammation in T2DM macaques, we performed a correlation analysis among the DEGs, differential metabolites, and differential microbes using spearman rank correlation. Four differential microbes in Lachnospiraceae family (Coprococcus, Lachnospira, Roseburia and Dorea longicatena) were significantly associated with three differential metabolites of PA, palmitoleic acid and oleic acid (|r|>0.5, adj p<0.05), suggesting the participation of Lachnospiraceae microbes in LCFAs accumulation in T2DM macaques (Figure 3K). In addition, the bacteria in class Coriobacteriia was also associated with the three LCFAs (|r|>0.5, adj p<0.05, Figures S1C and D).

Fecal microbiota transplantation (FMT) with high content PA food induce prediabetes in mice

To determine the causative effect of gut microbiota and PA on T2DM development, we collected feces from the spontaneous T2DM macaques and performed fecal microbiota transplantation (FMT) in antibiotic-pretreated mice. Mice were either administrated by FMT (FT), fed with high concentration PA diet (PA), or were combined FMT and PA diet (FTPA). A control group was used and they were fed with normal commercial food and lacked FMT (Figure 4A). FPG monitoring found that FPG levels in the FTPA group and FT group increased continuously from day 60 (Figure 4B), while the control group maintained stable FPG levels throughout the 120 days. The FTPA group showed the highest FPG of 6.7 mmol/L at day 120, which was significantly higher than the control group and indicated FTPA-treated mice had developed prediabetes (Figure 4C). The oral glucose tolerance test (OGTT) showed obvious glucose intolerance in the FTPA group (Figures 4D and E), and the FPI values and insulin tolerance test (ITT) were both significantly elevated in the FTPA group, indicating insulin resistance (Figures 4F and G). Meanwhile, body weight (BW), TC, and TG significantly increased in the FTPA group (Figures 4H-J). While mice in the FT group and PA group had significant weight gain and mild insulin resistance, they did not exhibit significant glucose intolerance and no significant elevation in TC and TG levels compared to controls. Histopathological changes in the pancreas and liver of mice were investigated using H-E staining. Hepatocytes focal necrosis with inflammatory cell infiltration was commonly observed in the FTPA mice, but not frequent in hepatocytes in FT and PA groups (Figure 4K). Furthermore, decreased volume and area in pancreatic islets and inflammatory cell infiltration were detected in the FTPA mice (Figure 4L), while such pathological changes were not found in the control group.

The FMT and high PA diet mice developed pre-T2DM characteristics.

(A) Experimental scheme of FMT and high PA diet treatment.

(B-H) Metabolic analysis, including the trend of FPG within 120 days (B), FPG (C, p=0.0003), OGTT (D), AUC of OGTT (E, p=0.028), FPI (F, p=0.007), ITT (G), and body weight change (H) on day 120.

(I-J) The contents of TC (I, p=0.005) and TG (J, p=0.041) in serum on day 120.

(K and L) Representative H-E staining images of liver (K) and pancreas (L). For all boxplots: centre lines, upper and lower bounds show median values, 25th and 75th quantiles; upper and lower whiskers show the largest and smallest non-outlier values. Significance was determined using one-way ANOVA. In d, g, and h: *p<0.05, **p<0.01. Data shown are from 4-6 individuals per group.

Specific structure of gut microbiota mediates the absorption of excess PA in the ileum

Given that the content of plasma PA significantly increased in the T2DM macaques (Figures 3F and G), we examined PA content in the feces, ileum, and serum in mice to compare the prediabetes mice and the spontaneous T2DM macaques. PA content in the serum and ileum of the FTPA group was significantly higher than the control group (p<0.05, Figures 5A and B), but its content in feces was significantly lower than the control group (p<0.05, Figure 5C). This suggested that the absorption of PA was significantly enhanced in the ileum leading to the increase of PA in serum. To verify this conclusion, expression of the Cd36 gene, a gene involved in the uptake and oxidation of LCFAs (48), was examined in the ileum of mice. Ileums showed a significantly upregulated expression of Cd36 in the FTPA group compared to the control group (p<0.05, Figure 5D). In contrast, the level of IL-17A, a protein inhibiting the expression of Cd36 (49), was significantly reduced in the ileum of the FTPA group (p<0.05, Figure 5E). Interestingly, the PA group’s content of PA, Cd36 expression and IL-17A was not significantly different from the control group. This indicated that without FMT, the ileum could not absorb PA effectively even when fed in high concentrations. Consequently, gut microbiota mediated the absorption of excess PA in the ileum.

The PA accumulation required the specific gut microbiota.

(A-C) Total PA contents in serum (A, p=0.013), ileum (B, p=0.016) and feces (C, p=0.014) on day 120.

(D) Quantitative RT-PCR for Cd36 transcripts in ileum on day 120 (p=0.049).

(E) The content of IL-17A in ileum on day 120 (p=0.027). For all boxplots: centre lines, upper and lower bounds show median values, 25th and 75th quantiles; upper and lower whiskers show the largest and smallest non-outlier values. Significance was determined using one-way ANOVA. Data shown are from 3-4 individual macaques per group.

(F) NMDS analysis (p=0.001, one-way ANOVA), ellipses represent the 95% confidence intervals.

(G and H) LEfSe analysis between FTPA and control groups (G), FT and control groups (H). Data shown are from 4 individuals per group.

(I) Specific gut microbiota structure promoted the absorption of excess PA by regulating the expression of IL-17A and Cd36, leading to the LCFAs accumulation and insulin resistance.

We then compared diversity and composition of microbiota communities between the FTPA, FT, PA and control groups, and a fifth group consisting of the microbiota transplants (TP) from T2DM macaques. Shannon and Simpson indices were lower in the FTPA, FT and PA groups than the control group (Figures S3A and B). NMDS analysis also indicated a distinct microbiota composition from the control group compared to other groups (NMDS1). In addition, microbiota composition of TP from T2DM macaques was distinct from other groups but a little closer to the FTPA and FT groups (NMDS2, Figure 5F). In particular, the Lachnospiraceae family showed the higher abundance in the TP, FTPA and FT groups than in the PA and control groups (Figure S3C). The abundances of three members of microbes in Lachnospiraceae (R. gnavus and Coprococcus sp. and Clostridium) in the FTPA and FT groups gradually increased over time after FMT (Figure S3D). The FTPA and FT groups shared many differential microbes compared to the control group, such as the significantly upregulated R. gnavus and Coprococcus sp., and the significantly downregulated Christensenellaceae, F16, Treponema sp. and Fibrobacter succinogenes (Figures 5G and H). It is noteworthy that these microbes were also differential microbes between T2DM macaques and controls, and their abundances changed in the same way as in macaques (Figure 1E). However, mice in the PA group did not share differential microbes with the spontaneous T2DM macaques. Correspondingly, the change of serum PA content in the PA group was not significantly higher than the control group (Figure S3E). Our results suggested that the transplanted microbiota from spontaneous T2DM macaques, especially the increased abundance of R. gnavus and Coprococcus sp. and decreased abundance of Treponema, F. succinogenes, Christensenellaceae and F16, promoted the absorption of excess PA by regulating the expression of IL-17A and Cd36, leading to the LCFAs accumulation and insulin resistance (Figures 5I).

Discussion

With a multi-omics technology, this study comprehensively characterizes the gut microbiota, metabolites and gene expression of spontaneous T2DM macaques. The gut microbiota diversity in T2DM macaques decreased. In particular, the abundance of bacteria R. gnavus and Erysipelotrichaceae were upregulated while the abundance of Christensenellaceae was downregulated in the T2DM macaques. Metabolome results demonstrated a decrease of microbiota-derived tryptophan and anti-inflammatory metabolites, indicating that T2DM macaques were prone to inflammation. Notably, the accumulation of acylcarnitine metabolites, the suggested biomarkers for human T2DM (46), indicated incomplete mitochondrial LCFA β oxidation in T2DM macaques. Transcriptome results identified many DEGs linked to insulin resistance, fatty acid β oxidation and inflammation, including IGF2BP2, LEPR, RAP1A, SESTRIN 3, and ITLN1 that have also been reported in human T2DM (37, 38, 39, 40, 41). Combining the multi-omics results we revealed the complex pathological mechanisms in the spontaneous T2DM macaques (Figure 6), which is comparable to T2DM humans. Firstly, expression of genes related to lipolysis, fatty acid oxidation, LCFA accumulation, inflammation, and insulin seretion are dysregulated, promoting the development of T2DM (Figure 6A). This is particularly related to lipid metabolism, where the increase of lipolysis and the downregulated expression of fatty acid metabolism-related genes HADHB and ACSM3 cause the accumulation of acylcarnitine and LCFAs, which results in incomplete LCFA oxidation in T2DM macaques (Figure 6B). We also suggest that the decrease of Lactobacillus sp. cause the reduction of Aryl hydrocarbon receptor (AhR) ligands (serotonin and indole-3-acetaldehyde) given that Lactobacillus is producer of AhR ligands (50). This then promotes the expansion of the PA producer Erysipelotrichacea (51, 52) and ultimately leads to the increase of PA levels. The increase in the abundance of Erysipelotrichacea and R. gnavus also promotes the development of T2DM by activation of inflammation (53, 54, 55) (Figure 6C). The accumulation of PA was reported to lead to insulin resistance by affecting genes in the insulin signaling pathway (16). Our study demonstrates that the significantly changed expression of RAP1A, SESTRIN3, and IRS1 in PA-mTORC1-Akt pathway causes insulin resistance in T2DM macaques. Moreover, the increase of PA can promote the development of T2DM by upregulating the NF-κB signaling pathway (Figure 6D).

Integration of multi-omics results.

(A) Insulin resistance, fatty acid oxidation disorders, LCFAs accumulation and inflammation occurred in spontaneous T2DM macaques.

(B) Incomplete mitochondrial LCFAs β oxidation. The expression levels of fatty acid metabolism-related genes HADHB and ACSM3 were downregulated in spontaneous T2DM macaques, which could lead to accumulation of acylcarnitine, including l-propionylcarnitine, hexanoyl-l-carnitine, (r)-butyrylcarnitine, and isovaleryl-l-carnitine.

(C) Gut inflammation. The decrease of Lactobacillus sp. likely caused the reduction of serotonin and indole-3-acetaldehyde, which promotes the expansion of PA producer Erysipelotrichacea and ultimately led to PA accumulation. Both Erysipelotrichacea and Ruminococcus gnavus promote the development of inflammation. Accumulation of PA and inflammation are important factors in the development of T2DM.

(D) Accumulation of PA promoted the development of insulin resistance. In the PA-mTORC1-Akt pathway, the changes of RAP1A, SESTRIN3, and IRS1 expression promoted the development of insulin resistance in spontaneous T2DM macaques. The increase of PA promoted the development of T2DM by up-regulating the NF-κB signaling pathway.

We validated the observations from the multi-omics analysis finding the significantly higher inflammatory cytokines IL-1β and LCFA accumulation, especially significant PA accumulation in the T2DM macaques. And found LCFA metabolites were significantly correlated with bacteria in the Lachnospiraceae family. Numerous studies imply an association of gut microbiota with T2DM development. By transplanting gut microbiota from healthy individuals to T2DM individuals, symptoms such as insulin resistance or inflammation were improved (6, 56). However, there is no evidence to date suggesting gut microbiota have directly causative effects on T2DM development. Our study confirmed the causative effect of gut microbiota and PA on T2DM development by transplanting fecal microbiota from spontaneous T2DM macaques to antibiotic pretreated mice. We successfully induced prediabetes in mice after combining FMT administration and high PA ingestion. However, when the treatments were administered on their own, the mice did not develop prediabetes. We determined, for the first time, that gut microbiota mediated the absorption of excess PA in the ileum by quantitative examining PA contents in feces, ileum, and serum, and analysis of the expression of Cd36 and IL-17A level in the ileum. Most notably, this then resulted in accumulation of PA in the serum and finally led to T2DM development. Without the transplanting gut microbiota, the ileum could not absorb the PA effectively even at a high concentration of ingested PA. Our study highlights the essential roles of gut microbiota in T2DM development, and this is probable also the reason that previous studies have failed to induce T2DM in macaques because they have only used a high-fat diet (28, 29).

We then determined the specific gut microbiota structure that related to T2DM development in the prediabetes mice and spontaneous T2DM macaques. We found that the increased abundance of R. gnavus and Coprococcus sp., and the decreased abundance of Treponema, F. succinogenes, Christensenellaceae, and F16, were involved in the T2DM development. R. gnavus can promote insulin resistance by regulating the content of tryptamine/phenethylamine (57). Moreover, R. gnavus is a mucin-degrading microbe that leads to an increase of inflammation (58, 59, 60). The intestinal mucous layer is an important barrier separating intestinal tissue from microbiota, and microbiota composition plays a major role in affecting the integrity of intestinal mucous layer (58). The increase of R. gnavus suggested a higher risk of damage to the integrity of the intestinal mucous layer. This is further supported by our results of the lower abundance of Lactobacillus sp., which are AhR ligand producers, and the decreased content of AhR ligands (tryptophan microbial metabolites) in spontaneous T2DM macaques. The deficiency of AhR ligand reduced the production of intestinal mucus and increased the risk of microbial invasion, which in turn affected the immune cell differentiation and cytokine production (61, 62). The cytokine IL-17A is a regulator of the fatty acid transporter Cd36 and lipid absorption can be promoted by reducing the inhibition of Cd36 expression by IL-17A (54). The decreased IL-17A level and increased Cd36 expression in prediabetes mice indicated that the specific gut microbiota promoted the absorption of excess PA by disrupting the integrity of the intestinal mucous layer and regulating the expression of IL-17A and Cd36. In addition, the beneficial bacteria decreased in abundance in T2DM macaques and prediabetes mice. These beneficial bacteria, such as F. succinogenes, Christensenellaceae, and F16, protect the mucosal barrier and improved insulin resistance (63, 64, 65). We inferred the collectively effects of these gut microbes determined the absorption of excess PA from ileum to serum, which might contribute to the development of T2DM.

In conclusion, spontaneous T2DM macaques that have never been treated with diabetes-related drugs provide a valuable model for our understanding of the pathological characteristics and pathogenesis of T2DM. This study characterized changes in gene expression, metabolites, and gut microbiota levels of spontaneous T2DM macaques using multi-omics techniques. We found gut abnormal microbiota, tryptophan metabolism and fatty acid β oxidation disorders, inflammation, and PA accumulation. We also successfully induced prediabetes in mice by transplanting fecal microbiota from T2DM macaques into antibiotic pretreated mice fed a high PA diet. Our study confirms the causative effect of gut microbiota composition change and PA in T2DM development. The microbiota composition, specifically higher abundance of R. gnavus and Coprococcus sp. and lower abundance of Treponema, F. succinogenes, Christensenellaceae and F16, promoted the absorption of excess PA thus led to the development of T2DM. This study provides new insights into the interaction of gut microbiota and metabolites in the development of T2DM, which expands our understanding of the pathogenesis of this metabolic disease and may provide novel insights for the treatment of T2DM in the future.

Materials and methods

The screening of spontaneous T2DM macaques

The experimental macaques used in this study were all from Greenhouse Biotechnology Co., LTD (Meishan, China). We obtained eight spontaneous T2DM macaques with FPG≥7 mmol/L and eight heathy control macaques with FPG≤6.1 mmol/L (three consecutive detections, each detection interval of one month) from a population of 1698 captive macaques. None of these 16 screened macaques received any medical treatment for diabetes.

Sample collection

Each experimental macaque was kept in a single cage, and was fasted for 12 hours but had free access to drinking water. We obtained serum, plasma and whole blood for the detection of physiological and biochemical parameters, and metabolome and transcriptome analysis. Fecal samples were collected within ten minutes after deposition. During the sampling process, fecal samples were loaded into a 50 mL sterile centrifuge tube were stored at −80℃. We followed animal welfare guidelines throughout the sample collection process, and all observations and samplings were approved by the Sichuan University’s Animal Care Committee.

Physiological and biochemical parameters

In this study, FGP, HbA1c and FPI were detected by hexokinase method, high performance liquid chromatography and electro-chemiluminescence method, respectively. TC, TG, HDL and LDL were detected by automatic biochemical analyzer. HOMA-IR is one of the criteria for T2DM, calculated as: (FPI×FPG)/22.5.

Feces 16S rRNA amplicon sequencing and analysis

Total DNA from fecal samples was extracted using the QIAamp Fast DNA Stool Mini Kit. The V3-V4 region of the 16S rRNA gene was amplified using the 341F/806R primer (341F: 5’-CCTAYGGGRBGCASCAG -3’, 806R: 5’-GGACTACNNGGGTATCTAAT -3’). The purification and quantification of the amplified products were performed and followed by the sequencing library preparation with TruSeq Nano DNA LT Library Prep Kit (Illumina, USA). The library sequencing was performed on Illumina MiSeq platform and 250 bp paired-end reads were generated. Raw sequencing reads were merged using FLASH (VI.2.7, http://ccb.jhu.edu/software/FLASH/) (66) and analyzed by QIIME2 (version 2020.11.1) pipeline with default parameters. Reads after denoising by DADA2 were clustered into OTUs at 99% similarity threshold. Taxonomy of OTUs was assigned based on Greengenes reference database (67). The QIIME2 diversity plugin was used to calculate alpha diversity (68). Principal Coordinate Analysis (PCoA) was determined by using the R package vegan. Differentially microbes were determined using linear discriminant analysis effect size (LEfSe) (69).

Feces shotgun metagenome sequencing and analysis

Total DNA from each sample was extracted using Tiangen DNA Stool Mini Kit (Tiangen Biotech Co., Ltd., China). Metagenome sequencing was performed using the Illumina NovaSeq 6000 platform with a paired-end sequencing length of 150 bp. Trimmomatic was used for removing the adapters and low-quality raw reads based on a four-base wide sliding window, with average quality per base>20 and minimum length 90 bp (70). The rhesus macaque potential sequences were removed using Bowtie2 (71) with the reference genome (assembly Mmul_10). Taxonomy of remaining reads was assigned using Kraken2 (72) with the option “--use-mpa-style”. De novo assembly of remaining reads was performed using MEGAHIT (73) with the option “-m 0.90 --min-contig-len 300”. Prodigal (74) was used for gene prediction. The construction of non-redundant gene catalogue was performed with CD-HIT (75) with the option “-c 0.95 -aS 0.90”. Quantification of the non-redundant genes in each sample was performed using Salmon (76). The amino acid sequences translated from non-redundant genes were aligned (--id 80% --query-cover 70% --evalue 1e-5) by DIAMOND (77) in the Carbohydrate-Active enZYmes (CAZy) database (78). The annotation metabolic pathway was performed using HUMANn3 (79).

RNA sequencing and DEG analysis

Total RNA was extracted using PAXgene Blood RNA kit. The cDNA Library was constructed following the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) manual, and index was added to each sample for sample differentiation. cBot Cluster Generation System was used to cluster the sequences with the same index. Illumina Hiseq 2500 sequencing platform was used to obtain the paired-end sequencing reads (PE150). NGS QC Toolkit v2.3.3 (80) was used to obtain high quality reads (high-quality paired reads with more than 90% of bases with Q-value≥20 were retained). Processed reads of each sample were mapped to the macaque reference genome using HISAT2 v2.1.0 (Kim, et al. 2015). Each alignment output file was assembled into a separate transcriptome using StringTie v1.3.6 (81), resulting in a transcript GTF file. To obtain the expression value of TPM (transcripts per million) and raw read counts of each gene, transcript GTF file was used as the reference annotation file. Differential expression analysis was performed using DESeq2 R package (82). DEGs were screened according to Foldchange value (FC) and p value corrected by FDR (Benjamini-Hochberg method was selected) (adj p<0.05, |log2FC|>1).

Weighted Gene Co-Expression Network Analysis (WGCNA) (83) was used to analyse the correlation between genes and phenotypes. We selected an appropriate “soft thresholding power” using the “picksoftthreshold” function in the WGCNA package (v1.61). Next, “blockwiseModules” function was used to construct the co-expression matrix with the option “checkMissingData=TRUE, power=16, TOMType=’unsigned’, minModuleSize=30, maxBlockSize=6000, mergeCutHeight=0.25”. The modules with high correlation to T2DM phenotype and p<0.05 were selected for downstream analysis, and genes in these modules were used for GO and KEGG functional enrichment analysis. Hub genes were identified based on the module eigengene-based connectivity (kME), |kME>0.8| as the cut-off criteria.

Instead of identifying the differentially expressed genes within a pathway between the two groups, Aggregate Fold Change (AFC) calculated the average multiple change for each gene and defined the pathway score as the average difference multiple for all genes in that pathway. A null hypothesis test was performed using the pathway scores of the gene expression dataset, and the significance of each pathway was estimated by p-value (84). STRING (v11.0) provides a tool for functional enrichment analysis based on AFC. GO and KEGG enrichment analysis was performed using g: Profiler (85), AFC enrichment analysis was performed using STRING (https://www.string-db.org/cgi/input?sessionId=bfRWW6asD8S4&input_page_show_search=on).

Untargeted metabolomics processing

The 200 μL homogenized fecal sample was mixed with 800 μL cold methanol/acetonitrile (1:1, v/v) to remove the protein. The mixture was centrifuged for 15 min (14000g, 4°C) followed by the drying of the supernatant in a vacuum centrifuge. For LC-MS analysis, the samples were re-dissolved in 100 μL acetonitrile/water (1:1, v/v) solvent. LC-MS/MS analysis was performed using an UHPLC (1290 Infinity LC, Agilent Technologies) coupled to a quadrupole time-of-flight (AB Sciex TripleTOF 6600) in Shanghai Applied Protein Technology Co., Ltd. The samples were separated by Agilent 1290 Infinity LC ultra-high performance liquid chromatography (UHPLC) HILIC column. In both ESI positive and negative modes, the mobile phase contained A=25 mM ammonium acetate and 25 mM ammonium hydroxide in water and B=acetonitrile. The gradient was 85% B for 1 min and was linearly reduced to 65% in 11 min, and then was reduced to 40% in 0.1 min and kept for 4 min, and then increased to 85% in 0.1 min, with a 5 min re-equilibration period employed. AB Triple TOF 6600 mass spectrometer was used to collect the primary and secondary spectra of per sample. After separation by UHPLC, the samples were analyzed by Triple TOF 6600 mass spectrometer (AB SCIEX). Positive and negative modes of electrospray ionization (ESI) were respectively detected.

Targeted medium-and long-chain fatty acid metabolomics processing

A total amount of 100 μL plasma per sample was taken in 2 mL glass centrifuge tubes and 1mL chloroform methanol solution was added. After 30min of ultrasound, the supernatant was taken and 2mL of 1% sulfuric acid-methanol solution was added. The mixed solution was placed in a water bath at 80℃, for 30min, and then 1mL of N-hexane and 5 mL of pure water were added in turn. Next, 500 μL of supernatant was absorbed, 25 μL of methyl N-nonaconate was added and mixed. The final sample was detected by GC-MS. All samples were separated by Agilent DB-WAX capillary column (30 m×0.25 mm ID×0.25 μm) gas chromatography. Agilent 7890/5975C gas-mass spectrometer was used for mass spectrometry. The chromatographic peak area and retention time were extracted by MSD ChemStation software. The content of medium-and long-chain fatty acids was calculated by drawing a standard curve.

Metabolomics statistical processing

The screening of significant changed metabolites was performed using univariate and multidimensional analysis. Student’s t test was applied to determine the significance of differences between two groups. The variable importance in the projection (VIP) value of each variable was obtained from the Orthogonal Partial Least Squares Discriminant Analysis (OPLS-DA) model, which was used to indicate its contribution to the classification. The screening criteria were p<0.05 and VIP>1.

Animal treatments

Specific-pathogen-free C57BL/6J male mice at 6 weeks of age (vendor: Chengdu dossy experimental animals Co., LTD) were randomly assigned to four groups (FTPA, FT, PA and Control). All the mice lived in cages with the same conditions, including 12h light and 12h dark cycles, temperature 22-25 °C and humidity 40-60%.

Diets

HPAD was prepared by adding 12% PA to conventional forage. Both the conventional forage and the HPAD were sterile and the fresh forage was renewed three times a week. On day 0, the diets of FTPA-treated mice and PA-treated mice were switched to HPAD, FT-treated mice and control mice were still fed with conventional forage. Drinking water was sterile and renewed twice a day.

Transplant preparation and use

After single cage feeding, FPG detection and fecal collection were performed, and fecal samples of seven T2DM macaques were mixed for the preparation of transplants. The appropriate volume of diluent was added to the fecal sample (i.e. add 4 ml diluent per gram of feces) and the preparation of diluent can be found in (86). Sodium L-ascorbic acid and L-cysteine hydrochloride monohydrate were added to all suspensions at final concentrations of 5% (w/v) and 0.1% (w/v), respectively (The sterile diluent of control group was added with the same amount of reagent). The mixture was homogenized and filtered with a 200-mesh sterile mesh screen to remove large particles from the feces, and the filtrate was passed through 400 and 800 sterile mesh screens to remove undigested food and smaller particulate matter. The filtrate was divided into 50 ml centrifuge tubes, centrifuged at 600×g for 5 min, and the precipitation was discarded. Finally, the fecal supernatant was divided into new centrifuge tubes in equal parts (400 μL per tube) and frozen at −80°C. For use, the transplant was quickly thawed in a 37℃ water bath.

Fecal microbiota transplantation (FMT)

After one week of changed feeding regime, the FTPA and FT mice groups were pre-treated with 1 g/L neomycin sulfate, 1 g/L ampicillin and 1 g/L metronidazole in the drinking water for 14 days, and the control group was not treated. For FMT treatment, the gavage with 400 μL transplant, which were thawed ahead of time, was performed for 14 days in FTPA and FT-treated mice. At the same time, the gavage in control and PA-treated mice were performed with sterile diluent.

Metabolic measurements

Throughout the experiment, body weight and feces were collected every month, FPG detection was every half month. OGTT was performed on day 110 and ITT was performed on day 115. For OGTT, after a 12h overnight fast, oral glucose gavage (1.2 g/kg of 12% dextrose solution) was performed and followed by blood sample collection from the tail vein at 0, 15, 30, 60, 90 and 120 min. For ITT, after a 6h fast, intraperitoneal insulin injection was performed (0.75 U/kg, human regular insulin), followed by the blood samples collection from the tail vein at 0, 20, 40, 60, 90 and 120 min. On day 120, after FPG detection, blood collection and issues were collected for quantitative RT-PCR, H-E staining and ELISA.

Isolation of tissue for quantitative RT-PCR and ELISA

Tissue was homogenized for RNA extraction, after adjusting the final concentration of RNA, and DNA reverse transcription was performed. Quantitative RT-PCR was performed on CFX Connect Real-Time PCR Detection System (Bio-rad, USA) using SYBR Green (Table 3-5). PA, IL-1β, IL-6, TNF-α and IL-17A were detected using Jiangsu Meimian ELISA kit and followed operator instructions.

Primers of RT-PCR

RT-PCR reaction components

RT-PCR cycle procotol

H-E staining

The isolated livers and pancreas were dehydrated and embedded after fixation with formalin for 24 h. Paraffin-embedded livers and pancreas specimens were cut at a thickness of 3 μm. All sections were stained with hematoxylin then eosin, and finally microscopy and image acquisition were performed.

Statistical analysis

In this study, one-way ANOVA was used to determine statistical significance for comparisons of more than three groups, and for comparisons of two groups, two-tailed t-test was used. p values are represented on figures as follows: ns, not significant, ∗p<0.05, ∗∗p<0.01.

Data availability

The raw data of transcriptomes, metagenomes, and 16S rRNA have been submitted to the China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences with BioProject accession no. PRJCA021499. We provide four links that reviewers can use to access the raw data. Transcriptomes data: https://ngdc.cncb.ac.cn/gsa/s/4xAVJ3Lk (GSA: CRA013604). 16S rRNA data (macaque): https://ngdc.cncb.ac.cn/gsa/s/bsDz9qrN (GSA: CRA013637). 16S rRNA data (mouse): https://ngdc.cncb.ac.cn/gsa/s/1YcuC464 (GSA: CRA013638). metagenomes data: https://ngdc.cncb.ac.cn/gsa/s/3k38T849 (GSA: CRA013607). The identified untargeted metabolites were listed in the Table S4. The identified targeted medium-and long-chain fatty acid metabolites were listed in the Table S5.

Acknowledgements

This work was supported by the Science and Technology Foundation of Sichuan Province (2021YJ0136) and the National Natural Science Foundation of China (No. 32171607). Special thanks to Sichuan Green-house Biotech Co., Ltd. for the sample collection and thanks to Prof. Jinchuan Xing and Dr. Megan Price for revising the manuscript.

Additional information

Author contributions

X.L., G.H., Q.H.L., S.Z.Y. and C.J. collected the samples; X.L. and S.Z.Y. performed the bioinformatics analyses; X.L., Y.C.X., K.S., C.J., J.X.L. and L.Z. performed the experiments; X.L. wrote the manuscript; J.L., Z.L.H., Z.X.F. and B.S.Y. revised the manuscript; J.L. conceived and designed the experiments.

Supplementary Figures

Figure S1 Metagenome analysis of microbiota.

(A) Differential microbes screened by metagenome analysis (p<0.05, two-tailed t-test). T2DM: type 2 diabetes mellitus

(B) The proportion of all metabolites.

(C) Correlation analysis between differential untargeted metabolites and differential microbes (Spearman’s Rho).

(D) Correlation analysis between differential targeted metabolites and differential microbes (Spearman’s Rho), *p<0.05. Data shown are from 5 individuals per group.

Figure S2 Pathway enrichment analysis by AFC and WGCNA.

(A) Total 26 differential pathways were enriched (FDR<0.05).

(B) The changes of gene expression in insulin resistance pathway compared to the control group, red color illustrating the up-regulation of genes and blue showed the down-regulation of genes in T2DM group, and the darker the color of the genes, the greater the |log10FC| value.

(C) WGCNA functional enrichment. Data shown are from 8 individuals per group.

Figure S3 The changes in gut microbiota in FTPA, FT, and PA-treated mice.

(A) Alpha diversity estimates (Shannon index) (p>0.05, one-way ANOVA, n=4).

(B) Alpha diversity estimates (Simpson index) between T2DM and control groups (p>0.05, one-way ANOVA, n=4).

(C) Family level taxonomy and relative abundance of five groups.

(D) The changes of members of Lachnospiraceae family from day −14 to day 120.

(E) LEfSe analysis between PA and control groups. Data shown are from 4 individuals per group.

Supplementary Tables

Table S1. Information on the macaque samples.

Table S2. List of 64 differential metabolites.

Table S3. List of genes in four modules significantly correlated with T2DM.

Table S4. List of all identified untargeted metabolites.

Table S5. List of all targeted medium-and long-chain fatty acid metabolites.