Meta-analysis To Define a Core Microbiota in the Swine Gut

1The mammalian gut microbiota is populated by an estimated 100 trillion bacterial cells that provide health benefits to the host through the production of short-chain fatty acids and vitamins from otherwise nondigestible food components, inhibition and prevention of colonization by pathogens, and development and maintenance of the host immune system. These properties also make the gut microbiota an attractive target for dietary interventions to improve the health and production of food-producing animals, including swine. Therefore, the characterization of the gut microbiota of swine has become an active area of research in recent years as the adaptation and availability of high-throughput sequencing methods have continued to expand.

2The majority of swine gut microbiota studies to date have focused on sequencing specific hypervariable regions within the 16S rRNA gene. The myriad of studies that have sequenced the 16S rRNA gene can be readily aggregated and compared in a meta-analysis that may reveal patterns among individual studies working on similar sample types. Recent meta-analyses of the avian gut microbiota and the indoor microbiota have reported that study-level variables, compared to biological factors such as sample type, have the greatest effect on the microbiota. In swine, researchers have examined the effect of a number of different feed additives on the swine gut microbiota, including antimicrobial agents, amylose/amylopectin, calcium-phosphorus, distillers’ dried grains with solubles, and resistant starch. The effect of challenge with Escherichia coli and Salmonella on the swine gut microbiota has also been investigated. All of those studies found various degrees of change in the swine gut microbiota as a result of these treatments. Aggregating and reanalyzing the data could reveal whether alterations observed in individual studies are still detectable when combined with data from similar studies or whether they would be lost due to technical variation between studies.

3Our objectives were to characterize the swine gut microbiota using 16S rRNA gene high-throughput sequencing data from gastrointestinal (GI) samples from 20 publically available studies and to determine if a core gut microbiota, defined here as being present in ≥90% of samples, exists in swine. We evaluated how technical variables, such as study design, country of origin, and hypervariable region sequenced, affect the observed swine gut microbiota.

4A total of 25,157,796 quality-filtered 16S rRNA gene sequences from 939 swine gastrointestinal samples were available for meta-analysis (Table 1). Overall, 18,743,484 (74.5%) of these sequences were at least 97% similar to a sequence in the SILVA database. These sequences represented 25,182 bacterial operational taxonomic units (OTUs), 35 phyla, and 887 genera prior to random subsampling. It was necessary to randomly subsample to a level of 1,000 sequences per sample to retain as many samples as possible; however, as few as 100 sequences per sample have been shown to be sufficient when large differences between sample groups are present. Not surprisingly, the majority of sequences available for meta-analysis were derived from swine feces, as these samples are the easiest to obtain from the GI tract and allow repeated sampling over time. Other regions of the swine GI tract, such as the ileal and colonic mucosa and digesta, have also been well represented among publically available sequences (Fig. 1). Conversely, the number of samples from the duodenum and jejunum (digesta and mucosa) was relatively small (n ≤ 5 samples) and belonged to a single study. Therefore, these samples were not included in the linear discriminant analysis (LDA) effect size (LEfSe) described below. Sequences were also available from three different continents and 10 countries, although approximately half of these were from studies conducted in the United States. The age of the pigs used in each study ranged from~3 weeks (preweaning piglets) to ~24 weeks (slaughter age).

5The idea of defining a core microbiota in the pig gut is intriguing as it may identify potential targets for dietary or therapeutic interventions in the swine production environment. Although a core microbiota may not exist in swine according to a strict definition that requires a particular taxon to be present in each sample, several phyla and genera were found in more than 90% of all GI samples. In the current meta-analysis, the Firmicutes and Bacteroidetes phyla accounted for nearly 85% of the total 16S rRNA gene sequences among all gastrointestinal locations, with Proteobacteria being the only other phylum common to all GI samples (Fig. 2A). No OTUs or genera were shared among all samples analyzed, but the genera Clostridium, Blautia, Lactobacillus, Prevotella, Ruminococcus, and Roseburia, the RC9 gut group, and Subdoligranulum were found in more than 90% of samples. Prevotella was also the most abundant genus among all genera that were identified (Fig. 2B; see also Table S1 in the supplemental material). In addition, a single OTU classified as an uncultured member of the Prevotella was present in at least 75% of all GI samples.

