Is There a Universal Endurance Microbiota?

ions of microbes sculpt the gut ecosystem, affecting physiology. Since endurance athletes' performance is often physiology-limited, understanding the composition and interactions within athletes' gut microbiota could improve performance. Individual studies describe differences in the relative abundance of bacterial taxa in endurance athletes, suggesting the existence of an “endurance microbiota”, yet the taxa identi ed are mostly non-overlapping. To narrow down the source of this variation, we created a bioinformatics work ow and reanalyzed fecal microbiota from four 16S rRNA gene sequence datasets associated with endurance athletes and controls, examining diversity, relative abundance, correlations, and association networks. There were no signi cant differences in alpha diversity among all datasets and only one out of four datasets showed a signi cant overall difference in bacterial community abundance. When bacteria were examined individually, there were no genera with signi cantly different relative abundance in all four datasets. Two genera were signi cantly different in two datasets (VeillonellaandRomboutsia). No changes in correlated abundances were consistent across datasets. A power analysis using the variance in relative abundance detected in each dataset indicated that much larger sample sizes will be necessary to detect a modest difference in relative abundance especially given the multitude of covariates. Our analysis con rms several challenges when comparing microbiota in general, and indicates that microbes consistently or

in correlated abundances were consistent across datasets. A power analysis using the variance in relative abundance detected in each dataset indicated that much larger sample sizes will be necessary to detect a modest difference in relative abundance especially given the multitude of covariates. Our analysis con rms several challenges when comparing microbiota in general, and indicates that microbes consistently or universally associated with human endurance remain elusive. Keywords:16S gene sequences; barcode; network;Romboutsia;Veillonella; work ow 1. Introduction Endurance has played an important role in human history from our origins in the African savannah, through historical times (i.e., the origin of the marathon) to the more recent fascination with ultra-endurance events. Success in endurance events is assumed to be a product of training, genetics, and psychological preparation to withstand extreme mental and physical challenges [1]. Physiologically, the factors limiting performance in endurance events have traditionally been divided into two broad categories—aerobic (i.e., VO2max) and anaerobic (i.e., lactic acid) [2]. However, what if one of the limiting factors to success in endurance events was not human at all? We are as microbial as we are human. Bacterial cells associated with the human body are at least equivalent in number to human cells, if not more abundant [3]. More than 1000 bacterial species may be found in the intestines of each person, and over 70% of the total human microbiome is contained in the gut [4]. Whereas human cells are slowly replaced with identical (or nearly identical) copies, the cells of the gut microbiota are constantly changing due to immigration, emigration and differential rates of division based on the dynamic gut environment, to which they are much more sensitive and responsive than human epithelial cells. If our fuel tanks are lined with bacteria, it is not surprising that gut microbiota may affect performance in fuel-limited sports such as endurance events. Furthermore, the dichotomy between human physiology and microbial metabolism is increasingly dif cult to disentangle, as more human genotype/microbiome interactions continue to be discovered [5]. Microorganisms2022,10, 2213.

is not surprising that gut microbiota may affect performance in fuel-limited sports such as endurance events. Furthermore, the dichotomy between human physiology and microbial metabolism is increasingly dif cult to disentangle, as more human genotype/microbiome interactions continue to be discovered [5]. Microorganisms2022,10, 2213.

Microorganisms2022,10, 2213 2 of 22 In some studies, exercise has been associated with increased gut microbial diversity, increasedBacteroidetes–Firmicutesratio, and proliferation of bacteria which can modu- late mucosal immunity and improve barrier functions [6–11]. All of these changes could contribute to increased performance, decreased in ammation, decreased gastrointesti- nal (GI) distress, and faster recovery times [6], as well as possible protection against GI infections [7–11] . However, factors other than intensive exercise could also contribute to these reported changes in the microbiota of endurance athletes. In fact, numerous factors have been suggested to affect the composition of the gut microbiota, including (but not limited to) method of birth, diet, sex, age, antibiotic use, geographic region, stress level, and disease history [12]. For example, dietary behaviors such as carbohydrate loading, a traditional prelude to an endurance event, could increasePrevotellaindependent of physical exertion [13]. Consumption of probiotics could also affect the microbiota irrespective of exercise level (for example, see [14]). A recent review found mixed effects of probiotics on performance, as 17 studies showed no signi cant results of probiotic consumption and seven showed an improvement in performance (summarized in [13]). Probiotics may also protect against some upper respiratory tract infections which in athletes could result in longer stretches of continuous, intensive training [15]. Interactions among bacteria may be important, as multi-strain probiotics are more likely than single strains to lead to improved endurance performance [13]. The mixed results regarding the effects of probiotics on athletic performance indicate that there may be an interaction between gut bacteria and endurance, but warrant further study. The concept of an “endurance microbiome” suggests there are assemblages of gut bacteria that are more common in athletes compared to sedentary controls, or gut bacteria that are enriched after an endurance event compared to before the event. The hypothet- ical endurance microbiota could be characterized by changes at multiple levels such as (1) increased overall diversity, (2) different levels of abundance, and (3) the emergence of bene cial bacterial associations. Below, we brie y summarize these three metrics. To date, studies on the effects of exercise on gut microbiota

