The Mouse Gastrointestinal Bacteria Catalogue enables translation between the mouse and human gut microbiotas via functional mapping

Lead contact

Further information and requests for resources, reagents and software should be directed to and will be fulfilled by the lead contact, Virginia A. Pedicord (vap33@cam.ac.uk).

Materials availability

Isolates of additional culturable species generated in this study are being deposited at the Leibniz Institute DSMZ-German Collection of Microorganism and Cell Cultures (DSMZ), and accession numbers are available under https://github.com/BenBeresfordJones/MGBC. Isolates generated in this study that have not been deposited at DSMZ will be made available without restriction upon request.

Faecal samples from mice

Mice were maintained under specific pathogen-free conditions at a Home Office-approved facility in accordance with the United Kingdom Animals (Scientific Procedures) Act of 1986. Faecal samples from 30 mice aged between 5 to 8 weeks were used to generate the Mouse Culture Collection. We used both male and female representatives from 10 different mouse colonies at the Wellcome Sanger Institute. Colony genotypes are detailed in Table S8. For shotgun metagenome sequencing, faecal samples were obtained from 48 8-week-old male C57BL/6N mice.

Bacterial culture

Fresh faecal samples were collected from mice into sterile 1.5mL Eppendorf tubes using aseptic technique. Sample processing and culturing were performed under anaerobic conditions (80% nitrogen, 10% carbon dioxide, 10% hydrogen) in a Don Whitley A95 anaerobic workstation. Faeces were homogenised in sterile, pre-reduced PBS (100mg/mL) and a 10-fold 1:10 dilution series performed. 200μL of each dilution was plated onto pre-reduced 140mm agar plates and incubated at 37°C. A range of agars were employed to maximise culturing yields and diversity (Table S8). After 2 days, individual colonies were picked and re-streaked onto fresh plates. This was then repeated until purity was achieved. Colonies were identified using 16S rRNA gene sequencing. Single bacterial colonies were scraped into 2mL screw cap tubes containing glass beads (acid-washed 425–600μm) and 500μL sterile PBS, and then physically lysed by 30 seconds shaking at speed 6.0 using a FastPrep Instrument (MPBio). After centrifugation at 14,000rpm for 5 minutes, 1μL of supernatant was taken to carry out a 16S PCR using the standard 7F and 1510R bacterial primers (Browne et al., 2016) and GoTaq Hot Start reagents (Promega). PCR products were sequenced by an external supplier (Eurofins Genomics) and mothur (Schloss et al., 2009) used to align the resulting sequences and create OTUs representing clusters of ≥97.8% sequence identity. For each OTU, a single sequence was taxonomically classified using NCBI BLAST (Johnson et al., 2008) and a single isolate selected for further culturing and whole-genome sequencing. 10mL of BHI or YCFA broth was inoculated for each new isolate identified and left to grow overnight. 500μL of the overnight culture was mixed with 500μL of 50% glycerol in a cryotube (performed in quadruplicate) and these were frozen at −80°C. The remaining overnight culture was centrifuged at 4,000rpm, and the cell pellet then washed with 10mL sterile PBS. For isolates that did not grow in broth, 2mL sterile PBS was pipetted onto mono-inoculated agar plates and colonies were dissolved using a bacterial scraper. Plate supernatants were then used in place of overnight cultures. Genomic DNA was extracted from the washed pellet using the ‘MasterPure Complete DNA and RNA Purification Kit’ (Lucigen). Genomic DNA was kept at 4°C until sequencing. While 276 strains were cultured to purity and sequenced, only 223 were subsequently recovered successfully from cryopreservation for banking at the Wellcome Sanger Institute and DSMZ.

Curation of public mouse gut-derived isolate genomes

