Higher physical activity levels are related to faecal microbiota diversity and composition in young adults

Design study and participants

A total of 92 (65 women) young healthy adults, aged 18–25 years, were included in the present cross-sectional study. This study was carried out within the framework of the ACTIBATE study [21], an exercise-based randomized controlled trial (Clinical Trials.gov ID: NCT02365129). All assessments were performed in Granada (Spain) between October and November in 2016. Inclusion criteria were: being engaged in less than 20 min of moderate-vigorous PA on less than 3 days/week, having a stable body weight over the last 3 months (< 3 kg change), not smoking, not taking any medication (including antibiotics in the last 3 months), not presenting any acute or chronic illness and not being pregnant. The study protocol and experimental design were applied in accordance with the last revised ethical guidelines of the Declaration of Helsinki. The study was approved by the Ethics Committee on Human Research of the University of Granada (no. 924) and the Servicio Andaluz de Salud (Centro de Granada, CEI-Granada); all participants signed informed consent.

Physical activity assessment

PA variables were objectively measured with one accelerometer on the non-dominant wrist (ActiGraph GT3X+, Pensacola, FL), during 7 consecutive days (24 h/day) [21]. Detailed information about how to wear the accelerometer was given to participants, including the instruction to remove it in daily water-based activities, such as washing dishes or showering.

The sampling frequency of 100 Hz was selected to store the raw accelerations of the accelerometers [22]. We exported and converted the raw accelerations to the “.csv” format using ActiLife v.6.13.3 software (ActiGraph, Pensacola, FL, US). Afterwards, the “ggir” [23] package in R software was used to process the raw “.csv” files. This processing consisted of: (i) local gravity data auto-calibration of accelerations according to the local gravitational acceleration [24], (ii) calculation of the Euclidean Norm of the raw accelerations Minus One G with negative values rounded to 0 (ENMO) calculated elsewhere [25], (iii) detection of non-wear time based on the raw acceleration of the three axes, (iv) determination of MAL detection of sustained functioning of the accelerometer by means of abnormal high accelerations incompatible with human movement (i.e., related to device malfunctioning), (v) imputation of non-wear time and abnormal high accelerations, (vi) identification of waking and sleeping time based on the automatized algorithm guided by the participants’ daily reports [26], and (vii) estimation of sedentary time and the time spent in light PA, moderate PA, vigorous PA, and moderate to vigorous PA using agespecific cut-points for a wrist-worn accelerometer, for Euclidean Norm Minus One (ENMO) [27]. We measured the mean ENMO (mg) during waking time, which is considered an overall indicator of the PA (overall PA). For the analyses we only included the participants who wore the accelerometers for ≥ 16 h/day during at least 4 days (including at least 1 weekend day).

Faecal microbiota analyses

Sequencing analysis

DNA extracted was amplified by polymerase chain reaction (PCR) by primer pairs – forward primer (5’CCTACGGGNGGCWGCAG3’) and reverse primer (5’GACTACHVGGGTATCTAATCC3’) – targeting the V3 and V4 hypervariable regions of the bacterial 16S rRNA gene. All PCRs were executed in 25 µL reaction volumes incorporating 12.5 µL of 2X KAPA HiFi Hotstart ready mix (KAPA Biosystems, Woburn, MA, USA), 5 µL of each forward and reverse primer (1 µM) and 2.5 µL of extracted DNA (10 ng) under the following cycling circumstances: (a) denaturation at 95°C for 3 min, (b) cycles of denaturation at 95°C for 30 s, (c) annealing at 55ºC for 30 s, (d) elongation at 72ºC for 30 s , (e) a final extension at 72°C for 5 min. To purify the PCR products from free primers and primer dimers we used AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA). Next, the index PCR attached dual indices and Illumina sequencing adapters using the Nextera XT Index Kit (Illumina, San Diego, CA, USA), on a thermal cycler using the requirements previously mentioned. After that, AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA) were used for purification of the pooled PCR products. The resultant amplicons were sequenced at MiSeq (Illumina, USA), using a paired-end (2 × 300 nt) Illumina MiSeq sequencing system (Illumina, San Diego, CA, USA).

Bioinformatics analysis