an endurance event compared to before the event. The hypothet- ical endurance microbiota could be characterized by changes at multiple levels such as (1) increased overall diversity, (2) different levels of abundance, and (3) the emergence of bene cial bacterial associations. Below, we brie y summarize these three metrics. To date, studies on the effects of exercise on gut microbiota diversity have exhibited mixed results. Some investigations describe increased microbial diversity associated with exercise [6,16–20], while others report no signi cant change in diversity [21–24]. Lensu and Pekkala [25] have conducted a recent and thorough literature review of the effect of exercise on the gut microbiota concluding that exercise has a “bene cial effect on the gut microbiota” and is associated with “healthy gut microbiota”. However, when the relative abundance of individual bacterial genera is examined, there is very little consistency among the studies they reviewed. Deciphering these inconsistent results requires homogenizing the bioinformatics pipeline and subsequent statistical analyses. Therefore, we reanalyzed the raw data from three studies (four comparison groups) that all utilize 16S barcoding data collected on the Illumina platform (at least in part) yet have come to largely different conclusions regarding the changes in relative abundance associated with endurance exercise. First, in a comparison of runners before and after the Boston Marathon, Scheiman et al. [26] found onlyVeillonellaexhibited a signi cant increase in abundance after the marathon. The authors then showed thatVeillonellaimproves endurance performance by metabolizing lactate that has crossed the blood–gut barrier into short-chain fatty acids that can improve muscle performance (in mice) [26]. Second, when Zhao et al. [23] examined Chinese half-marathon runners, they found 12 genera whose abundance increased after the event, not includingVeillonella. Third, Peterson et al. [18] compared competitive cyclists with variable training intensities and found some of the same genera as Zhao et al. [23], yet identi ed several new genera with differential abundance based on 16S and transcriptome data. These studies identi ed a diversity of candidate genera that may comprise endurance-associated microbiota, but due to methodological differences in their bioinformatic work ows, it is dif cult to make direct

training intensities and found some of the same genera as Zhao et al. [23], yet identi ed several new genera with differential abundance based on 16S and transcriptome data. These studies identi ed a diversity of candidate genera that may comprise endurance-associated microbiota, but due to methodological differences in their bioinformatic work ows, it is dif cult to make direct comparisons among the studies. Endurance-associated microbiota could also manifest in associations among bacteria, or groups of bacteria. For example, several studies have highlighted a putative trade-off

Microorganisms2022,10, 2213 3 of 22 between two common groups of gut bacteria implicated in endurance athletes—Prevotella andBacteroides[27,28]. Other investigations have found strong correlations among bacterial lineages and in response to environmental stresses [29]. Many exercise-related studies have looked for inter-bacterial associations.Prevotellahas been implicated in athletic performance and associated withStreptococcus,Enterococcus,Desulfovibrio,Lachnospiraceae,Succinivibrio, Oscillospira,Xylanibacter, andButyrivibrio[20]. Nevertheless, Gorvitovskaia et al. [28] report no consistent bacterial correlations among the four studies they reviewed. Potentially, bacterial interactions associated with endurance are broader than pairwise correlations, perhaps represented better as a connectivity network of the entire microbiota community. Active people who exercise regularly have been reported to have more complex gut bacterial networks than sedentary controls [17]; however, the connectivity networks underlying endurance-associated microbiota are largely unexplored [30], and no meta-analysis has compared networks using the same methodology across datasets. We hypothesize that if a “universal” endurance-associated assemblage of microbes exists, it should manifest regardless of geographic location, type of sport, or speci c endurance event. Our goal is to reanalyze several relevant studies using a single bioinfor- matics pipeline and consistent downstream statistical analyses to determine whether there are repeated changes in diversity, relative abundance, or associations between genera in response to intensive endurance exercise. The variance in bacterial abundances were used in a power analysis to determine the necessary sample sizes to detect a modest difference. 2. Materials and Methods 2.1. Datasets After searching the literature for relevant studies (Table S1), we selected and rean- alyzed four gut microbiota datasets involving endurance athletes from three previously published studies (Table). All studies utilized Illumina sequencing of the V3-V4 region of the 16S rRNA gene (amplicon sequencing) [18,23,26]. Below is a brief description of each dataset. Table 1. Datasets used in this work and number of microbial genera detected in our bioinformat- ics pipeline. Event and Location Treatment Group Sample Size Sampling Frequency Reference No. of Genera Detected Boston Marathon, Boston, MA, USA Runners Before 15 Multiple samples taken before the event [26] 221 Boston Marathon, Boston, MA, USA Runners After 15 (paired with above) Multiple samples taken after the event [26] 233 Boston Marathon, Boston, MA, USA