To build a comprehensive genome collection, we curated 319 publicly available mouse-derived isolate genomes from the ENA, including the genomes of previously published mouse gut isolate collections (Lagkouvardos et al., 2019, 2016; Liu et al., 2020). As genome assemblies for the isolates of the miBC (Lagkouvardos et al., 2016) had not been made available, we assembled these from raw reads according to our standard genome processing pipeline. Only genomes that passed our quality control thresholds (n=288; “Quality control of genome assemblies”) were included in the analyses for this study. The metadata for these public mouse isolate genomes are included in Table S1.

Curation of public human gut-derived bacterial genomes

To perform taxonomic and functional analyses between mouse and human gut bacterial species, we curated 204,939 non-redundant human gut microbial genomes from the Unified Human Gastrointestinal Genome (UHGG) catalogue (Almeida et al., 2021). We applied the same quality control and taxonomy assignment pipelines to the UHGG genomes as with the MGBC. In total, 100,456 non-redundant, high-quality human gut microbial genomes, representing 3,006 species, were curated for comparison with the MGBC. Metadata for these genomes are provided in Table S3.

Mouse gut metagenome cohort curation

To create a comprehensive mouse gut-derived metagenome catalogue for MAG synthesis and quantification of species global abundance and prevalence profiles, we utilised the Advanced Search functionality of the ENA (Harrison et al., 2021) to identify all WGS raw read samples with an NCBI Taxonomy metadata value of “mouse gut metagenome” (taxid: 410661; last accessed February 2021). In total, 8,418 samples were identified which were then manually assessed according to our exclusion criteria. We additionally performed a review of the available literature to identify further studies that our ENA search might have overlooked due to constraints with manual metadata entry. Samples were technically excluded if they were unpublished, i.e., where no publication listed the study accession number, or if they represented 16S rRNA amplicon sequencing datasets. Samples were biologically excluded if mice had been exposed to antibiotics, had an active gastrointestinal infection, or had received faecal microbiota transplantation derived from a non-mouse host. Furthermore, we additionally excluded samples that would disrupt the validity of species prevalence and abundance analyses, e.g., following enrichment for viral particles, or deriving from an ex-germfree mouse reconstituted with a simplified microbiota. Details of included and excluded studies are provided in Table S2. Of the included metagenomes, 64 were kept aside as independent samples for read classification analyses. Metadata for these samples are included in Table S2.

These publicly available samples were supplemented with 48 newly sequenced faecal samples from 8-week-old male C57BL/6N mice. In total, 2,913 metagenome samples were used to generate MAGs; data for these samples are supplied in Table S2. 2,446 samples yielded ≥1 MAG and were used to generate abundance and prevalence profiles of mouse gut microbial species. The metadata for these samples are provided in Table S2.

Whole-genome sequencing and assembly

Bacterial genomic DNA was sequenced using the Illumina Hi-Seq Ten platform at the Wellcome Sanger Institute with library fragment sizes of 200–300 bp, a read length of 150 bp and a target read depth of 100x. Annotated assemblies were produced using a previously described pipeline (Page et al., 2016). Briefly, multiple assemblies were generated from sequence reads using Velvet v1.2 (Zerbino and Birney, 2008) and VelvetOptimiser v2.2.5. An assembly improvement step was applied to the assembly with the best N50 (Page et al., 2016), and contigs were scaffolded using SSPACE v2.1.1 (Boetzer et al., 2011) and sequence gaps filled using GapFiller (Boetzer and Pirovano, 2012). Automated annotation was performed using Prokka v1.14.5 (Seemann, 2014).

Faecal sample collection and shotgun metagenomic sequencing

Faecal samples were collected directly from C57BL/6N mice using aseptic techniques, and immediately stored at −80°C until DNA extraction. DNA was extracted using the ‘FastDNA Spin Kit for Soil’ (MPBio) according to manufacturer’s instructions and stored at −20°C until metagenomic sequencing. DNA samples were quantified using a Qubit 4 Fluorometer (Thermo Fisher), and samples with ≥100 ng DNA material proceeded to paired-end (2×150 bp) shotgun metagenomic sequencing on the HiSeq 4000 platform.