We analysed the FASTQ files with the “dada2” [28] package in R software, obtaining 11,659,014 paired ends with an average of 126,728 ± 33,395 reads per sample. The cut-off of 10,000 reads was surpassed for all samples. Samples were resampled to a minimum sequencing depth of 30,982 reads using the “phyloseq” [29] package in R software, returning 11,158 phylotypes.

The “Classifier” function from the Ribosomal Database Project (RDP) was used to assign the taxonomic affiliation of phylotypes, based on the naïve Bayesian classification [30] with a pseudo-bootstrap threshold of 80%. We obtained a total of 209 genera belonging to 16 different phyla. The “seqmatch” [31] function from RDP was performed to define the discriminatory power of each sequence read with the purpose of annotating species assignments; we executed annotation according to previously published criteria [32]. Microbial communities were analysed at different taxonomic levels (phylum to genus), calculating relative abundances, expressed as percentages.

We performed the analyses with those bacteria with more than 0.5% on average in their relative abundance.

Next, alpha and beta diversities were estimated based on the identified microbial communities. Alpha diversity takes into account the number of different phylotypes and relative abundances of a single sample [33], whereas beta diversity shows differences in microbial community composition between individuals, which is the degree to which samples differ from one another [34]. Alpha diversity was assessed based on Chao richness, Shannon, inverse Simpson and evenness Camargo indexes with the “microbiome” [35] package in R software. Chao richness estimates diversity according to the number of different phylotypes in the community [36]; that is, higher Chao richness indicates higher diversity in the community. Shannon diversity increases as both the richness and the evenness of the community increase [37]; the inverse of Simpson diversity is calculated from the classical Simpson diversity and indicates richness in a community with uniform evenness [38], and evenness Camargo determines the equitability of phylotype frequencies in a community [39]. Beta diversity was measured quantitatively using permutational multivariate analysis of variance (PERMANOVA) based on Bray-Curtis dissimilarity.

Anthropometric and body composition measurements

Participants’ weight and height were measured, without shoes and wearing the standard clothes, using a SECA scale and stadiometer (model 799, Electronic Column Scale, Hamburg, Germany). Body mass index (BMI) was calculated as weight (kg)/height (m²). Body composition was evaluated by dual energy X-ray absorptiometry (DEXA, HOLOGIC, Discovery Wi, Marlborough, MA). The lean mass index (LMI) and fat mass index (FMI) were calculated as lean body mass and fat body mass, respectively, in kg, divided by height in m². The fat mass percentage was determined as the body fat mass divided by the total body mass and multiplied by 100.

Cardiometabolic profile