detected in our bioinformat- ics pipeline. Event and Location Treatment Group Sample Size Sampling Frequency Reference No. of Genera Detected Boston Marathon, Boston, MA, USA Runners Before 15 Multiple samples taken before the event [26] 221 Boston Marathon, Boston, MA, USA Runners After 15 (paired with above) Multiple samples taken after the event [26] 233 Boston Marathon, Boston, MA, USA Sedentary Controls 10 Multiple samples taken from controls [26] 228 Chongqing International Half Marathon, Chongqing, China Runners Before 20 runners 1 Once before the event [23] 194 Chongqing International Half Marathon, Chongqing, China Runners After 20 (paired with above) Once after the event [23] 197 Competitive Cyclists, USA 2 Low (6–10 h/wk) 8 One time point [18] 115 Competitive Cyclists, USA 2 Medium (11–15 h/wk) 17 One time point [18] 133 Competitive Cyclists, USA 2 High (16–20+ h/wk) 8 One time point [18] 115 1 One before-event sample (BEF09) could not be used because it had an inconsistently formatted fastq le (Zhao attempted personal communication, unrequited). 2 Petersen et al. [18] reports very high consistency in relative abundances estimated from both 16S amplicon sequencing and whole-genome shotgun sequencing. We reanalyzed the 16S results for direct comparison with other endurance datasets.

Microorganisms2022,10, 2213 4 of 22 2.1.1. Boston Marathon Study Scheiman et al. [26] recruited 15 elite athletes running in the 2015 Boston Marathon, along with 10 sedentary controls. They conducted amplicon sequencing on 209 fecal samples taken daily from participants up to one week before to one week after the marathon using Earth Microbiome Project primers targeting the v4 region of the 16S rRNA gene [using primers 515F (Caporaso) and 806R (Caporaso); and-standards/16s/ onE. coli16S rRNA REFSEQ NR_024570) were sequenced using 150 bp Illumina paired-end reads and were processed with the DADA2 pipeline and phyloseq. Generalized linear mixed-effect models and leave-one-out cross validation were used to determine signi cant associations. According to their supplemental data, some samples were “rerun” on the Illumina sequencer. In these cases, we only used the data from the reruns (i.e., SG10 and SG27). 2.1.2. Chongqing Half-Marathon Study Zhao and colleagues [23] recruited 20 amateur athletes who were running in the 2016 Chongqing International Half Marathon. A total of 40 fecal samples were collected—each runner was sampled the morning before the race and again after the race. Zhao et al. used very similar reverse primers to Scheiman et al.; however, they added CC to the 3 0 end and used T instead of the ambiguity W nine base pairs from the 3 0 end. Their forward primer lands 176 bp further upstream than the Scheiman et al. primers generating a ~465 bp amplicon (based onE. coli16S rRNA REFSEQ NR_024570). Because of the larger amplicon length, they collected 250 bp paired-end reads generated on an Illumina HiSeq. Two samples from runner nine were eliminated from our reanalysis because of ambiguous labeling of the data. Each participant was given the same kind of food during the period between the rst and second sample collection. 2.1.3. Competitive Cyclist Study Petersen et al. [18]studied 33 competitive cyclists categorized into four non-overlapping training groups based on their average training time per week: 6–10 h, 11–15 h, 16–20 h and 20+ h per week. A total of 33 samples were collected, one from each cyclist. They collected both whole-genome shotgun sequence data

the rst and second sample collection. 2.1.3. Competitive Cyclist Study Petersen et al. [18]studied 33 competitive cyclists categorized into four non-overlapping training groups based on their average training time per week: 6–10 h, 11–15 h, 16–20 h and 20+ h per week. A total of 33 samples were collected, one from each cyclist. They collected both whole-genome shotgun sequence data and 16S rRNA gene amplicon sequence data. We reanalyzed the latter to compare to the two studies (Petersen et al. 2017, Additional File 1). Petersen et al. used 150 bp paired-end Illumina reads from 16S rRNA amplicons of the v4 hypervariable region with primers comparable to Scheiman et al. [24]. We grouped the cyclists into low (6–10 h/wk;n= 8), medium (11–15 h/wk;n= 17), and high (16–20+ h/wk; n= 8) categories for diversity analyses, then focused on the two most extreme training groups (with balanced sample sizes) when searching for a universal endurance microbiota (e.g., “low” vs. “high” training). 2.1.4. Sample Collection, Storage and DNA Extraction Fecal sample collection and storage and affect the estimates of relative abundance from 16S amplicon sequencing [31]. Fecal samples from all three studies reanalyzed herein were self-collected; however, Scheiman et al. [24] and Zhao et al. [23] used polypropylene tubes for collection and stored samples at 4 C short-term, while Petersen et al. [18] used polyethylene tubes and stored samples with “frozen freezer packs” short-term. All three studies 80 C for long-term sample storage. Although different DNA extraction methods were used across all three studies, Rintala et al. demonstrate that the impact of different DNA extraction methods on relative abundance estimates from 16S amplicon sequencing is relatively minor [31].