Workflow for metagenome-assembled genome (MAG) synthesis

MAGs were generated using a custom in-house pipeline that leveraged MetaWRAP v1.2.3 (Uritskiy et al., 2018) for single sample assembly, binning and bin refinement. First, metagenomes were quality controlled using KneadData v0.7.3 with default settings. Host reads were removed from samples using the GRCm39 reference genome and Bowtie2 v2.3.5 (Langmead et al., 2019). In addition, reads were aligned to the phi-X174 genome and removed. MetaSPAdes v3.10.1 (Nurk et al., 2017) was used for the assembly of paired-end samples. In line with other reports, runtime with MetaSPAdes was excessively long (>48 hours) for some samples (Kim et al., 2020). In these cases, and for samples with only unpaired reads, MEGAHIT v1.1.1-2-g02102e1 (Li et al., 2016) was used for assembly. To generate genome bins, we utilised the ‘binning’ module from MetaWRAP to run MetaBAT2 v2.9.1 (Kang et al., 2019), MaxBin 2.0 v2.2.4 (Wu et al., 2016) and CONCOCT v0.4.0 (Alneberg et al., 2014) in parallel on each sample. These bins were then consolidated and refined using the ‘bin_refinement’ module.

Quality control of genome assemblies

Both isolate genomes and MAGs were subjected to stringent genome quality criteria to be included in our analyses. Completeness, contamination, and genome fragmentation were estimated using CheckM v1.1.2 (Parks et al., 2015). Genome assemblies with ≥90% completeness, <5% contamination, maximum contig count ≤500 (Browne et al., 2016), maximum genome size ≤8 Mb (Zou et al., 2019), N50 ≥10,000 kb (Nayfach et al., 2019) and mean contig length ≥5 kb (Zou et al., 2019), were defined as high-quality genomes, in line with guidelines and previous studies (Bowers et al., 2017; Parks et al., 2017; Pasolli et al., 2019). Any isolate genome that did not meet these thresholds was excluded from analysis. For MAGs, we additionally defined medium-plus quality genomes as those assemblies with ≥50% completeness, <5% contamination and a quality score ≥50, where

(Parks et al., 2017). This definition exceeds the medium-quality thresholds as defined by MIMAGs (Bowers et al., 2017). Only high-quality genomes were included in our analyses, unless otherwise indicated. Data for MAG yields from included samples are provided in Table S2.

Additional quality metrics were calculated to complement the contamination scores generated with CheckM. For the non-redundant, high-quality genomes of the MGBC genome chimerism was quantified using GUNC v1.0.4 (Orakov et al., 2021) and strain heterogeneity was assessed using CMseq v1.0.3 (Pasolli et al., 2019). Quality data for these genomes are provided in Table S3.

Taxonomic classification of genomes and species clustering

To remove redundancy from our collection of high-quality genomes, we used dRep v2.5.4 (Olm et al., 2017) to remove conspecific genomes that shared ≥99.9% ANI (options: -pa 0.999 --SkipSecondary). We taxonomically classified our genomes using the ‘classify_wf’ workflow from GTDB-Tk v1.3-r95 (Chaumeil et al., 2019). To generate species clusters for genomes that could not be assigned to a species-level taxonomy using GTDB-Tk, we clustered these genomes at ≥95% ANI using a two-step genomic distance analysis implemented by dRep (options: -comp 50 -con 5 -pa 0.9 -sa 0.95 -nc 0.6). Previously calculated quality data from CheckM were supplied for each genome with the --genomeInfo flag to reduce computation time. Genomes that shared ≥95% ANI at 0.6 alignment fraction were considered the same species (Jain et al., 2018; Nielsen et al., 2014; Varghese et al., 2015).

To determine genome representatives for each species cluster, we ranked each genome according to a modified quality score

and used the highest scoring genome from each species as the representative.

Determining prior cultured status of isolates