6Gut locations and feces yielded discrete core microbiotas. In the fecal microbiota, more than 99% of fecal samples (12 studies) contained sequences from Prevotella, Clostridium, Alloprevotella, Ruminococcus, and the RC9 gut group. Certain genera were present in all samples from a particular GI tract location (n = >50 samples). These included Prevotella and Blautia in samples of colonic digesta (n = 91; 5 studies); Anaerovibrio, Clostridium, Phascolarctobacterium, Ruminococcus, Sarcina, and Streptococcus in cecal digesta (n = 52; 4 studies); and Clostridium in ileal mucosal samples (n = 91; 3 studies). Clostridium spp. were also present in all but one sample of ileal digesta (n = 81; 5 studies).

7In commercial swine production, pigs are typically fed a cereal grain-based diet that is relatively high in carbohydrate content. Meanwhile, Prevotella spp. have been reported to be associated with dietary carbohydrates in humans, producing polysaccharidases such as glucanase, mannanase, and xylanase. Members of Prevotella are also relatively abundant in the rumen and/or gut microbiota of other animals such as cattle, chimpanzees, goats, and sheep, attesting to the ubiquity of this genus in the mammalian gut. Clostridium, Blautia, and Ruminococcus are all members of the Clostridiales order, and, similarly to Prevotella, are widely found in the mammalian gut. These genera can produce butyrate, a short-chain fatty acid (SCFA), most often from acetate (also an SCFA) via the butyryl-coenzyme A (CoA):acetate CoA-transferase pathway. Prevotella spp., which produce acetate in the gut, are thereby able to provide a source of energy for butyrate-producing bacteria. Importantly, butyrate decreases inflammation in the gut of the host, and cells in the intestinal epithelium can use it as an energy source.

8The RC9 gut group, which belongs to the Rikenellaceae family, was identified in 90.5% of all GI samples and was among the most abundant genera in the lower GI tract (Fig. 2B; Table S1). This finding is noteworthy as its presence in the swine gut microbiota was reported only recently. The Rikenellaceae family itself is a relatively new taxonomic classification within the Bacteroidales order, with only three genera currently described. Similarly to Prevotella, also a member of the Bacteroidales, some bacteria in the Rikenellaceae family are acetate producers, suggesting a high level of functional redundancy for this metabolite in the swine gut microbiome.

9We determined which variables most strongly affected the structure of the swine gut microbiota using a permutational multivariate analysis of variance (PERMANOVA) test of the unweighted and weighted UniFrac distances and Bray-Curtis dissimilarities (Table 2). Among all distance metrics, each metadata category tested was significantly associated with the microbial community structure. However, this finding is likely a result of the heterogeneity of dispersion within each group in each metadata category, as shown by the permutational analysis of multivariate dispersions (PERMDISP; P < 0.05). Overall, “study” was the most important factor affecting the swine gut microbiota, as it explained the greatest variation in the data set (Fig. 3; also see Fig. S1 in the supplemental material). The GI location that was sampled was also strongly associated with the structure of the microbiota as assessed using weighted UniFrac distances (Fig. 3A); however, this finding was not replicated in the results obtained with respect to the unweighted UniFrac distances and Bray-Curtis dissimilarities (Fig. 3B and C). Therefore, researchers could reasonably expect that the taxa that are relatively abundant in their swine gut samples are similar to those in samples taken from the same GI location of different pigs compared with those found in samples obtained from different GI locations.

10We evaluated the effect of age on the swine gut microbiota using samples from the 19 studies (of 20) that provided specific information for pig age (n = 920; Table 2; Fig. S2). Although the effect of the age of the pigs sampled was not as large as the effects of study and GI sampling location, it did explain more variation in the data set than the remaining three metadata categories for the weighted UniFrac distances. It should be noted that pig age, like several of the other metadata categories, is also associated with “study,” as many studies sampled from pigs at only one particular time point. In addition, the 3-week-old category included both preweaned and postweaned pigs, and weaning is known to result in significant changes to the swine gut microbiota.