Microorganisms2022,10, 2213 5 of 22 2.2. Target Genera Among the Boston Marathon, Chongqing half marathon, and competitive cyclist studies, there were 1, 12, and 6 bacterial genera identi ed as having signi cantly different abundances between treatment groups, respectively (Table S2). Of these, no single genus was identi ed in all three studies and only three genera were found in two of these studies. We used the 16 unique genera from all three studies as our “target genera” in a hypothesis- testing framework (using alpha < 0.05 as a cutoff). Subsequently, we expanded our analyses to all remaining genera, and corrected for multiple tests, since the comparisons were not based on a priori hypotheses (i.e., using the Benjamini–Hochberg false-discovery-rate correction—see below for details). The target “genera” from Zhao et al. [23] included both individual species (e.g.,Pre- votella corporiswhich we treated asPrevotella) and genera. The authors emphasized the signi cantly differential abundance in the family Coriobacteriaceae before and after the half marathon, but we only includedCollinsellasince we assume this genus was driving the signi cant result based on their Figure 2B. We did not include “unclassi ed Porphyromon- adaceae”. Finally, they report “Phaseolus vulgaris” in their Figure 2 results [23], but they discussRomboutsialater in their ndings. We assume this is a technical mistake (Phasaeolus vulgarisis a species of legume, not a bacteria).Phaseolussensu Zhao et al. [23] is hereafter treated asRomboutsia. 2.3. Microbiome Assembly (Bioinformatics Pipeline) Although there are many published tools for measuring bacterial abundance from 16S rRNA gene amplicon sequencing using the Illumina platform (e.g., QIIME [32]; MOTHUR [33]; DADA2 [34]), we developed a simple work ow in Geneious Prime 2021.2.2 by adapting their Amplicon Metagenomics tutorial [(https://www.geneious.com/ tutorials/metagenomic-analysis/ studies was downloaded from usegalaxy.org in the form of fastq les and imported into Geneious Prime. We used Illumina paired-end, inward pointing reads with sequences inter- laced within each fastq le. The minimum quality (q) cutoff was empirically determined to be 13 in a pilot study, and the merge rate was set at “very high” to maximize reads mapped, while reducing incorrectly mapped reads. After using BLAST to Genbank (Release 242)

and imported into Geneious Prime. We used Illumina paired-end, inward pointing reads with sequences inter- laced within each fastq le. The minimum quality (q) cutoff was empirically determined to be 13 in a pilot study, and the merge rate was set at “very high” to maximize reads mapped, while reducing incorrectly mapped reads. After using BLAST to Genbank (Release 242) for OTU clustering to construct a reduced 16S rRNA sequence database per dataset, we mapped reads back onto the reduced database to classify the reads to genus based on 90% minimum overlap identity per Geneious' Amplicon Metagenomics tutorial and their Sequence Classi er tutorial. Although 95% 16S rRNA gene sequence identity has often been used as a cut-off for bacterial genus-level operational taxonomic units [35–38], our more liberal cut-off was intended to maximize the number of classi ed reads while accounting for the sequence variation within genera. Using this cutoff, we were not trying to detect bacterial species, only genera for consistency in making comparisons across studies. All subsequent analyses were based on the relative abundance which is the proportion of reads mapped to a genus compared to total reads mapped per sample (per Gloor et al. [39]). Clas- si cations were pruned to genus, with all higher taxonomic-level BLAST results removed. Our Geneious work ow is available from FigShare (DOI: 10.6084/m9. gshare.c.6036347). 2.4. Diversity We measured the alpha diversity for each dataset and each treatment group as the number of unique genera identi ed. We then used Simpson and Shannon indices to compare diversity considering each genus' relative abundance in the vegan package (v. 2.6-2) in R Studio (2022.02.3). For the competitive-cyclist study [18], diversity measures among the three independent treatment groups were compared using an ANOVA in R. All other diversity comparisons for the Boston marathon and Chongqing half-marathoners were carried out usingt-tests on independent samples (athletes vs. controls) or paired samples (before vs. after) in R.

This study analyzes gut microbiota in endurance athletes to identify potential universal microbial associations.