We inferred the cultured status of our isolates using a two-step approach. First, we compared our isolate genomes to the publicly available mouse-derived isolate genomes of the MGBC. Isolates were considered to be the same species if (1) they were designated as the same species by GTDB-Tk, or (2) they shared ≥95% ANI across an alignment fraction of ≥0.6, in the case of a non-species-level classification. Next, we searched our isolates that were not represented in the public mouse gut isolate collections against NCBI RefSeq release 205 (O’Leary et al., 2016) using Mash v2.2.2 (Ondov et al., 2016), after which the most similar RefSeq genome to each isolate was then compared using FastANI v1.3 (Jain et al., 2018). As RefSeq excludes metagenome-derived genomes, an isolate was designated as “previously uncultured” if it shared <95% ANI with the closest related genome.

Benchmarking binner performance for production of high-quality MAGs

We compared high-quality bins generated from a subset of 2,303 publicly available mouse gut metagenomes by MetaBAT2, MaxBin 2.0 and CONCOCT, as well as the consolidated bins generated by MetaWRAP refinement. Binner performance was compared across four metrics:

• i.

  CheckM estimates of genome quality (completeness, contamination, quality score) across all high-quality bins

• ii.

  Taxonomic coverage of high-quality bins (i.e., the number of species represented)

• iii.

  Average quality score of high-quality bins on a per species basis

• iv.

  Core genome conservation by high-quality bins

Akkermansia muciniphila (Am), Bifidobacterium globosum (Bg), Ligilactobacillus murinus (Lm), and Lactobacillus johnsonii (Lj) were selected for comparing the conservation of the core genome between binners as these species ranked among the top 10 most commonly binned species across all four binning tools, and each have ≥50 isolate genomes with which to build a robust baseline core genome. In addition, these species represent three phyla, reducing the potential for any taxonomic bias. Isolate genomes were compiled from RefSeq and the MGBC, and high-quality isolate genomes that were designated as the correct species by GTDB-Tk were annotated using Prokka v1.14.5 with default settings. For each species, Panaroo v1.2.4 (Tonkin-Hill et al., 2020) was used to build an ‘isolate-only’ core genome and an ‘isolate+bins’ core genome for each binner using the following options: −clean-mode strict --core_threshold 0.99. The number of isolate genomes used for each species was as follows: Am, 136; Bg, 62; Lm, 58; Lj, 54. For each species, a standardised number of subsampled bins was used to build the ‘isolate+bins’ core genome for each binner: Am, 90; Bg, 35; Lm, 150; Lj, 60. To quantify core genome conservation, 100 iterations of bin subsampling and core genome analysis were performed for each binner, and the core genome size distribution was calculated as a percentage of the ‘isolate-only’ core genome.

Comparison to other mouse microbiota resources

The MAGs of the iMGMC and co-abundance gene groups (CAGs) from Xiao et al. (2015) (MGCv1) were accessed and processed according to the quality control and taxonomic assignment protocols above. 8,509 MAGs of the iMGMC were defined as high-quality genomes and used for benchmarking the MGBC. Comparison of this resource to the MGBC consisted of three stages: taxonomic coverage, genome quality, and metagenome coverage. Shared and unique species were identified between the two collections using GTDB-Tk and dRep as performed for the MGBC alone. For each shared species, the quality score was compared between the representative genome of each collection.

