Literature search and data retrieval
We performed a literature search on the internet (Google Scholar, Web of Science) using keywords including “gut microbiota”, “animal microbiome”, “microbiome 16S” and similar to retrieve relevant studies. This search yielded 222 articles published between 2014 and 2020. The materials and methods of these articles were analysed to ascertain whether the study met the following criteria: i) all wild and captive samples were processed using identical procedures, ii) compared wild and captive animals were phylogenetically closely related (members of the same species or species complex), iii) captive individuals were born in captivity, or no information was provided about the origin of the captive animals; i.e., wild animals brought into captivity and sampled some time later were excluded, iv) captive animals that underwent a deliberate selection process (e.g. inbred mice or domestic animals) were also excluded for considering them genetically not comparable to the wild counterparts, and v) only datasets with sample sizes over 12 individuals were considered for analysis. Raw data were extracted from the databases and repositories indicated in the articles (accession numbers listed in the “Bioinformatic resources”).
Bioinformatic sequencing data processing
Datafiles from the different studies were i) stored at the University of Copenhagen’s Electronic Research Data Repository (ERDA), ii) assigned a unique study identifier and iii) re-processed in the Danish National Supercomputer for Life Sciences ‘Computerome2’ using a new bioinformatic pipeline we developed for processing data with different characteristics, including sequencing mode, read length and 16S rRNA gene fragment. The entire code can be found in the “Bioinformatic resources''. In short, for each individual dataset, we quality-filtered (mean phred score of q=25) and (if necessary) trimmed and merged the paired-end reads based on the sequence overlap using AdapterRemoval2 16. Primers (if present) were trimmed using Cutadapt 17, and reads were dereplicated with USEARCH Derep 18 using a relative minimum copy number threshold of 0.01% of the total sequencing depth. Reads were then converted into zero-ratio OTUs using the denoising algorithm UNOISE3 19, and USEARCH was used to map the reads back to the OTUs and create an OTU table. HS-Blast 20 was used to assign taxonomy against the non-redundant Silva 132 database 21, and taxonomic assignments were filtered using different identity thresholds for each taxonomic level: 97% for genus-level taxonomy, 95% for family-level taxonomy, 92% for order-level taxonomy and 90% for higher taxonomic levels 22. To minimise the impact of incorrectly assigned taxa, taxonomic annotations below these identity thresholds were converted into unclassified. This pipeline yielded relative read abundances assigned to different taxa for each individual dataset analysed.
Data quality filtering
Individual data files generated by the aforementioned pipeline were aggregated by study and host species into genus-level abundance tables. The two datasets of Sarcophilus harrisii retrieved from two different studies were processed independently. We then proceeded to quality-filter the genus-level abundance tables of each species through filtering individuals by minimum sequencing depth, minimum diversity coverage and taxonomic annotation. Only individual datasets with more than 1000 reads and diversity coverage values over 99% were retained, and final genus-level abundance tables that contained at least five animals in each contrasting group were considered for analysis. Read counts assigned to taxonomic groups not belonging to Bacteria or not present in the LTPs132_SSU release of the SILVA Living Tree (https://www.arb-silva.de/projects/living-tree) used for measuring the phylogenetic relationships among bacteria were removed to ensure accurate measurements of phylogenetic diversities. In the cases where one group (either wild or captive) outnumbered the other, samples were randomly selected to ensure even sample sizes.
Diversity analyses were carried out in the R statistical environment v.3.6.3 23 based on the Hill numbers framework. The operations explained below were conducted using the R packages hilldiv 24, vegan 25, ape 26, phytools 27, meta 28, dmetar 29 and dendextend 30.
The Hill numbers framework encompasses the group of diversity measures that quantify diversity in units of equivalent numbers of equally abundant taxa 31,32 — in our context bacteria genera. Hill numbers provide a general statistical framework that is sufficiently robust and flexible to address a wide range of scientific questions that molecular ecologists regularly try to answer through measurement, estimation, partitioning and comparison of diversities 33.
Diversity has multiple components, which can be analysed separately through decomposition of diversity or jointly through specific diversity metrics. The three main components of diversity are richness, evenness and regularity 34. Richness is the most intuitive component of diversity, as it only counts whether taxa are present or absent in the system. Evenness quantifies how unbiased a system is in terms or relative abundances of the taxa that comprise it, and the relative weight of abundant and rare taxa can be modulated through the parameter q, also known as the order of diversity. Finally, regularity quantifies the degree of co-relationship among the different taxa, which in our study referred to phylogenetic distances. Hill numbers-based metrics can account for one, some or all components of diversity within a unified statistical framework. To obtain a complete vision of the gut microbiome differences between wild and captive animals, we conducted all our diversity and compositional analyses based on two contrasting metrics: the so-called dR, which only accounted for richness (i.e. whether bacteria taxa were present or not, q=0) and the dRER, which considered Richness + Evenness of order of diversity 1 + Regularity (i.e. whether bacteria taxa were present or not, their relative abundance values that proportionately weighs rare and abundant taxa and their phylogenetic relationships - this metric is also known as the phylogenetic Shannon diversity). Further explanations about the use of Hill numbers in molecularly characterised systems can be found elsewhere 11,34,35.
RER metric requires a phylogenetic tree to compute the dissimilarity values among taxa. As our datasets contained different fragments of the 16S rRNA gene, we were unable to generate a phylogenetic tree from our DNA sequence data. Therefore, we relied on the SILVA Living Tree, and used the LTPs132_SSU release as the reference phylogenetic tree for computing homogeneity values.
Alpha and beta diversity metrics
We computed individual-based diversity metrics using the function hilldiv::hill_div, and obtained average alpha diversity metrics per species, as well as wild and captive populations per species. We also computed beta diversity values for all pairwise combinations of the 644 analysed individuals using the function hilldiv::pair_dis, from which a dissimilarity matrix was derived using the hilldiv::beta_dis function.
We used a Kruskal-Wallis (K-W) test as implemented in the function hilldiv::div_test to ascertain whether the mean diversity values varied across analysed host species, and a PERMANOVA (PMV) test using vegan::adonis function based on the pairwise dissimilarity matrix to test whether host species were compositionally distinct 25. In both cases, statistical significance level was set at a FDR-adjusted p-value of 0.05.
Wild-captive diversity differences meta-analysis
Average alpha diversity metrics of wild and captive populations per species were used to conduct a random-effects-model meta-analysis with raw effect sizes using the function meta::metacont. We used the Sidik-Jonkman estimator for the between-study variance and the Knapp-Hartung-Sidik-Jonkman adjustment method. The overall effect was calculated using Hedge's g (SMD) and its 95% confidence interval and p-value. An identical analysis was performed for the entire dataset and two representative subsets of five species, containing only datasets derived from primates and cetartiodactylans.
Higgin’s & Thompson’s I2 test, Tau-squared T2 and Cochran’s Q were used for quantifying the heterogeneity between the included species. Due to the high heterogeneity found in the study
we evaluated whether the between-study heterogeneity was caused by outliers with extreme effect sizes, which could be distorting our overall effect. We defined an outlier if the species’s confidence interval did not overlap with the confidence interval of the pooled effect using dmetar::find.outliers function.The function detected three outliers in dR metric (GOGO, PEMA and TUTR) and seven in dRER (RHBR, PYNE, PEMA, TUTR, MOCH, CENI and AIME). Even when these outliers were excluded from the analysis the I2 heterogeneity value was substantial for dR (I2 from 79.3 to 71.4%) and moderate for dRER (I2 from 86.9% to 54.2%) and significant for both (Cochran’s Q, p-value <0.001). We performed a sensitivity analysis removing the outliers from the meta-analysis, yet the results of the random effects model did not change (R: SMD=0.345, p-value= 0.075; REH: SMD=0.015, p-value= 0.928). We also performed Graphic Display of Heterogeneity (GOSH) plots to explore the patterns of effect sizes and heterogeneity in our data. We used three supervised machine learning (k-means, DBSCAN and the Gaussian Mixture Model) algorithms to detect clusters in the GOSH plot data and determine which studies contribute the most to them automatically using dmetar::gosh.diagnostics function. This analysis detected five datasets which might potentially contribute to the cluster imbalance, yet even removing these datasets from our study the results of the meta-analysis did not change. Funnel plot and Egger’s regression test were used to evaluate publication bias. Egger’s test was not significant (intercept=0.773, 95%=-2.34, CI=-3.89, t=0.486, p-value=0.632), which means that there is no substantial asymmetry in the Funnel plot that could be caused by publication bias. In light of the above, we decided to include all datasets.
Wild-captive compositional differences
We performed two complementary analyses to study compositional differences between wild and captive animals within each species. On the one hand, we quantified compositional differences between wild and captive animals within each species by means of Sørensen‐type overlap-complement derived from the beta diversity values between wild and captive populations using the function hilldiv::beta_dis. On the other hand, we tested for compositional differences using PERMANOVA (vegan::adonis) based on the pairwise dissimilarity matrix of individuals within each host species.
Contrast of compositional differences with pre-defined models
We defined four models that described four different scenarios of gut microbiota variation between wild and captive animals. Each model was characterised as a different combination of percentage of taxa detected only in wild individuals, in both populations and only in captive individuals, respectively. Model A (5%, 95%, 5%) assumed barely no difference between the gut microbiota of wild and captive animals. Model B (75%, 20%, 5%) described a scenario in which the gut microbiota of captive animals is a subset of that of wild counterparts. Model C (33%, 34%, 33%) defined a situation in which captive animals recruit a proportional set of microorganisms that is different from that of wild counterparts, yet maintain a considerable overlap. Finally, Model D (45%, 10%, 45%) described a situation in which the gut microbiota of wild and captive animals are totally different. Subsequently, we measured the percentage of taxa detected only in wild individuals, only in captive individuals and in both populations for each host species, and calculated the Euclidean distances between the observed data and the four models. The model with the lowest distance to the observed data was selected as the most representative for each dataset.
Compositional similarities among populations within and between host species
We split the abundance table of each species into two subtables only containing either wild or captive individuals, and these subtables were converted into incidence data, which yielded one incidence-based microbiota profile representing each origin (wild and captive) per dataset. Pairwise dissimilarities of all dataset*origin combinations were computed using hilldiv::pair_dis. For each of these 50 datasets, we identified the closest match and calculated the percentage of cases in which the closest match was not the wild or captive counterpart of each species.
Correlation and topological differences between wild- and captive-data derived host comparisons
The profiles generated in the previous step were employed to compute pairwise dissimilarities of the microbiota profiles among species, based solely on wild (25 datasets) and captive (25 datasets) individuals. The correlation between the two dissimilarity matrices was computed in terms of Pearson Product-Moment Correlation test using the function stats::cor.test. For each pairwise comparison, we also calculated the evolutionary distance in terms of millions of years of divergence time. Pairwise dissimilarity matrices derived from wild and captive individuals were subsequently used to generate compositional cladograms based on hierarchical clustering (UPGMA method).
Charts and figures
All charts and figures in the manuscript were originally generated either in R or in Mac OS X Numbers based on csv files created with R scripts. These were subsequently modified in Adobe Illustrator to achieve the desired layout without distorting the dimensions of the quantitative elements. Fig. 1d and 1e were created using the function meta::forest, Fig. 2a and 2b using the function hilldiv::dis_nmds and Fig. 2c using the function dendextend::tanglegram.
All raw DNA sequences are available in NCBI SRA. Accession numbers are provided in the SI. Bioinformatic procedures and relevant metadata are included in the SI. The studies included in the meta-analysis are listed in the Supplementary Dataset. Bioinformatic resources, including scripts and data files have been archived in Zenodo under the doi: [the files are currently in Github (https://github.com/ostaizka/Wild_Captive_16S) but will be transferred to Zenodo, frozen and assigned a doi when the manuscript is accepted for publication and the analyses are considered definitive]. Bioinformatic Code A (613 lines) contains the pipeline used for the re-generation of gut microbiota profiles from raw DNA sequence data. Bioinformatic Code B (914 lines) contains the pipeline employed for the diversity analyses.