11Genera that were enriched in specific GI samples were identified using linear discriminant analysis (LDA) effect size (LEfSe) (Fig. 4). Among the more notable genera, Lactobacillus was more abundant in the gastric mucosa, Prevotella, Helicobacter, and Campylobacter in the colonic mucosa, Clostridium in the ileal mucosa, Alloprevotella in the cecal mucosa, Bacteroides in the colon, and Treponema in fecal samples. These differences in the relative abundances of certain genera along the gastrointestinal tract likely reflect the physiological conditions at each section. For example, the pH is significantly lower in the stomach than at any other GI location and Lactobacillus spp. are relatively acid tolerant, possibly explaining their enrichment in this location. The lactobacilli are also able to attach to the epithelial and mucosal layers, where they can form biofilm-like communities, thereby aiding their persistence. Lactobacillus spp. were present among all types of GI samples; therefore, the population of lactobacilli in the stomach may serve as a source of these bacteria for other gut locations.

12Comparisons of alpha-diversity measures also revealed significant differences among samples from different GI locations, including feces (Table 3). In particular, there was a clear demarcation between upper and lower GI samples, with significantly higher diversity, richness, and evenness in the lower GI samples (P < 0.05). Within the lower GI tract, fecal and cecal mucosa samples were the most rich in terms of the number of OTUs and phylogenetic diversity (PD). These observations are similar to those reported in the longitudinal swine gut microbiota study conducted by Looft et al.. The only significant difference among upper GI tract sections was observed between the samples of ileal digesta and mucosa, where the mucosal samples had greater phylogenetic diversity as determined using PD whole tree. These results are expected, given that transit is faster in the upper GI tract than in the lower GI tract, thereby limiting adhesion and colonization of bacteria in the upper GI tract. Also, competition for nutrients between the host and bacteria is greater in the upper GI tract.

13Although it is reassuring that GI sampling location is a strong predictor of the swine gut microbial community structure, study-level effects are also significant, as evidenced above. The importance of technical variation among studies has been documented in recent meta-analyses of the avian gut microbiota and the indoor microbiota, where study was the most important driver of the microbiota. In the present meta-analysis, a wide range of ages were included for the pigs, as well as different diets, pig breeds, and administered treatments (Table 1). In addition, at least seven different DNA extraction methods were used and six separate 16S rRNA gene hypervariable regions were sequenced on three different sequencing platforms. The use of standard or uniform DNA extraction, 16S rRNA gene library preparation, and sequencing protocols would likely help reduce some of this study-to-study variation. With that said, eliminating all sources of technical variation in a single study may prove difficult as even technical replicates from different sequencing centers can differ significantly.

14The choice of which 16S rRNA gene hypervariable region(s) to sequence is an important consideration for any 16S-based microbial community study as even so-called “universal primers” differ greatly in their amplification of certain bacterial groups. To determine the effect of different 16S hypervariable region sequences on the characterization of the swine gut microbiota, we analyzed only fecal samples to eliminate one of the largest sources of variability among samples (i.e., the GI sample type). When unweighted and weighted UniFrac distances and Bray-Curtis dissimilarities were plotted using PCoA, samples clustered by which hypervariable region was sequenced, although the amount of variation explained was smaller for the weighted UniFrac distances and there was more overlap between samples sequenced with different hypervariable regions (Fig. 5). In the fecal samples, the study effect was still stronger than the effect of the hypervariable region sequenced with PERMANOVA (adonis R2 values of 0.31, 0.28, and 0.36 for weighted and unweighted UniFrac distances and Bray-Curtis dissimilarities, respectively). However, this is to be expected given that each study sequenced only one specific 16S rRNA gene hypervariable region. Fecal samples that were sequenced using primers that targeted the V4 region did have more phylogenetic diversity (PD whole tree) and a higher Simpson’s reciprocal value than those that were sequenced using primers that targeted the other 16S hypervariable regions analyzed, although the richness and the Shannon index data were not significantly different (Table 4).