Fasting blood samples were collected for assessment of the cardio-metabolic profile. Serum glucose, total cholesterol, high density lipoprotein-cholesterol (HDL-C) and triglycerides were measured following standard methods using an AU5832 automated analyser (Beckman Coulter Inc., Brea CA, USA). Low-density lipoprotein cholesterol (LDL-C) was estimated as: [total cholesterol – HDL-C – (triglycerides/5)], in mg/dL. Serum insulin was measured using the Access Ultrasensitive Insulin chemiluminescent immunoassay kit (Beckman Coulter Inc., Brea CA, USA). The homeostatic model assessment for insulin resistance (HOMA-IR) index was calculated as (insulin (µU/mL) × glucose (mmol/L)/22.5.

Dietary assessment

Dietary intake was registered using three non-consecutive 24-hour recalls, 2 weekdays and a weekend day. These 24-hour recalls were performed in the laboratory via face-to-face interviews with dietitians. To improve the accuracy of food quantification, we used coloured photographs of different portion sizes of food during the interviews [40]. All 24-hour recalls were analysed for total energy (kcal), fat, proteins, carbohydrates, and fibre intake (g) by EvalFINUT software, which is based on the United States Department of Agriculture (USDA) and “Base de Datos Española de Composición de Alimentos” (BEDCA) databases.

Statistical analysis

This is a secondary study derived from the ACTIBATE trial [41]; therefore, there is not a sample size calculation for this study. Data normality was explored using the D’Agostino & Pearson omnibus, visual histograms and Q-Q plots (data not shown). None of the variables followed a normal distribution; therefore data were presented as median ± interquartile range and non-parametric tests were used for all analyses. Moreover, no sex interaction was detected (all P > 0.05), so both sexes were pooled together. Spearman correlations were performed to investigate the correlation between the PA variables and faecal microbiota diversity, using the “psych” [42] and “corrplot” [43] packages in R software. Since faecal microbiota diversity can be modified by several factors including sex [44], BMI [45] and dietary intake [46], we repeated the aforementioned correlations adjusted for sex, BMI and dietary intake in separate models (data not shown). Moreover, we repeated this analysis by adjusting for accelerometer non-wear time and glucose levels in separate models (data not shown) as possible confounders of PA variables. Overall PA and the time spent in vigorous PA were computed as tertiles according to number of participants with SPSS (SPSS v. 22.0, IBM SPSS Statistics, IBM Corp. Armonk, NY), because they were the only variables with a significant correlation with faecal microbiota diversity. The tertiles’ values for overall PA were low (13.45–29.44 mg), intermediate (30.02–35.41 mg), and high (35.49–67.10 mg), whereas for the time spent in vigorous PA the values were low (0.02–0.83 min/day), intermediate (0.87–2.67 min/day), and high (2.75–14.40 min/day). Tertiles of overall PA and the time spent in vigorous PA were compared using one-way PERMANOVA with 9,999 permutations for significance testing with the Paleontological Statistics (Past3) software [47] for the calculation of beta diversity. Kruskal-Wallis tests were performed to investigate whether there were significant differences in body composition, dietary intake and cardiometabolic profile as well as faecal microbiota alpha diversity and composition outcomes across tertiles of overall PA and the time spent in vigorous PA. Analysis of covariance was used to compare the relative abundance of genera across tertiles of the time spent in vigorous PA adjusted for protein intake with the data transformed by Blom’s formula. All P values were corrected by Benjamini and Hochberg multiple testing to control the false discovery rate (FDR, shown as q-values) [48]. The level of significance was set at P < 0.05 and q < 0.05. R software (V.3.6.0; http://www.r-project.org) and GraphPad Prism version 8.0.0 for Windows (GraphPad Software, San Diego, California, USA, (http://www.graphpad.com) were used for the statistical analysis and graphical plots.

Characteristics of participants

A total of 92 participants had data from analysis of faecal microbiota diversity and composition, but only 88 participants (24 men, age = 22.0 ± 2.0; and 64 women, age = 21.6 ± 2.0) had valid PA measurements (as they wore the accelerometer for < 16 h/day during at least 4 days), who were finally included in the analyses. Table 1 shows the descriptive characteristics of the included participants (age 21.7 (19.8–23.9) years and BMI 23.6 (21.6–28.1 kg/m²)), of whom 72.7% were women. We performed tertiles of overall PA and the time spent in vigorous PA and we observed that, generally, body composition, dietary intake and cardiometabolic profile were similar across them (Table S1), with the exception of protein intake and glucose levels (P = 0.018 and P = 0.003, respectively; Table S1).

Relationship between physical activity and faecal microbiota diversity

Overall PA and the time spent in vigorous PA were positively correlated with alpha diversity indexes, more specifically with Shannon and inverse Simpson diversity indexes (all rho ≥ 0.23, P ≤ 0.034; Fig. 1). Only the time spent in vigorous PA was positively correlated with the Chao richness index (rho = 0.24, P = 0.023; Fig. 1). However, we did not observe any significant correlation between other PA variables and alpha diversity indexes (all P > 0.05; Fig. 1). The results were similar when sex, BMI, energy and macronutrient intake, as well as accelerometer non-wear time and glucose levels, were included as confounders in separate models (data not shown). Moreover, we found that individuals with high time spent in vigorous PA had a higher Chao richness, Shannon and inverse Simpson diversity indexes than individuals with low and intermediate time spent in vigorous PA (all P ≤ 0.038; data not shown). However, there were no differences across tertiles of overall PA and the time spent in vigorous PA in the beta diversity at any taxonomic levels (all P ≥ 0.060; Table 2).

FIG. 1

Spearman correlation of physical activity variables with faecal microbiota diversity. Boxes only represent the statistically significant (P < 0.05) correlations and the values within the boxes show the Spearman correlation coefficient. Blue boxes indicate a positive correlation whereas red squares indicate a negative correlation between physical activity variables and faecal microbiota diversity indexes. mg: milli-gravitational units; PA: physical activity.

Relationship between physical activity variables and faecal microbiota composition

We analysed the differences across tertiles of overall PA and the time spent in vigorous PA on faecal microbiota composition at all taxonomic levels. There were no significant differences across tertiles of overall PA on the relative abundance of all bacteria at the different taxonomic levels (all P > 0.05; Fig. 2). Similarly, we observed no differences across tertiles of time spent in vigorous PA on the relative abundance of bacteria at the phylum taxonomic level (all P ≥ 0.318; Fig. 3A). However, we observed that individuals with low time spent in vigorous PA had a higher relative abundance of the Gammaproteobacteria class (Proteobacteria phylum) than individuals with intermediate time spent in vigorous PA (q = 0.021, FDR = q-value; Fig. 3B). Moreover, individuals with high time spent in vigorous PA had higher relative abundance of unclassified Firmicutes class (Firmicutes phylum) and Porphyromonadaceae family (Bacteroidetes phylum) than individuals with low time spent in vigorous PA (q = 0.027, and q = 0.031, respectively; Fig. 3B and D). Finally, we found that individuals with intermediate time spent in vigorous PA had a higher relative abundance of the Alistipes genus (Bacteroidetes phylum) than individuals with low time spent in vigorous PA (q = 0.015; Fig. 3E). Interestingly, these same participants presented differences in protein intake (q = 0.011; Table S1). Thus, we repeated the analyses after adjusting for the protein intake and the differences in the relative abundance of the Alistipes genus disappeared (P = 0.080; Table S2).

FIG. 2

Faecal microbiota composition according to tertiles of overall physical activity (PA) levels (L: low, 13.45–29.44 mg; I: intermediate, 30.02–35.41 mg; H: high, 35.49–67.10 mg). Panels indicate relative abundance of the faecal microbiota at phylum (A), class (B), order (C), family (D) and genus (E) taxonomic levels according to tertiles of overall PA. Stacked bar represents percentage abundance. Kruskal-Wallis test was used to test for each pairwise comparison, correcting for multiple comparisons FDR (q < 0.05) (GraphPad Prism 8.00).

FIG. 3

Faecal microbiota composition according to tertiles of the time spent in vigorous physical activity (PA) (L: low, 0.02–0.83 min/day; I: intermediate, 0.87–2.67 min/day; H: high, 2.75–14.40 min/day). Panels indicate relative abundance of the faecal microbiota at phylum (A), class (B), order (C), family (D) and genus (E) taxonomic levels according to tertiles of the time spent in vigorous PA. Stacked bar represents percentage abundance. Symbols * and ◦ mean statistical significance differences between low and high time spent in vigorous PA, and symbol • represents statistical significance differences between low and intermediate time spent in vigorous PA. Kruskal-Wallis test was used to test for each pairwise comparison, correcting for multiple comparisons FDR (q < 0.05) (GraphPad Prism 8.00).

DISCUSSION

In the present study, overall PA and the time spent in vigorous PA were found to be positively correlated with alpha diversity indexes in young adults. Moreover, there were differences across the tertiles of time spent in vigorous PA in the relative abundance of the Gammaproteobacteria class (Proteobacteria phylum), Porphyromonadaceaefamily and Alistipes genus (both Bacteroidetes phylum). These findings indicate that PA may play a role in faecal microbiota diversity and composition in young adults, although further studies are needed to confirm these findings.

Our results showing the positive correlation between PA and alpha diversity agree with recent findings [14, 49, 50]. However, the mechanisms by which PA may promote higher faecal microbiota diversity are unknown. A possible explanation could be the changes in the gastrointestinal tract due to intrinsic adaptations of performing PA [51]. Interestingly, from an ecological perspective, microbial diversity may be a key factor in allowing an ecosystem to continue operating properly [52]. In fact, greater species diversity has been associated with a healthy phenotype’s host [53]. This is due to the potential effects that the bacteria can exert via metabolites, such as short-chain fatty acids and neurotransmitters locally and extra-intestinal tissues in the host [54].

Our data showed that the participants with low time spent in vigorous PA had higher relative abundance of the Gammaproteobacteria class than individuals with higher time spent in vigorous PA. The relative abundance of the Gammaproteobacteria class (Proteobacteria phylum) has been reported to be increased in obese mice [55] and individuals with obesity [56], and disease states such as metabolic diseases and intestinal inflammation [57]. In fact, many common human pathogens, known as sulphur producers [58], are found in the Gammaproteobacteria class, for example, Escherichia, Shigella, and Yersinia genera [58]. In agreement with our findings, sedentary women and participants with low cardiorespiratory fitness [59] had a higher relative abundance of the Gammaproteobacteria class than active women and participants with high cardiorespiratory fitness, respectively. Similarly, a very recent study performed in > 8,000 individuals using accelerometers observed that PA levels were associated differently with faecal microbiota composition, suggesting that the higher the PA level is, the higher is the diversity [60]. Moreover, several studies have shown that exercise seems to decrease the relative abundance of the Gammaproteobacteria class [61, 62]. Thus, our results suggest that performing less than 1 min/day of vigorous PA could be related to having a higher relative abundance of the Gammaproteobacteria class, bacteria considered health-detrimental.

In contrast, we observed that individuals with high and intermediate time spent in vigorous PA had a higher relative abundance of the Porphyromonadaceae family and Alistipes genus (both Bacteroidetes phylum) than individuals with lower time spent in vigorous PA. Accordingly, in a cross-sectional study of professional martial arts athletes, the relative abundance of the Porphyromonadaceae family was higher in the higher-level athletes in comparison with the lowerlevel athletes [63]. Moreover, regular swimming training [64] and voluntary wheel running [65], both in mice, were able to increase the relative abundance of the Porphyromonadaceae family. In fact, it has recently been found that lean individuals had a significantly higher relative abundance of the Porphyromonadaceae and Rikenellaceae families than individuals with obesity [45]. Of note is the fact that the Alistipes genus belongs to the Rikenellaceae family. In resistance-trained mice, the relative abundance of the Alistipes genus was positively correlated with resistance performance [66]. In humans, the relative abundance of the Alistipes genus is increased after consuming an animal-based diet intake, rich in protein, for 5 days [64]. Certain species that belong to the Alistipes genus are involved in amino acid metabolism; specifically, they can hydrolyse tryptophan to indole [67]. Since tryptophan is an essential amino acid that cannot be produced by animal cells, humans rely on dietary intake, mainly proteins, for incorporating it into the organism [68]. In our study, the individuals with intermediate time spent in vigorous PA had higher protein intake than individuals with low time spent in vigorous PA. In fact, when the protein intake was included as a confounder, the differences in the relative abundance of the Alistipes genus between these individuals disappeared. Considering the relationship between the Alistipes genus and protein metabolism [67], and the results observed in the present study, it seems possible that these differences were explained by protein intake. Therefore, our data suggest that spending time on vigorous PA, in the range 3–14 min/day, could be related to having a higher relative abundance of Porphyromonadaceae family bacteria, whereas the protein intake seems to modulate the relative abundance of the Alistipes genus in individuals with intermediate time spent in vigorous PA. Even so, the possible effect of time spent in vigorous PA on the relative abundance of the Gammaproteo-bacteria class, Porphyromonadaceae family and Alistipes genus deserves further analysis.

Limitations and strengths

A limitation to consider in the current study is that it followed a cross-sectional design, which prevents a causal interpretation of our results. Well-designed randomized controlled trials should be carried out to elucidate the role of PA in faecal microbiota diversity and composition. In addition, we do not know whether our findings apply to older people or individuals presenting any metabolic disease. As for strengths of this study, we sequenced the microbiota composition using the latest technology (Illumina platform) and annotations were made with RDP to the genus taxon level. Moreover, PA was objectively measured by accelerometry during 7 consecutive days (24 h/day) [21], and we used a cut-point-free approach to assess overall PA since PA intensities estimated from cut-points might be biased by poor calibration studies [69].