To assess metagenomic read classification performance, custom Kraken2 (Wood et al., 2019) databases were built for all high-quality genomes (MGBC, iMGMC) and medium-plus quality genomes (mq MGBC, mq iMGMC) of the MGBC and iMGMC. As CAGs cannot be defined as high-quality due to failure to meet minimum genome fragmentation criteria (Bowers et al., 2017), a custom Kraken2 database for the 239 CAGs defined as medium-plus quality was built for comparison (MGCv1). In addition, custom Kraken2 databases were built for the post-qc isolate genomes of the miBC (Lagkouvardos et al., 2016) (n=43), the mGMB (Liu et al., 2020) (n=100), all publicly available mouse-derived isolates (Public, n=288), the Mouse Culture Collection (MCC, n=276), all mouse isolates (MCC+Public, n=564), and representative genomes of the high-quality species of the Unified Human Gastrointestinal Genome (Human, n=3,006). The standard Kraken2 database for all NCBI genomes (accessed 2^nd December 2020) was also included. 64 independent, post-qc mouse gut metagenome samples that had not been included in the generation of the MGBC were analysed with the different Kraken2 databases and percentage read classification was utilised as a proxy for database efficiency. The metadata for these samples are included in Table S2.

Construction of phylogenetic trees

Maximum likelihood phylogenetic trees were built de novo from protein sequence alignments of 120 core bacterial genes generated by the GTDB-Tk ‘align’ module. The phylogenetic trees of the Mouse Culture Collection and the MGBC representative genomes were built using FastTree v2.1.10 (Price et al., 2010) with default settings (BLOSUM45 matrix; JTT+CAT model). IQ-TREE v1.6.10 (Nguyen et al., 2015) was used with default settings to build a phylogenetic tree of representative genomes for 3,006 human-derived and 1,094 mouse-derived bacterial species. LG+F+R10 was identified as the best fit protein substitution model based on the Bayesian information criterion (Kalyaanamoorthy et al., 2017). Trees were visualised using Interactive Tree Of Life v5.6.3 (Letunic and Bork, 2021).

Metagenome classification and analysis

For analysis of mouse gut shotgun metagenome samples, taxonomic classification was performed using Kraken2 v2.0.8 (Wood et al., 2019) and Bracken v2.5.2 (Lu et al., 2017). To enable species-resolved metagenomic analyses, we built a custom Kraken2/Bracken database (options: -k 31 -l 150) with the 26,640 high-quality genomes of the MGBC using a custom GTDB taxonomy (Parks et al., 2020, 2018). Only post-qc metagenomes that were of sufficient read depth to generate MAGs were used in metagenomic analyses (n=2,446). The resulting Bracken outputs were compiled and analysed with R v4.0.2 (R Core Team, 2020). To calculate prevalence, a threshold of 0.01% assigned classified reads was used to define presence of a species in a sample.

Due to the compositional nature metagenome analyses, we determined the Aitchison distances between samples (Aitchison, 1992). We performed Bayesian-multiplicative treatment of count zeros (Martín-Fernández et al., 2015) using the zCompositions v1.3.4 R package (Palarea-Albaladejo and Martín-Fernández, 2015) and transformed data using a center log-ratio transformation. Finally, the Euclidean distances of samples were determined using the vegan R package v2.5-6 (Oksanen et al., 2014). To assess the ability of metadata variables to explain variance in microbial communities of the mouse microbiome, we applied the ‘adonis’ function from vegan to calculate the Permutational Multivariate Analysis of Variance of the Aitchison distance matrices using 999 permutations. The PERMANOVA summary statistics are provided in Table S5.

For institute analyses, samples from “control” C57BL/6 mice were curated that were 1) faecal samples, 2) from wildtype mice, 3) not exposed to a wild mouse microbiota, and 4) fed only a regular chow diet. Only institutes that were represented by ≥10 samples were compared. Center log-ratio transformation of the data was performed and a heatmap generated using the pheatmap R package v1.0.12 for the top 20 most abundant species of these samples against institute. For laboratory vs wild analyses, 65 samples from ‘wild’ gut microbiotas were compared against 1,065 samples from control ‘laboratory’ microbiotas. Hybrid microbiotas, where wild and laboratory mice were crossed, were excluded from these analyses.

Taxonomic comparison of the mouse and human gut microbiotas

For taxonomic comparison, species were considered shared between the UHGG and the MGBC if they were annotated as the same species by GTDB-Tk, or, if they could not be assigned at a species-level, the representative genomes shared ≥95% ANI across a minimum alignment fraction of 0.6.