15Comparing the region V1 to V3 and region V4 sequences from fecal samples using LEfSe, the region V1 to V3 sequences contained a higher proportion of the Firmicutes phylum and Blautia genus (LDA score [log10] = >4.0). Sequences classified at the phylum level as Spirochaetae and at the genus level as Bacteroides, Lactobacillus, and Treponema were more abundant in fecal samples that were sequenced using the V4 region. Another notable difference at the genus level among the different hypervariable regions was the nearly complete absence of Bifidobacterium in fecal samples from studies using the region V1 to V3 and region V3 16S primer sets (1 sample of 288; Table S2). In contrast, 74.3% of fecal samples that were sequenced using the V4 hypervariable region and 93% of fecal samples using region V3 to V4 had at least one Bifidobacterium 16S rRNA gene sequence. The generic 27F primer sequence (5′-AGAGTTTGATCCTGGCTCAG-3′) which is often used to amplify regions V1 to V3 has been previously reported to underestimate the relative abundance of bifidobacteria unless modifications are made to the primer sequence. Although the relative abundance of Bifidobacterium in fecal samples remained low (<0.35%) in the present meta-analysis, this may be an important consideration for researchers interested in this particular genus.

16Primers targeting the V4 region, however, poorly amplify sequences from Spirochaeta and Desulfitibacter, as these genera were identified in only 12.1% and 0.5% of fecal samples that were sequenced with V4 primers, compared to 91.5% and 77.8% of fecal samples analyzed with region V1 to V3 primers, respectively (Table S2). Although only 14 fecal samples from a single study were available for analysis of region V3 to V4, Treponema spp. were completely absent from these samples, an unexpected finding given the ubiquity of this genus among fecal samples sequenced using other 16S primers; therefore, those results may reflect the primer set used rather than the hypervariable region sequenced.

17Antimicrobials are frequently used in swine production for the purpose of growth promotion, disease prevention, or disease treatment. However, their usage in pigs may help maintain a reservoir of antimicrobial-resistant bacteria and resistance determinants. The development of alternatives to antimicrobials in swine requires a better understanding of how the microbiota is affected by these agents. As 16S rRNA gene sequences were available from six separate studies that investigated the effect of antimicrobial use on the swine fecal microbiota, we combined these studies to determine if any pattern would emerge that was associated with the use of any type of antimicrobial agent in swine. Samples once again clustered by study rather than by antimicrobial use, as nonantimicrobial control samples were not separated from antimicrobial treatment samples within each study (Fig. S3). Also, no genera were significantly different between pigs that received antimicrobials and those that did not (LDA [log10] score, <4.0). This result is perhaps not surprising given that there is a strong study effect and that antimicrobials with different mechanisms of action and dosages were used in each study. Each study also reported findings that are not necessarily in agreement with each other. For example, the Bacteroidetes phylum has been reported to be relatively less abundant at specific time points in pigs fed tylosin and ASP250. However, carbadox treatment has been shown to actually increase the proportion of Bacteroidetes during the early phase of administration.

18These findings emphasize the need for additional research into the effects of antimicrobials on the swine gut microbiota to better predict the response of the microbiota to these agents and to aid in the development of alternatives. The use of larger numbers of animals, longitudinal sampling, and the reporting of pig mass and other production attributes in each trial would all improve our understanding of the effect that antimicrobial agents have on gut microbiota of swine. It is also possible that any beneficial health or production effects observed in the pig are a result of changes in areas of the GI tract that preclude repeated sampling such as the ileum or cecum or that directly impact the host animal rather than the microbiota.

19The results from our meta-analysis demonstrate that after “study” effects, GI tract location is the strongest predictor of the swine gut microbiota composition. In particular, the relatively abundant taxa are similar among samples from the same GI location, independently of other experimental variables. However, technical variation between individual studies may mask other biological differences among swine gut samples. For example, it was not possible to identify any consistent patterns related to antimicrobial use when six studies that included an antimicrobial treatment were analyzed together. The specificity of 16S rRNA gene primers and the hypervariable region that is sequenced also have a large effect on the analysis of the microbiota, as demonstrated in fecal samples from studies using different 16S regions. Nonetheless, a portion of the gut microbiota is shared among at least 90% of GI samples regardless of experimental variables, and this includes Clostridium, Blautia, Lactobacillus, Prevotella, Ruminococcus, Roseburia, the RC9 gut group, and Subdoligranulum at the genus level. These genera represent bacteria that are well adapted to the swine gut and may serve as potential markers of a typical swine gut microbiota. The inclusion of comprehensive metadata to accompany sequencing data sets would enhance future meta-analyses and allow researchers to compare their results directly with those from similar studies.

