Study design

Healthy young adults between 18 and 35 years old were recruited. Exclusion criteria included antibiotic treatment within the previous 6 months, current prescribed pharmaceutical drug utilization, or active acute or chronic diseases. All participants were verbally interviewed on their dietary habits and CRF was determined using a VO 2 peak cycle test. Participants were then provided a stool collection kit with instructions and were asked to provide a sample within a week following their lab visit.

Nutritional data collection

On the day of VO2peak testing, nutritional data, including supplements, was collected by means of a 24-h dietary recall interview and assessed by a research nutritionist using FoodWorks nutrient analysis software (version 16.0). Food items described by participants that were not available in the software were manually added as needed. A sample copy of a completed questionnaire is available in Additional file 1. On average, over 100 food categories per participant were produced by the FoodWorks software. A manual screening was applied to select categories of interest based on a priori interest and existing literature showing a significant interaction between those categories and intestinal microbiota. The selected 24 food category data are available in the uploaded metadata mapping file.

Cardiorespiratory fitness testing

Participants initially completed a physical activity readiness questionnaire (PAR-Q) to rule out any contraindications to vigorous exercise. A continuous incremental ramp maximal exercise test on an electronically braked cycle ergometer (Lode Excalibur, the Netherlands) was used to determine VO 2 peak and peak power output (Wpeak). Expired gas was collected continuously by a metabolic cart (Parvomedics TrueOne 2400, Salt Lake City, Utah, USA) calibrated with gases of known concentration. The test started at 50 W and increased by 30 W/min. The test was terminated upon volitional exhaustion or when revolutions per minute fell below 50. VO 2 peak was defined as the highest 30-s average for VO 2 (in ml/kg/min). Criteria for achieving VO 2 peak were the following: (i) respiratory exchange ratio >1.15; (ii) plateau in VO 2 ; (iii) reaching age-predicted HRpeak (220-age); and/or (iv) volitional exhaustion. Following VO 2 peak assessment, participants were categorized to either low (LOW), average (AVG), or high (HI) fitness based on their sex and age according to a modified Heyward normal VO 2 max reference chart (Additional file 2).

Stool collection and storage

Participants were provided with a home stool collection kit including a sterile 120 ml polypropylene container (Starplex, Etobicoke, Ontario), sterile tongue depressor and gloves, and an ice box. Participants were instructed to avoid alcohol for 3 days prior to stool collection. Stool samples were immediately stored in the participant’s freezer overnight and transported on ice to the lab and stored in −80 °C until further analysis. Frozen portions from the inner area of the samples were scrapped using sterile razor blades for DNA extraction and SCFA analysis.

SCFA analysis

SCFAs (acetic, propionic, heptanoic, valeric, caproic, and butyric acid) were analyzed from the feces by gas chromatography (GC) as described previously [14]. In brief, ~50 mg of stool was homogenized with isopropyl alcohol, containing 2-ethylbutyric acid at 0.01 % v/v as internal standard, at 30 Hz for 13 min using metal beads. Homogenates were centrifuged twice, and the cleared supernatant was injected to Trace 1300 Gas Chromatograph, equipped with Flame-ionization detector, with AI1310 auto sampler (Thermo Fisher Scientific) in splitless mode. Data was processed using Chromeleon 7 software. An aliquot of 50 mg of stool was freeze dried to measure the dry weight, and measurements are expressed as mass % (g of SCFA per g of dry weight stool).

High-throughput sequencing

DNA was extracted from feces using QIAmp DNA Stool Mini Kit (Qiagen) according to the manufacturer’s instructions following 3 × 30 s of homogenization using metal beads on a Retsch MixerMill MM 400 homogenizer. Sequencing libraries were prepared according to the Illumina MiSeq system instructions. In brief, the V3 and V4 region of the 16S bacterial rRNA gene was amplified using recommended primers [15] (IDT, Vancouver, Canada): Forward 5′ TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG, and Reverse 5′ GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC, which create amplicons of ~460 bp. Amplicons were cleaned using AMPure XP bead step, and then, adapters and dual-index barcodes (Nextera XT) were attached to the amplicons to facilitate multiplex sequencing. After another clean-up step, libraries were validated on an agarose gel, quantified, normalized, and sent to The Applied Genomic Core (TAGC) facility at the University of Alberta (Edmonton, Canada) for sequencing using the Illumina MiSeq platform. The resulting ~16,000,000 paired-end reads were merged using PEAR software [16] and screened to exclude sequences containing one or more base calls with a Phred score <20. The average read per sample was ~350,000 with a min/max of ~165,000/452,000 reads. Rarefaction curves demonstrated that sufficient sampling depth had been reached amongst all samples (Additional file 3).

Bioinformatics

