The study species
The Damaraland mole-rat (Fukomys damarensis) is a social subterranean rodent that live in cooperatively breeding family groups which can breed in captivity in artificial tunnel systems or can be studied in the wild by trapping individuals in their natural burrow systems [33–36]. Despite being a relatively small rodent (adult body mass 90 to 200 g), Damaraland mole-rats can reach ages of more than 10 years in the wild and likely more than 15 years in captivity [37]. They are strictly herbivores and feed on geophytes with their diet often dominated by a single species, the gemsbok cucumber (Acanthosicyos naudinianus) [38]. The gemsbok cucumber has tubers that are high in fibres, but low in protein and starch content and contain the animals’ entire requirement for water [33,39]. To locate the tubers and expand the tunnel system the animals dig with their large frontal teeth and push up the sand to the surface. Energy requirements of digging behaviour is high, and has been measured for captive Damaraland mole-rats to about five times that of resting metabolic rate [40]. To gain sufficient amount of energy from the tubers, the Damaraland mole-rats are believed to efficiently ferment fibres in their guts [38], which in turn suggests that a healthy and stable gut microbiota of wild Damaraland mole-rat is crucial for the host’s health and fitness.
Sample and data collection
The samples used in this study were collected from 53 captive and 59 wild non-breeding Damaraland mole-rats. For each of the two populations, individuals (51 females and 58 males) from multiple social groups (14 wild, 20 captive) were sampled. The wild animals were captured as a part of a long-term population study of the Damaraland mole-rat population at the Kalahari Research Centre (-26.977439, 21.832659), South Africa, and the captive animals were part of the captive population at the laboratory facility at the Kalahari Research Centre within the reserve. The captive population was founded by wild caught individuals captured around the Kalahari Research Centre in 2013 (-26.938854, 21.691686; -26.890933, 22.079785; -27.112075, 22.061217) and all but three of the sampled individuals of the captive population were F1 and F2 generation from wild caught individuals. All individuals within both locations were pit-tagged to allow individual identification.
The individuals from the wild population were housed in a separate laboratory from the captive population during captures until release back to their burrow system in the wild after a maximum of 7 days in the laboratory. Careful measures were taken to avoid transmission of bacteria between the two populations and contamination of samples. The captive Damaraland mole-rats were fed with sweet potato (Ipomoea batatas) while wild groups were provided their natural diet during group captures and while temporarily housed inside the laboratory. In contrast to the natural diet in the wild, sweet potato is richer in starches and protein but poorer in fibres [39]. Captive groups in this study were provided daily with sand from the nearby area to promote digging behaviour and, while living in captivity, remained exposed to the soil microbiome of their natural habitat.
The fecal samples were collected by placing animals inside a sterilised plastic box provided with paper and a small piece of food. The boxes were checked frequently until defecation. Subsequently, the animals were released back to their group members. For wild-caught animals, the fecal samples were the first fecal pellets after capture. Samples were collected and placed into a 1.5 ml sterile tube and then stored in a minus 80°C freezer on site until transported on dry ice to the laboratory at Linnaeus University, Kalmar, Sweden.
Library Preparation and Sequencing
The 16S library preparation, sequencing protocol and bioinformatic pipeline used in this study has previously been described in Bensch et al. [41] where detailed information on the workflow and pipeline can be found. Fecal samples of captive and wild Damaraland mole-rats were randomised on three 96-well plates, and for this analysis we used 56 samples from wild caught animals collected within the time range of the 53 samples from animals in the captive population, between the 6th of September and 9th of November 2019. On each plate we included four negative control samples by excluding the sample added the first step of extraction and one mock community standard (25ml ZymoBIOMICS Microbial Community DNA Standard) and DNA was extracted using the DNeasy PowerSoil Pro Kit (Qiagen). We amplified the DNA using the primers 341F (5’-CCTACGGGNGGCWGCAG-3’) and 805R (5’-GACTACHVGGGTATCTAATCC-3’) targeting the hypervariable V3-V4 region of the 16S rRNA gene and including adapter sequences for Illumina n5/n7 index primers [42,43] using 25 ml reactions. PCR products were purified using AMPure XP magnetic beads and were used as templates for a second PCR adding a unique combination of Illumina n5/n7 index primers to each sample using 50 ml reactions. PCR-products were purified, DNA concentrations were determined using a Qubit fluorometer (Thermo Fisher Scientific) and equimolar amounts of each sample library were pooled together per 96-well plate into pools with final concentration 4 ng/ml. Pools were 300-bp paired end sequenced following standard Illumina sequencing protocols on an Illumina MiSeq platform at the Swedish National Genomics Infrastructure (NGI) at SciLifeLab in Uppsala, Sweden.
Bioinformatics and sequencing filtering
For bioinformatics, we followed Bensch et al. [41] and processed the raw reads from FastQ inputs using the Ampliseq workflow v1.2.0dev (https://nf-co.re/ampliseq/1.2.0, [44]) which uses Cutadapt v.2.8 [45] and the implementation of DADA2 v.1.10.0 [46] in QIIME2 v2019.10.0 [47] to create Amplicon sequencing variants (ASVs) tables. Quality of the sample reads was checked with FastQC v0.11.8 [48] and MultiQC v1.9 [49], and taxonomy was assigned against the SILVA database v.132 [50].
Quality check and filtering of NGS data
All analyses post Ampliseq were conducted in R version 4.1.2 [51], using functions within the packages tidyverse, vegan and phyloseq [52–54]. To increase the number of reads per sample, we combined reads of samples on plates (plate 2 and 3) that had been sequenced twice.
We filtered away 210 ASVs identified as contaminants by the decontam-package v1.8.0, using a threshold of 0.5 and plate number as batch argument [55], identifying a total 8,709,371 reads and 5,042 unique amplicon sequence variants (ASVs) in the 109 fecal samples from 53 captive and 56 wild Damaraland mole-rats. Mean number of sequences per sample was 79,902 (SD = 40,183), and samples from wild animals had significantly larger library sizes than samples from captive animals (mean captive = 65,054 +- 33,006 SD, mean wild = 93,955 +- 41,559 SD, LMM p = < 0.001, Table S1).
Statistical analysis and measures of diversity
To test for differences in beta diversity between wild and captive animals, we performed a principal component analysis (PCA) on centred log ratio (CLR) transformed counts of ASVs using the rda function in the vegan package [53]. We performed a Permutational Multivariate Analyses of Variance (PERMANOVA) with the adonis2 function in vegan [53] on a Euclidean distance matrix of CLR-transformed counts with population and library size as factors with plate number as strata argument to explore the marginal effects explained by population. To test for differences in dispersion in beta diversity between the two populations, we ran multivariate homogeneity of groups dispersion test with betadisper function in vegan [53].
We explored the taxa driving the microbiome community composition in the directions of the two populations by focusing on the ASVs with the 2% highest and lowest loading scores on the first PC axis (N = 101 unique ASVs in each direction) which clearly separated the two populations. Further, we tested for differently abundant ASVs between the two populations with an Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC) with the ANCOMBC package [56], correcting for multiple testing with Benjamini-Houchberg false discovery rate correction [57].
Sample ASV richness was estimated with the breakaway package [58] which provides standard errors and correction for incomplete sampling. This method uses the complete dataset without rarefying, which has been common practice in alpha diversity estimates of microbiome data but can result in the false impression of unequal richness [59]. We used the betta_random function [58] to test the hypothesis that alpha diversity did not differ between the two populations while controlling for variation between plates using plate number as a random factor and taking the uncertainty and error of the diversity estimate into account. To explore the taxa where ASV richness differed between the two populations, we counted the number of unique ASVs within each sample group of the most common phyla. We tested for differences in the number of unique ASVs per phylum between the two populations with Wilcox signed rank tests as the number of ASVs were non-randomly distributed for most phyla, using Bonferroni correction of p-values for multiple testing (padj).
We tested for differences in library sizes and the ratio of relative abundances of Firmicutes/Bacteroidetes (F/B-ratio) between the two populations with linear mixed models using the lme4 package [60], with plate number as random factor and population as fixed factor. Difference in body mass index (fatness) between captive and wild animals was tested with a linear model fitting population as fixed factor. For all statistical tests, we defined p < 0.05 as our threshold of statistical significance.