20A total of 20 studies were included in the meta-analysis and are described in Table 1. Studies were identified through a literature search of the Short Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra), NCBI PubMed (http://www.ncbi.nlm.nih.gov/pubmed), and Google Scholar (http://scholar.google.com). To be included in the meta-analysis, all studies were required (i) to have been sampled from the swine GI tract, (ii) to have used high-throughput sequencing of the 16S rRNA gene, (iii) to have associated metadata and quality score files, and (iv) to have provided data that were publically available before 31 March 2016. Sequence files from each study were downloaded from the SRA or European Nucleotide Archive (ENA; http://www.ebi.ac.uk/ena).

21All downloaded 16S rRNA gene sequences were processed using the QIIME software package (v. 1.9.1). When paired-end sequence files were provided from Illumina MiSeq or HiSeq 2000-based studies, reads were joined using fastq-join with an allowed maximum difference of 15% and a minimum overlap of 35 bp. Sequences derived from the 454 FLX platform were converted from the fastq file format to fasta and qual files, and sequences were subsequently quality filtered using the split_libraries.py command with a Phred score cutoff of 25, a minimum read length of 200 bp, and a maximum homopolymer run of 6 bp. Illumina MiSeq and HiSeq 2000 reads were quality filtered using the split_libraries_fastq.py command with a Phred score cutoff of 25, reads were truncated following three consecutive base calls with a Phred score of <25, and reads with <75% of the original read length following truncation were discarded. The maximum length of each read was set based on the expected read length, as reads spanning regions V1 to V3, for example, are significantly longer than those corresponding to the V4 region.

22The quality-filtered 16S rRNA gene sequences were clustered into operational taxonomic units (OTUs) at 97% similarity using a closed-reference approach, the SILVA release 119 database, and USEARCH v. 6.1.544, with reverse strand matching enabled. In this method, reads that failed to match one of the 16S rRNA gene sequences in the SILVA database at 97% similarity were discarded. Depending upon the hypervariable region targeted, some 16S rRNA gene primers also amplify archaeal 16S rRNA gene sequences. Therefore, all sequences classified as Archaea, as well as chloroplast, mitochondrial, and eukaryotic sequences, were removed prior to downstream analysis. All samples were randomly subsampled to a level of 1,000 sequences per sample to account for uneven sequencing depth among samples and studies.

23Samples were categorized based on sample type, treatment, diet, country of origin, age, sequencing platform, and 16S rRNA hypervariable region sequenced, where available. Only information about sample type, country of origin, sequencing platform, and hypervariable region sequenced was available for each sample. Linear discriminant analysis effect size (LEfSe) was used to determine which genera were significantly enriched in each sample type. Genera that were relatively more abundant in a particular sample group were identified by LEfSe using the Kruskal-Wallis test (P < 0.05), and the effect size of each of these genera was estimated using linear discriminant analysis. An LDA score (log10) of 4.0 was used as the cutoff for identifying differentially abundant genera.

24The between-sample (beta) diversity was assessed using the unweighted and weighted UniFrac distances and Bray-Curtis dissimilarities. Principal-coordinate analysis (PCoA) was used to visualize these distances using Emperor. Permutational multivariate analysis of variance (PERMANOVA) using the adonis function with 9,999 permutations was implemented in QIIME to analyze the unweighted and weighted UniFrac distances and the Bray-Curtis dissimilarities for each gastrointestinal location, country of origin, hypervariable region, sequencing platform used, and study. Pig age was also assessed but with one study removed due to a lack of metadata. Permutational analysis of multivariate dispersions (PERMDISP) was used to test the homogeneity of dispersion for each metadata category as implemented in QIIME using the betadisper function of vegan. The within-sample (alpha) diversity, richness, and evenness were calculated within QIIME using the Shannon index, phylogenetic diversity (PD whole tree), Simpson reciprocal index, and equitability (evenness) index. These metrics, as well as the metric of hypervariable region sequenced for fecal samples only, were compared among gastrointestinal sample types using a two-way ANOVA in R (v. 3.2.5) with hypervariable region and sample type as the independent factors, followed by Tukey’s honestly significant difference (HSD) post hoc pairwise comparison test (agricolae package). All results were considered significant at P values of <0.05.