Pangenome synthesis and functional annotation

To generate species pangenomes for functional annotation, we first concatenated 76,937,350 pre-clustered (100% sequence identity) proteins derived from the 100,456 non-redundant, high-quality genomes of the UHGG with 67,768,723 predicted proteins from the 26,640 non-redundant, high-quality genomes of the MGBC, and performed protein clustering using the ‘linclust’ function from MMseqs2 v10-6d92c (-c 0.8 --cov-mode 1 --cluster-mode 2 --kmer-per-seq 80) (Steinegger and Söding, 2018, 2017). Proteins were clustered at 100%, 90%, 80% and 50% sequence identity.

Next, we generated species pangenomes by concatenating all non-redundant (90% sequence identity) protein-coding sequences from member genomes. Pangenomes were then functionally annotated using both InterProScan v5.39-77.0 (Jones et al., 2014) and EggNOG-mapper v2.0.1 (Huerta-Cepas et al., 2017).

To assess global functional overlap between host gut microbiotas, we calculated the number of InterPro (IPR) protein families and KEGG orthology (KO) groups that were shared in total between human and mouse pangenomes.

Taxonomic and functional distance analyses

Taxonomic distances were calculated from the human-mouse representative genome phylogeny (“Construction of phylogenetic trees”). Phylogenetic branch lengths between the representative genomes of each human and mouse species were calculated using the ‘cophenetic.phylo’ function from ape v5.5 (Paradis and Schliep, 2019).

To maximise the resolution of functional comparisons, functional annotations from multiple schemes – including InterPro, KEGG (Kanehisa et al., 2017), MetaCyc (Caspi et al., 2016), CAZy (Lombard et al., 2014), Reactome (Jassal et al., 2020), and Gene Ontology (Ashburner et al., 2000) – were considered together (all functions) in the context of each species pangenome. To facilitate functional distance analyses, pangenome-feature matrices were constructed where each functional annotation was scored according to the fraction of genomes per pangenome that were annotated with that feature. Inter-pangenomic functional distances were then calculated using the Jaccard Index.

A global comparison of the interspecies taxonomic and functional relationships was performed using a Mantel test. The taxonomic and functional distance matrices were ordinated using the ‘cmdscale’ function, and the two most dominant principal coordinates were visualised using the R packages ggpubr v0.4.0 and ggplot2 v3.3.5.

Functional comparison of Phocaeicola species

To identify the functional pathways that underpin the divergence between human and mouse Phocaeicola dorei (PdH, PdM) and the relative similarity of human P. dorei and mouse P. vulgatus (PvM), InterProScan v5.39-77.0 and Genome Properties v2.0.1 (Richardson et al., 2019) were run on all genomes for these species (PdH, n=1,954; PdM, n=15; PvM, n=177). The tabular outputs were concatenated and the proportion of ‘complete pathway’-encoding genomes for each species was calculated for each property. The proportion of property-encoding genomes was then compared between each species pair using two-proportion z-tests with Yates’ continuity correction. The data from these analyses are provided in Table S10.

Drug metabolism analyses