Bioinformatics analyses on the demultiplexed paired reads were conducted using QIIME 1.8.0 software suites [17]. Reads were clustered at 97 % identity using the uclust method into operational taxonomic units (OTUs) then aligned to the most recent available version (2013/08) of Greengenes bacterial database [18]. Singleton and doubletons were removed, and the produced OTU table was normalized using phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) [19] to adjust for different 16S rRNA gene copy numbers. Instead of rarefying our OTU table to the lowest sample depth [20], uneven variance as a result of differential sample sequencing depth was stabilized using the cumulative sum scaling (CSS) method [21] of “metagenomeSeq” package in R. Alpha diversity indexes, rarefaction curves, OTU tables, and distance metrics were also generated using QIIME.

Statistical analysis

All statistical analyses were performed using R [22] version 3.2.0 unless stated otherwise.

The groups’ age and VO 2 peak data were tested for normality using Shapiro-Wilk test, and a one-way analysis of variance (ANOVA) with Tukey’s multiple-comparison test was used to compare mean differences amongst groups. Kruskal-Wallis non-parametric test was used for comparing BMI as this dataset failed normality tests even after several transformation attempts. For comparison of dietary intake amongst groups, a permutational multivariate ANOVA (PERMANOVA) with 999 random permutations was used. Due to the inherent high variability of dietary data, we further searched for dietary patterns amongst groups by looking at a principal component analysis (PCA) plot of participants’ dietary scores using the ggbiplot package [23]. To facilitate comparisons with previous work, we first compared average alpha diversity amongst the three fitness categories using a one-way ANOVA, followed by a Tukey’s multiple comparison. To simultaneously evaluate the role of CRF alongside other potential predictors of alpha diversity (sex, age, BMI, and dietary components), we performed a multiple regression analysis. Given our comparatively low sample size (n = 39), and the general rule that multiple regressions should include at least 10 observations per predictor variable [24], we first screened potential predictors using a Spearman correlation matrix. Those that showed a significant correlation with alpha diversity were retained for entry in the multiple regression model. Multicolinearity was checked using the variable inflation factor (VIF) index with a maximum cutoff score of 10.

Microbial communities in fecal samples were ordinated using the Bray-Curtis and weighted and unweighted UniFrac distance metrics. Principal coordinate analysis (PCoA) based on the Bray-Curtis dissimilarity metric was conducted using the cmdscale function in the base “stats” package in R, while PCoA based on the weighted and unweighted unifrac distances was made using EMPeror tool [25]. Microbial communities were analyzed using two complementary multivariate approaches: (1) constrained ordination and (2) generalized linear models (GLM). For the constrained ordination approach, redundancy analysis (RDA) was used, which focuses on assemblage composition differences in relation to predictors of interest (VO 2 peak, sex, age, BMI, and dietary components). This was implemented using the “vegan” package [26] version 2.2-1 in R. Abundance data at each taxonomical resolution (phyla, class, order, family, and genus) were first Hellinger-transformed [27] to accommodate counts data with large occurrences of low and zero abundance. Variable selection in RDA was implemented using the ordistep function of vegan using both forward and backward stepwise inclusion. Predictors selected by this method at each classification level are presented in Additional file 4. To identify genera that significantly contributed to total variance, we evaluated Spearman correlations between genus abundance and the first two RDA axes. OTUs with a significant correlation coefficient (evaluated at Bonferroni adjusted alpha level) were displayed on the RDA plots with type II scaling. To evaluate the association of genus abundance with explanatory variables, we implemented multiple negative binomial GLMs using the “mvabund” package [28]. This multiple GLM method utilizes a series of univariate F tests of the effects of predictor variables on the abundance of each taxon. Regression assumptions were assessed using residual diagnostics. Taxa that made up less than 0.1 % of the total count and occurring in less than 75 % of samples were first removed (cf. [29]). Fifty taxa met the inclusion criteria and were included in the model. The default implementation of the multi-GLM method adjusts P values to account for multiple tests. Classification of relative abundance data according to the previously described enterotypes [30] was carried out using the Calinski-Harabasz (CH) index as described online (http://enterotype.embl.de/enterotypes.html).

Bacterial phylogeny is sufficiently linked to their functional capabilities and can be used to computationally predict the functional composition of the community metagenome [19]. The normalized genus abundance OTU table was used to predict the microbiome’s metagenomic functions using PICRUSt’s extended ancestral-state reconstruction approach. A new abundance matrix of predicted functional categories based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was created. We constructed a biplot from the output of a PCA of functional category data and visually assessed clustering patterns based on CRF groupings. Next, to isolate the influence of specific predictor variables, an RDA was also performed using these functional categories as response variables and the same variables and selection methods described above.

Similarly, to determine the role of our exploratory variables in explaining variance in fecal SCFAs, an RDA was conducted using SCFA abundance data as the response variables (cf. [31]).