Host Species and Captivity Distinguish the Microbiome Compositions of a Diverse Zoo-Resident Non-Human Primate Population

: Vast numbers of microorganisms inhabit the mammalian gastrointestinal tract in a complex community referred to as the gut microbiome. An individual’s microbiome may be impacted by genetics, diet, and various environmental factors, and has been associated with many health states and diseases, though speciﬁc explanations are lacking. While these communities are well-studied in human populations, non-human primates (NHPs), in particular zoo-resident or captive NHPs, offer distinct advantages to increasing our understanding of factors that inﬂuence gut microbiome composition. Here, we characterize the gut microbiome composition of a phylogenetically diverse cohort of NHPs residing in the same urban zoo. We show that despite overlapping and controlled environmental contexts, gut microbiomes are still distinguished between NHP host species. However, when comparing the zoo cohort to wild NHPs, we show that captivity status strongly distinguishes zoo-resident NHPs from their wild counterparts, regardless of host phylogeny. Microbial orders unique to captive NHPs include taxa commonly present in human gut microbiomes. Together, these results demonstrate that differences between NHP species are strongly associated with gut microbiome composition and diversity, suggesting that species-speciﬁc approaches should be considered when investigating environmental factors’ inﬂuence on gut microbiome composition.


Introduction
Most macroorganisms are home to diverse and complex networks of bacteria, archaea, viruses, and fungi that play vital roles in the day-to-day functioning of their hosts [1,2]. Next generation DNA sequencing (NGS) has shown that different regions of the bodyincluding but not limited to the skin, urogenital system, oral cavity, and gastrointestinal tract-harbor their own compositionally and functionally unique communities, commonly called microbiomes. Of these, the gastrointestinal (gut) microbiome is the most well-studied, due largely to its potential influence on overall host health [1][2][3][4]. The gut microbiome is a vital component of the digestive system and has been linked to immune system function and neurological development through the "gut-brain axis" [5,6]. Importantly, dysbiosis, or the negative perturbation of the gut microbiome, is associated with both acute gastrointestinal distress and chronic diseases, including diabetes, obesity, and some  [7,8]. However, the mechanisms that structure microbial communities in diverse host species, and the ways in which dysbiosis affects host health remain unclear.
Studying the gut microbiomes of non-human primates (NHPs), particularly captive NHPs, provides a unique opportunity to analyze the links between environment, gut microbiome composition, and host factors like species relatedness, genetics, physiology, and health [2,9,10]. Many otherwise healthy captive NHPs experience chronic and unexplained gastrointestinal distress not observed in their wild counterparts [9,[11][12][13][14]. Understanding and mitigating the underlying causes of this distress is crucial to maintaining healthy populations of NHPs, both permanently in institutions like zoos and temporarily during conservation efforts like rehabilitation and relocation. As our closest living relatives, NHPs are also excellent model organisms for the study of the human gut microbiome both in healthy and diseased states. In particular, the comparison of captive NHPs and wild NHPs may assist with determining how industrialization impacts the human gut microbiome [10,15].
In wild NHP populations, host phylogeny is overwhelmingly the strongest driver of microbial community composition [15][16][17][18]. Longitudinal sampling has shown that gut microbiome composition is highly heritable in wild NHPs [19], and phylogenetic analyses of both host and microbial DNA provides strong evidence of coevolution between wild NHP species and their gut microbiomes [20]. Additionally, Clayton and colleagues demonstrated that despite substantial differences in the gut microbiomes of wild red-shanked doucs (Pygathrix nemaeus) and mantled howler monkeys (Alouatta palliata), in captivity the gut microbiomes of both species converged towards a more "human" composition [15]. Notably, "semi-captive" individuals (housed in a sanctuary) possessed gut microbiomes that were compositionally located between fully captive and wild individuals, providing evidence that captivity was driving this change. Our study aims to build on these findings by comparing a phylogenetically broad cross-section of captive, specifically zoo-resident, NHP species.
This project aimed to assess the impact of captivity on the NHP gut microbiome using fecal samples collected from NHPs in a single urban zoo. This cohort was composed of a phylogenetically diverse array of multiple individuals from eight different species [21], but unlike in the wild, all species reside in the same building. Here, we leverage this unique dataset to assess the broader influence of captivity and historical antibiotic use on the diversity and composition of the NHP gut microbiome. We show that even when housed in the same building and sharing many diet components, zoo-resident NHP microbiomes are strongly linked to host phylogeny.

Materials and Methods
Study location, individuals, and sample material. The NHP collection at the Como Zoo in Saint Paul, MN, USA, an American Zoological Association (AZA) accredited zoo, served as the primary study population. Research protocols were reviewed by the University of Minnesota Animal Care and Use Committee and prepared in consultation with the Como Zoo. A subset of samples was described in a previous publication [15]. Fecal samples (n = 129) were collected from mid-July to mid-August, 2009. All fecal samples were collected from freshly voided feces using a sterile spatula, with care to collect from the top of the sample to avoid ground contamination. After collection, samples were immediately frozen, transported on dry ice, and stored at −80 • C until processing.
Samples in this cohort were collected from 33 individuals representing eight primate species (Table S1). The species' phylogenetic relatedness is described in references [21] and [15]. The species sampled include western lowland gorilla (Gorilla gorilla gorilla) (individuals n = 3, samples n = 15), Sumatran orangutan (Pongo abelii) (individuals n = 4, samples n = 22), DeBrazza's monkey (Cercopithecus neglectus) (individuals n = 2, samples n = 11), black-handed spider monkey (Ateles geoffroyi) (individuals n = 5, samples n = 25), blue-eyed black lemur (Eulemur macaco flavifrons) (individuals n = 2, samples n = 10), white-faced saki (Pithecia pithecia) (individuals n = 5, samples n = 16), Geoffroy's tamarin (Saguinus geoffroyi) (individuals n = 2, samples n = 5), and emperor tamarin (Saguinus imperator subgrisescens) (individuals n = 10, samples n = 25). All sampled individuals were housed by species while on exhibit and at night, with the exception of the Geoffroy's tamarins and white-faced saki, who were all housed together during and prior to the period of sample collection ( Figure S1), along with a sloth species (since sloths are not primates, and there is only one, this individual was excluded from this analysis). Animals were not moved between exhibits, and water sources kept separate though derived from the same building water supply; enrichment and physical items were sanitized before entering or moving between animal areas.
Publicly available samples from several species of wild NHPs and from humans living in the United States were incorporated into some analyses. Wild mantled howler monkey (Alouatta palliata) (n = 7) and red shanked douc (Pygathrix nemaeus) (n = 7) samples were originally published by Clayton and colleagues [15]. Wild central chimpanzee (Pan troglodytes troglodytes) (n = 6), and western lowland gorilla (n = 10) samples were originally published by Hicks et al. [22]. United States human (Homo sapiens) (n = 11) samples were randomly chosen from those originally published by the Human Microbiome Project Consortium [3].
DNA extraction and 16S rDNA profiling. Microbial DNA extractions and sequencing were carried out by the University of Minnesota Genomics Center (UMGC). Extractions were carried out according to the Earth Microbiome Project (EMP) protocols [23]; 16S rRNA gene V4 region amplicons were generated according to EMP protocols on all samples and sequenced on an Illumina MiSeq with v3 chemistry and 2 × 300 bp reads. Resulting DNA reads were processed using two common microbiome pipelines in parallel: closed and open reference ( Figure S2).
Closed reference processing consisted of first using SHI7 to stitch paired end reads and filter based on quality scores, removing reads with average quality < 33 and trimming reads to quality score > 30 [24]. One sample was dropped based on quality thresholds. These quality-controlled, preprocessed reads were aligned to the GreenGenes database version 13_8 using the optimal aligner BURST at 97% identity in CAPITALIST mode [25]. Unaligned reads were removed at this stage and were not incorporated into any downstream analysis.
During open reference processing, FASTQ files were pre-processed with SHI7 for adapter removal. The resulting FASTQ files were imported into QIIME2 v2019.4 for further processing [26]. Forward and reverse reads were truncated to the point where mean quality score consistently fell below 25 (forward reads remained at 250 bp; reverse reads were truncated at 233 bp). Then, all sequences were denoised de novo using DADA2 within QIIME2 [27,28]. This pipeline dropped the same sample as the closed reference analysis. All other samples possessed adequate sampling depth and were retained for analysis. The amplicon sequence variants (ASVs) generated by DADA2 processing were then assigned taxonomy using both the SILVA database version 132 [29], and the GreenGenes database version 13_8, through the QIIME2 taxonomic classifier pipeline. All ASVs were retained for downstream analysis.
Visualization and statistical analyses. Unless otherwise stated, data analysis and visualization of all samples from both pipelines was conducted using QIIME2 v2019.4 or R v4.1.2, including the R packages ggplot2, reshape2, and vegan [30]. Significance tests were conducted with an alpha of 0.05, and false discovery rate correction for multiple hypotheses where appropriate. Outputs of the open and closed reference processing pipelines were compared using a Procrustes analysis of unweighted UniFrac distance matrices [31,32]. The statistical significance of the difference between these unweighted UniFrac matrices was further assessed using a Mantel test. Database alignment was preliminarily compared using the top 20 most commonly identified genera in each cohort.
For diversity analyses, sequences were rarefied to the highest possible depth that included all remaining post-QC samples and ranged from 25,000 to 35,000 reads per sample, depending on the subset. Statistical significance in alpha diversity was assessed using a Kruskal-Wallis test. Beta diversity was primarily assessed using unweighted UniFrac, weighted UniFrac, and Bray-Curtis distances and tested for statistical significance using a PERMANOVA test. Finally, differential abundance testing was conducted using ANCOM analyses [26,33].

Pipeline and Database Comparison
To assess the potential influence of computational method on our analysis, we compared two common but distinct community profiling pipelines. 16S rDNA amplicon sequences were processed in parallel by closed reference and open reference pipelines, and two common reference databases. Despite fundamental differences in the ways each pipeline functions, a Procrustes analyses showed that the resulting two community profiles were highly similar ( Figure S3; Mantel test p < 0.001). Unless stated, all following analyses are derived from the open reference (DADA2 in QIIME2) pipeline; this method was chosen because it allows for the retention of all detected ASVs regardless of successful taxonomic classification, and some prior studies have shown it to be more robust to errors [27,28,34,35].
We also compared DADA2 s output with two reference databases: GreenGenes (version 13_8) and SILVA (version 132). Preliminary comparisons showed few substantial differences in taxonomic assignment at the order level. However, SILVA enabled greater genus-level specificity than GreenGenes. Therefore, moving forward we prioritized the use of DADA2 with SILVA as it was more specific, more recently updated, and features more current microbial nomenclature.

Host Species Are Distinguished by Microbiome Composition in Diverse Zoo-Resident Primates
To investigate patterns of microbiome composition within this diverse cohort, we first assessed broad trends in gut microbiome composition. Through beta diversity analyses, we found host phylogeny was overwhelmingly the strongest driver of microbial community composition ( Figure 1). Not only did samples cluster most strongly by host species, but more closely related hosts-e.g., the two tamarin species or phylogenetic groupings of apes, non-ape catarrhines, platyrrhines, and prosimians-had more similar microbial communities than more distantly related hosts. This pattern held even in the case of the white-faced sakis and Geoffroy's tamarins, who were housed in the same exhibit space.
We next compared individual primates' alpha diversity across the cohort, using an averaged value from the multiple samples collected for a given individual. Diversity patterns for each species were similar across metrics measuring evenness, richness, and phylogenetic diversity (Figure 3a-d); in general, the tamarin species and white-faced sakis had the lowest alpha diversity values. Additionally, while our cohort lacked sufficient prosimian and non-ape catarrhines for a statistical comparison (only 2 individuals in each group), platyrrhines as a group had significantly lower alpha diversity than the apes (Figure 3e; p = 0.0002, Faith's phylogenetic diversity; p = 0.0001, observed features; p = 0.004, Shannon index; p = 0.07, Pielou's evenness). While many factors may contribute to these observations, we compared the diets of these primates at the time of sample collection. We found that two ape species received, on average, 44% more types of produce than the platyrrhines (14 vs. 9.7 different produce foods, respectively). While diet may play a contributing role in determining species' microbial diversity, looking at all food sources across the NHP cohort suggests that it is not the major driver of high or low diversity in these individuals ( Figure S4). Bray-Curtis distances generated from DADA2 de novo ASV clustering. Samples were collected from eight primate species, including one lemur (blue-eyed black lemur) two apes (western lowland gorilla and Sumatran orangutan), one non-ape catarrhine (De Brazza's monkey), and four platyrrhines (black-handed spider monkey, white-faced saki, emperor tamarin, and Geoffroy's tamarin). In (a-c), primate species is indicated by color. The same microbiome profiles based on unweighted UniFrac, weighted UniFrac, and Bray-Curtis distances are plotted again in (d-f), respectively, with coloring by phylogenetic group: ape, non-ape catarrhine, platyrrhine, and prosimian. In all plots, sequential samples from each unique individual are connected by lines in the order as indicated by the arrowhead.
In order to determine if any specific taxa were responsible for this pattern, we calculated the average taxonomic composition of each sampled individual's gut microbiome (Figure 2a). Even at the order level, multiple taxa were found at significantly different abundances across host species. An ANCOM analysis revealed that 15 microbial orders, shown in Figure 2b, differed significantly between host species (W statistic ≥ 52): Izimaplasmatales, Verrucomicrobiales, Bifidobacteriales, Paracaedibacterales, Bradymonadales, Methanobacteriales, Selenomonadales, Lactobacillales, Coxiellales, Figure 1. Within a zoo-resident NHP cohort, host species and phylogeny were the strongest driver of gut microbiome composition. Principal coordinate analysis plots based on (a) unweighted UniFrac, (b) weighted UniFrac, and (c) Bray-Curtis distances generated from DADA2 de novo ASV clustering. Samples were collected from eight primate species, including one lemur (blue-eyed black lemur) two apes (western lowland gorilla and Sumatran orangutan), one non-ape catarrhine (De Brazza's monkey), and four platyrrhines (black-handed spider monkey, white-faced saki, emperor tamarin, and Geoffroy's tamarin). In (a-c), primate species is indicated by color. The same microbiome profiles based on unweighted UniFrac, weighted UniFrac, and Bray-Curtis distances are plotted again in (d-f), respectively, with coloring by phylogenetic group: ape, non-ape catarrhine, platyrrhine, and prosimian. In all plots, sequential samples from each unique individual are connected by lines in the order as indicated by the arrowhead. Enterobacteriales, class Mollicutes order RF39, Aeromonadales, Erysipelotrichales, class Kiritimatiellae order WCHB1-41, and Desulfovibrionales.

Figure 2. Fifteen microbial orders differ significantly between host species in the zoo-resident
NHPs. (a) Taxa bar plots were generated using DADA2 de novo ASV clustering and taxonomic alignment to the SILVA database. Relative abundance is square root adjusted for visual clarity. Additionally, microbiome composition is shown at the order level to best visualize differences between species while maintaining clarity. (b) ANCOM analysis revealed 15 microbial orders that differed significantly in relative abundance between primate host species. These discriminating taxa were then plotted in a heatmap to illustrate these differences in abundance.

Figure 2. Fifteen microbial orders differ significantly between host species in the zoo-resident
NHPs. (a) Taxa bar plots were generated using DADA2 de novo ASV clustering and taxonomic alignment to the SILVA database. Relative abundance is square root adjusted for visual clarity. Additionally, microbiome composition is shown at the order level to best visualize differences between species while maintaining clarity. (b) ANCOM analysis revealed 15 microbial orders that differed significantly in relative abundance between primate host species. These discriminating taxa were then plotted in a heatmap to illustrate these differences in abundance.
Shannon index; p = 0.07, Pielou's evenness). While many factors may contribute to these observations, we compared the diets of these primates at the time of sample collection. We found that two ape species received, on average, 44% more types of produce than the platyrrhines (14 vs. 9.7 different produce foods, respectively). While diet may play a contributing role in determining species' microbial diversity, looking at all food sources across the NHP cohort suggests that it is not the major driver of high or low diversity in these individuals ( Figure S4).

Zoo-Resident Primates Have Distinct Microbiomes Compared to Wild Counterparts
To further investigate the influence of captivity on the NHP gut microbiome, we compared the zoo-resident NHP profiles to a set of publicly available samples from four wild NHP species (mantled howler monkey, red-shanked douc, central chimpanzee, and western lowland gorilla) and humans living in the United States. As the sequences from wild gorilla and chimpanzee samples were only available without their associated quality scores, it was not possible to process them through the open reference pipeline. Therefore, closed reference processing was used for this analysis. With the addition of these samples, host species was no longer the strongest association with gut microbial community composition ( Figure 4). Instead, wild/zoo status overwhelmingly drove gut microbiome composition, even when comparing wild and zoo-resident individuals from the same species (western lowland gorilla). scores, it was not possible to process them through the open reference pipeline. Therefore, closed reference processing was used for this analysis. With the addition of these samples, host species was no longer the strongest association with gut microbial community composition ( Figure 4). Instead, wild/zoo status overwhelmingly drove gut microbiome composition, even when comparing wild and zoo-resident individuals from the same species (western lowland gorilla). Color indicates species while shape indicates wild/zoo status (zoo = •, wild = ▲, and human = ◾). All samples grouped together by wild/zoo status regardless of how closely related host species were. This distinction between wild and zoo-resident samples is present even within individuals from the same species (western lowland gorilla). All ellipses have a 95% confidence interval.
To determine the taxa responsible for this pattern, we used an ANCOM analysis comparing relative abundance at the order level between wild and zoo-resident NHPs. This analysis revealed that 13 taxa had significantly different abundances depending on host captivity status ( Figure 5). Taxa found to be significantly more abundant (W statistic ≥ 53) in wild individuals were: Streptophyta, Rickettsiales, class Verruco-5 order WCHB1-41, Figure 4. With the addition of wild NHP and human samples, zoo/wild status defines clustering by gut microbiome composition. Principal coordinate analysis plot based on unweighted UniFrac distances between zoo-resident, wild, and human samples, generated by closed reference analysis. Color indicates species while shape indicates wild/zoo status (zoo = •, wild = , and human = ). All samples grouped together by wild/zoo status regardless of how closely related host species were. This distinction between wild and zoo-resident samples is present even within individuals from the same species (western lowland gorilla). All ellipses have a 95% confidence interval.
To determine the taxa responsible for this pattern, we used an ANCOM analysis comparing relative abundance at the order level between wild and zoo-resident NHPs. This analysis revealed that 13 taxa had significantly different abundances depending on host captivity status ( Figure 5). Taxa found to be significantly more abundant (W statistic ≥ 53) in wild individuals were: Streptophyta, Rickettsiales, class Verruco-5 order WCHB1-41, and class Lentisphaeria order Z20. Conversely, the taxa found to be more common in zoo-resident individuals (W statistic ≥ 52) were the microbial orders Actinomycetales, Bacteroidales, Spirochaetales, Erysipelotrichales, Lactobacillales, Desulfovibrionales, Aeromonadales, Bifidobacteriales, and Methanobacteriales. and class Lentisphaeria order Z20. Conversely, the taxa found to be more common in zooresident individuals (W statistic ≥ 52) were the microbial orders Actinomycetales, Bacteroidales, Spirochaetales, Erysipelotrichales, Lactobacillales, Desulfovibrionales, Aeromonadales, Bifidobacteriales, and Methanobacteriales.

Figure 5. Thirteen microbial orders differ significantly between zoo-resident and wild NHPs.
ANCOM analysis revealed 13 microbial orders that differed significantly (W ≥ 52) in relative abundance between wild and zoo-resident NHP hosts.

Historic Antibiotic Usage Is Linked to Lower Gut Microbiome Diversity
We next assessed the impact of historic antibiotic usage on gut microbiome diversity in the zoo cohort. Despite widespread use of antibiotics in captive NHP populations, their long-term influence is understudied. Of the 33 individuals sampled in the zoo-resident cohort, 70% (n = 24) had received antibiotics at some point prior to sample collection, 48% (n = 16) had been prescribed antibiotics multiple times, and 30% (n = 10) had received at least one course in the previous 9 months (Table S2); however, no individuals were receiving antibiotics during the study sample collection period. In contrast, direct antibiotic exposure in wild primates is nearly nonexistent, with the exception of certain rare conservation-driven interventions [36]. Therefore, we hypothesized that antibiotic usage was contributing to the observed differences in gut microbiome composition between wild and zoo-resident primates.
To eliminate confounding effects of host phylogeny (Figure 1), we reviewed the antibiotic usage patterns within individual zoo-resident species in our cohort. This data showed that the emperor tamarins were the only group where multiple individuals both had (n = 3) and had not (n = 7) been prescribed antibiotics prior to sample collection; they also comprise the largest single-species group of zoo-resident individuals (n = 10). Though small sample size is still a limitation, across four alpha diversity metrics there was a consistent pattern of lower alpha diversity in individuals with a history of antibiotic usage (Figure 6a-d; Faith's phylogenetic diversity, p = 0.12; observed features, p = 0.12; Pielou's Evenness, p = 0.07; Shannon, p = 0.07), consistent with long-term gut microbiome diversity suppression associated with previous antibiotic exposure. Furthermore, while a single observation, the individual with the greatest number of previous antibiotic treatments (four prescriptions, versus one each for the other two individuals) possessed the lowest microbiome diversity of the whole emperor tamarin population by three of the four metrics calculated ( Figure 6). Neither sex (p = 0.18, Mann-Whitney U test) nor age (p = 0.88, Spearman's correlation) of the individuals was significantly correlated with alpha diversity in this cohort.

Historic Antibiotic Usage Is Linked to Lower Gut Microbiome Diversity
We next assessed the impact of historic antibiotic usage on gut microbiome diversity in the zoo cohort. Despite widespread use of antibiotics in captive NHP populations, their long-term influence is understudied. Of the 33 individuals sampled in the zoo-resident cohort, 70% (n = 24) had received antibiotics at some point prior to sample collection, 48% (n = 16) had been prescribed antibiotics multiple times, and 30% (n = 10) had received at least one course in the previous 9 months (Table S2); however, no individuals were receiving antibiotics during the study sample collection period. In contrast, direct antibiotic exposure in wild primates is nearly nonexistent, with the exception of certain rare conservation-driven interventions [36]. Therefore, we hypothesized that antibiotic usage was contributing to the observed differences in gut microbiome composition between wild and zoo-resident primates.
To eliminate confounding effects of host phylogeny (Figure 1), we reviewed the antibiotic usage patterns within individual zoo-resident species in our cohort. This data showed that the emperor tamarins were the only group where multiple individuals both had (n = 3) and had not (n = 7) been prescribed antibiotics prior to sample collection; they also comprise the largest single-species group of zoo-resident individuals (n = 10). Though small sample size is still a limitation, across four alpha diversity metrics there was a consistent pattern of lower alpha diversity in individuals with a history of antibiotic usage (Figure 6a-d; Faith's phylogenetic diversity, p = 0.12; observed features, p = 0.12; Pielou's Evenness, p = 0.07; Shannon, p = 0.07), consistent with long-term gut microbiome diversity suppression associated with previous antibiotic exposure. Furthermore, while a single observation, the individual with the greatest number of previous antibiotic treatments (four prescriptions, versus one each for the other two individuals) possessed the lowest microbiome diversity of the whole emperor tamarin population by three of the four metrics calculated ( Figure 6). Neither sex (p = 0.18, Mann-Whitney U test) nor age (p = 0.88, Spearman's correlation) of the individuals was significantly correlated with alpha diversity in this cohort. Figure 6. Historic antibiotic usage may lead to lower microbiome diversity in zoo-resident emperor tamarins. Diversity analysis of antibiotic usage in emperor tamarins, generated using DADA2 de novo ASV clustering. Only emperor tamarins were chosen for several reasons: to reduce confounding effects of inter-species differences, it is the largest number of individuals in a single species (n = 10), and this was the only species where multiple individuals both had (n = 3, orange) and had not (n = 7, blue) taken antibiotics. Of the three previously treated individuals (orange shapes), two (•) had received a single treatment, while one (▲) had received four. Despite the small samples size, in each of the four metrics-(a) Faith's phylogenetic diversity, (b) observed features, (c) Pielou's evenness index, and (d) Shannon's diversity index-previous antibiotic use trends toward lower alpha diversity.

Discussion
A detailed understanding of the role the gut microbiome plays in host health could lead to the development of a suite of clinical applications for NHPs, humans, and other animals [12,37]. As the gut microbiome can be quickly and noninvasively sampled through the collection of fresh feces, assessing gut microbiome composition could become a staple early diagnostic test of host health. Additionally, targeted alterations in the gut microbiome could provide more patient-specific treatment than current therapeutics allow for some chronic diseases [6,32,[38][39][40].
Determining specific mechanisms by which the microbiome affects host health to the level of detail required for animal health interventions or clinical applications is a current challenge in the field; however promising techniques have emerged. One of these includes building an understanding of how environment and other exposures impact the host's Figure 6. Historic antibiotic usage may lead to lower microbiome diversity in zoo-resident emperor tamarins. Diversity analysis of antibiotic usage in emperor tamarins, generated using DADA2 de novo ASV clustering. Only emperor tamarins were chosen for several reasons: to reduce confounding effects of inter-species differences, it is the largest number of individuals in a single species (n = 10), and this was the only species where multiple individuals both had (n = 3, orange) and had not (n = 7, blue) taken antibiotics. Of the three previously treated individuals (orange shapes), two (•) had received a single treatment, while one ( ) had received four. Despite the small samples size, in each of the four metrics-(a) Faith's phylogenetic diversity, (b) observed features, (c) Pielou's evenness index, and (d) Shannon's diversity index-previous antibiotic use trends toward lower alpha diversity.

Discussion
A detailed understanding of the role the gut microbiome plays in host health could lead to the development of a suite of clinical applications for NHPs, humans, and other animals [12,37]. As the gut microbiome can be quickly and noninvasively sampled through the collection of fresh feces, assessing gut microbiome composition could become a staple early diagnostic test of host health. Additionally, targeted alterations in the gut microbiome could provide more patient-specific treatment than current therapeutics allow for some chronic diseases [6,32,[38][39][40].
Determining specific mechanisms by which the microbiome affects host health to the level of detail required for animal health interventions or clinical applications is a current challenge in the field; however promising techniques have emerged. One of these includes building an understanding of how environment and other exposures impact the host's microbiome composition. Here, we compared the microbiomes of a diverse NHP cohort sharing many environmental factors, but showed that despite this, composition was conserved by host species. A powerful relationship between host species and gut microbiome composition is consistent with prior studies of wild NHPs [15,16,18]. While the gut microbiome is not directly heritable like ones' genome, there is mounting evidence that the composition of the gut microbiome has co-evolved with its host species, resulting in the clear phylogenetic patterns seen in this and other studies [18,20,32].
For a variety of NHP species, the simplified diets they consume in captivity have been linked to a reduction in alpha diversity and the alteration of community composition in the gut microbiome [15]. Frankel and colleagues examined the gut microbiomes of both folivorous and generalist NHP species, and found that across both groups captive individuals experienced lowered alpha diversity and altered gut microbial composition [16]. Notably, this pattern was far more pronounced in folivorous species, whose diets were more substantially altered in captivity than generalist species'. As folivorous species are unable to self-produce the necessary enzymes for plant digestion, they rely disproportionately on their gut microbiomes to fulfill nutritional requirements, leaving them more vulnerable to dysbiosis. While the diets of the black howler monkey (Alouatta caraya) and the mantled guereza (Colobus guereza), two folivorous species, consist of an average of~70% and~61% foliage, respectively, in the wild, in captivity their diets only consisted of~52% and~31% foliage [16]. This substantial reduction in foliage corresponds to a reduction in complex carbohydrates like fiber and secondary plant metabolites, the metabolism of which has been linked to some microbial species more prevalent in the gut microbiomes of wild folivores. Therefore, we wanted to know, in this unnatural context, whether host species was still a strong enough factor to support distinct microbiomes. We show here that host species consistently had a strong correlation with zoo-resident NHP gut microbiome composition, despite potentially confounding factors including highly similar diets, close proximity, and medication usage. While more work needs to be done on the underlying causes of this pattern, differences in diet and gut physiology among primate species are likely playing a role [9,12,18], such as in the difference between apes and platyrrhine primates. Additionally, quantitative dietary intake data, while challenging to obtain, would enable more robust testing for the influence of specific foods and/or nutrients.
Several microbial orders were particularly associated with specific primate species. For example, both of the tamarin species sampled in this cohort (emperor tamarins and Geoffroy's tamarins) had significantly higher abundances of Bifidobacteriales, a bacterial order typically associated with positive health effects in the human microbiome, and significantly lower abundances of Desulfovibrionales, a microbial order linked to diseases such as IBS and periodontitis in humans [41,42]. Lactobacillales, a microbial order commonly found in human probiotics, was significantly more abundant in white-faced sakis than any other species. Several rare taxa were also only present in one or two host species, including Paracaedibacterales (in blue-eyed black lemurs) and Coxiellales (in orangutans and western lowland gorillas).
While the pattern of discriminating taxa revealed in this analysis is intriguing, more work at the genus and species level, ideally coupled with metagenomic sequencing, is needed to fully characterize and compare functional differences that these taxa represent. Additionally, more work is needed to study the roles that these taxa play in zoo-resident and wild NHP gut microbiomes. While several of these taxa, including Bifidobacteriales, Methanobacteriales, and Lactobacillales are actively studied in humans and mice, very little work has investigated their roles in the gut microbiomes of non-model mammalian species. While one might assume that their influence on host health will be similar, the physiological differences between human and NHP digestive tracts are substantial enough that significant differences in microbe-host associations could occur.
Captivity status has a stronger influence on the gut microbiome than host phylogeny. Despite the powerful influence host phylogeny appears to have on the gut micro-biomes of zoo-resident NHPs, it was supplanted by captivity status when the zoo cohort was compared to wild NHPs and humans (Figure 4). Our ability to find wild counterparts for each of the zoo-resident NHP species sampled was limited to publicly available amplicon sequences; therefore, we focused on acquiring samples from a diverse group of NHPs from several different research studies [3,15,22]. Despite the extensive primate phylogenetic diversity of both the zoo cohort and the wild NHP samples, and the fact that both groups included samples from western lowland gorillas, all samples clustered strongly by captivity status regardless of host phylogeny.
While captivity clearly delineates these groups of samples, our findings still support the influence of host phylogeny on gut microbiome composition in NHPs [15][16][17][18]43]; within the zoo and wild cohorts, samples still cluster by host species (Figure 4). Interestingly, while previous studies suggest that NHP gut microbiomes converge toward a more human composition when going from the wild to captivity, we still observed strong separation between wild, zoo-resident, and human samples. However, our analysis includes comparatively few wild NHP and human samples, and no sanctuary or rehabilitation samples, and our study was not designed to address this question.
By comparing wild and zoo-resident NHP taxonomic profiles, we found 13 discriminating taxa, four of which were more common in wild NHPs and nine of which were more common in zoo-resident NHPs ( Figure 5). The four microbial orders that were significantly more abundant in wild NHPs were Streptophyta, Rickettsiales, the Verruco-5 order WCHB1-41, and the Lentisphaeria order Z20. Two of these taxa (WCHB1-41 and Z20) are not well described, and thus far have only been identified to the order level, making it difficult to speculate about their functional role in the gut microbiome. The increased prevalence of Streptophyta in wild NHP populations is particularly intriguing, because these sequences are actually derived not from gut bacteria but from the chloroplasts of plants consumed by the host [44]. The significantly higher abundances of chloroplasts in the wild NHP samples is an indication that their diets include much more plant matter than the zoo-resident NHPs. This is consistent with previous data and studies that link simplified diets of captive NHPs to changes in the gut microbiome [12,22,32,45].
Interestingly, many of the taxa more commonly found within the zoo-resident population are typical members of the human gut microbiome. Taxa from the order Bacteroidales are the most abundant Gram-negative bacteria of the human gut microbiome [3,4]. A subset of the taxa enriched in zoo-resident NHPs, including Actinomycetales. Desulfovibrionales, Erysipelotrichales, and Spirochaetales, are associated with dysbiosis and disease in the human gut microbiome [42,[46][47][48]. However, Spirochaetales has also been found to be negatively correlated with host obesity in captive cynomolgus monkeys (Macaca fascicularis), making its exact role in the NHP gut microbiome unclear, and suggesting that other human-study-derived health associations may not transfer to NHPs [47]. Several other discriminating taxa, including Lactobacillales, Bifidobacteriales, and Methanobacteriales, are associated with health benefits in humans, although their roles in the NHP gut microbiome again remain unknown [41,49,50]. Further research at higher resolution and over larger, multi-center cohorts is needed to determine functional and health consequences of the compositional patterns revealed here.
Antibiotic usage impacts gut microbiome diversity. While many factors are likely contributing to captivity's influence on the NHP gut microbiome, our dataset enabled us to ask particular questions about diet and historic antibiotic usage. Despite their widespread use in human and captive NHP populations, antibiotic use in NHPs is not well-studied. One study looking at the immediate effect of the antibiotic cephalosporin on the gut microbiomes of wild western lowland gorillas found patterns of lowered stability and a change in the relative abundance of certain key taxa [36]. However, as these wild individuals had never been given antibiotics before, it is unclear how directly applicable these results are to captive NHPs, some of whom have received multiple courses of antibiotics.
We chose to assess the impact of antibiotics within just emperor tamarins, to eliminate the confounding effect of host phylogeny as they were the largest group of one zoo-resident species (n = 10) and the only species where multiple individuals both had (n = 3) and had not (n = 7) been administered antibiotics in the past. We found that across four different metrics, individuals who had taken antibiotics had lower gut microbial diversity than those who had not, and one individual with the most antibiotic treatments (four, versus one each for others) had the lowest microbiome diversity among all of the emperor tamarins in three of the four diversity metrics ( Figure 6). These findings are in agreement with the limited research conducted on antibiotic use in primates thus far, as several prior studies have found a link between antibiotics and lowered gut microbial diversity [36,51,52]. It should be noted that the vast majority of prior work on this subject focused on the short-term effects of antibiotic usage (i.e., days to weeks post treatment), making this study one of the first to explore the long-term effects of antibiotics on the primate gut microbiome. While this analysis is limited by its sample size and the amount and complexity of associated animal data, our findings warrant further study including additional primate species.

Conclusions
Zoo-resident NHPs provide a unique opportunity to investigate host-microbiome relationships. Not only are many environmental factors monitored and recorded, but several could be readily altered in efforts to improve host health. Additionally, captive NHPs have an inherent "control" group, as wild individuals exhibit fewer GI issues and possess distinct gut microbiomes compared to captive individuals. A better understanding of the ways in which the environment, host health, and gut microbiome composition intersect are vital to mitigating chronic GI issues in captive NHPs, as current treatments have been largely unsuccessful [9,12]. These GI issues pose a risk not only to permanently captive individuals, but also to individuals who are transiently held in captivity for conservation efforts such as rehabilitation and relocation. As anthropogenic habitat destruction makes these invasive conservation efforts become more frequent and necessary, keeping affected NHPs (as well as other threatened species) healthy during captivity will be imperative. Should a re-released NHP acquire a chronic GI illness while in captivity, it will not only reduce that specific individual's chances of survival, but could also risk introducing novel pathogens into already vulnerable wild populations. Furthermore, maintaining a primate's "wild" microbiome composition during rehabilitation could be a critical measure of readiness for release back into the wild. A deeper and more robust understanding of the microbiome features, especially metabolic functions, that distinguish captive from wild primates will be transformative in applying this knowledge to rehabilitation and successful wild release practices.
Understanding the role of gut microbiome dysbiosis in chronic health issues seen in captive NHPs can also shed light onto similar phenomena observed in humans living in highly industrialized societies. In humans, gut microbiome dysbiosis has been linked to a suite of illnesses; however, it is much more difficult to untangle the effects of specific environmental factors in humans than it is in captive NHPs. Importantly, there are many key parallels between the environmental conditions experienced by captive NHPs and humans in highly industrialized societies. Both groups typically have diets that are higher in fats and proteins and lower in fiber and complex carbohydrates [16,45,53,54]. In addition, both groups are regularly exposed to antibiotics, and it is not uncommon for individuals to undergo substantial geographic relocations multiple times during their lifetimes, which has been shown to affect gut microbiome composition and function [51,52,54,55]. However, the long-term study of specific environmental factors like diet and environment in humans is challenging. As NHPs are the closest living relatives to humans, determining the effects of diverse exposures and environments on gut microbiome composition and stability can shed light on the causes and consequences of similar patterns observed in humans [2,10].
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/d14090715/s1, Figure S1: Zoo cohort NHP exhibit layout; Figure  S2: Simplified workflow for data analysis in QIIME2, showing both closed and open reference pipelines. Figure S3: The outputs of both closed and open reference processing were statistically similar. Figure S4: Number of different food sources and microbiome alpha diversity patterns in the