Bile acids drive the newborn’s gut microbiota maturation
Ethic statement
All animal experiments were performed in compliance with the German animal protection law (TierSchG) and approved by the local animal welfare committees, the Landesamt für Natur, Umwelt und Verbraucherschutz, North Rhine Westfalia (84–02.04.2016.A207 and 84–02.04.2015.A293). All C57BL6J wildtype mice were bred locally and held under specific pathogen-free conditions at the Institute of Laboratory Animal Science at RWTH Aachen University Hospital. The day of birth was considered day 0, i.e., animals screened at day 1 were approximately 24 h old and verified to have ingested breast milk (abdominal milk spot). Mice were weaned at PND21.
In vivo study design
To monitor microbiota and host metabolic development throughout the neonatal period into adulthood, intestinal and hepatic tissues were obtained from C57BL/6 J mice in two different approaches. First, total small intestinal, colon, and liver tissues were obtained from C57BL/6 J mice aged 1, 7, 14, 21, 28, and 56 days (n = 5 per timepoint). In a separate set of experiments, similar tissues were obtained from mice aged 0, 6, 12, 18, and 24 h (n = 10 per timepoint). To rule out potential litter and cage effects, we obtained tissues from a single animal of one litter for a given age group and repeated this by selecting a new animal from the same litter for every other age group (Fig. 1a). Every animal was only examined once at the indicated age (PND). This means that in total 30 animals from 5 litters were used for monitoring the microbiota development between age 1 and 56 days and 46 animals from 10 litters to monitor the microbiota development within the first 24 h. For oral bile acid administration, 7–14 (UDCA and Control, n = 14;, GCA, n = 11; TCA, n = 10; βTMCA, n = 7) PND7 animals received the indicated bile acid (Sigma-Aldrich, Biomol) at a concentration of 70 µg/g body weight or PBS daily for 3 days by oral gavage (average 5 µL) with 100% succession47. Body weight and food/water consumption was monitored daily. To determine the effect of bile acid administration, the intestine was aseptically cut into 10–20 parts and alternately assigned to two collections for microbial profiling and assessment of metabolites. Samples from an additional group of adult 8–12 week old animals (n = 5) were processed in parallel to provide an adult microbiota control. Tissues were transferred into sterile micro-centrifuge tubes and stored at −80 °C before analysis. Liver tissues of 7-, 14-, 21- and 56-day-old germ-free mice were obtained from the Institute for Laboratory Animal Science at Hannover Medical School. Tissues were collected and stored in sterile tubes stored at −80 °C until further analysis.
DNA isolation and generation of sequence data
Total metagenomic DNA was isolated from snap frozen small intestinal and colonic tissues by repeated-bead-beating (RBB) combined with chemical lysis plus a column-based purification method48. Approximately 200 mg tissue was added to a 2.0 mL screw-cap tube containing 0.5 g of 0.1 mm zirconia beads (Biospec Products, Bartlesville, OK, USA) and 1 mL ASL lysis buffer from the QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany). Samples were incubated at 15 min. at 95 °C and subsequently two successive rounds of bead-beating were employed using a FastPrep®−24 instrument (MP Biomedicals, Inc., USA) at 5.5 ms−1 (3 × 1 min. for each round). To minimize physical damage of DNA in the RBB method, the lysate fraction produced from the first round of bead beating was drawn after centrifugation at full speed (~160,000 × g) for 5 min at 4 °C. A second round of bead beating was performed upon adding 300 μL fresh lysate buffer after which supernatants were pooled. To precipitate nucleic acids, 260 μL 10 M ammonium acetate was added to the lysate tubes, mixed and incubated on ice for 5 min. After centrifugation at 4 °C for 15 min. at full speed, supernatants were transferred to a new 1.5 mL tube to which an equal volume of isopropanol was added. Next, samples were incubated on ice for 30 min., centrifuged at RT for 15. min. at full speed and supernatants were removed by decanting. The nucleic acid pellet was washed with 500 μL 70% EtOH and dried under vacuum for 3 min. The nucleic acid pellet was dissolved in 100 μL TE (Tris-EDTA) buffer. Two microliter DNase-free RNase (1 mg/mL) was added and samples were incubated at 37 °C for 15 min. Next, 15 μL proteinase K and 200 μL AL buffer were added, samples were vortexed for 15 s. and incubated at 70 °C for 10 min. After adding 200 μL ethanol (96–100%) and vortexing, the lysate was transferred to QIAamp spin columns (Qiagen, Hilden, Germany). The DNA was finally purified using the QIAamp DNA Stool Mini Kit according to the manufacturer’s instructions and eluated in 200 μL AE Buffer. For each DNA isolation batch, additional isolation was performed on PCR-grade water as a negative control. Generation of amplicon libraries and sequencing was performed according to a previously published protocol49. Briefly, the V4 region of the 16S rRNA gene was PCR amplified from each DNA sample in duplicate in 50 μL volumes containing 10 pmol of both primers (5′-GTGCCAGCMGCCGCGGTAA*-3′ [515 F] and barcoded 5′-GGACTACHVGGGTWTCTAAT*-3′ [806 R]), 5 μL Accuprime buffer II, 0.2 μL Accuprime Hifi polymerase (Thermo Fisher Scientific, Waltham, WA, USA) and 2 μL DNA. After an initial denaturation step at 94 °C for 3 min., amplification was carried out for 30 s. at 94 °C, 45 s at 50 °C and 1 min. at 72 °C. Amplification was carried out in 35 cycles for PND1-PND56; the samples with a low yield of DNA required 40 cycles (0–24 h) The PCR program ended with a final post-PCR incubation step of 10 min at 72 °C to promote complete synthesis of all PCR products. Pooled amplicons from the duplicate reactions were purified using AMPure XP purification (Agencourt, Massachusetts, USA) according to the manufacturer’s instructions and subsequently quantified by Quant-iT PicoGreen dsDNA reagent kit (Invitrogen, New York, USA). Amplicons were mixed in equimolar concentrations, to ensure equal representation of each sample, and sequenced on an Illumina MiSeq instrument using the V3 reagent kit (2 × 250 cycles). All V4 16S rDNA bacterial sequences generated in this study have been submitted to the Qiita and ENA databases under accession No. 10719 and ERP116798, respectively.
Sequence analysis
Data demultiplexing, length and quality filtering, pairing of reads and clustering of reads into Operational Taxonomic Units (OTUs) at 97% sequence identity was performed using the online Integrated Microbial Next Generation Sequencing (IMNGS, www.imngs.org) platform using default settings50. Removal of primers and technical reads resulted in fragments of approximately 250 bases. Sequencing was performed from both the 3′ and 5′ side resulting in sufficient resolution. IMNGS is a UPARSE based analysis pipeline51. Pairing, quality filtering and OTU clustering (97% identity) was done by USEARCH 8.052. The analysis was based on OTUs rather than amplicon sequence variants (ASVs) since we aimed at aggregating taxa at a higher level and wanted to avoid overestimation of prokaryotic diversity due to Intragenomic heterogeneity of 16S rRNA genes53. Chimera filtering was performed by UCHIME (with RDP set 15 as a reference database54. Taxonomic classification was done by RDP classifier version 2.11 training set 15.8 Sequence alignment was performed by MUSCLE and treeing by Fasttree55,56. A total of 21,372,397 V4 reads were generated over two runs. After trimming, quality filtering, removal of potential chimeric reads, de-multiplexing and removal of low abundant operational taxonomic units (OTUs), 15,692,587 sequences belonging to 478 OTUs were retained for downstream analysis. Negative controls were evaluated based on their number of sequences and composition compared to other samples. We used for each batch, sampling blank controls, DNA blank extraction controls and no-template amplification controls, and monitored the lack of contaminant bacterial DNA load herein by a gel-based principle. Moreover, we compared the acquired OTUs composition from the negative controls to our low abundant microbial samples to ensure that our findings were not driven by potential contaminant taxa. Subsequently, samples of the negative controls and with low sequencing depth (less than 6,749 sequences/sample) were excluded from subsequent analysis. For the remaining samples, the number of sequences per sample ranged from 6749 to 239,395 (median 87,227).
Richness, diversity, taxonomy, and enterotype analyses
Data normalization, diversity, taxonomical binning and group comparisons were performed using the Rhea package version 1.657. In order to not discard informative data, normalization in Rhea was performed by dividing OTU counts per sample for their total count (sample depth) followed by multiplying all of the obtained relative abundance for the lowest sample depth (6749 reads/sample). Alpha- (observed species and Shannon index) and the generalized Unifrac beta-diversity index were calculated using the Rhea package58. Additional beta-diversity indices (weighted Unifrac, unweighted Unifrac, and Bray–Curtis distance) were calculated using the R package (version 3.6.1.) Phyloseq package version 1.30.059. Ordination of samples according to their microbial composition expressed as Hellinger transformed genus abundance data or beta-diversity indices was visualized using Principal Component Analysis (PCA) and Principal Coordinate Analysis (PCoA), respectively. All ordinations were constructed using the R package Phyloseq and included 95% confidence ellipses. Dirichlet multinomial mixtures models (DMM) were used to calculate genus-level enterotypes60. When including samples from all time-points, Laplace approximation revealed an optimal number of three clusters.
Downstream microbial analyses and presentation
Smoothing of the kinetic for dominant taxa (Fig. 1e and f, as well as Supplementary Fig. 1d and e) was generated using the geom_smooth function of the ggplot package 3.2.1. with default settings. The lines reflect the mean values of the relative abundance. The appearance and disappearance of OTUs of the dominant phyla (Actinobacteria, Proteobacteria, Firmicutes, Bacteroidetes) with age was visualized in Sankey-plots (SankeyMATIC.com). For readability of Sankey-plots, only OTUs present in >10% of all samples per timepoint and with a prevalence of >20% in the entire dataset, were included. Ecosystem specific functional metagenome predictions were created by the novel PICRUSt-iMGMC workflow (using PICRUSt version 1.1.3) with the de novo picked OTUs and using mouse metagenome-assembled genomes linked to 16S rRNA genes61. The derived KEGG orthologs were mapped into multiple pathways or modules. Differentially changed KEGG modules were identified using the pathways enrichment analyses from MicrobiomeAnalyst with default settings62. The identification of lactobacillus-related OTUs were performed using EZbiocloud and used to construct a phylogentic tree of lactobacilli by MEGA7 (MUSCLE) for alignment and iTOL v4 for the final annotations (itol.embl.de) with default settings.
Liver metabolomics and bile acid analyses
The metabolome analyses were carried out with the AbsoluteIDQ® p180 Kit (Biocrates Life Science AG, Austria. The kit allowed identification and quantification of 188 metabolites from 5 compound classes (acyl carnitines, amino acids, glycerophospho-, and sphingolipids, biogenic amines, and hexoses). The kit used flow injection tandem mass spectrometry (FIA) for the non-polar metabolites and LC-MS/MS for the more polar compounds. The integrated MetIDQ Software (version Boron 2623) streamlined the data analysis by automated calculation of metabolite concentrations. Quantification of analytes utilized stable isotope-labeled or chemically homologous internal standards (IS). Controls were included for 3 different concentration levels. For calibration, the kit contained a calibrator mix at 7 different concentrations. The measurements were carried out with an ABI Sciex API5500 Q-TRAP mass spectrometer via Electrospray ionization (ESI) by Multi Reaction Monitoring (MRM) mode for high specifity and sensitivity. 158 MRM pairs were measured in positive ion mode (13 IS) and 2 MRM pairs were measured in negative mode (1 IS). The following additional chemicals for LC-MS were used: water, Millipore; PITC, Fluka; pyridine, Fluka (p.a.); methanol, Merck; Lichrosolv for LC/MS; acetonitrile, Merck; formic acid, Sigma Aldrich. Metabolites were extracted from liver samples by adding H2O/acetonitrile (1:1,v:v) per mg sample followed by homogenization with a tissue disruptor (10 min, 30 Hz, 4 steel balls). The samples were centrifuged (1400×g, 2 min) and the supernatant analyzed. The targeted analysis was performed by adding 10 μL of extracted liver sample to the AbsoluteIDQ® p180 Kit (Biocrates Life Science AG, Innsbruck, Austria), following the vendor’s instructions63. For bile acid measurements the MS-based Bile Acids Kit (Biocrates Life Sciences AG, Innsbruck, Austria), a 96-well plate format assay, was used following the manufacturer´s instructions with normalization based on mass64. The following settings were used: turbo spray for ion source, 20 for Curtian Gas, medium for CAD Gas, 40 psi for ion source gas 1 and 50 psi for ion source gas 2. For the bile acid measurements we used an ion spray voltage of −4500 V and a temperature of 600 °C; for the p180 kit we used an ion spray voltage of 5500 V and a temperature of 500 °C. All metabolomic data generated in this study have been submitted to the Metabolomics Workbench and have been assigned the Study ID ST001388, ST001389, ST001396, ST001397, the Project ID PR000952 and the Project DOI 10.21228/M8N397.
The contribution of each metabolite to metabolomic variation was derived from age-constrained redundancy analysis (RDA) based on all metabolites stratified in group levels with additional scaling by normalization based on z-scores (Phyloseq package). Moreover, PCA was used to illustrate changes in the composition of the different metabolic groups during the postnatal period for both PC1 and PC2, and PC2 and PC3. For the bile acids the 2nd and 3rd component were chosen for plotting, since the first component was mainly driven by the strong separation observed between the first and subsequent time points (PND1 vs. PND7–56).
Multi-omics analyses
Regularized canonical correlation analyses (rCCA) were performed (Mixomics package 6.10.8)65 to unravel specific correlations between bile acids and OTUs with a minimal presence of 20% in all samples. Samples were excluded from the microbiota-data if they were not measured in the bile acid analyses (i.e. PND1 of litter 1 and 2). Prior to rCCA a hyperbolic sine transformation was used on OTU-counts and a log-transformation for the bile acids. For the estimation of regularization (penalization) parameters λ1 and λ2, the cross-validation procedure (CV) method was used. We used a λ1 = 0.0001, λ2 = 1 with a CV-score = 0.4779644 and 2 components. OTUs with a correlation between −0.3 and 0.3 on the first 2 components were filtered out to optimize the rCCA.
The Spearman’s rank correlation coefficient was calculated between bile acids (weight corrected) and bacterial genera (relative abundances) with a minimal presence of 20% in all samples. Benjamini and Hochberg FDR correction was performed to correct for multiple testing (p > 0.05). For the heatmap only significant correlations with adjusted p-value of More
