Study participants and data collection
The study participants were recruited from 2014 to 2016 during their annual health check-ups at the University of Tokyo Hospital. The individuals included both male and female Japanese individuals aged from 20 to 75 years. The exclusion criteria were as follows: established diagnosis of diabetes, routine use of medications for diabetes and/or intestinal diseases, use of antibiotics within 2 weeks before sample collection and loss of 3 kg of body weight in the 3 months before sample collection. Written consent was obtained from the participants after a thorough explanation of the nature of the study at their health-checkups.
To normalize the participants’ clinical characteristics, we planned to recruit around 100 healthy individuals, 100 individuals with obesity (BMI ≥25, based on the Japanese definition) and 100 individuals with a prediabetic condition (FBG ≥110 mg dl−1 and/or HbA1c ≥6.0%) on the basis of their clinical data, and stopped recruiting when the number of participants almost reached the goal. The sample size was determined on the basis of previous metagenomics studies showing microbial signatures of diabetic patients5,6. We enrolled 112, 100 and 101 individuals for the normal, obese and prediabetic groups, respectively. The participants were provided with instructions to fast overnight before their visits, and all clinical information and blood samples were collected in the morning during their hospital visit. Blood samples were immediately centrifuged to collect plasma and then stored at −80 °C until the sample preparation and analysis. The participants were also instructed to collect faecal samples in the morning and were provided with guidance on how to collect and preserve faecal samples, along with a kit comprising a sampling tube and an ice pack. The faecal samples were then transported to the hospital either by refrigerated shipping or by the participants themselves. In both scenarios, the samples were delivered in a chilled state within 24 h after collection and stored at −80 °C until sample preparation and analysis. Consequently, 256 participants collected their faeces in the morning on the day of their hospital visit. As for the remaining participants, they collected their faeces in the morning between 2 days before and 7 days after their hospital visit, with the exception of 5 individuals who collected their faeces in the morning more than 7 days after their hospital visit, 2 individuals who reported collecting their faeces in the evening 1 day before their hospital visit, and 5 individuals who did not provide faecal samples. Moreover, two individuals withdrew from the study after enrolment. Thus, 306 individuals who underwent physical examination, laboratory tests, faecal sampling for faecal 16S rRNA pyrosequencing and metabolomic analyses, and plasma sampling for plasma metabolomic analyses were included for the analysis. Owing to the limited samples, faecal metagenomics data were available for 290 individuals; CAGE analysis data for 298 individuals; and plasma cytokine and insulin data for 282 individuals. The number of samples included in each analysis is described in the figure legends. The clinical study was approved by the institutional review board of RIKEN and The University of Tokyo and performed in accordance with the institutes’ guidelines.
Although we determined the criteria for enrolment, these criteria were not necessarily appropriate for subsequent analyses. For example, those in the prediabetes group were significantly leaner than those in the obese group (27.3 kg m−2 versus 25.2 kg m−2, P < 0.0001). Moreover, owing to the nature of the study participants (that is, those participated in regular health checkups), the blood glucose and HbA1c of the prediabetes group were significantly but only marginally higher than those of the obese group (FBG, 106 mg dl−1 versus 94 mg dl−1, P < 0.0001; and HbA1c, 6.2% versus 5.6%, P < 0.0001). We therefore reasoned that, in these subclinical conditions of diabetes, many metabolic traits may be overlapping between prediabetes and obesity groups and they do not necessarily capture their distinct features in metabolic and clinical continuums. This hinders us from distinguishing microbial and metabolomic characteristics directly related to human metabolic dysfunctions. We therefore considered that individual indices representing participants’ clinical conditions (that is, IR and MetS, as described below) may offer a better interpretation of the participants’ metabolic traits and data. Nevertheless, we observed consistent results even with the clinical criteria of obesity and prediabetes (Extended Data Fig. 2d).
Phenotypic outcomes
IR is defined as HOMA-IR ≥2.5, as has been set for the Japanese population13. Similarly, normal IS was defined as HOMA-IR ≤1.6. HOMA-IR is calculated using the following formula: fasting plasma insulin (μU ml−1) × fasting plasma glucose (mg dl−1)/405. HOMA-IR values could be calculated for 282 individuals only, owing to the limited data of plasma insulin in some participants. MetS is diagnosed according to the Japanese criteria44, which require an abdominal circumference of ≥85 cm for male and ≥90 cm for female individuals and at least two out of the following three clinical abnormalities: (1) dyslipidaemia, defined as triglyceride levels of ≥150 mg dl−1 and/or HDL-C levels of <40 mg dl−1; (2) elevated blood pressure, defined as systolic blood pressure of ≥130 mmHg and/or diastolic blood pressure of ≥85 mmHg; and (3) impaired fasting glucose, defined as FBG levels of ≥110 mg dl−1. Individuals who meet the criteria of abdominal circumference but only one clinical abnormality were defined as pre-MetS, as reported previously45.
Measurement of plasma cytokines
Plasma cytokines were measured using Human Adipokine Magnetic Bead Panel 2 (Millipore, HADK2MAG-61K) and Human Obesity Premixed Magnetic Luminex Performance Assay Kit (R&D, FCSTM08) according to the manufacturers’ instructions. Measurements below the lower detection limits were considered to be zero, and those above the upper detection limits were considered to be the highest values of analysed cytokines.
Preparation for faecal samples
Aliquots (5 g) of faeces were blended with 30 ml methanol and filtrated with 100 µm of mesh filter to remove food residue after vigorous vortexing. The filtrate was centrifuged at 15,000g for 10 min at 4 °C and the supernatant (methanol extract) was used for metabolomics analysis. DNA of the faecal microbiome was extracted from the pellet.
Extraction and measurement for hydrophilic metabolites of faecal and plasma samples
We followed the extraction process and gas chromatography-tandem MS (GC–MS/MS) measurement methods for water-soluble metabolites described previously46 with some modifications. In brief, a 10 μl aliquot of plasma was added to 150 μl methanol, 125 μl Milli-Q water, 15 μl internal standard solution (1 mM 2-isopropylmalic acid) and 60 μl CHCl3. For faecal samples, a 25 μl aliquot of methanol extract was added to 125 μl methanol, 150 μl Milli-Q water containing internal standard (100 μM 2-isopropylmalic acid) and 60 μl CHCl3. The solution was shaken at 1,200 rpm for 30 min at 37 °C. After centrifugation at 16,000g for 5 min at room temperature, 250 μl of the supernatant was transferred to a new tube and 200 μl of Milli-Q water was added. After mixing, the solution was centrifuged at 16,000g for 5 min at room temperature, and 250 μl of the supernatant was transferred to a new tube. The samples were evaporated dry using a vacuum evaporator for 20 min at 40 °C and lyophilized using a freeze dryer. Dried extracts were derivatized with 40 μl of 20 mg ml−1 methoxyamine hydrochloride (Sigma-Aldrich) dissolved in pyridine and shaken at 1,200 rpm for 90 min at 30 °C. The solution was then mixed with 20 μl of N-methyl-N-trimethylsilyl-trifluoroacetamide (MSTFA, GL Science) and incubated for 30 min at 37 °C with shaking at 1,200 rpm. After derivatization, the samples were centrifuged at 16,000g for 5 min at room temperature, and the supernatant was transferred to a glass vial. The analysis was performed using a GC–MS/MS platform on the Shimadzu GCMS-TQ8030 triple quadrupole mass spectrometer (Shimadzu) with a capillary column (BPX5, SGE Analytical Science). The GC oven was programmed as follows: 60 °C held for 2 min, increased to 330 °C (15 °C min−1), and finally 330 °C held for 3.45 min. GC was operated in constant linear velocity mode set to 39 cm s−1. The detector and injector temperatures were 200 °C and 250 °C, respectively. Injection volume was set at 1 μl with a split ratio of 1:30.
We followed the SCFA extraction and GC–MS/MS measurement methods as previously described47 with some modifications. A 90 μl aliquot of plasma was added to 10 μl Milli-Q water containing internal standards (2 mM [1,2-13C2]acetate, 2 mM [2H7]butyrate and 2 mM crotonate). For faecal samples, a 25 μl aliquot of methanol extract was added to 10 μl Milli-Q water containing internal standards and then centrifugally concentrated at 40 °C and reconstituted with 100 μl of Milli-Q water. Then, 50 μl of hydrochloric acid (HCl) and 200 μl of diethyl ether were added to the solution and mixed well. After centrifugation at 3,000g for 10 min, 80 μl of the organic layer was transferred to a glass vial and 16 μl N-tert-butyldimethylsilyl-N-trifluoroacetamide (MTBSTFA, Sigma-Aldrich) was added to derivatize the samples. The vials were incubated at 80 °C for 20 min and allowed to stand for 48 h before injection. The analysis was performed using a Shimadzu GCMS-TQ8030 triple quadrupole mass spectrometer with a capillary column (BPX5). The GC oven was programmed as follows: 60 °C held for 3 min, increased to 130 °C (8 °C min−1), increased to 330 °C (30 °C min−1) and finally 330 °C held for 3 min. The detector and injector temperatures were 230 °C and 250 °C, respectively. GC was operated in constant linear velocity mode set to 40 cm s−1. Injection volume was set at 1 μl with a split ratio of 1:30. The data were processed and concentration was calculated by LabSolutions Insight (Shimadzu).
Overall, 195 and 100 metabolites in the faecal and plasma samples, respectively, were detected by our GC–MS/MS platform and passed quality control. The values below the limit of detection were replaced with zero. Consequently, 110 faecal and 88 plasma metabolites that were detected (that is, above zero) in more than 75% of participants were included in subsequent analyses, for which they were combined into a common analysis pipeline and defined as hydrophilic metabolites.
Lipidomics of faecal and plasma samples
The lipidomics analysis was performed according to a previously reported study20. Methanol, isopropanol, chloroform and acetonitrile of liquid chromatography (LC)–MS grade were purchased from Wako. Ammonium acetate and EDTA were purchased from Wako and Dojindo, respectively. Milli-Q water was purchased from Millipore (Merck). EquiSPLASH was purchased from Avanti Polar Lipids. Palmitic acid-d3 and stearic acid-d3 were purchased from Olbracht Serdary Research Laboratories.
For plasma lipid extraction, an aliquot of 20 μl of human plasma sample was added to 200 μl of methanol containing 5 μl of EquiSPLASH, 10 μM palmitic acid-d3 and 10 μM stearic acid-d3, and vortexed for 10 s. Then, 100 μl of chloroform was added and vortexed for 10 s. After incubation for 2 h at room temperature, the solvent tube was centrifuged at 2,000g for 10 min at 20 °C. A total of 200 μl of supernatant was transferred to an LC–MS vial (Agilent Technologies). For faecal lipid extraction, 50 μl of the methanol extract was added to 145 μl of methanol containing 5 μl of EquiSPLASH, 10 μM palmitic acid-d3 and 10 μM stearic acid-d3 in a 2 ml glass tube, and vortexed for 10 s. Then, 100 μl of chloroform was added and vortexed for 10 s. After incubation for 1 h at room temperature, 20 μl of water was added and vortexed for 10 s. After 10 min incubation at room temperature, the solvent was centrifuged at 2,000g for 10 min at 4 °C, and the supernatant was transferred to the LC–MS vial. All of the samples were divided into four batches for plasma analyses and five batches for faecal analyses, with 70–80 and 55–60 samples per batch after randomization, respectively. For each batch, a series of samples was prepared, and subsequent LC–MS/MS measurements were performed. A quality control sample was prepared by mixing the same volume of plasma from the first batch subjects. A procedure blank was prepared by using the same volume of water instead of a biological sample. The blank sample was analysed at the beginning and the end of each analysis batch, and the quality-control sample was injected every ten study samples.
The LC system consisted of a Waters Acquity UPLC system. Lipids were separated on an Acquity UPLC Peptide BEH C18 column (50 × 2.1 mm; 1.7 µm) (Waters). The column was maintained at 45 °C at a flow rate of 0.3 ml min−1. The mobile phases consisted of (A) 1:1:3 (v/v/v) acetonitrile:methanol:water with ammonium acetate (5 mM) and 10 nM EDTA; and (B) 100% isopropanol with ammonium acetate (5 mM) and 10 nM EDTA. A sample volume of 0.5−3 µl, depending biological samples, was used for the injection. The separation was conducted under the following gradient: 0 min, 0% B; 1 min, 0% B; 5 min, 40% B; 7.5 min, 64% B; 12 min, 64% B; 12.5 min, 82.5% B; 19 min, 85% B; 20 min, 95% B; 20.1 min, 0% B; and 25 min, 0% B. The sample temperature was maintained at 4 °C.
MS detection of lipids was performed on a quadrupole/time-of-flight mass spectrometer TripleTOF 6600 (SCIEX). All analyses were performed in high-resolution mode in MS1 (~35,000 full width at half-maximum) and the high sensitivity mode (~20,000 full width at half-maximum) in MS2. Data-dependent MS/MS acquisition (DDA) was used. The parameters were MS1 and MS2 mass ranges, m/z 70–1,250; MS1 accumulation time, 250 ms; MS2 accumulation time, 100 ms; collision energy, +40/−42 eV; collision energy spread, 15 eV; cycle time, 1,300 ms; curtain gas, 30; ion source gas 1, 40(+)/50(−); ion source gas 2, 80(+)/50(−); temperature, 250 °C(+)/300 °C(−); ion spray voltage floating, +5.5/−4.5 kV; declustering potential, 80 V. The other DDA parameters were dependent product ion scan number, 16; intensity threshold, 100 cps; exclusion time of precursor ion, 0 s; mass tolerance, 20 ppm; ignore peaks, within m/z 200; and dynamic background subtraction, true. The mass calibration was automatically performed using an APCI positive/negative calibration solution through a calibration delivery system.
MS-DIAL (v.4.48)20,48 was used with the following parameters: (data collection) retention time begin, 1.0 min; retention time end, 18 min; MS1 and MS2 mass range begin, 0 Da; MS1 and MS2 mass range end, 2,000 Da; MS1 tolerance, 0.01 Da; MS2 tolerance, 0.025 Da; (peak detection) minimum peak height, 3,000 amplitude; mass slice width, 0.1 Da; smoothing method, linear weighted moving average; smoothing level, 3 scan; minimum peak width, 5 scan; exclusion mass list, none; (identification) retention time tolerance, 1.5 min; MS1 accurate mass tolerance, 0.01 Da; MS2 accurate mass tolerance, 0.05 Da; identification score cut off, 70%; all lipid subclasses were used as the search space; (alignment) retention time tolerance 0.15 min; MS1 tolerance, 0.015 Da. The default values were used for other parameters. In faecal lipidomics, a total of 48,790 and 20,367 chromatographic peaks were detected in positive- and negative-ion mode data, respectively. Of these, 2,654 unique lipid molecules were annotated and semi-quantified in the MS-DIAL software program and used for further statistical analyses. Likewise, in plasma lipidomics, 1,469 and 2,167 chromatographic peaks were detected in positive- and negative-ion mode data, respectively, and 635 unique lipid molecules were annotated and semi-quantified. The semi-quantitative value of lipids was calculated by the internal standards according to the previous study20. The abbreviations of lipids are listed in Supplementary Table 27. Details of lipid subclass characterization follow the previous study20.
Co-abundance clustering of metabolites
To generate co-abundance clusters, 110 hydrophilic metabolites and 2,654 lipid metabolites detected in more than 75% of participants were included. These metabolites were clustered based on their co-abundance using the R package WGCNA49 (v.1.72-1). The following parameters were used for the analysis. For hydrophilic metabolites, soft thresholding β = 12, minimum cluster size = 3, deep split = 4, cut height = 0.9999, PAM clustering = F. For lipid metabolites, soft thresholding β = 14, minimum cluster size = 20, deep split = 4, cut height = 0.999, PAM clustering = F. As soft thresholding of WGCNA was not able to cluster all of the metabolites, the remaining metabolites that did not fit the criteria were subsequently clustered on the basis of biweight midcorrelation. The following parameters were used for the secondary clustering. For hydrophilic metabolites, minimum cluster size = 3, deep split = 4, cut height = 0.9999, PAM clustering = F. For lipid metabolites, minimum cluster size = 6, deep split = 4, cut height = 0.999, PAM clustering = F. The clusters with biweight midcorrelation above 0.8 were merged. The first principal component (PC1) of each cluster was calculated using the moduleEigengenes command of WGCNA and used as the representative value of the cluster for further analyses. The representative classes of the clusters were described in Supplementary Tables 2 and 3. KEGG pathway enrichment analysis of CAGs was performed on MetaboAnalyst (v.5.0)50 using 84 metabolite sets based on the KEGG pathway. Hypergeometric test and false-discovery rate (FDR)-adjusted P values were used to test significance. The enrichment ratio was calculated as the ratio of actual metabolite number to the expected value in each pathway.
Reanalysis of publicly available metabolomic data
To validate the associations between clinical markers and faecal metabolites, we used the metabolomic data of TwinsUK17 and HMP2 (ref. 18). The metabolome data of the TwinsUK cohort included 1,116 metabolites including 36 carbohydrates. The median (interquartile range) of age and BMI were 65 years (60–71 years) and 25.4 (22.8–28.8), and the proportion of males was 6.6%. As reported previously17, the metabolite levels were scaled by run-day medians. The data were then log-transformed and scaled. For regression analyses, we filtered out the metabolites detected in less than 50% of participants; as a result, 759 metabolites including 29 carbohydrates were used for further analyses. The record of BMI and HOMA-IR were used for phenotypic outcomes. For BMI, we retrieved 786 samples measured on the same day of faecal collection. For HOMA-IR, plasma glucose and insulin obtained in the same year of the faecal collection were used for the following calculation: plasma glucose (mM) × insulin (pM)/6.945/22.5. We identified 550 individuals who underwent both faecal collection and glucose and insulin measurement in the same year and included them in the analysis. The HMP2 data were obtained from the Inflammatory Bowel Disease Multi’omics Database (https://ibdmdb.org/). Among the 26 out of 106 samples from non-IBD control, BMI data were available for 20 samples. We further excluded four individuals aged <10 years. As HMP2 is a longitudinal study, only the first faecal sampling for metabolomics was used for the current analysis to avoid redundancy. The intensity of fructose, glucose and/or galactose was log-transformed and scaled.
DNA extraction from faecal samples
DNA extraction was performed according to a protocol described previously47 with slight modifications. Before DNA extraction, the faecal pellet was washed once with PBS and suspended in a 10 mM Tris-HCl/20 mM EDTA buffer (pH 8.0). Lysozyme (Sigma-Aldrich), achromopeptidase (Wako) and proteinase K (Merck) were subsequently added to the samples for cell lysis. DNA was recovered by a phenol–chloroform extraction method. To purify the extracted DNA, RNA was digested with RNase (Nippon Gene). DNA was then precipitated in a solution containing polyethylene glycol 6000 (Hampton Research). The DNA concentration was quantified using Quant-iT PicoGreen (Thermo Fisher Scientific).
16S rRNA gene sequencing and taxonomic assignment
The hypervariable V1–V2 region of the 16S rRNA gene was amplified by PCR using barcoded primers. PCR amplicons were purified using AMPure XP magnetic purification beads (Beckman Coulter), and quantified using the Quant-iT PicoGreen dsDNA Assay Kit (Life Technologies Japan). Equal amounts of each PCR amplicon were mixed and then sequenced using the MiSeq (Illumina) system.
On the basis of sample-specific barcodes, reads were assigned to each sample using bcl2fastq. Next, the reads lacking both forward and reverse primer sequences were removed using BLAST and parasail followed by trimming of both primer sequences. Data were further denoised by removing reads with average quality values of <25 and possible chimeric sequences. Reads with BLAST match lengths of <90% with the representative sequence in the 16S databases (described below) were considered to be chimeras and were removed. The filter-passed reads were used for further analysis. The 16S database was constructed from three publicly available databases: the Ribosomal Database Project (RDP; v.10.27), CORE (http://microbiome.osu.edu/) and a reference genome sequence database obtained from the NCBI FTP site (ftp://ftp.ncbi.nih.gov/genbank/, December 2011).
Operational taxonomic unit (OTU) clustering and UniFrac analysis from the filter-passed reads, 3,000 high-quality reads per sample were randomly chosen. All reads (the number of samples × 3,000) were then sorted according to their average quality value and grouped into OTUs using UCLUST (http://www.drive5.com/) with a sequence-identity threshold of 97%. The representative sequences of the generated OTUs were processed for homology search against the databases mentioned above using the GLSEARCH program for taxonomic assignments. For assignment at the phylum, genus and species levels, sequence similarity thresholds of 70%, 94% and 97% were applied, respectively.
Shotgun metagenomic sequencing
Metagenome shotgun libraries (insert size of 500 bp) were prepared using the TruSeq Nano DNA kit (Illumina) and sequenced on the Illumina NovaSeq platform. After quality filtering, reads mapped to the human genome (HG19) or the phiX bacteriophage genome were removed. For each individual, the filter-passed NovaSeq reads were assembled using MEGAHIT (v.1.2.4). Prodigal (v.2.6.3) was used to predict protein-coding genes (≥100 bp) in the contigs (≥500 bp) and singletons (≥300 bp). Finally, 6,458,217 non-redundant genes were identified in the 290 samples by clustering the predicted genes using CD-HIT with a 95% nucleotide identity and 90% length coverage cut-off. Functional assignment of the non-redundant genes was performed using DIAMOND (e-value ≤ 1.0 × 10−5) against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (release 2019-10-07) to obtain the KEGG orthologues. The genes with the best hit correlating to eukaryotic genes were excluded from further analysis.
Quantification of annotated genes in human gut microbiomes
For taxonomic assignment of metagenomic reads, 1 million filter-passed reads were processed for mOTU analysis (v.2.5.1)51 to obtain the relative abundance at the species level. To functionally annotate the predicted genes, 1 million filter-passed metagenomic reads per individual were mapped to the combined reference gene set consisting of non-redundant genes identified in this study, JPGM52 and IGC53 using Bowtie2 with a 95% identity cut-off. Multi-mapped reads, that is, the reads that mapped to multiple genes with identical scores, were normalized to the proportion of the number of other reads that uniquely mapped to these genes, according to a strategy outlined in a previous report52. The proportion of KEGG orthologues was calculated from the number of reads mapped to them. For the enrichment analysis of KEGG pathways, the significantly and positively (negatively) associated KEGG orthologue gave +1 (−1) for all of the upstream pathways linked to the KEGG orthologue, and the points were summarized as the ratio to the number of KEGG orthologues in the pathway. For the KEGG-orthologue-level analyses of PTS, those including ‘phosphotransferase system (PTS)’ in the KEGG pathway (02060) were selected for the following correlation analyses. In the analyses of glucosidases, ‘glycoside hydrolases’ defined in the CAZy database on the basis of EC numbers54 were selected. We further selected those included in ‘starch and sucrose metabolism’ in the KEGG pathway (00500). We defined intracellular glucosidase by their substrate described in the KEGG pathway map; those cleave phosphorylated carbohydrates were recognized as intracellular, and the rest of the genes were recognized to possess extracellular enzymatic activities. The pathways were further summarized into carbohydrate metabolism (09101), amino acid metabolism (09105), lipid metabolism (09103) and membrane transport (09131) on the basis of the KEGG Orthology database.
Comparison of KEGG organism genomes
The list of KEGG organisms used for this genome analysis is listed in Supplementary Table 28. All KEGG organisms from genera Alistipes, Bacteroides, Flavonifractor, Blautia, Dorea and Coprococcus, which showed the top three positive or negative correlations with faecal carbohydrates in Fig. 2d, were selected for this analysis. The lists of genes involving the ‘starch and sucrose metabolism’ pathway (00500) in these KEGG organisms were extracted using the R package KEGGREST (v.1.32.0). The representative protein sequences of Blautia hydrogenotrophica strain 2789STDY5608857 (taxonomy ID 53443), Dorea longicatena strain 2789STDY5608851 (taxonomy ID 88431) and Dorea formicigenerans strain ATCC 27755 (taxonomy ID 411461) were downloaded from the NCBI Datasets (https://www.ncbi.nlm.nih.gov/datasets/genomes/). KEGG annotation of these protein sequence files was performed using BlastKOALA (https://www.kegg.jp/blastkoala/) with ‘Bacteria’ used as the taxonomy group. The presence of KEGG orthologues relating to extracellular glycoside hydrolases in starch and sucrose metabolism pathways shown in grey in Extended Data Fig. 6f was summarized.
RNA extraction from PBMC
Blood samples were collected in Vacutainer CPT tubes (Becton Dickinson) and mixed with the anticoagulant by gently inverting the tubes 8 to 10 times. After centrifugation of the blood for 30 min at 1,500g, PBMCs were isolated as a diffuse layer above the gel. The plasma was removed, and the PBMCs were collected in conical tubes with 500 μl RNAlater (Thermo Fisher Scientific). The conical tubes were centrifuged at 1,000g at room temperature for 3 min to pellet the cells and the supernatant was discarded. The RNA was then isolated using the Maxwell 16 LEV simplyRNA Blood Kit (Promega) according to the manufacturer’s instructions. The quality of the RNA was assessed using Bioanalyzer (Agilent), as recommended by the manufacturer. The RNAs were quantified using the GloMax plate reader (Promega) and Quant-iT RiboGreen RNA Assay Kit (Thermo Fisher Scientific).
CAGE analysis
The CAGE libraries were constructed according to the dual-index nanoCAGE protocol, a template-switching-based variation of the standard CAGE protocol designed for low quantities of RNA55,56. cDNA libraries were prepared with RNA extracted from PBMC samples and sequenced on the Illumina HiSeq 2000 (50 bp paired-end). The sequenced reads were processed with the MOIRAI pipeline57: low quality and rDNA reads were first removed, then the remaining reads were mapped to the human genome version hg38 patch 1 using BWA v.0.5.9 (r16). The mapped reads were overlapped with the FANTOM5 robust promoter set (http://fantom.gsc.riken.jp/5/datafiles/latest/extra/CAGE_peaks/) and mapped to the nearest GENCODE v.27 annotations within (500 bp)58,59. The mapped reads falling under each FANTOM5 CAGE cluster were summed to produce the raw expression counts. Expression counts were converted to counts per million (CPM), and CAGE clusters expressed in less than 100 samples with at least 1 CPM and greater than 1 sample with at least 10 CPM were removed from further analysis. For each sample, the richness index was calculated using the R package vegan’s rarefy function with a subsample size of 100 on the filtered raw counts. Samples with a read library size of less than 1,000,000 and a number of unique CAGE clusters of <11,000 and richness less than 44 were removed as outliers, with the thresholds selected from visual inspection of the respective distributions. Cell type specificities of promoters of interest were determined using the FANTOM5 hg38 human promoterome view.11 in the ZENBU genome browser (https://fantom.gsc.riken.jp/zenbu/). Top-hit cells for analysed promoters were described. For cell-type gene set enrichment analysis of genes significantly associated with IR, annotated genes were analysed using Enrichr60,61 and the Human Gene Atlas database60, and the results of cell types with Padj < 0.05 were selected. The Enrichr combined score is defined as: c = log[p] z, where c is the combined score, p is the P value based on Fisher’s exact test and z is the z-score60.
Metabolite measurement in bacterial culture
The following strains were used for this culture analysis: A. indistinctus (JCM16068), A. finegoldii (JCM16770), Alistipes putredinis (JCM 16772), B. thetaiotaomicron (JCM 5827), Bacteroides xylanisolvens (JCM 15633), Bacteroides ovatus (JCM 5824), Bacteroides caccae (JCM 9498), Parabacteroides merdae (JCM 9497), Parabacteroides distasonis (JCM 5825), D. formicigenerans (JCM 31256), D. longicatena (JCM 11232), B. hydrogenotrophica (JCM 14656), Blautia producta (BP, JCM 1471), Coprococcus comes (JCM 31264), Faecalibacterium prausnitzii (JCM 31915), Flavonifractor plautii (JCM 32125), Clostridium spiroforme (JCM1432), Coriobacterium glomerans (JCM 10262), Roseburia hominis (JCM 17582), Adlercreutzia equolifaciens subsp. equolifaciens (JCM 14793), Eggerthella lenta (JCM 9979) and Collinsella aerofaciens (JCM 10188). All strains were obtained from RIKEN BioResource Research Center. All of the strains were cultivated in EG medium (JCM Medium No. 14) supplemented with 5% Fildes extract prepared by pepsin-digested horse blood instead of horse blood itself. For measurement of metabolites in bacterial culture, 60 μl of the bacterial culture grown in the EG medium was inoculated into 3 ml of the experiment medium (EG medium) and cultivated for 24 h. The samples were centrifuged, and the cell-free supernatant was collected for analysis. GC–MS was performed to measure hydrophilic metabolites as described above. We identified 261 metabolites by the analysis and used 198 metabolites observed in at least 30% of samples for the following analyses.
Animal experiments
C57BL6/N male mice aged 6 weeks were purchased from CLEA Japan. They were randomly assigned to either the control or treatment group and housed in a conventional animal facility of Yokohama City University. The mice were fed Quick Fat (CLEA Japan) for 3 weeks before bacterial administration and continued to be fed for 3 weeks during bacterial challenges. A. indistinctus (JCM16068), A. finegoldii (JCM16770), B. thetaiotaomicron (JCM 5827), B. xylanisolvens (JCM 15633), P. merdae (JCM 9497), F. prausnitzii (JCM 31915) and C. spiroforme (JCM1432) were used to broadly compare the efficacy of bacterial administration on the animal model. These strains were cultivated in EG medium overnight, and the concentration was adjusted to 2.5 × 108 CFU per ml by PBS. The bacteria and PBS, a negative control, were orally administered to the mice at a dose of 200 μl per mouse. The bacteria and PBS as the vehicle control were provided 3 times a week for 3 or 4 weeks. Body mass was measured before oral gavage. Postprandial blood glucose measurement and insulin tolerance test were performed 3 weeks after the initiation of bacterial challenges. After the insulin tolerance test, the mice were subjected to 5 h fasting before insulin injection, and 0.85 U kg−1 human regular insulin (Eli Lilly) was subsequently administered intraperitoneally. The intraperitoneal glucose tolerance test was performed 4 weeks after the initiation of bacterial challenges. The mice were subjected to 5 h fasting before glucose infusion, and 2.0 g per kg glucose (Nacalai Tesque) was administered intraperitoneally. In both experiments, the blood glucose was collected from the tail vein and serially measured using GLUCOCARD G Black (Arkray). For the necropsy, the mice were euthanized by isoflurane (MSD), and the fat mass of perigonadal and mesenteric fats was measured. Blood was drawn through cardiac puncture after the anaesthesia. HDL-C (Wako), triglycerides (Wako) and adiponectin (Otsuka) were measured in accordance with the manufacturers’ instructions. The Yokohama City University animal facility is maintained under a 12 h–12 h light–dark cycle at 24 ± 1.5 °C and 55 ± 10% humidity.
To assess the metabolism, dietary intake and locomotor activity of mice, C57BL6/N male mice at the age of 6 weeks were purchased from CLEA Japan and were maintained in a vinyl isolator of SPF animal facility at RIKEN Yokohama branch. Using the same experimental protocol in the conventional condition, the mice were fed Quick Fat (CLEA Japan) for 3 weeks before bacterial administration and continued to be fed the diet during bacterial challenges and metabolic measurement. We gave three oral gavages of A. indistinctus or PBS (vehicle control) every other day and then placed the mice individually in acrylic cages. We further gave one oral gavage 2 days after the start of individual housing. Their metabolic activity, dietary intake and physical activity were subsequently monitored. There was no significant difference in body mass at the start of metabolic measurement (mean ± s.d. of body mass were 25.7 ± 2.6 g and 26.1 ± 1.4 g in the vehicle and A. indistinctus groups, respectively). Oxygen and carbon dioxide concentration was measured using the ARCO-2000 system, an open-circuit metabolic gas analysis system with a mass spectrometer (Arco Systems). VO2, VCO2, energy expenditure, fat oxidation rate, carbohydrate oxidation rate and respiratory quotient were calculated within the system. Dietary intake and physical activity were simultaneously monitored through ACTIMO-100M and MFD-100M (Shinfactory). The differences in diurnal variation were tested using two-way mixed ANOVA, and P values for interactions between time and group were reported. The RIKEN animal facility is maintained under a 12 h–12 h light–dark cycle at 23 ± 2 °C and 50 ± 10% humidity. The sample size was determined on the basis of our preliminary experiments. Bacterial administration and body mass measurements were performed by an independent researcher who was not involved in the grouping and outcome assessments. All experimental procedures were approved by the Institutional Animal Care and Use Committee of the RIKEN and Yokohama City University and performed in accordance with the institutes’ guidelines.
Western blot analysis of phosphorylated AKT
To analyse phosphorylation of AKT (p-AKT) at Ser473, the mice administered with A. indistinctus, A. finegoldii and PBS (vehicle control) 3 times a week for 4 weeks were subjected to 6 h fasting before the insulin injection, and 0.85 U kg−1 human regular insulin (Eli Lilly) was subsequently administered from the inferior vena cava. The liver, epididymal fat (eWAT) and gastrocnemius muscle were subsequently collected 5 min after the insulin injection, weighed and snap-frozen by liquid nitrogen. To prepare the lysates for western blotting, the tissues were homogenized in buffer A (25 mM Tris-HCl, pH 7.4, 10 mM sodium orthovanadate, 10 mM sodium pyrophosphate, 100 mM sodium fluoride, 10 mM EDTA, 10 mM EGTA and 1 mM phenylmethylsulfonyl fluoride). Thereafter, the lysates were resolved on 10% SDS–PAGE. Phosphorylated or total protein of AKT was isolated by immunoblotting using specific antibodies after the tissue lysates were resolved by SDS–PAGE and transferred to a Hybond-P PVDF transfer membrane (Amersham Biosciences). Bound antibodies were detected with HRP-conjugated secondary antibodies using ECL detection reagents (Amersham Biosciences). Rabbit polyclonal antibodies directed against AKT and p-AKT (Ser473) were purchased from Cell Signaling Technology. Precision Plus Protein All Blue Standards (Bio-Rad) were used for the molecular mass markers.
Hyperinsulinaemic–euglycaemic clamp test
The protocol has been published elsewhere62,63. Mice administered with A. indistinctus or PBS (vehicle control) for 5 to 6 weeks were used for the experiment. Jugular vein catheterization was performed 1 day before the clamp test. In brief, a mouse was anaesthetized with isoflurane (MSD), and the right jugular vein was exposed. A double-channel catheter was subsequently inserted to the vein. The next day, the mice were subjected to 4 h fasting before the clamp test. Human regular insulin (Eli Lilly) was intravenously administered at 7.5 mU kg−1 min−1, and the blood glucose levels were monitored every 5 min for 120 min. 50% glucose solution containing 6,6-d2-glucose (Sigma-Aldrich) was simultaneously infused to keep blood glucose levels around 100 to 120 mg dl−1. To separate the plasma, approximately 25 μl of blood was also drawn from tail vein at 0, 90, 105 and 120 min, placed into a tube containing 2 μl of heparin (Mochida Pharmaceutical) and centrifuged at 12,000g at 4 °C for 5 min. The plasma levels of glucose and 6,6-d2-glucose were measured using GC–MS. In brief, a 5 µl aliquot of plasma was extracted and derivatized with methoxyamine hydrochloride (Sigma-Aldrich) and N-methyl-N-(trimethylsilyl)trifluoroacetamide (GL Sciences), as previously described46. The analysis was performed using a GC–MS/MS platform on a Shimadzu GCMS-TQ8040 triple quadrupole mass spectrometer (Shimadzu) with a capillary column (BPX5) (SGE Analytical Science/Trajan Scientific and Medical). The programme of GC–MS/MS analysis was previously described46 with minor modifications. We integrated each derivative peak to obtain mass isotopomers of glucose for the following ions: m/z 319.1, 320.1 and 321.1. The glucose infusion rate was determined as the infusion rate at 90, 105 and 120 min. The rate of glucose disappearance was determined on the basis of the plasma levels of 6,6-d2-glucose and total glucose using a non-steady-state equation as described previously63,64 and considered as the whole-body glucose disposal after insulin stimulation. Hepatic glucose production was determined as the subtraction of glucose disappearance rate and glucose infusion rate.
Analysis of triglyceride contents in the liver
For the necropsy, the mice were anesthetized using isoflurane (MSD), and the left half of liver was dissected, weighed and frozen in liquid nitrogen. The extraction of triglyceride contents from the liver tissue has been reported elsewhere62,64. In brief, the samples were homogenized in buffer A (25 mM Tris-HCl at pH 7.4, 10 mM sodium orthovanadate, 10 mM sodium pyrophosphate, 100 mM sodium fluoride, 10 mM EDTA, 10 mM EGTA and 1 mM phenylmethylsulfonyl fluoride) and mixed with chloroform/methanol (2:1, v/v). The mixture was shaken for 15 min, centrifuged and the organic layer was collected. The extraction step was repeated three times. The collected samples were evaporated and resuspended in 1% Triton X-100/ethanol. The triglyceride content was assessed using Triglyceride E-test Wako (Wako) according to the manufacturer’s instructions.
Statistical methods and comparisons
For general statistical comparisons, two-sided Wilcoxon rank-sum tests were used for two-group comparisons, Kruskal–Wallis tests followed by Dunn’s post hoc analysis were used for comparisons of more than two groups, and Fisher’s exact tests were used for comparison of categorical variables. For general correlation analyses, Spearman’s rank correlation in the function corr.test of the R package psych v.2.1.6 was used. For partial correlation analyses, partial Spearman’s rank correlation in the function pcor.test of the R package ppcor v.1.1 was used. To predict the metabolite levels and their CAGs (Fig. 1b,d and Extended Data Figs. 2c,d and 3a), rank-based regression analyses were performed using the function rfit of the R package Rfit (v.0.24.2)65. For the ordinal independent variables (that is, IR, MetS, and original categories with obese and prediabetes), IS, no MetS, and healthy categories were considered as the references, respectively, and the coefficients and P values for other categories were calculated against these reference categories. For the analyses involving generalized linear models (GLM) such as Fig. 2b and Extended Data Figs. 5i and 6c, the dependent variables were assumed to follow a Gamma distribution and arcsine square root transformation was applied to the relative-abundance values of microbiota and KEGG orthologues. To enhance comparability, the standardized coefficient was also calculated by standard deviations of dependent and independent variables using the function lm.beta of the R package QuantPsyc v.1.5 in Extended Data Fig. 5i. In the reanalysis of TwinsUK data, we fitted generalized linear mixed-effects models with age, sex, zygosity and BMI as fixed effects and sample collection year as a random effect using the function glmer of R package lme4 v.1.1-27.1 to estimate the associations between HOMA-IR and faecal carbohydrate metabolites (Extended Data Fig. 3b,c). Similarly, in the reanalysis of HMP2 data, we fitted a generalized linear mixed-effects model with consent age and sex as fixed effects and sample collection site as a random effect to estimate the associations between BMI and faecal fructose, glucose and/or galactose (Extended Data Fig. 3d). To analyse the associations between the participants’ clusters and clinical markers in Fig. 2b, the clusters were reordered before regression analyses according to their proportion of individuals with IR, where cluster C showing the lowest proportion of IR was set as the reference. To calculate the KEGG pathway enrichment associated with the participant clusters (Fig. 2e and Extended Data Fig. 6a,b), the KEGG orthologues were compared between each cluster and the remaining three clusters using a two-sided Wilcoxon rank-sum test, and significant (Padj < 0.05) KEGG orthologues were summarized into the pathway level (Supplementary Table 16). For comparison of metabolites in bacterial cultures (Extended Data Fig. 8), one-way ANOVA followed by Tukey’s post hoc test was performed, followed by multiple testing corrections based on the Benjamini–Hochberg procedure. For comparisons of time-series data such as insulin tolerance test, two-way repeated-measures ANOVA was used and the between-group difference was analysed by estimated marginal means. P < 0.05 was considered to be significant. To analyse the body mass change in animal experiments, ANCOVA analysis was performed to adjust baseline body mass (that is, body mass change as a dependent variable and group and baseline body mass as independent variables). We also validated the assumption of this ANCOVA model, that is, homogeneity of regression slopes, homogeneity of variances and normality of residuals. For multiple-testing corrections, P values were corrected using the Benjamini–Hochberg procedure using the R function p.adjust. Padj < 0.05 was used as a cut-off unless otherwise specified. All data were collected using Microsoft Excel 2016. All statistical and graphical analyses were conducted using R v.4.1.1 using R studio v.1.4.1717, unless otherwise specified.
ROC curve analysis of omics datasets
To analyse ROC curves of omics datasets, the datasets of faecal metabolomics, including hydrophilic and lipid metabolites, faecal 16S rRNA gene sequencing at the genus level, faecal metagenome consisting of KEGG orthologues and clinical metadata, were included. We first selected feature variables in each dataset, that is, the best explaining variables in the given dataset, using the minimum redundancy maximum relevance (mRMR) algorithm15. The function mRMR.classic of the R package mRMRe v.2.1.2.1 was used for the calculation. The datasets were square-root-transformed before mRMR calculation. We selected 5 to 50 variables in 5 increments as the maximum number of genera was 50. Using the selected variables, we next established random-forest models using the R package caret v.6.0-88 to classify the individuals into IR or not. Specifically, the results of mRMR were split into train and test datasets in a 3:1 ratio. The generated random-forest models were evaluated using a tenfold cross-validation method and applied to the test datasets to obtain probability scores. The accuracy of each classification model was described by the AUC of ROC curves using the R package pROC v.1.17.0.1.
Construction of microorganism–metabolite networks
To construct the co-abundance networks of genus-level bacteria, we selected 28 genus-level microorganisms that were observed in more than 40% of the participants and calculated the correlations using the R package CCREPE (compositionality corrected by renormalization and permutation)66 v.1.28.0 with Spearman’s correlations and the default settings. Interactions with Padj < 0.05 were selected for further analysis. Bacteria that exhibited a positive correlation with one another were determined to be members of an independent co-abundance microbial group, except for the interaction between Bacteroides and Robinsoniella. We decided to categorize Robinsoniella into the Blautia and Dorea group owing to its stronger correlation with Blautia in comparison to Bacteroides, both of which showed the highest centrality within their respective networks. Those weakly associated with each other or negatively associated with the members of other CAGs were classified as miscellaneous (Extended Data Fig. 5f). To characterize the microbial profiles of the study participants, the individuals were clustered on the basis of the abundance of 28 genera, which includes 20 genera in co-abundance microbial groups identified with CCREPE and 8 unclustered genera, using the ward.D function of the R package pheatmap v.1.0.12. Four distinct clusters of participants were determined, and the proportion of IR was compared using Fisher’s exact tests. Microorganism–metabolite networks were constructed on the basis of the correlations between the 28 genera observed in at least 40% of samples and the faecal metabolites, including all hydrophilic metabolites (n = 110) and bacteria-related lipid metabolites (n = 259). Bacteria-related metabolites were defined according to previous reports20,21. The following classes were selected: DGDG, PE-Cer, MGDG O, FAHFA, Cer-AS, Cer-BDS, NAGly, NAGlySer, PI-Cer, SL, AcylCer, bile acids, DGDG O and AAHFA. Positive and negative Spearman’s correlations with Padj < 0.05 were separately depicted in the networks. The networks were visualized using Cytoscape (v.3.7.0)67.
Construction of cross-omics networks
To construct and visualize a correlation-based network of omics data, we first analysed IR-associated host signatures using plasma cytokines, plasma metabolites and CAGE promoter expression data. We identified the significant host markers through the following models: (1) GLM with a gamma distribution: HOMA-IR as a dependent variable and host markers, age and sex as independent variables; (2) logistic regression model: IR (HOMA-IR ≥ 2.5 = 1, HOMA-IR ≤ 1.6 = 0) as a dependent variable and significant host markers in the model 1, age and sex as independent variables. In both models, host markers with Padj < 0.05 were considered to be significant. We finally identified 6, 21 and 36 significant associations from plasma cytokines, plasma metabolites and CAGE promoter expression data, respectively (Supplementary Tables 19–21). In terms of bacteria, 20 genera with significant interactions between each other, which were identified with CCREPE as shown in Extended Data Fig. 5f, were included. In terms of faecal metabolites, 15 carbohydrates associated with IR in the CAG analysis as shown in Fig. 1b were included. Pairwise partial Spearman’s rank correlations adjusted by age, sex, BMI and FBG between all given factors were calculated with the R package ppcor v.1.1. The correlations with Padj < 0.05 were selected for visualization. The size of nodes was determined as the ratio of median abundance in IR over IS. As the median values of genera Robinsoniella and Rothia were zero, these elements were removed from the visualization. The width of lines was determined as the absolute value of partial Spearman’s coefficient. The networks were visualized using Cytoscape v.3.7.0. as in the microorganism–metabolite networks described above.
Explained variance of plasma cytokines by omics data
To assess the explained variance of ten plasma cytokines, we established random-forest models using the R package caret v.6.0-88 to predict the plasma cytokine levels using 15 IR-associated faecal carbohydrates identified in Fig. 1b; 20 genera with significant interactions with each other that were identified in Fig. 2a; 21 IR-associated plasma hydrophilic metabolites (Supplementary Table 20); or 36 IR-associated CAGE promoters (Supplementary Table 21). Plasma cytokines were log10-transformed and scaled before the regression analyses. The data were split into train and test datasets at a 4:1 ratio. The generated random-forest models were evaluated using a tenfold cross-validation method and applied to the test datasets to obtain predictions. The explained variance shown as R2 was calculated as its definition: 1 − sum(test − predict)2/sum(test − mean(test))2. The negative values were considered as zero.
Causal mediation analysis
To infer the effects of plasma cytokines on in silico causal relationships between faecal carbohydrates and IR markers (HOMA-IR, BMI, triglycerides and HDL-C), we performed causal mediation analysis using the R package mediation (v.4.5.0)38. As previously reported68, we first screened significant associations (Padj < 0.05) between 15 IR-associated faecal carbohydrates and four IR markers, and significant associations between ten plasma cytokines and four IR markers. Age and sex were included as independent variables in both models. We then performed causal mediation analyses with the following models: (1) Mediator models: cytokine ~ metabolite + age + sex; (2) outcome models: IR marker ~ metabolite + age + sex + cytokine. In both models, faecal carbohydrate and plasma cytokine values were scaled before the analyses, and GLM with Gaussian distribution was used. A nonparametric bootstrap procedure was used to calculate the significance, followed by multiple testing corrections using the R function p.adjust. Average causal mediation effects and average direct effects with Padj values from representative models are reported in Extended Data Fig. 7d, whereas all of the results including the total effects and proportion mediated are reported in Supplementary Table 23.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.