The protein sequences of genes validated for drug metabolism by the human gut microbiota from four independent studies (Table S7) were accessed via UniProt (The UniProt Consortium, 2021). Each gene was queried against a UHGG-MGBC combined protein catalogue (pre-clustered at 100% sequence identity) using BLAST+ v2.7.1 (Camacho et al., 2009). The ‘hm_blast’ module of the MGBC-Toolkit (https://github.com/BenBeresfordJones/MGBC-Toolkit) was utilised for these analyses. Hits with ≥95% sequence identity were considered functionally equivalent, and the genomes of origin for each hit were identified. For bt_4096, additional hits with sequence identity down to 50% were considered.

Butyrate synthesis analyses

The ‘feature_search’ module of the MGBC-Toolkit was applied to search bacterial pangenomes for the IPR family identifiers of the terminal pathways of butyrate synthesis: IPR023990 (BCOAT), IPR011245 (BUK), IPR014079 (PTB). Only species encoding both BUK and PTB were considered as butyrate producers using the PTB/BUK terminal pathway. The fraction of genomes per pangenome that were predicted to encode terminal pathway genes was calculated and based on these data a threshold of 70% was used to define butyrate-producing species.

Validating butyrate terminal pathway predictions using isolate cultures

From the functional annotations of their genomes, isolates from the Mouse Culture Collection were identified that encoded either the BCOAT pathway, the PTB/BUK pathway or neither. For each pathway, three isolates that ranked in the top 30 most abundant predicted butyrate-producing species according to the global mouse gut metagenome catalogue curated above (“Metagenome classification and analysis”) were selected for culturing. Isolates for known butyrate-producing species of the human gut microbiota (Agathobacter rectalis, Coprococcus eutactus_A) were included as positive controls, and an isolate from the MCC that lacked predicted butyrate terminal pathway genes (Lachnospiraceae_NOV MGBC000113) was included as a negative control.

Under anaerobic conditions and using pre-reduced reagents, isolates were streaked onto YCFA (Duncan et al., 2002) agar and single colonies picked into 10 mL YCFA broth. Broth cultures were incubated at 37°C for 48 hours. Culture turbidity and pH were measured to confirm growth, and 16S rRNA gene sequencing was performed to check for contamination. Bacterial cultures were centrifuged at 3,600rpm for 5 minutes to pellet cells. The supernatant was pipetted off and sterile-filtered before being immediately stored at -80°C until analysis. These experiments were performed in triplicate.

Short chain fatty acid quantitation by GC-MS

Bacterial media samples (100 μL) were extracted using 400 μL methanol containing acetate-d3, propionate-d5, butyrate-d7 and valerate-d9 as internal standards (Cambridge Isotope Laboratories). After vortexing, samples were incubated at -80°C for >1 h to promote protein precipitation, then centrifuged for 20 min at 20,000g at 4°C. 100 μL of the resulting supernatant was added to 100 μL of 100 mM borate buffer (pH 10). Subsequently, 400 μL of 100 mM pentafluorobenzyl bromide (Thermo Scientific) diluted in acetonitrile (Fisher) and 400 μl of cyclohexane (Acros Organics) were added and reaction vials were sealed. Samples were derivatised by heating to 65°C for 1 h with agitation, then cooled to room temperature and centrifuged at 2,000g for 2 min to promote phase separation. 100 μL of the cyclohexane (upper) phase was transferred to a fresh autosampler vial and diluted 1:100 with cyclohexane prior to analysis. Gas chromatography-mass spectrometry (GC-MS) was performed using an Agilent 7890A gas chromatograph and Agilent 5975C MS detector operating in negative chemical ionisation mode. A 1μL splitless injection was made onto a VF-1701ms column (30 m × 0.25 mm, 0.25 μm; Agilent Technologies). Helium (1.2 mL/min) was the carrier gas and methane (2 mL/min) was used as the chemical ionisation reagent gas. For SCFA quantification, the peak areas of acetate (m/z 59) and propionate (m/z 73) were normalized to acetate-d3 (m/z 62) and propionate-d5 (m/z 78) internal standards respectively; the C4 compounds butyrate and isobutyrate (m/z 87) were normalized to butyrate-d7 (m/z 94) and the C5 compounds 2-methylbutyrate, valerate and isovalerate (m/z 101) were normalized to valerate-d9 (m/z 110). Calibrator and quality control (QC) samples were prepared in borate buffer and derivatised using the same procedure covering the range 0.05–125 mM. All data analyses were performed with Agilent MassHunter quantitative analysis software (version 10.1, Agilent Technologies) and QC samples were confirmed to be within ±15% accuracy. Butyrate production data are provided in Table S9.