How being synanthropic affects the gut bacteriome and mycobiome: comparison of two mouse species with contrasting ecologies

Background The vertebrate gastrointestinal tract is colonised by microbiota that have a major effect on the host’s health, physiology and phenotype. Once introduced into captivity, however, the gut microbial composition of free-living individuals can change dramatically. At present, little is known about gut microbial changes associated with adaptation to a synanthropic lifestyle in commensal species, compared with their non-commensal counterparts. Here, we compare the taxonomic composition and diversity of bacterial and fungal communities across three gut sections in synanthropic house mouse (Mus musculus) and a closely related non-synanthropic mound-building mouse (Mus spicilegus). Results Using Illumina sequencing of bacterial 16S rRNA amplicons, we found higher bacterial diversity in M. spicilegus and detected 11 bacterial operational taxonomic units with significantly different proportions. Notably, abundance of Oscillospira, which is typically higher in lean or outdoor pasturing animals, was more abundant in non-commensal M. spicilegus. ITS2-based barcoding revealed low diversity and high uniformity of gut fungi in both species, with the genus Kazachstania clearly dominant. Conclusions Though differences in gut bacteria observed in the two species can be associated with their close association with humans, changes due to a move from commensalism to captivity would appear to have caused larger shifts in microbiota.

from their wild ancestors (e.g. see [7][8][9][10]). Despite intensive research in this field, several aspects of phenotypic change due to a commensal lifestyle remain unexplored. In particular, little is known about the effect of a commensal lifestyle on the composition and function of microbial communities harboured within commensal bodies [11] and putative phenotype changes that may be induced by variation in microbial populations.
Animal-associated microbial composition has an important influence on both host health and physiology [12][13][14], including digestive capacity [15] and interactions with the host's central nervous [16] and immune systems [17][18][19][20] and, as such, can be considered a component of the animal's phenotype. At the same time, animal-associated microbiota exhibit considerable plasticity due to changing physiological states and environmental conditions, including diet variation, stress or parasite infection [16,[21][22][23]. Consequently, there is a high expectation that remodulation of a host's microbiota due to a commensal lifestyle may parallel widely observed changes in microbial populations after translocation of free-living individuals into captivity [24][25][26][27].
Murine rodents represent a suitable model group for exploring the effect of commensalism on associated microbiota as the switch between a commensal vs. noncommensal lifestyle has taken place repeatedly during the course of their evolution [1,28]. To date, however, most research on murine gut microbiota has focused on captive murine species [12,15,21,[29][30][31][32], commensal populations of murine rodents [33][34][35] or noncommensal murine species that are phylogenetically distant to commensal taxa [11,36]. Consequently, there is little comparative data available aimed directly at the assessment of commensalism on microbial structure. This lack of knowledge is even greater as regards eukaryotic microbiota, with the major component in noncommensal animals, the fungal mycobiome, having only been studied exceptionally [37], and never characterised in wild rodents.
To gain an insight into microbial changes due to commensalism in murine rodents, we studied the gut microbiota of two closely related mouse species, the mound-building mouse (Mus spicilegus; hereafter MS) and the house mouse (Mus musculus; hereafter MM). Together with two other mouse species, MM and MS form a monophyletic group with a genetic distance between them of less than 1% [38]. MS are found in the steppes or in agricultural landscapes and use large mounds built from harvested plants as communal shelters and reproduction sites; importantly, the MS lifecycle is independent of human buildings [39]. Unlike the noncommensal MS, MM benefit from a tight commensal association with humans established ca. 8000 years ago that has enabled its cosmopolitan spread [1].
Consequently, the MM lifecycle is dependent on human infrastructure over most of its current distributional range, including our sample sites.
In our study, we use culture-independent gut microbiota profiling based on high-throughput amplicone sequencing, focusing specifically on the gut bacteriome (hereafter GB), which was characterised using 16S rRNA profiling. Moreover, for the first time, our research describes the gut mycobiome (hereafter GM) structure of a free-living mouse population using ITS2 sequencing. As gut microbial variation in different gut sections from the same individual may exceed interindividual variation of microbial communities sampled within single gut compartments [35], we assessed whether there are any species-specific patterns in gut microbial variation across three gut sections (colon, caecum and ileum).
According to both OTU prevalence-based (i.e. Jaccard) and relative abundance-based (i.e. Bray-Curtis)   (Fig. 3). PERMANOVA analysis indicated that the effect of species and sample location was highly significant for both dissimilarity types, whereas differences between gut sections were only supported by PERMANOVA for relative abundance-based dissimilarity. The effect of sex, despite being marginally significant, explained only < 5% of GB variation ( Table 2). Using negative binomial GLMMs, we detected 11 OTUs (represented by 11% of all reads in our dataset) whose abundances varied between MM and MS. Remarkably, four Oscillospira, two OTUs from the family Clostidiaceae and one Odoribacter OTU were more abundant in non-commensal MS than commensal MM. The two species also differed in relative abundance of two Helicobacter OTUs, one being more abundant in MM and the other in MS (Fig. 4). According to phylogenetic placement analysis, these OTUs represent phylogenetically distant Helicobacter species, where Campylobacter jejuni (GenBank accession: EU127548.1) is considered as a root (Fig. S1). Furthermore, OTUlevel analysis identified 55 OTUs (representing 71% of all reads) whose abundances varied between gut sections (Fig. 5).

Gut mycobiome
Unlike the GB, GM profiles were highly homogenous across all murine samples, with fungi of the genus Kazachstania consistently dominant (representing 97% of all reads and 28 OTUs). According to species-level assignment, most samples comprised K. pintolopesii, with the GM of a single MM individual being dominated by K. heterogenica. Other fungal taxa were represented by a low proportion of reads ( Fig. 6). Alongside the murine samples, we also undertook GM profiling (using the same methodology) of ten faecal microbiota samples from a passerine bird (barn swallow [Hirundo rustica], see [40] for details). The GM composition of these samples was much more diverse than the murine GM and covered a broader range of phylogenetically distant fungal taxa. Moreover, Kazachstania was almost absent in these non-murine samples (Fig. S2). As such, we conclude that the high homogeneity and low diversity of the murine GM profiles did not arise as an artefact of wet laboratory procedures.

Discussion
Changes in gut microbiota due to association with humans have been described for a range of vertebrate species [11, 24-27, 34, 41, 42] living in wild vs. in captivity, i.e. zoological gardens [27,42], domesticated animals [26,41] and laboratory animals [34,43]. Overlaying wild vs. captivity microbiome changes for different species potentially allows for the identification of general patterns of microbiome change due to association with humans. On the other hand, the opportunity to establish altered (co)adaptations between microbiota and their hosts in the context of human-altered living environments and lifestyle of the host is limited in captivity. A change in the microbiota of the first generation of animals bred in captivity, for example, may simply represent an unstable, transient state that does not yet show shifts typical for human associated lifestyles. Here, we contrast the gut microbiome of two closely-related mouse species adapted to ecological niches differing greatly in tightness of association with humans (the commensal house mouse [MM] and the non-commensal mound-building mouse [MS]), to get insight on how can commensal association with humans affect bacterial and fungal communities residing in the gut.
In our study, MS displayed a higher gut bacteriome diversity than MM. This may parallel the commonly observed decrease in microbial diversity after introduction of wild animals to captivity [24,25,42,44], ascribed to  reduced variation either in environmental factors affecting GB richness (such as diet composition) or in environmental bacteria colonising gut. While GB diversity between the two mouse species showed significant differences, GB diversity variation between different gut sections was even more dramatic, with ileal GB being less diverse than that from the colon and caecum. These within-gut diversity changes are consistent with most previous studies on different mammalian taxa, including rodents [32,35,36], and can be explained by spatial gradients in immune function, acidity level, nutrient concentrations and various host-derived secretions within the gut [30,45,46]. In addition to changes in GB diversity, we also observed differences in GB composition between the two species. Subsequent differential-abundance analysis identified 11 OTUs (represented by 11% of high-quality reads) whose proportions varied between the two species, with the greatest difference being the increase in abundance of four Oscillospira OTUs in MS. Interestingly, in parallel to our findings, a previous study on cattle showed that individuals relying on outdoor pasture hosted more abundant Oscillospira populations than indoor-bred individuals [47]. According to a few independent studies on human subjects, an increase in Oscillospira abundance was associated with reduced body weight [48]. Moreover, a comparative study on several vertebrate species exposed to a fasting regime revealed consistent enrichment of their GB by Oscillospira [44]. Together, these observations suggest that variation in the abundance of Oscillospira OTUs between MS and MM may reflect differences in nutritional conditions between the commensal and non-commensal niches. Along with Oscillospira, two OTUs from the family Clostridiaceae and one Odoribacter OTU were more abundant in non-commensal MS than commensal MM. All these taxa (including Oscillospira) are anaerobes, capable of fermenting complex plant polysaccharides and producing bioactive short-chained fatty acids (SCFA) with anti-inflammatory capacity [49,50]. As deficit of SCFA is associated with increased emergence of metabolic disorders and impaired cognitive and immune functions [51], an experimental study linking differences in the abundance of specific SCFA-producing bacteria to Helicobacter, a common inhabitant of murine GB [36,52], can induce pathological states under certain circumstances [29,53]. Some Helicobacter species are transferred across generations by social contact between community members and, consequently, exhibit phylogeographic codivergence with their hosts [54,55]. This does not appear to be the case regarding the differences in abundance of Helicobacter OTUs in MM and MS, however, as they represent phylogenetically distant lineages and group with Helicobacter clades of low host specificity. Particularly, Helicobacter OTUs abundant in MS exhibited relatedness to H. canicola and H. bilis, whereas those abundant in MM clustered with H. apodemus.
In this study, considerable differences in GB composition were observed between gut sections, irrespective of host species identity, with the abundance of 55 OTUs (representing 71% of all reads) varying between gut sections. Variation in GB taxonomic content between gut sections was generally congruent with most previously published data on mammalian species [35,36]. In particular, the distal parts of the gut (i.e. the colon and caecum) were enriched with OTUs of the phylum Bacteroidetes (taxa Bacteroides, Prevotella, S24-7, Rikenellaceae), class Clostridia (genera Ruminococcus, Oscillospira and unassigned Lachnospiraceae) and the genus Helicobacter (from phylum Proteobacteria). On the other hand, the ileal GB was characterised by an increased abundance of Candidatus Arthromitus, Lactobacillus (both of the phylum Firmicutes), Mycoplasma (phylum Tenericutes), Aggregatibacter (phylum Proteobacteria) and Mucispirillum (phylum Defferibacteres). GB composition also varied between sample locations, implying that the individuals sampled were exposed to different environmental bacterial pools, despite being only 10 km apart, which is consistent with previous studies on wild MM microbiota [33,43]. Finally, males and females exhibited only slight variation in GB, accounting for~4% of total variation in GB composition.
Surprisingly, both MM and MS exhibited highly uniform gut mycobiome structure, with the genus Kazachstania (including two species, K. pintolopesii and K. heterogenica) as the dominant component, representing ca. 75-100% of high-quality reads. Kazachstania, a species complex of Ascomycetous yeasts, has been isolated from a number of captive species, including K. heterogenica and K. pintolopesii from rodents, K. slooffiae from pigs and horses and K. bovina from cows and birds [56]. Although little is known about its functions, recent experimental studies have shown the importance of Kazachstania for the development of a healthy porcine microbiome, where Kazachstania populations support growth of SCFA-producing bacteria [57,58]. On the  [59]. Similarly, mixed infection of Kazachstania and Escherichia coli was reported as a causal agent of a fatal disease in captive bred primates [60]. Interestingly, the GM of free-living mice from our populations differed from GM profiles previously reported for captive MM, with captive individuals exhibiting much higher diversity and presence of other prevailing yeast genera, such as Candida or Saccharomyces [23,31].

Conclusions
Massive changes in gut bacteriome composition, explaining~20% of total GB variation, have been documented in several previous studies following introduction of MM from commensal populations to breeding facilities [34]. These changes had a dramatic effect on host physiology, immune function and fitness [43,61], with important implications on the usage of captive MM colonised by breeding facility-specific gastrointestinal microbiota as models for biomedical research. Similarly, the gut mycobiome profiles reported from captive MM appear to show pronounced differences to the commensal population GM described in this study [23,31]. The fact that such massive changes in GB and associated host phenotype occurred between two human-associated lifestyles (commensal and captive) leads to the obvious assumption that transition from a wild-living to a commensal niche will have an even more dramatic effect. Nevertheless, our study comparing gut microbiota differences between commensal MM and non-commensal MS suggests relatively low differences in GB content and diversity between the two species, accounting for < 10% of total GB variation. Moreover, mouse GM was surprisingly uniform and consistently dominated by Kazachstania, both in commensal MM and non-commensal MS. Consequently, we suggest that translocation from a commensal population to captivity has a comparatively higher effect on mouse gut microbiota than gut microbiota changes associated with adaptation to a commensal niche. Further observation research (1) covering more commensal vs. non-commensal species pairs or (2) focusing on temporal dynamics of gut microbiota changes as well as (3) experimental studies aimed at direct microbiota manipulations would help to uncover details behind variation in gut microbiota vs. host physiology associated with commensal life style.

Sample collection
The mice used in this study were live-trapped in east Slovakia between the 7th and 9th of December 2015, with 4 individuals of MM trapped in a farm granary at Drienovec (N 48°36.683′, E 20°56.333′) and 8 individuals of MS obtained from field habitat near Drienovec and Čečejovce (N 48°33.881′, E 21°4.137′). In each case, the traps were controlled twice per day, with trapped individuals taken back to the laboratory and caged separately in clean cages with sterile bedding to avoid microbial cross-infection. All mice were sacrificed by cervical dislocation and dissected within 12-h of trapping. Each of the three gut sections (ileum, caecum and colon) was placed on a sterile Petri dish. They were cut longitudinally. The content was gently washed out with sterile physiological saline solution and the tissue samples were placed separately into sterile cryotubes, rapidly frozen in liquid nitrogen and stored at − 80°C. Whole compartments were taken from the caecum and colon, and the distal part only from the ileum (1 cm). Sampled species are not under legislative protection. Ethical statement regarding sample collection and experimental procedures is provided in "Declarations" section.

Gut bacteriome genotyping
Metagenomic DNA from gut tissue samples was extracted in April 2016 in a laminar flow cabinet using the PowerSoil DNA isolation kit (MO BIO Laboratories Inc., USA). We failed to collect a caecum sample from one MS individual; consequently, the final dataset included 35 samples of metagenomic DNA (Table S2).
Primers flanking the V3-V4 variable region on bacterial 16S rRNA gene (i.e., S-D-Bact-0341-b-S-17 [CCTACGGGNGGCWGCAG] and S-D-Bact-0785-a-A-21 [GACTACHVGGGTATCTAATCC]) were used during the polymerase chain reaction (PCR) step [62], both forward and reverse primers being tagged with 10bp barcodes for sample demultiplexing. For the PCR, we used 8.2 μl of KAPA HIFI Hot Start Ready Mix (Kapa Biosystems, USA), 0.56 μM of each primer and 6.2 μl of DNA template. PCR conditions were as follows: initial denaturation at 98°C for 5 min, followed by 30 cycles each of 98°C (15 s), 55°C (20 s) and 72°C (40 s) and a final extension at 72°C (5 min). The PCR products, together with negative controls (PCR products for blank DNA isolates), were run on 1.5% agarose gel, the concentration of the PCR product being assessed based on gel band intensity using GenoSoft software (VWR International, Belgium). The samples were subsequently pooled at equimolar concentration and run on 1.5% agarose gel, with bands of appropriate size excised from the gel and purified using the High Pure PCR product Purification Kit (Roche, Switzerland), according to the manufacturer's instructions. Sequencing adaptors were ligated using TruSeq nano DNA library preparation kits (Illumina, USA) and the resulting amplicon libraries sequenced on a single MiSeq run (Illumina, USA) using v3 chemistry and 2 × 300 bp paired-end reads at Centre de Biologie pour la Gestion des Populations (CBGP, Montferrier sur Lez Cedex, France). Technical PCR duplicates were sequenced for individual samples.

Gut mycobiome genotyping
The ITS2 region was amplified using universal ITS3 (GCATCGATGAAGAACGCAGC) and ITS4 (TCCTCC GCTTATTGATATGC) primers [63] flanked by oligonucleotides compatible with Nextera adaptors (Illumina, USA). For the PCR reaction, we used 5 μl KAPA HIFI Hot Start Ready Mix (Kapa Biosystems, USA), 0.2 μM of each primer and 4.6 μl of DNA template. PCR conditions were as follows: initial denaturation at 95°C for 3 min, followed by 30 cycles each of 95°C (30 s), 53°C (30 s) and 72°C (30 s) and a final extension at 72°C (5 min). Nextera sequencing adaptors were appended to the resulting PCR products during the second PCR round using 10 μl of KAPA HIFI Hot Start Ready Mix (Kapa Biosystems, USA), 2 μM of each primer and 6 μl of PCR product from the first PCR. PCR conditions were as follows: initial denaturation at 95°C for 3 min, followed by 15 cycles each of 95°C (30 s), 55°C (30 s) and 72°C (30 s) and a final extension at 72°C (5 min). PCR products of the second PCR were run on 1.5% agarose gel and pooled at equimolar concentration. The final library was cleaned up using Agencourt AmpureXP beads (Beckman Coulter Life Sciences). Products of the desired size (250-700 bp) were extracted by PipinPrep (Sage Science Inc., USA) and sequenced on Illumina Miseq (v3 kit, 300 bp paired-end reads) at the Central European Institute of Technology (CEITEC), Masaryk University, Brno (Czech Republic). Technical PCR duplicates were sequenced for individual samples as in the GB analysis. We failed to amplify and sequence ITS2 in two DNA samples.

Bioinformatics analysis
Skewer [64] was used for both sample demultiplexing and detection and trimming of gene-specific primers. In the next step, reads of low quality (expected error rate per paired-end read > 1) were eliminated. Dada2 [65] was used for denoising of quality-filtered reads and quantification of 16S rRNA and ITS2 haplotypes (hereafter operational taxonomic units; OTUs) in each sample. Next, UCHIME [66] was employed for detection of chimeric OTUs. The gold.fna (available at: https:// drive5.com/uchime/gold.fa) and UNITE databases [67] were used as references for chimeras filtering from the GB and GM datasets, respectively. After elimination of chimeric haplotypes, taxonomy for non-chimeric OTUs was assigned at 80% posterior confidence by RDP Classifier [68]. The GreenGenes database version gg.13.8 [69] was used for bacterial OTU annotation, and the UNITE database [67] for fungal OTU annotation. Bacterial OTU sequences were aligned using PyNast [70] and their phylogenetic tree constructed using FastTree [71]. We did not conduct phylogenetic reconstruction for ITS2 OTUs as ITS2 evolution is driven to a large extent by insertions and deletions, which complicates identification of homologous positions for phylogenetically disparate taxa. OTU tables (i.e. OTU read counts in individual samples), OTU sequences, their taxonomic annotations and phylogeny along with sample metadata were merged into separate GB and GM databases using the package phyloseq [72] in R (R Core Team 2015).

Statistical analysis
The GM database comprised 589,406 high-quality sequences grouped in 153 non-chimeric OTUs. The number of ITS2 reads per sample ranged between 1826 and 32,002 (median = 17,896). Only basic descriptive tools were used in GM analysis as the sample structure was highly homogenous (detailed below). The GB database comprised 1,088,282 high quality sequences grouped in 1634 non-chimeric OTUs. As the number of bacterial reads varied between samples (range = 1415 -58,106, median = 29,188), we rarefied the OTU table to achieve even sequencing depth per sample and used the downsampled database for further statistical analysis unless otherwise stated. The Shannon index was used as an alpha diversity measure, with per-sample Shannon diversities included as a response variable in generalised linear mixed effect models (GLMM), while sex, sampling location, gut section and species identity were included as alpha diversity predictors. We also tested for potential effects of species vs. gut section interaction. Individual identity was included a random term in order to account for pseudo-replication as multiple gut sections were analysed for each individual. Nonsignificant predictors were eliminated from the initial model in a step-wise manner in order to obtain the minimal adequate model [73]. Marginal effects for significant predictors (i.e. effect of predictor in questions controlled for the effects of other significant predictors) are reported. Dissimilarity based on OTU prevalence (binary Jaccard index) and relative abundance (Bray-Curtis index) was used to study divergence in microbiota composition between samples. First, we visualised the between-sample divergence pattern using Principal Coordinate Analysis (PCoA), then applied PERMANOVA (adonis2 function from the vegan R package) to test whether predictors already included in the alpha diversity analysis drove variation in microbial profile composition. To account for within-individual covariance, individual identity was included as a constraint for permutations (i.e. strata). We reported marginal probability values, i.e. significance of the variable in question controlled for the effect of all other variables included in the PERMANOVA model. To identify OTUs exhibiting variation in abundance between different gut sections and between the two species, we employed mixed models (from R package BhGLM [74]) assuming negative binomial error distribution. These analyses were only conducted for a subset of dominant OTUs (represented by > 0.1% reads and present in > 5 samples, n = 119). The marginal effect of the predictor in question (i.e. effect of species identity statistically controlled for systematic variation between gut sections and vice versa) was tested using likelihood-ratio tests, with qvalue [75] used as a multiple testing correction method. As in the previous analyses, individual identity was included as a random effect. Tukey post-hoc tests were employed after likelihood ratio testing to identify specific gut section pairs where the OTU in question exhibited significant variation in abundance. In specific cases, phylogenetic placement was used to provide a more detailed insight into OTU taxonomy. A set of reference 16S rRNA sequences exhibiting > 97% similarity with the OTUs in question was extracted from the Silva database v. 132 [76]. Sequences were aligned using mafft [77] and a phylogenetic tree constructed using RAxML [78]. Bootstrap analysis based on 1000 replications was used to estimate support for the tree nodes.