Processes shaping gut microbiota diversity in allopatric populations of the endemic lizard Podarcis lilfordi from Menorcan islets (Balearic Islands)

Processes shaping gut microbiota diversity in allopatric populations of the endemic lizard... Abstract Compositional variation of the gut microbiota across host allopatric populations can reflect both adaptation and stochasticity since the time of separation. Major factors shaping this variation include the host phylogeographic and demographic history, the microbiota inheritance, environmental inputs and dispersal of bacteria. Here we explored the impact of these factors in driving gut community diversity in seven allopatric populations of the omnivorous lizard Podarcis lilfordi from the Menorcan coastal islets, all descending from an ancestral mainland population. Using 16S rRNA Illumina sequencing, we showed that ‘islet’ and ‘age’ (time since islet separation from mainland) were the only significant variables in microbial community clustering, suggesting a partial islet-restricted diversification following these lizards phylogeography. Despite a significant variation, islets/populations were characterized by a remarkably low bacterial uniqueness (2.4% of total OTUs) and a minor differential enrichment of taxa, indicating a negligible impact of local inputs and important host common constraints. Overall, the extant pattern of similarity/dissimilarity among islets is compatible with partial retention of the ancestral mainland microbial pool, with differences among islets potentially explained by a differential loss of bacteria following population fragmentation and bottlenecks (i.e. ecological drift). While more quantitative data are needed to validate this hypothesis, this study unveils the importance of considering both neutral and niche-driven processes in driving contemporary patterns of gut metacommunity diversity. island, drift, dispersal, inheritance, population, allopatry INTRODUCTION A multitude of factors shape the gut microbiota composition and challenge its stability over ecological and evolutionary times. Most critical among them are the host constraints (e.g. gut morphology, physiology and immune system) (Rawls et al.2006; Ley et al.2008; Benson et al.2010; Moeller et al.2016) and environmental inputs (e.g. diet and parasites) (Muegge et al.2011; David et al.2014; Delsuc et al.2014; Aivelo and Norberg 2016; Kohl et al.2017). Although diet and host genotype have been consistently found to be the most influential determinants, large part of animal intraspecific microbial variation remains unexplained by adaptive variables. Hence, the role of neutral processes (drift and dispersal) is gaining more attention (Costello et al.2012; Lankau, Hong and Mackie 2012; Nemergut et al.2013; Zeng et al.2015; Burns et al.2016; Martinson, Douglas and Jaenike 2017; Shoemaker, Locey and Lennon 2017). Population-level studies of gut microbiotas dynamics might represent the right scale to explore the influence of both neutral and selective forces on metacommunity structure. While substantial interindividual microbial variation might exist, populations do typically carry discrete gut metacommunities (Yatsunenko et al.2012; Song et al.2013), resulting from a restricted, although not exclusive circulation of gut bacteria favored by parental transmission (Funkhouser and Bordenstein 2013), preferential dispersal by shared ecological niche (e.g. common diet), reproductive strategy (including parental care) and social behavior (Nemergut et al.2013; Song et al.2013; Ren et al.2016). How gut metacommunities respond to restrictions in ‘free’ circulation among host members, for instance in the event of population split, remains to date poorly understood/investigated in natural systems (Lankau, Hong and Mackie 2012), owing to the concomitant effect of adaptive variables (ecological traits). Allopatric but ecologically similar populations found in islands might provide excellent natural systems to explore the impact of geographic barriers and host demographic history on gut microbiota diversity. In the event of fragmentation of an original insular population following isolation of some specimens into novel islets, distinct factors can act over time to shape the microbiota compositional structure (for simplicity, hereafter use of microbial/microbiota will refer to bacteria only). These include the conservation of taxa present in the original pool by stable transmission through host generations (i.e. microbiota inheritance) (Funkhouser and Bordenstein 2013; Zeng et al.2015), selective and stochastic changes in taxa relative abundances (the latter known as ecological drift) (Chase and Myers 2011; Lankau, Hong and Mackie 2012), acquisition of novel taxa from local pools (i.e. local exposure) and dispersal of taxa among populations (Nemergut et al.2013). Of these, ecological drift during and following population isolation and local exposure should drive microbiota divergence among populations (Lankau, Hong and Mackie 2012), whereas microbiota inheritance and dispersal of common taxa across islands should maintain/increase microbial similarity among populations (Zeng et al.2015). Disentangling the contribution of each of these processes on extant patterns of similarity/dissimilarity of gut microbiota would provide a better understanding of the compositional stability and dynamics of these microbial communities over spatial and temporal changes. We explored this scenario using allopatric populations of the lizard species Podarcis lilfordi inhabiting coastal islets around Menorca, one of the Balearic Islands (Fig. 1). The species originated during the re-flooding of the Mediterranean at the end of the Messinian salinity crisis, about 5.33 Mya (Krijgsman et al.1999), with a first cladogenic event (2.88–2.66 Ma) leading to the separation between P. lilfordi (currently present in coastal islets of Mallorca and Menorca, and Cabrera) and P. pityusensis (Ibiza, Formentera and adjacent coastal islets) (Brown et al.2008; Terrasa et al.2009). To date, endemic populations of P. lilfordi inhabit 16 of the 19 Menorca islets (Pérez-Mellado 1989; Terrasa et al.2009) (Fig. 1). According to mtDNA data, these Menorca populations were largely panmictic during the Pleistocene (Terrasa et al.2009); the following sea level rise during the past 10 000 years (Holocene) (Vecchi et al.2016) led to coastal fragmentation and progressive isolation of islet populations up to date (Pretus, Marques and Pérez-Mellado 2004; Brown et al.2008; Vecchi et al.2016). Archaeological data indicate that the mainland stock became extinct ∼2000 years ago subsequent to the introduction of competitive species and foreign predators (Alcover, Moyá-Solà and Pons-Moyá 1981; Sanders 1984; Pretus, Marques and Pérez-Mellado 2004), while islets provided valuable refugia. Extant populations have been geographically and reproductively isolated for the past 2000 years, although time since separation most likely ranges between 4000 and 9000 years, according to bathymetry with respect to mainland and Holocene sea level curves (Vecchi et al.2016). While morphologically diverse, suggesting incipient speciation (Salvador 2009), mtDNA data show low genetic differentiation among P. lilfordi populations, in accordance with their recent segregation (Terrasa et al.2009). The presence of additional vertebrates (birds, rabbits and humans) is restricted to few islets (Pérez-Mellado 1989; Salvador 2009). Several studies indicate that a dietary shift occurred with variable intensity in the different populations, from an ancestral insectivorous niche to a current omnivorous niche, including the uptake of a large proportion of arthropods and vegetable matter, with seasonal variations (Saez and Traveset 1995; Traveset and Riera 2005; Perez-Cembranos, Leon and Perez-Mellado 2016). Figure 1. View largeDownload slide Map of Menorca and sampled islets (circled). Population size according to Pérez-Mellado et al. (2008) is provided below each islet (no estimates available for Ravl, although expected population size is likely below 1000 individuals, personal observations). Bars below each islet correspond to a 50-mt scale. Islet abbreviations: AdGr (Addaia Gran), Aire (Aire), Colm (Colom), PrCv (Porros de Cavalleria), Ravl (Ravells), Rei (Rei) and Tosq (Tosqueta). Figure 1. View largeDownload slide Map of Menorca and sampled islets (circled). Population size according to Pérez-Mellado et al. (2008) is provided below each islet (no estimates available for Ravl, although expected population size is likely below 1000 individuals, personal observations). Bars below each islet correspond to a 50-mt scale. Islet abbreviations: AdGr (Addaia Gran), Aire (Aire), Colm (Colom), PrCv (Porros de Cavalleria), Ravl (Ravells), Rei (Rei) and Tosq (Tosqueta). Overall, two aspects make these populations a rare system for studying natural processes shaping the vertebrate gut microbiota: their genetic and trophic similarity, which allowed us to detract the differential effect of these two major factors affecting the microbiota, and their distribution into several allopatric populations descending from a single ancestral population. Here we aimed to address some key questions: Was the geographic isolation and subsequent exposure to distinct local pools of bacteria sufficient to drive significant divergence among enteric communities in these populations? Do their gut microbial communities still carry an identifiable signature of their past history, i.e. a core component that has been putatively inherited from the ancestral mainland population? What is the actual impact of ecological drift and dispersal in explaining the current microbiota variation? To explore these questions, we sampled whole intestines from P. lilfordi specimens from seven Menorca islets and performed comparative analyses of their gut microbial community structure. We specifically partitioned their gut microbiota into shared and variable components and investigated forces responsible for shaping contemporary patterns of similarity/dissimilarity. MATERIALS AND METHODS Sample collection and morphometric body measures A total of 72 individuals from 7 of the 16 populated islets surrounding the major Menorca island (3 to 18 individuals per population/islet) were sampled in summer 2015, between 27 July and 1 August (Fig. 1 and Table S1, Supporting Information). Summer represents the season when trophic resources are most limited for these lizard populations. Sampled islets were chosen to vary along a gradient of bathymetry (from deep to shallow channels between islets and Menorca) and distance to the coastline (together providing a proxy for time since population isolation), and for their different exposure to the open sea (including both islets found in open sea and within protected bay areas). Specimens were collected using food traps containing fruits or fruit juices. For each specimen, the following morphometric measurements were taken: number of ventral-row scales (VS), weight (W), snout-vent length (SVL), head width (HW) and head length (HL). Of the 72 specimens, 33 adult males were chosen for microbiota sequencing (five for population, except for Rei, for which three specimens were sampled), while the remaining specimens were immediately released (Table S1, Supporting Information). Adulthood stage was assigned when SVL ≥ 51 mm for females and SVL ≥ 55 mm for males, according to the mean values for the smallest subspecies as described in Salvador (2009). Sex was assigned primarily according to number of VS (<27 for males) and validated by body size and femoral pores (males are larger than females and show pores with visible lipophilic compounds (Salvador 2009) (Table S1, Supporting Information). Statistical correlations were assessed through the cor.test function in the R Stats package version 3.2.2. Specimen sampling and manipulation were carried out in accordance with the ethics guidelines and recommendations of the Species Protection Service (Department of Agriculture, Environment and Territory, Government of the Balearic Islands), under permit No. CEP 31/2015 given to LB and JLP. Microbiota DNA extractions and 16S rRNA Illumina sequencing Sampled specimens for the microbiota study were kept in individual plastic containers for up to 6 h, and then transferred to –20°C until death. The full digestive tracts were dissected with sterile instruments and preserved in RNAlater. Samples were immediately stored at 4°C for 24 h and then kept at –20°C until further processing. Sex was then carefully reassessed, revealing the inclusion of three adult females in the final sample set (Table S1, Supporting Information). Under a sterile hood, the entire digestive tracts were wiped from RNAlater using a Kimwipes filter paper (Sigma-Aldrich) and lengths measured with a Vernier Caliper (Table S1, Supporting Information). Stomachs and intestines were separated according to their different tissue characteristics (stomachs had visibly ticker walls). Stomach contents were removed for examination under a stereoscope. The intact intestines were transferred to tubes with a 800-μL lysis buffer (see Baldo et al.2015), chopped using sterile scissors for facilitating the lysis process, vigorously shacked and incubated at 70°C for 30 min with gentle agitation. After addition of 0.3 g of 0.1-mm zirconia-silica beads, tube contents were homogenized with a Precellys Evolution bead beating instrument (Bertin Technologies) at 5500 rpm for 2 × 45 s, and centrifuged at 4°C for 5 min at 16 000 × g. Supernatant was stored in a separate tube. Additional 300 μL of lysis buffer was added to the pellet, and samples were homogenized again in the Precellys instrument for 2 × 30 s. Stored supernatant was placed back into the tube with beads, 20 μL of proteinase K (to a final of 10 mg/mL) was added and samples incubated at 55°C for 1 h at 450 rpm. Extracted DNA was purified according to Baldo et al. (2015). Briefly, residues were eliminated with ammonium acetate and DNA precipitated with isopropanol and cleaned with standard ethanol washes followed by Qiagen columns (DNeasy Blood and Tissue kit). The region V3-V4 (∼460 bp) of 16S rRNA was amplified using the following non-barcoded primers including TruthSeq adapters: S-D-Bact-0341-b-S-17 (5΄-CTACACGACGCTCTTCCGATCT CCTACGGGNGGCWGCAG-3΄) and S-D-Bact-0785-a-A-21 (5΄-CAGACGTGT GCTCTTCCGATCT GACTACHVGGGTATCTAATCC -3΄). Amplicons from individual specimens were barcoded and pooled according to Baldo et al. (2017). The final library was sequenced on Illumina MiSeq v3 (600-cycle cartridge, 300 paired-end reads) at the Centre for Genomic Regulation in Barcelona (Spain). Sequence analyses After demultiplexing, paired-end reads were merged and quality-filtered using the scripts make.contigs and screen.seqs (maxambig = 0, maxhomop = 8, maxlength = 490) implemented in mothur v.1.25.0 (Schloss et al.2009). Chimeras were identified using uchime_ref in vsearch v1.4.4 (VSEARCH GitHub repository) against the Chimera Slayer reference database (‘gold’) in the Broad Microbiome Utilities and discarded. Sequences were input into the QIIME pipeline (Caporaso et al.2010), using Macqiime 1.9.1., and clustered into operational taxonomic units (OTUs) at 97% similarity through the command pick_open_reference_otus.py. Taxonomy was assigned with the UCLUST method against the Greengenes gg_13_8_otus with RDP classifier (80% confidence). Mitochondrial and chloroplast-derived OTUs were discarded, retaining only bacterial OTUs with ≥10 total sequences. Sequences were aligned with PyNAST default parameters (Caporaso et al.2010) using the reference database ‘core_set_aligned.fasta’ from Greengenes website and default lanemask. Alignments were further trimmed to ensure overlapping ends, and a phylogenetic tree was built with FastTree (Price, Dehal and Arkin 2009). Alpha diversity measures (phylogenetic diversity (PD) whole tree and Shannon Index) were estimated on multiple rarefied biom tables to an even depth of 38 160 reads (10 iterations). Statistical differences in alpha diversity across islets were calculated with a non-parametric two-sample t-test (999 permutations). Beta diversity and PERMANOVA analyses Unweighted Unifrac and binary Jaccard distance matrices were built on the rarefied OTU table (38 160 sequences). Statistical correlations between microbiota distances and categorical variables (‘islet’, ‘age’, presence of ‘rats’ and stable colonies of ‘seagulls’, and ‘exposure’ to open sea, see Table S1, Supporting Information) were assessed by permutational multivariate analysis of variance (PERMANOVA, Anderson 2001; Anderson, Gorley and Clarke 2008) for PRIMER v7 (Clarke and Gorley 2015). All factors were considered fixed save for ‘islet’, which was defined as a random factor and included in all designs as nested within the levels of fixed factors. Categorical ‘age’ provides an approximation of each islet isolation time and was largely assigned according to bathymetry (as depth of the shallowest channel between islets and main island) choosing a threshold of ≥9 m to discriminate between old and more recently separated islets, in accordance with the speed of sea level rising, which slowed down after 7500 yr BP (Vecchi et al.2016) (Table S1, Supporting Information). Bathymetry and distance from mainland were also combined to define different age categories, grouping Aire and PrCv (old) against all remaining islets (recent) (see the Results section). To summarize microbial metacommunity diversity by islet, OTU counts, with and without previous rarefaction, were binned by islet and the data were used to calculate ‘uniqueness’ by islet (i.e. taxa and OTUs occurring only in one islet to the exclusion of all other islets). To test whether microbiota distances correlated with the geographic distances among islets, the original OTU table was filtered for OTUs above 50 total counts, reads were binned by islet and rarefied to the minimum depth per islet. Jaccard and unweighted Unifrac distance matrices were built on this OTU table and tested against geographic (i.e. linear) distances across islets with a Mantel test. Islet distances from mainland were instead approximated by the ‘age’ category and tested with PERMANOVA. Taxa and OTUs enrichment among islets To statistically compare differences in taxa and OTU relative abundances among islets, we used the script group_significance.py in QIIME with the non-parametric Kruskal–Wallis test. In all cases, a false discovery rate (FDR) correction was applied to assess significant differences. Differentially enriched OTUs (i.e. indicator OTUs) were estimated on the original OTU table rarefied to 38 160 reads. For differentially enriched taxa (i.e. indicator taxa), OTU counts were binned to each taxon level. Results were filtered for OTU/taxa with total sum of means across islets greater than 10. Core analyses We used the non-rarefied OTU table to estimate core OTUs and taxa (from phylum to species, estimated by binning OTU counts at each taxon level, as above) at different depths (shared by 100%, 90% and 80% of total specimens), retaining those taxa and OTUs with representation in all islets. Core OTUs were searched with batch BLASTN against the local nt database (cut-off e-value = 30), retrieving the first two best hits. Outputs were parsed using in-house perl scripts. Metadata associated with each BLASTN hit was extracted from Genbank records with Biopython scripts. The 134 core OTUs present in at least 90% of the specimens were used to explore a potential signature of within-population diversification. For this purpose, all reads associated with core OTUs were pooled and reclustered at 98% and 99% of similarity (now referred as to 98% and 99% OTUs or more generally as to sub-OTUs), picking the most abundant as representative and retaining OTUs with at least 50 reads total and presence in at least two specimens. A binary Jaccard distance matrix was built on each sub-OTU abundance table, rarefied to the minimum depth coverage across samples (13 600 and 11 000 for 98% and 99% OTUs, respectively), and significance of clustering was tested with PERMANOVA, as for 97% OTUs (see above). The 23 core OTUs that occurred in all specimens (representing the 100% core) were analyzed individually for their sub-OTU contribution to the islet-based structuring. Specifically, reads from each OTU were reclustered at 99% similarity and significance tested with permutational analysis of variance using adonis (999 permutations) with the script compare_categories.py in QIIME. RESULTS Morphometrics and stomach content The large majority of the captured specimens were adult males, putatively because they were less shy in approaching food sources, although a skewed sex ratio cannot be ruled out (Table S1, Supporting Information). Males were on average significantly larger than females (weight, W was 8.2 g ± 2.17 for males, 5.4 g ± 1.55 for females) and had longer SVL (65.2 mm ± 4.20 for males, 58.6 mm ± 3.68 for females) and larger HW/SVL and HL/SVL values (Table S1, Supporting Information). Considering males only, Addaia Gran (AdGr) hosted the smallest lizards, which showed among the largest HL and HW/SVL ratios, while Aire hosted the biggest and longest specimens. These results are compatible with previous findings (Salvador 2009). Overall, SVL positively correlated with W (Pearson's r = 0.677, P < 0.0001), HW (r = 0.6735581, P < 0.0001) and HL (r = 0.863, P < 0.0001), but not with intestine length (r = 0.132, P = 0.333). Stomach content analyses indicated widespread and abundant presence of arthropods, particularly insects (especially ants) in most specimens and all islets, compatible with a dominant role of insectivory during the dry season. Sequence summary We characterized the gut microbiota of 33 total specimens (30 males and 3 females), with three to five representatives per population (Table S1, Supporting Information). After quality filtering, removal of the three OTUs found in the blank control (water only) and OTUs with less than 10 reads total, we obtained a total of 2 014 276 high-quality reads and an average of 60 081 reads per specimen (ranging between 38 163 and 106 206), corresponding to 6380 OTUs. These were assigned to 17 phyla, 32 classes (6378 annotated OTUs), 58 orders (6369), 99 families (4975) and 147 genera (2321). Rarefaction curves based on the number of observed OTUs indicated that we captured most of the sample diversity, although we did not reach saturation at the chosen rarefied depth of 38 160 sequences (Fig. S1, Supporting Information). The taxonomic composition of the P. lilfordi gut microbiota is highly conserved Podarcis lilfordi showed the typical vertebrate phyla composition, with three phyla representing altogether 97% of the microbiota: Firmicutes (average 65% across islets), Bacteroidetes (18%) and Proteobacteria (14%) (Fig. 2). The proportion of Firmicutes was highly comparable across islets, while the representation of both Bacteroidetes and Proteobacteria varied to a greater extent. Firmicutes and Bacteroidetes were represented by six and five orders, respectively, while being dominated by a single order in all islets: Clostridiales (Firmicutes) and Bacteroidales (Bacteroidetes) (Fig. 2). Two families within both phyla showed the highest relative abundances: Lachnospiraceae and Ruminococcaceae within Firmicutes, and Porphyromonadaceae and Bacteroidaceae within Bacteroidetes (Fig. 2). Patterns of family relative abundance within islet were remarkably conserved, especially for Firmicutes. On the other hand, Proteobacteria were more diverse, including a larger number of represented orders (24) and families (33). Among the most represented taxa within this phylum were four orders (Desulfobrionales, Enterobacteriales, Campylobacteriales and Pseudomonadales) and three families (Helicobacteraceae, Desulfovibrionaceae and Enterobacteriaceae), present in all islets (Fig. 2). Proteobacteria varied widely in their quantitative profiles across islets at all taxonomic levels. Nonetheless, alpha diversity estimates for this phylum were not significantly different among islets (two-sample t-test P > 0.05, both Shannon and PD whole tree for all islet pairwise). Figure 2. View largeDownload slide Order and family-level taxonomic composition of the three major phyla, Firmicutes, Bacteroidetes and Proteobacteria, as proportions of OTUs. Values were estimated on the rarefied dataset and averaged by islet. Legends list only the 10 most abundant taxa. Figure 2. View largeDownload slide Order and family-level taxonomic composition of the three major phyla, Firmicutes, Bacteroidetes and Proteobacteria, as proportions of OTUs. Values were estimated on the rarefied dataset and averaged by islet. Legends list only the 10 most abundant taxa. As a whole, gut communities did not show significantly distinct taxonomic diversity and richness across islets, according to both Shannon and PD-whole tree (Fig. S2, Supplementary; P > 0.05 at both indexes for any islet pairwise), although a visual trend indicates Porros de Cavalleria (PrCv) as the islet hosting the largest diversity. Diversity and richness did not correlate with intestinal length (Pearson's r = –0.092, P = 0.49 for PD whole tree and –0.194, P = 0.15 for Shannon). Overall, the major taxonomic structure of the P. lilfordi gut microbiota was remarkably conserved across islets, with Proteobacteria content representing the most variable quantitative trait. Islet as the major determinant in the P. lilfordi microbiota clustering According to both unweighted Unifrac and Jaccard distances on 97% OTUs, ‘islet’ was the main significant factor accounting for 13.5% of the community variation, as measured by partition into components of variation (PERMANOVA, P = 0.0001, Table 1). There were no significant effects of age (but see below), exposure or presence of rats (P > 0.05). The presence of seagull colonies was barely significant, accounting for only 3.5% of microbiota variation. The islet Rei was grouped with the oldest islets according to bathymetry, but it presents several unique features, including its location within a natural harbor at short distance from the mainland and in a context of a largely anthropized environment (the islet hosted a paleochristian basilica in the 6th century and a hospital in the 19th century). When both bathymetry and distance from the mainland were used as grouping factors, thus excluding Rei as an ‘ancient’ islet and retaining only Aire and PrCv, ‘age’ became a significant factor (P = 0.0001, Table 1). The only interaction that could be tested (age by exposure) was not significant. Microbiota distances across islets did not correlate with geographic distances across islets (Mantel test, r = 0.00588, P = 0.980 for Jaccard and r = –0.02849, P = 0.906 for Unifrac), but perhaps indirectly with their distances from the mainland, here approximated by the ‘age’ category (see the Methods section). The ‘islet’ and ‘age’ effects are also illustrated by a principal coordinate analysis based on unweighted Unifrac distances, depicting the specimens clustering according to the islet origin (color-coded, Fig. 3). The two oldest islets, Aire and PrCv, clearly diverged from the remaining islets along PC2. Figure 3. View largeDownload slide Principal coordinate analysis according to unweighted Unifrac distances. Distances were calculated for OTUs at 97% similarity using the rarefied dataset (38 160 sequences). Dots correspond to single specimens. Specimens within each islet are connected to the median coordinates for each islet. Figure 3. View largeDownload slide Principal coordinate analysis according to unweighted Unifrac distances. Distances were calculated for OTUs at 97% similarity using the rarefied dataset (38 160 sequences). Dots correspond to single specimens. Specimens within each islet are connected to the median coordinates for each islet. Table 1. Results of PERMANOVA analyses, separately for the whole dataset and core OTUs after reclustering at 98 and 99% similarity.     Whole dataset  Core (90% of specimens)      Unweighed Unifrac  Jaccard 97%  Jaccard sub-OTU 98%  Jaccard sub-OTU 99%  Model  Term  P-value  Variation  P-value  Variation  P-value  Variation  P-value  Variation        estimate    estimate    estimate    estimate  Islet  Islet  0.0001  0.028  0.0001  0.035  0.0001  0.029  0.0001  0.031  Age + islet(A)  Age  N.S. [0.0001]  0.004 [0.011]  N.S. [0.0001]  0.006 [0.012]  N.S. [0.0001]  0.005 [0.012]  N.S.  0.003 [0.008]    Islet(A)  0.0001  0.024 [0.022]  0.0001  0.029 [0.028]  0.0001  0.024 [0.023]  0.0001  0.028 [0.027]  Exposure + islet(E)  Exposure  N.S.  0.000  N.S.  0.000  N.S.  0.000  N.S.  0.000    Islet(E)  0.0001  0.028  0.0001  0.034  0.0001  0.029  0.0001  0.031  Age + exposure + islet(A, E)  Age  N.S.  0.002 [0.010]  N.S.  0.005 [0.012]  N.S. [0.0286]  0.002 [0.012]  N.S.  0.001 [0.008]    Exposure  N.S.  0.003 [0.000]  N.S.  0.003 [0.000]  N.S.  0.003 [0.001]  N.S.  0.002 [0.001]    A*E  N.S.  0.003  N.S.  0.002  N.S.  0.003  N.S.  0.003    Islet(A, E)  0.0001  0.022  0.0001  0.027  0.0001 [N.T.]  0.022  0.0001  0.026  Rats + islet(R)  Rats  N.S.  0.001  N.S.  0.001  N.S.  0.002  N.S.  0.002    Islet(R)  0.0001  0.027  0.0001  0.034  0.0001  0.028  0.0001  0.030  Seagulls + islet(S)  Seagulls  0.0465  0.007  0.0388  0.008  0.0299  0.008  0.0340  0.007    Islet(S)  0.0001  0.024  0.0001  0.030  0.0001  0.025  0.0001  0.027    Residual    0.178    0.201    0.161    0.175      Whole dataset  Core (90% of specimens)      Unweighed Unifrac  Jaccard 97%  Jaccard sub-OTU 98%  Jaccard sub-OTU 99%  Model  Term  P-value  Variation  P-value  Variation  P-value  Variation  P-value  Variation        estimate    estimate    estimate    estimate  Islet  Islet  0.0001  0.028  0.0001  0.035  0.0001  0.029  0.0001  0.031  Age + islet(A)  Age  N.S. [0.0001]  0.004 [0.011]  N.S. [0.0001]  0.006 [0.012]  N.S. [0.0001]  0.005 [0.012]  N.S.  0.003 [0.008]    Islet(A)  0.0001  0.024 [0.022]  0.0001  0.029 [0.028]  0.0001  0.024 [0.023]  0.0001  0.028 [0.027]  Exposure + islet(E)  Exposure  N.S.  0.000  N.S.  0.000  N.S.  0.000  N.S.  0.000    Islet(E)  0.0001  0.028  0.0001  0.034  0.0001  0.029  0.0001  0.031  Age + exposure + islet(A, E)  Age  N.S.  0.002 [0.010]  N.S.  0.005 [0.012]  N.S. [0.0286]  0.002 [0.012]  N.S.  0.001 [0.008]    Exposure  N.S.  0.003 [0.000]  N.S.  0.003 [0.000]  N.S.  0.003 [0.001]  N.S.  0.002 [0.001]    A*E  N.S.  0.003  N.S.  0.002  N.S.  0.003  N.S.  0.003    Islet(A, E)  0.0001  0.022  0.0001  0.027  0.0001 [N.T.]  0.022  0.0001  0.026  Rats + islet(R)  Rats  N.S.  0.001  N.S.  0.001  N.S.  0.002  N.S.  0.002    Islet(R)  0.0001  0.027  0.0001  0.034  0.0001  0.028  0.0001  0.030  Seagulls + islet(S)  Seagulls  0.0465  0.007  0.0388  0.008  0.0299  0.008  0.0340  0.007    Islet(S)  0.0001  0.024  0.0001  0.030  0.0001  0.025  0.0001  0.027    Residual    0.178    0.201    0.161    0.175  N.S.: not significant (P > 0.05). N.T.: no test performed. Rarefaction depth used: whole dataset (38 160 sequences); core sub-OTU 98% (13 600), sub-OTU 99% (11 000). In brackets are estimates for ‘age’ when excluding Rei from the ‘old’ category (see Table S1, Supporting Information). Residual are common to all models. View Large Despite a main significant islet effect, a large fraction of the P. lilfordi gut microbiota was common to all populations: according to Jaccard pairwise distances across islets, mean distance between any two islets was only 0.38, with the Aire and PrCv representing the closest pair (mean Jaccard = 0.28). Hereafter, we investigated both shared and variable microbiota components driving the above pattern. Low microbiota uniqueness and differential enrichment among islets To explore the variable component responsible for the observed clustering and distances, we considered both taxa/OTU that are islet-specific in terms of presence/absence (i.e. uniqueness) and those that are differentially enriched. Taxa and OTUs unique to one islet can either represent novel gut acquisitions from the local pool, contamination by allochthonous bacteria, or simply loss of these taxa in all other islets. Total number of OTUs per islet, without rarefaction, ranged between 3013 and 4816 for a comparable sampling effort (five individuals per islet, except for Rei, with three individuals sampled) (Table 2). A total of 1093 OTUs (17% of total OTUs) were shared among all islets (although not all specimens). Uniqueness by islet was remarkably low, accounting for only between 1.1% and 5.2% of an islet's OTU content (average 2.4%) and showing limited distribution (1–2 individuals). This suggests that we either failed to sample them in most specimens or, more likely, that they represent part of the allochthonous component. Of the 666 total unique OTUs detected (Table 2), none corresponded to unique taxa when checked against the unfiltered OTU table including all reads. In all cases, most unique OTUs occurred in Aire, with Tosqueta (Tosq) hosting the lowest uniqueness. Overall, at least 95% of each islet OTU content was shared between two or more islets. The same analyses performed on the rarefied OTU table returned highly comparable results. Table 2. OTU uniqueness and shared component.   Number of OTUs      Islet  Totala  Unique (%)  Shared %b    Tosq  3287  37 (1.13)  98.87    Rei  3013  61 (2.02)  97.98    Ravl  3370  71 (2.11)  97.89    Colm  4334  82 (1.89)  98.11    PrCv  4816  87 (1.81)  98.19    AdGr  4133  110 (2.66)  97.34    Aire  4219  218 (5.17)  94.83    Average  3881.7  95.1 (2.40)  97.60      Total  Unique  Sharedb  Commonc  90% of specimens  All islets  6380  666  5714  1093  134    Number of OTUs      Islet  Totala  Unique (%)  Shared %b    Tosq  3287  37 (1.13)  98.87    Rei  3013  61 (2.02)  97.98    Ravl  3370  71 (2.11)  97.89    Colm  4334  82 (1.89)  98.11    PrCv  4816  87 (1.81)  98.19    AdGr  4133  110 (2.66)  97.34    Aire  4219  218 (5.17)  94.83    Average  3881.7  95.1 (2.40)  97.60      Total  Unique  Sharedb  Commonc  90% of specimens  All islets  6380  666  5714  1093  134  aof OTUs was estimated from a non-rarefied dataset, considering OTUs > 10 total sequences. bby two or more islets. cin all islets, but not all specimens. View Large As for presence/absence, differences in relative abundances of taxa/OTUs among islets were also minor. According to enrichment analyses, only a few taxa (Table S2, Supporting Information) and OTUs (total of 30, Table 3) showed a significant differential representation among islets (kw, FDR P < 0.05). Among indicator taxa (Table S2, Supporting Information), the phylum Synergistetes and its only recognized family Synergistaceae were a signature of Aire and PrCv, with only a few reads detected in AdGr and Colom (Colm). Aire had a particularly high representation of this family in terms of reads (376) and number of specimens (all). The islet of Rei showed a particular enrichment in Oceanospirillales and Halomonadaceae, while Ravells (Ravl) was enriched in the genus Turicibacter. Table 3. Differentially enriched OTUs among islets (Kruskal Wallis, FDR P < 0.05).   Mean of reads    Finest taxonomy  BLASTN Hits  OTU_IDa  Aire  AdGr  Rei  Tosq  Ravl  Colm  PrCv  Specificity  Phylum  Order  Family  Genus/species  Accession number  Identity (%)  Hostb  N.Ref.OTU895  516.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Peptococcaceae    GQ502520  91.8  Termite  721569  26.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    CU926622  98.9    175037  22.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Christensenellaceae    FJ362815  97.7  Human  NC.Ref.OTU185722  20.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      FJ535004  96.8    NC.Ref.OTU27689  20.6  0  0  0  0  0  0  Aire  Bacteroidetes  Bacteroidales  Bacteroidaceae  Bacteroides  HQ745181  96.7  Human  NC.Ref.OTU224538  18.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Peptococcaceae    GQ502520  90.8  Termite  NC.Ref.OTU27099  17.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      LC028698  97.3  Buffalo  NC.Ref.OTU212111  15  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    GQ294881  96.4  Mouse  566434  14  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      FJ535004  97.1    NC.Ref.OTU43839  10.4  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    AB218341  92.3  Dugong  NC.Ref.OTU59527  0  0  0  0  0  0  10.6  PrCv  Firmicutes  Clostridiales      AF371819  96.8  Swine  688934  0  0  0  0  25.2  0  0  Ravl  Proteobacteria  Enterobacteriales  Enterobacteriaceae  Serratia  KY124195  98.3    NC.Ref.OTU106014  0  0  0  0  22.8  0  0  Ravl  Proteobacteria  Enterobacteriales  Enterobacteriaceae    FR848151  98.1  Millipedes  182970  128.6  0  0  0  0  0  664.2  Aire, PrCv  Verrucomicrobia  Verrucomicrobiales  Verrucomicrobiaceae  Akkermansia  LT629973  99.1  Python  N.Ref.OTU217  46.4  0  0  0  0  0  66.2  Aire, PrCv  Firmicutes  Clostridiales      KU506346  94.4  Pig  N.Ref.OTU651  28  0  0  0  0  0  73.4  Aire, PrCv  Bacteroidetes  Bacteroidales  [Barnesiellaceae]    HQ792288  91.5  Human  NC.Ref.OTU124139  108  0  0  0  0  0  6  Aire, PrCv  Firmicutes  Clostridiales      AB606338  95.7  Mouse  NC.Ref.OTU110653  67.6  0  0  0  0  0  0.4  Aire, PrCv  Firmicutes  Clostridiales      JQ248100  97.3  Chicken  NC.Ref.OTU229167  2.2  0  0  0  0  0  16.8  Aire, PrCv  Firmicutes  Clostridiales  Ruminococcaceae    HM123952  95.5  Human  3970651  1041.6  0  0  0  0  23  0  Aire, Colm  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  JX457157  99.3  Cockroach  NC.Ref.OTU12308  15.6  0  0  0  0  0.2  0  Aire, Colm  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  JX457157  96.7  Cockroach  NC.Ref.OTU231957  7.4  0  12  0  0  0  0  Aire, Rei  Firmicutes  Clostridiales  Ruminococcaceae    CU919465  94.4    NC.Ref.OTU25986  0  0  0  0  100.4  0  3.2  Ravl, PrCv  Firmicutes  Clostridiales      AB606272  94.3  Mouse  347529  25  6.2  0  0  329.2  5.2  26.6  Several  Firmicutes  Turicibacterales  Turicibacteraceae  Turicibacter  LT624967  99.1    727945  8.4  2.6  0  0  11  1.4  2.2  Several  Firmicutes  Clostridiales  Eubacteriaceae  Anaerofustis  LT223622  99.3  Human  NC.Ref.OTU117227  1.4  0  0  0  0.2  0  16.4  Several  Firmicutes  Clostridiales  Lachnospiraceae  Coprococcus  EU779145  96.6  Sheep  NC.Ref.OTU63406  7  42.2  0  0  0  0  54.2  Several  Bacteroidetes  Bacteroidales  Bacteroidaceae  Bacteroides  GU105959  92.8  Human  N.Ref.OTU255  10  0  0  0.2  0  14.2  7.2  Several  Firmicutes  Clostridiales  [Mogibacteriaceae]  Anaerovorax  JQ084279  94.8  Mouse  N.Ref.OTU41  274.4  134  3.7  0  0  0  86.2  Several  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  HQ767608  94.3  Human  N.Ref.OTU544  251.6  0.2  46.7  0  0  0  51.6  Several  Bacteroidetes  Bacteroidales  Porphyromonadaceae  P. distasonis  HQ797049  97.4  Human    Mean of reads    Finest taxonomy  BLASTN Hits  OTU_IDa  Aire  AdGr  Rei  Tosq  Ravl  Colm  PrCv  Specificity  Phylum  Order  Family  Genus/species  Accession number  Identity (%)  Hostb  N.Ref.OTU895  516.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Peptococcaceae    GQ502520  91.8  Termite  721569  26.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    CU926622  98.9    175037  22.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Christensenellaceae    FJ362815  97.7  Human  NC.Ref.OTU185722  20.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      FJ535004  96.8    NC.Ref.OTU27689  20.6  0  0  0  0  0  0  Aire  Bacteroidetes  Bacteroidales  Bacteroidaceae  Bacteroides  HQ745181  96.7  Human  NC.Ref.OTU224538  18.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Peptococcaceae    GQ502520  90.8  Termite  NC.Ref.OTU27099  17.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      LC028698  97.3  Buffalo  NC.Ref.OTU212111  15  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    GQ294881  96.4  Mouse  566434  14  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      FJ535004  97.1    NC.Ref.OTU43839  10.4  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    AB218341  92.3  Dugong  NC.Ref.OTU59527  0  0  0  0  0  0  10.6  PrCv  Firmicutes  Clostridiales      AF371819  96.8  Swine  688934  0  0  0  0  25.2  0  0  Ravl  Proteobacteria  Enterobacteriales  Enterobacteriaceae  Serratia  KY124195  98.3    NC.Ref.OTU106014  0  0  0  0  22.8  0  0  Ravl  Proteobacteria  Enterobacteriales  Enterobacteriaceae    FR848151  98.1  Millipedes  182970  128.6  0  0  0  0  0  664.2  Aire, PrCv  Verrucomicrobia  Verrucomicrobiales  Verrucomicrobiaceae  Akkermansia  LT629973  99.1  Python  N.Ref.OTU217  46.4  0  0  0  0  0  66.2  Aire, PrCv  Firmicutes  Clostridiales      KU506346  94.4  Pig  N.Ref.OTU651  28  0  0  0  0  0  73.4  Aire, PrCv  Bacteroidetes  Bacteroidales  [Barnesiellaceae]    HQ792288  91.5  Human  NC.Ref.OTU124139  108  0  0  0  0  0  6  Aire, PrCv  Firmicutes  Clostridiales      AB606338  95.7  Mouse  NC.Ref.OTU110653  67.6  0  0  0  0  0  0.4  Aire, PrCv  Firmicutes  Clostridiales      JQ248100  97.3  Chicken  NC.Ref.OTU229167  2.2  0  0  0  0  0  16.8  Aire, PrCv  Firmicutes  Clostridiales  Ruminococcaceae    HM123952  95.5  Human  3970651  1041.6  0  0  0  0  23  0  Aire, Colm  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  JX457157  99.3  Cockroach  NC.Ref.OTU12308  15.6  0  0  0  0  0.2  0  Aire, Colm  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  JX457157  96.7  Cockroach  NC.Ref.OTU231957  7.4  0  12  0  0  0  0  Aire, Rei  Firmicutes  Clostridiales  Ruminococcaceae    CU919465  94.4    NC.Ref.OTU25986  0  0  0  0  100.4  0  3.2  Ravl, PrCv  Firmicutes  Clostridiales      AB606272  94.3  Mouse  347529  25  6.2  0  0  329.2  5.2  26.6  Several  Firmicutes  Turicibacterales  Turicibacteraceae  Turicibacter  LT624967  99.1    727945  8.4  2.6  0  0  11  1.4  2.2  Several  Firmicutes  Clostridiales  Eubacteriaceae  Anaerofustis  LT223622  99.3  Human  NC.Ref.OTU117227  1.4  0  0  0  0.2  0  16.4  Several  Firmicutes  Clostridiales  Lachnospiraceae  Coprococcus  EU779145  96.6  Sheep  NC.Ref.OTU63406  7  42.2  0  0  0  0  54.2  Several  Bacteroidetes  Bacteroidales  Bacteroidaceae  Bacteroides  GU105959  92.8  Human  N.Ref.OTU255  10  0  0  0.2  0  14.2  7.2  Several  Firmicutes  Clostridiales  [Mogibacteriaceae]  Anaerovorax  JQ084279  94.8  Mouse  N.Ref.OTU41  274.4  134  3.7  0  0  0  86.2  Several  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  HQ767608  94.3  Human  N.Ref.OTU544  251.6  0.2  46.7  0  0  0  51.6  Several  Bacteroidetes  Bacteroidales  Porphyromonadaceae  P. distasonis  HQ797049  97.4  Human  Taxonomy in brackets is not officially recognized (under revision). aOTUs were filtered for total sum of means > 10. bPresence in gut or feces samples. Blanks correspond to putative environmental hits. View Large Among indicator OTUs (Table 3), the majority belonged to the order Clostridiales (Firmicutes). Except for six putative environmental OTUs, all others had best hits to gut/feces of animals, typically vertebrates. Of them, 13 were single islet-specific, 10 being uniquely found in Aire. One Aire-specific OTU (N.Ref.OTU895) was present in high abundance (516 mean reads) in most specimens and had best BLASTN hit with low identity (91.8%) to the termite gut. Six OTUs were specific to Aire and PrCv, three of them with large mean read abundance in both islets (Table 3). OTU-182970, corresponding to a member of the genus Akkermansia, was particularly abundant in PrCv and had the best BLASTN hit to the reticulated python gut. This is overall the variable component responsible of the islet-based compositional structure of the P. lilfordi microbiota. Hereafter, we focused on the large shared component, the core microbiota. Core composition of the P. lilfordi gut microbiota We considered core taxa those with presence in at least 90% of specimens and all islets (calculated on the non-rarefied dataset). The core comprised 7 phyla, 12 classes, 13 orders, 25 families, 27 genera and 7 species (for species we considered presence in at least 80% of specimens) (Table 4). Notably, the only represented order within Cyanobacteria was YS2, a recognized member of the novel non-photosynthetic phylum Melainabacteria (Di Rienzi et al.2013) that occurred at very low abundances. Among core genera, most abundant were Bacteroides (average 10% of total reads per specimen) and Parabacteroides (4.6%), followed by Oscillospira, Dorea, Ruminococcus, Desulfovibrio and Coprococcus (all above 1%) (Fig. S3, Supporting Information). Table 4. List of core taxa (90%) by their finest taxonomic resolution, with number of core OTUs per taxon shown in parentheses (total = 134). Phylum  Class  Order  Family  Genus  Species1  Firmicutes (118)  Erysipelotrichi (5)  Erysipelotrichales (5)  Erysipelotrichaceae (5)  Coprobacillus (1)            PSB-M-3            Holdemania      Clostridia (113)  Clostridiales (113)  Lachnospiraceae (37)  Ruminococcus (8)  [Ruminococcus] gnavus          Blautia (2)  Blautia producta          Dorea (10)            Coprococcus            Roseburia          Ruminococcaceae (51)  Oscillospira (24)            Anaerotruncus          Eubacteriaceae (1)  Pseudoramibacter (1)            Anaerofustis          Clostridiaceae  Clostridium          Dehalobacteriaceae  Dehalobacterium          Christensenellaceae (4)  Christensenella (1)          [Mogibacteriaceae] (3)  Anaerovorax          Peptostreptococcaceae (1)            Peptococcaceae (3)            Veillonellaceae            [Tissierellaceae]        Bacilli  Lactobacillales  Lactobacillaceae  Lactobacillus          Leuconostocaceae          Bacillales  Bacillaceae  Bacillus    Actinobacteria (6)  Coriobacteriia (6)  Coriobacteriales (6)  Coriobacteriaceae (6)  Adlercreutzia (1)    Proteobacteria (6)  Deltaproteobacteria (6)  Desulfovibrionales (6)  Desulfovibrionaceae (6)  Desulfovibrio (5)            Bilophila      Gammaproteobacteria  Enterobacteriales  Enterobacteriaceae        Alphaproteobacteria  Campylobacterales  Helicobacteraceae        Epsilonproteobacteria          Bacteroidetes (4)  Bacteroidia (4)  Bacteroidales (4)  Bacteroidaceae (3)  Bacteroides (3)  Bacteroides fragilis            Bacteroides ovatus        Porphyromonadaceae (1)  Parabacteroides (1)  Parabacteroides distasonis            Parabacteroides gordonii        Rikenellaceae  Rikenella          [Odoribacteraceae]  Odoribacter    Verrucomicrobia  Verrucomicrobiae  Verrucomicrobiales  Verrucomicrobiaceae  Akkermansia  Akkermansia muciniphila  Tenericutes  Mollicutes  Anaeroplasmatales  Anaeroplasmataceae  Anaeroplasma        RF39        Cyanobacteria (Melainabacteria)  4C0d-2  YS2        Phylum  Class  Order  Family  Genus  Species1  Firmicutes (118)  Erysipelotrichi (5)  Erysipelotrichales (5)  Erysipelotrichaceae (5)  Coprobacillus (1)            PSB-M-3            Holdemania      Clostridia (113)  Clostridiales (113)  Lachnospiraceae (37)  Ruminococcus (8)  [Ruminococcus] gnavus          Blautia (2)  Blautia producta          Dorea (10)            Coprococcus            Roseburia          Ruminococcaceae (51)  Oscillospira (24)            Anaerotruncus          Eubacteriaceae (1)  Pseudoramibacter (1)            Anaerofustis          Clostridiaceae  Clostridium          Dehalobacteriaceae  Dehalobacterium          Christensenellaceae (4)  Christensenella (1)          [Mogibacteriaceae] (3)  Anaerovorax          Peptostreptococcaceae (1)            Peptococcaceae (3)            Veillonellaceae            [Tissierellaceae]        Bacilli  Lactobacillales  Lactobacillaceae  Lactobacillus          Leuconostocaceae          Bacillales  Bacillaceae  Bacillus    Actinobacteria (6)  Coriobacteriia (6)  Coriobacteriales (6)  Coriobacteriaceae (6)  Adlercreutzia (1)    Proteobacteria (6)  Deltaproteobacteria (6)  Desulfovibrionales (6)  Desulfovibrionaceae (6)  Desulfovibrio (5)            Bilophila      Gammaproteobacteria  Enterobacteriales  Enterobacteriaceae        Alphaproteobacteria  Campylobacterales  Helicobacteraceae        Epsilonproteobacteria          Bacteroidetes (4)  Bacteroidia (4)  Bacteroidales (4)  Bacteroidaceae (3)  Bacteroides (3)  Bacteroides fragilis            Bacteroides ovatus        Porphyromonadaceae (1)  Parabacteroides (1)  Parabacteroides distasonis            Parabacteroides gordonii        Rikenellaceae  Rikenella          [Odoribacteraceae]  Odoribacter    Verrucomicrobia  Verrucomicrobiae  Verrucomicrobiales  Verrucomicrobiaceae  Akkermansia  Akkermansia muciniphila  Tenericutes  Mollicutes  Anaeroplasmatales  Anaeroplasmataceae  Anaeroplasma        RF39        Cyanobacteria (Melainabacteria)  4C0d-2  YS2        Taxonomy in brackets is not officially recognized (under revision). Not all core OTUs were identified below order level. aPresence in 80% of the specimens and all islets. View Large At OTU level, 23 OTUs were present in 100% of the specimens, 134 at 90% and 303 at 80%. Four phyla hosted the 134 core OTUs at 90% prevalence, Firmicutes, Proteobacteria, Actinobacteria and Bacteroidetes (Table 4), with the majority of core OTUs (118) belonging to the phylum Firmicutes and, within this, to the order Clostridiales (113), families Ruminococcaceae (51), Lachnospiraceae (37), and genera Oscillosphira (24), Ruminococcus (8) (both Ruminococcaceae) and Dorea (10, Lachnospiraceae) (Table 4). Despite the fact that Bacteroidetes and Proteobacteria were among the most abundant phyla, they hosted only four and six core OTUs, respectively. All Proteobacteria OTUs belonged to the family Desulfovibronaceae, five of which were classified to the genus Desulfovibrio. All six Actinobacteria core OTUs were Coriobacteriaceae. Overall, these 134 core OTUs accounted for only 2.1% of the total number of OTUs detected, with individual core OTUs typically occurring at low abundances, ranging on average across all specimens between 0.01% and 2.7% of total reads per specimen (Table S3, Supporting Information). The most abundant core OTU was N.Ref.OTU202 (2.7%, classified as Bacteroides), followed by OTU-113 542 (1.6%, Desulfovibrio), N.Ref.OTU980 (1.6%, Erysipelotrichaceae) and OTU-194 286 (1.5%, Parabacteroides gordonii). Altogether, the 134 core OTUs contributed 34.8% of the total reads across all samples and 35.1% average reads per specimen. According to best BLASTN hits and taxonomic classification (Table S3, Supporting Information), core OTUs (i) typically belonged to bacteria known as anaerobic or strictly anaerobic, suggesting a limitation in their habitat range and potential for dispersal; (ii) were not classified as pathogens at present abundance; and (iii) in most cases matched to stool/gut from humans and other vertebrates, rarely invertebrates (representing the main food source). These features suggest large specificity to the vertebrate gut, making these OTUs good candidates to investigate microbiota heritability across host generations and the effect of long-term geographic isolation. Sub-OTU resolution recovers islet signal To address the above scenario, we exploited the sub-OTU diversity of the core (sequence variants within each OTU as inferred by re-clustering of sequences at 98 and 99% similarity). If core OTUs represent a conserved component inherited from the ancestral microbial pool and then ‘confined’ in each islet after population splits, we would expect the sub-OTU variants to carry a signature of their islet-restricted diversification, including local drift and mutation/recombination events. On the contrary, if acquisition of this core occurred by more recent dispersal across islets, the sub-OTU variants should present a stochastic distribution or correlate with geographic distances, but not with islet structure. For this purpose, we used individual and pooled core OTUs, which carry no information for discrimination among specimens/islets, and reclustered all reads at 98 and 99% similarity (the latter approximating the haplotype level). PERMANOVA analyses on all previously tested categorical variables indicated ‘islet’ as the main significant factor (Table 1, P = 0.0001), explaining 15% of the core microbiota sub-OTU variance, a proportion comparable to the one explained by clustering of the whole dataset at 97% similarity. According to Jaccard distance matrices built on de novo 98% and 99% OTUs, distances were significantly shorter within islets than between islets (two-sided Student's two-sample t-test, Bonferroni-corrected P < 0.001). Moreover, mean Jaccard distances across islets mirrored the 97% OTU-based distances: in all cases, the closest islet pair was Aire-PrCv (mean Jaccard = 0.25 and 0.30 at 98% and 99% clustering, respectively). Again, microbiota distances did not correlate with geographic distances (Mantel r = –0.08838, P = 0.752 for 99% clustering and r = –0.10154, P = 0.674 for 98% clustering). To understand the individual contribution of core OTUs to this pattern, the 23 OTUs with presence in 100% of the specimens were individually analyzed (99% OTU clustering, Table S4, Supporting Information). Intra-OTU variation analysis revealed that 18 of the 23 OTUs carried a significant signal by islet (adonis, P < 0.05), explaining between 14% and 23% of the total sub-OTU variation. The five OTUs with no significance had their best BLASTN hits to putative environmental locations, human skin and a novel bacterium from the human gut (i.e. Emergencia timonensis). DISCUSSION Here we provide the first characterization of the P. lilfordi gut microbiota, exploring its compositional variation in adult males from seven isolated populations inhabiting coastal islets around Menorca. Taking advantage of the peculiar phylogeographic history of these populations and the expected minor diversifying role of two crucial variables in microbiota studies (host genotype and diet), we performed a preliminary investigation of the forces driving current patterns of microbiota similarity/dissimilarity among populations, including inheritance, ecological drift, dispersal and local environmental inputs. Geographic separation and age drive the diversity and structuring of the Podarcis microbiota Compatible with a significant spatial and temporal effect, islet was the most significant factor in structuring the Podarcis microbiota, accounting for 15% of the observed microbiota variation (Table 1). Aire in particular, with its geographic remoteness and morphological diversity (it hosts a unique melanic form) (Salvador 2009), presented the most distinct microbial community as indicated by Unifrac distances (Fig. 3) as well as by taxa/OTU uniqueness and enrichment (Tables 2 and 3). Among the fixed factors, within which islet was nested, only ‘age’ was significant when considering Aire and PrCv as the oldest islets. Most of the microbiota variation remained unexplained by all other factors tested, including geographic distances among islets. Consistent with a minor although significant islet effect, microbiota beta-diversity distances across populations were relatively short (average pairwise Jaccard distance = 0.38) and partly associated with time since islet/population isolation. Although the phylogeographic history of the Menorcan populations is far from being fully resolved (and out of the scope of this study), Aire and PrCv represent the two first islets/populations to separate from mainland about 8000–9000 years (Pérez-Mellado 1989), in agreement with both bathymetry and curves of sea level rise in the Mediterranean during Holocene (Vecchi et al.2016). Despite being the most geographically distant islets within our sample set (Fig. 1), their gut microbial composition showed the shortest pairwise Jaccard distances and presented several unique commonalities (Table 3 and Table S2, Supporting Information). Retention of ancestral taxa in these two islets or convergence driven by similar environmental opportunities (open sea and scarce vegetation) could both be invoked. However, none of the ecological variables tested with PERMANOVA were strictly associated with these two islets (Table S1, Supporting Information), largely excluding convergence by niche selection. In terms of gut community similarities across islets, OTU uniqueness per islet/population was remarkably low, with the large majority of OTUs (98% on average, Table 2) occurring in more than one islet. Quantitative differences in common taxa and OTUs representation across islets were also minor (Table 3 and Table S2, Supporting Information) and associated with bacteria belonging to core taxa (i.e. order Clostridiales and family Ruminococcaceae, Table 4), suggesting minor functional community differences among islets. While, to some extent, this large microbial membership and structural similarity might result from common host genotypic constraints and lack of major dietary differences (all populations were largely omnivorous), important differences in ecological niches exist among islets: some show rich vegetation and large anthropic impact (Colm and Rei), while others are characterized by more natural environments, with scarce vegetation and large exposure to the open sea (Aire and PrCv). Overall, this suggests that the striking compositional similarity observed is, at least in part, independent from the environmental inputs/pressures. Ancestral retention versus dispersal effect To explore the origin of the shared taxonomic content among populations, we looked at the nature of core OTUs (defined as presence in 90% of the specimens and all islets). Core OTUs could either represent a conserved inherited component since divergence from the mainland population or an independent acquisition after isolation through taxa dispersal across islets followed by gut selective retention (common niche selection). Excluding recent host immigration events (at least during the past 2000 years), as supported by bathymetric data and the degree of morphological and pigmentation diversity among these populations (Salvador 2009), dispersal should be largely restricted to other biotic or abiotic vectors able to cross the sea barrier. Invertebrates consumed with the diet represent the most likely biotic vectors, particularly flying insects (i.e. all populations are predominantly insectivorous, especially during the summer season) and occasional mollusks and crustaceans captured along the shores. According to best BLASTN hits and taxonomic classification, however, the majority of core OTUs typically matched gut/stool samples of vertebrates, with rare invertebrate hits (Table S3, Supporting Information). This suggests that despite invertebrates (mostly insects) represent the largest component of these lizards’ diet, they are not a major source of gut bacteria. This is in accordance with previous findings showing that invertebrate communities did not significantly contribute to the lizard gut microbiota (Ren et al.2016). Likewise, vertebrates (including rats, rabbits, seagull colonies and humans) are improbable vectors of pervasive dispersal of shared OTUs, as their presence is largely restricted to a few islets (Table S1, Supporting Information; except for Colm and Rei, none of the islets hosts stable human settlements). Aside from animal vectors, recent studies have shown that the phyllosphere microbiome might play a substantial role in shaping the gut microbiota in both lizards (Kohl et al.2017) and desert woodrats (Kohl and Dearing 2012); in this respect, P. lilfordi populations are known to consume vegetable matter (Salvador 2009; Perez-Cembranos, Leon and Perez-Mellado 2016), a finding that is also confirmed by stomach content analysis from this study, revealing plant material and flowers in the majority of samples, and by the presence of several fiber-degrader taxa in the gut (see below). While our current lack of microbial samples from the surrounding environment does not allow us to test the impact of the phyllosphere in these lizards’ microbiota, both best BLASTN matches and large differences in vegetation cover among islets make this scenario improbable. A second important observation is that the 98% and 99% OTU-based structure of the core presented a significant clustering by islet and age (Table 1), which is compatible with an islet-restricted diversification of core OTUs (putatively driven by loss of variants, mutation and recombination events). This hypothesis clearly needs to be corroborated by strain-level analyses. Moreover, these data might suffer from undersampling in some populations/individuals, despite the high sequence depth coverage per specimen obtained (38 160 sequences after rarefaction). Finally, whether any of these sub-OTU variants are under selection or exist in a transitory neutral state remains to be investigated. Indirect evidence of gut microbiome heritability Clearly, the presence of shared OTUs, that are unlikely the result of dispersal, indirectly stands for the inheritance of a microbial component over relatively large temporal and spatial scales. Up to date, levels and mechanisms of gut microbiota inheritance in lizards, and reptiles in general, remain unclear (Colston 2017). Gut microbiotas from Anolis lizards showed a remarkably high intraspecific variation, with only 10% on average of shared OTUs between any two conspecifics in nature (Ren et al.2016), suggesting low microbiota heritability in these lizards. Although their large generalist diet was proposed as a major diversifying factor, this is unlike our findings for P. lilfordi, also a generalist species, where any two conspecifics within an islet share on average 31% of their OTUs and populations are allopatric. Methodological differences, such as the use of full intestines in our study versus fecal samples for Anolis, could be partly responsible for the observed pattern. In contrast, for viviparous lizards of genera Liolaemus and Phymaturus, recent evidence indicated around 34% of microbiome heritability in captivity (Kohl et al.2017), an estimate that nicely matches the percentage of core microbial sequences observed in Podarcis (34.8%). While the Podarcis core summed up to more than one-third of the microbiota, individual core OTUs typically occurred at low abundances (mean = 0.26 ± 0.4% over total reads per specimen, Table S3, Supporting Information), suggesting that, while being more vulnerable to stochastic loss (Chase and Myers 2011), their persistence over time is sustained by important selective forces and putative high inheritability. This is compatible with the general understanding that core gut microbiotas in vertebrates are rather stable in the absence of differential adaptive pressures (such as type of metabolism) (Ley et al.2008; Lozupone et al.2012). In this respect, the overall taxonomic composition of the P. lilfordi gut microbiota mirrored that found in other vertebrates (Fig. 2), with the majority of core families and genera detected (Table 4) representing stable residents of human and mouse guts (Xiao et al.2015), and being consistently found in other lizards (Ren et al.2016; Kohl et al.2017). Core taxa were particularly enriched in bacteria known as active fiber degraders, including the families Ruminococcaceae and Lachnospiraceae (Fig. 2, both Clostridiales) (Biddle et al.2013) and the genera Bacteroides, Blautia, Akkermansia, Oscillospira and Desulfovibrio (Fig. S3, Supporting Informtion). Bacteroides, for instance, which was the most conspicuous genus in our samples, also makes up a very large proportion of the mammalian gut where it is an active degrader of plant material (Ley et al.2008; Xiao et al.2015). Oscillospira, also highly abundant in humans, where it shows high heritability and is associated with leanness, is largely enriched in plant-fed lizards (Kohl et al.2017). Overall, the P. lilfordi core composition is compatible with major conserved functions in vertebrates and the consumption of plant material, suggesting common selective pressures for its maintenance across host generations. Nonetheless, considering that parental care in reptiles is rare and not documented in the genus Podarcis (oviparous), means of bacteria recycling at each new host generation remain unclear. Coprophagy, cannibalism or other types of social behavior have been proposed as mechanisms for bacteria transmission (Colston and Jackson 2016; Perez-Cembranos, Leon and Perez-Mellado 2016; Colston 2017); however, they might suffer from much stochasticity, compromising their ability to preserve a metacommunity structure over time. Covering the eggs with a bacterial layer by means of fecal contamination during oviposition through the single urogenital opening present in lizards (the cloaca) might represent a more reliable mechanism for a stable microbiota transmission. Overall, although our results provided important indirect evidence for gut microbiome heritability in P. lilfordi, confirming previous data for lizards (Kohl et al.2017), this topic definitely demands more empirical investigation. Islet microbiotas as nested subsets of the ancestral mainland microbial pool: a hypothesis To summarize, our main findings indicated that (i) there is a large taxonomic microbiota conservation among all seven allopatric populations studied, down to OTU level, with a reduced impact of local exposure (i.e. uniqueness is negligible); (ii) the minor variation detected is explained by ‘islet’ and ‘age’, supporting a role of the host phylogeographic history in shaping microbiota composition; (iii) dispersal of OTUs among islets appears to be low; and (iv) core variation at 98%–99%-OTU level carries a significant ‘islet’ and ‘age’ signature, compatible with an islet-restricted evolution. Altogether, these findings suggest that persistence of taxa from the mainland stock and ecological drift might have played a substantial role in shaping current pattern of microbiota similarity/dissimilarity among P. lilfordi populations. This intriguing hypothesis clearly requires more extensive quantitative data, targeting additional specimens, islets and microbial markers, as well as ad hoc statistical tools to test drift in host-associated microbiota across populations, a field still in its infancy. Although overlooked in current literature, host population demographic history is likely to represent a crucial factor in explaining current gut metacommunity diversity, as host drift could, to some extent, cause microbiota ecological drift (Orrock and Watling 2010; Chase and Myers 2011). Up to date, only few studies have attempted to address this topic. In the Galápagos archipelago, processes of ecological drift were proposed to explain intraspecific differences in marine iguanas (Lankau, Hong and Mackie 2012). Likewise, in two Anolis lizard species, the complex pattern of similarity/dissimilarity observed, only partially attributed to geography and with no differences among ‘ecomorphs’ (Ren et al.2016), could be potentially explained by chance variation driven by drift and dispersal. In our study system, population densities varied greatly across islets (Fig. 1), with Aire presenting one of the highest population densities of lizards/vertebrates worldwide (Pérez-Mellado et al.2008). Although past population demography for these islets is poorly known (Pérez-Mellado et al.2008), the largest populations (Aire, PrCv and Colm) carried slightly higher diversities (PD whole tree) compared to smaller islets (Tosq, Ravl and Rei), though differences were not significant (Fig. S2, Supporting Information). The small population found in Tosqueta, in particular, which showed the strongest genetic drift among all 16 islets according to unpublished microsatellite data (personal observations), carried the least diverse microbial community and lowest uniqueness, despite the extensive vegetation cover and ample opportunities for local inputs in this islet. These observations, although preliminary, point to an interesting correlation between population demographic history and drift effect in gut metacommunities, which could potentially explain a large part of the extant variation, a topic that should definitely warrant further investigation. CONCLUSIONS Overall, the minor dietary and genotypic differences among our sampled populations (suggesting weak diversifying selection) and the geographically restricted areas to which they are confined, reducing the exposure to microbial pools from other vertebrates, have most likely contributed to the large compositional stability of P. lilfordi gut microbiota, ultimately revealing the strength of drift. We envision that, owing to a progressive within-islet microbiota evolution, in a context of no host immigration events, the core OTU signature of these populations’ past host legacy will gradually disappear, while resemblance at higher taxonomic levels will persist for a longer time. Integrating our data with a fine-scale resolution of P. lilfordi phylogeographic history and a bacterial strain-level characterization that targets core OTUs (under development) might offer a way to disentangle the dynamics of this microbiota component over such short evolutionary time scales. Finally, extending the sampling to all remaining populations from Balearic Islands, including the closely related species P. pityusensis, will provide important comparative data for a most integrative understanding of the gut microbiota diversification in these lizards. Data accessibility Raw sequences are available in the Sequence Read Archive (SRA) database at NCBI under BioProject ID PRJNA429123PRJNA429123. SUPPLEMENTARY DATA Supplementary data are available at FEMSEC online. Acknowledgements We thank Sam Pons, Jordi García-Darocas and Montserrat Farré for aid in sampling, and the Animal Biology department (University of Barcelona) for use of equipment. FUNDING This work was supported by funding from the Spanish Ministerio de Economía y Competitividad to LB (RYC-2012-11443). Conflicts of Interest. None declared. REFERENCES Aivelo T, Norberg A. Parasite-microbiota interactions potentially affect intestinal communities inwildmammals. J Anim Ecol  2017, doi: 10.1111/1365-2656.12708. Alcover JA, Moyá-Solà S, Pons-Moyá J. Les quimeres del passat. Els vertebrats fòssils del Plio-Quaternari de les Balears i Pitiüses . Palma, Espanya: Editorial Moll, 1981. Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecol  2001; 26: 32– 46. Anderson MJ, Gorley RN, Clarke KR. PERMANOVA+ for PRIMER: Guide to Software and Statistical Methods . Plymouth, UK: PRIMER-E, 2008. Baldo L, Pretus JL, Riera JL et al.   Convergence of gut microbiotas in the adaptive radiations of African cichlid fishes. ISME J  2017; 11: 1975– 87. Google Scholar CrossRef Search ADS PubMed  Baldo L, Riera JL, Tooming-Klunderud A et al.   Gut microbiota dynamics during dietary shift in eastern african cichlid fishes. PLoS One  2015; 10: e0127462. Google Scholar CrossRef Search ADS PubMed  Benson AK, Kelly SA, Legge R et al.   Individuality in gut microbiota composition is a complex polygenic trait shaped by multiple environmental and host genetic factors. P Natl Acad Sci USA  2010; 107: 18933– 8. Google Scholar CrossRef Search ADS   Biddle A, Stewart L, Blanchard J et al.   Untangling the genetic basis of fibrolytic specialization by lachnospiraceae and ruminococcaceae in diverse gut communities. Divers Distrib  2013; 5: 627– 40. Google Scholar CrossRef Search ADS   Brown RP, Terrasa B, Perez-Mellado V et al.   Bayesian estimation of post-Messinian divergence times in Balearic Island lizards. Mol Phyl Evol  2008; 48: 350– 8. Google Scholar CrossRef Search ADS   Burns AR, Stephens WZ, Stagaman K et al.   Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. ISME J  2016; 10: 655– 64. Google Scholar CrossRef Search ADS PubMed  Caporaso JG, Kuczynski J, Stombaugh J et al.   QIIME allows analysis of high-throughput community sequencing data. Nat Methods  2010; 7: 335– 6. Google Scholar CrossRef Search ADS PubMed  Chase JM, Myers JA. Disentangling the importance of ecological niches from stochastic processes across scales. Philos T Roy Soc B  2011; 366: 2351– 63. Google Scholar CrossRef Search ADS   Clarke KR, Gorley RN. PRIMER v7: User Manual/Tutorial . Plymouth, UK: PRIMER-E, 2015. Colston TJ. Gut microbiome transmission in lizards. Mol Ecol  2017; 26: 972– 4. Google Scholar CrossRef Search ADS PubMed  Colston TJ, Jackson CR. Microbiome evolution along divergent branches of the vertebrate tree of life: what is known and unknown. Mol Ecol  2016; 25: 3776– 800. Google Scholar CrossRef Search ADS PubMed  Costello EK, Stagaman K, Dethlefsen L et al.   The application of ecological theory toward an understanding of the human microbiome. Science  2012; 336: 1255– 62. Google Scholar CrossRef Search ADS PubMed  David LA, Maurice CF, Carmody RN et al.   Diet rapidly and reproducibly alters the human gut microbiome. Nature  2014; 505: 559– 63. Google Scholar CrossRef Search ADS PubMed  Delsuc F, Metcalf JL, Wegener Parfrey L et al.   Convergence of gut microbiomes in myrmecophagous mammals. Mol Ecol  2014; 23: 1301– 17. Google Scholar CrossRef Search ADS PubMed  Di Rienzi SC, Sharon I, Wrighton KC et al.   The human gut and groundwater harbor non-photosynthetic bacteria belonging to a new candidate phylum sibling to Cyanobacteria. eLife  2013; 2: e01102. Google Scholar CrossRef Search ADS PubMed  Funkhouser LJ, Bordenstein SR. Mom knows best: the universality of maternal microbial transmission. PLoS Biol  2013; 11: e1001631. Google Scholar CrossRef Search ADS PubMed  Kohl KD, Brun A, Magallanes M et al.   Gut microbial ecology of lizards: insights into diversity in the wild, effects of captivity, variation across gut regions and transmission. Mol Ecol  2017; 26: 1175– 89. Google Scholar CrossRef Search ADS PubMed  Kohl KD, Dearing MD. Experience matters: prior exposure to plant toxins enhances diversity of gut microbes in herbivores. Ecol Lett  2012; 15: 1008– 15. Google Scholar CrossRef Search ADS PubMed  Krijgsman W, Hilgen FJ, Raffi Y et al.   Chronology, causes and progression of the Messinian salinity crisis. Nature  1999; 400: 652– 5. Google Scholar CrossRef Search ADS   Lankau EW, Hong PY, Mackie RI. Ecological drift and local exposures drive enteric bacterial community differences within species of Galapagos iguanas. Mol Ecol  2012; 21: 1779– 88. Google Scholar CrossRef Search ADS PubMed  Ley RE, Lozupone CA, Hamady M et al.   Worlds within worlds: evolution of the vertebrate gut microbiota. Nat Rev Microbiol  2008; 6: 776– 88. Google Scholar CrossRef Search ADS PubMed  Lozupone CA, Stombaugh JI, Gordon JI et al.   Diversity, stability and resilience of the human gut microbiota. Nature  2012; 489: 220– 30. Google Scholar CrossRef Search ADS PubMed  Martinson VG, Douglas AE, Jaenike J. Community structure of the gut microbiota in sympatric species of wild Drosophila. Ecol Lett  2017; 20: 629– 39. Google Scholar CrossRef Search ADS PubMed  Moeller AH, Caro-Quintero A, Mjungu D et al.   Cospeciation of gut microbiota with hominids. Science  2016; 353: 380– 2. Google Scholar CrossRef Search ADS PubMed  Muegge BD, Kuczynski J, Knights D et al.   Diet drives convergence in gut microbiome functions across mammalian phylogeny and within humans. Science  2011; 332: 970– 4. Google Scholar CrossRef Search ADS PubMed  Nemergut DR, Schmidt SK, Fukami T et al.   Patterns and processes of microbial community assembly. Microbiol Mol Biol R  2013; 77: 342– 56. Google Scholar CrossRef Search ADS   Orrock JL, Watling JI. Local community size mediates ecological drift and competition in metacommunities. Proc Biol Sci  2010; 277: 2185– 91. Google Scholar CrossRef Search ADS PubMed  Perez-Cembranos A, Leon A, Perez-Mellado V. Omnivory of an insular lizard: sources of variation in the diet of Podarcis lilfordi (Squamata, Lacertidae). PLoS One  2016; 11: e0148947. Google Scholar CrossRef Search ADS PubMed  Pérez-Mellado V. Estudio ecológico de la lagartija balear Podarcis lilfordi (Günther, 1874) en Menorca. Rev de Menorca  1989; 53: 455– 511. Pérez-Mellado V, Hernández-Estévez JA, García-Díez T et al.   Population density in Podarcis lilfordi (Squamata, Lacertidae), a lizard species endemic to small islets in the Balearic Islands (Spain). Amphibia-Reptilia  2008; 29: 49– 60. Google Scholar CrossRef Search ADS   Pretus JL, Marqués R, Pérez-Mellado V. Holocene sea level rise 975 and range of fragmentation of Podarcis lilfordi on Minorcan islets: a vicariance scenario reviewed through a mtDNA tree. In Pérez-Mellado V, Riera N, Perera A, (eds). The Biology of lacertid lizards: Evolutionary and Ecological Perspectives . Institut Menorquí d'Estudis . 2004; 8: 279– 291. Price MN, Dehal PS, Arkin AP. FastTree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol Biol Evol  2009; 26: 1641– 50. Google Scholar CrossRef Search ADS PubMed  Rawls JF, Mahowald MA, Ley RE et al.   Reciprocal gut microbiota transplants from zebrafish and mice to germ-free recipients reveal host habitat selection. Cell  2006; 127: 423– 33. Google Scholar CrossRef Search ADS PubMed  Ren T, Kahrl AF, Wu M et al.   Does adaptive radiation of a host lineage promote ecological diversity of its bacterial communities? A test using gut microbiota of Anolis lizards. Mol Ecol  2016; 25: 4793– 804. Google Scholar CrossRef Search ADS PubMed  Saez ET, Traveset A. Fruit and nectar feeding by Podarcis lilfordi (Lacertidae) on Cabrera archipelago (Balearic islands). Herpetol Rev  1995; 26: 121– 3. Salvador A. Lagartija balear – Podarcis lilfordi (Günther, 1874). In Salvador A, Marco A., (eds). Museo Nacional de Ciencias Naturales, Madrid . Enciclopedia Virtual de los Vertebrados Españoles, 2009. http://www.vertebradosibericos.org. Sanders ED. Evidence Concerning Late Survival and Extinction of Endemic Amphibia and Reptilia from the Bronze and Iron Age settlement of Torralba d’En Salort (Alayor, Menorca) . Palma de Mallorca: Monografies Cientìfiques, Palma de Mallorca. 1984. Schloss PD, Westcott SL, Ryabin T et al.   Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microb  2009; 75: 7537– 41. Google Scholar CrossRef Search ADS   Shoemaker W, Locey KJ, Lennon JT. A macroecological theory of microbial biodiversity. Nat Ecol Evol  2017; 1: 107. Google Scholar CrossRef Search ADS PubMed  Song SJ, Lauber C, Costello EK et al.   Cohabiting family members share microbiota with one another and with their dogs. eLife  2013; 2: e00458. Google Scholar PubMed  Terrasa B, Perez-Mellado V, Brown RP et al.   Foundations for conservation of intraspecific genetic diversity revealed by analysis of phylogeographical structure in the endangered endemic lizard Podarcis lilfordi. Divers Distrib  2009; 15: 207– 21. Google Scholar CrossRef Search ADS   Traveset A, Riera N. Disruption of a Plant-Lizard seed dispersal system and its ecological effects on a threatened endemic plant in the Balearic islands. Conserv Biol  2005; 19: 421– 31. Google Scholar CrossRef Search ADS   Vecchi M, Marriner N, Morhange C et al.   Multiproxy assessment of Holocene relative sea-level changes in the western Mediterranean: Sea-level variability and improvements in the definition of the isostatic signal. Earth-Sci Rev  2016; 155: 172– 97. Google Scholar CrossRef Search ADS   Xiao L, Feng Q, Liang S et al.   A catalog of the mouse gut metagenome. Nat Biotechnol  2015; 33: 1103– 8. Google Scholar CrossRef Search ADS PubMed  Yatsunenko T, Rey FE, Manary MJ et al.   Human gut microbiome viewed across age and geography. Nature  2012; 486: 222– 7. Google Scholar PubMed  Zeng Q, Sukumaran J, Wu S et al.   Neutral models of microbiome evolution. PLoS Comput Biol  2015; 11: e1004365. Google Scholar CrossRef Search ADS PubMed  © FEMS 2017. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png FEMS Microbiology Ecology Oxford University Press

Processes shaping gut microbiota diversity in allopatric populations of the endemic lizard Podarcis lilfordi from Menorcan islets (Balearic Islands)

Loading next page...
 
/lp/ou_press/processes-shaping-gut-microbiota-diversity-in-allopatric-populations-eJ9MYq2Eo9
Publisher
Oxford University Press
Copyright
© FEMS 2017. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com
ISSN
0168-6496
eISSN
1574-6941
D.O.I.
10.1093/femsec/fix186
Publisher site
See Article on Publisher Site

Abstract

Abstract Compositional variation of the gut microbiota across host allopatric populations can reflect both adaptation and stochasticity since the time of separation. Major factors shaping this variation include the host phylogeographic and demographic history, the microbiota inheritance, environmental inputs and dispersal of bacteria. Here we explored the impact of these factors in driving gut community diversity in seven allopatric populations of the omnivorous lizard Podarcis lilfordi from the Menorcan coastal islets, all descending from an ancestral mainland population. Using 16S rRNA Illumina sequencing, we showed that ‘islet’ and ‘age’ (time since islet separation from mainland) were the only significant variables in microbial community clustering, suggesting a partial islet-restricted diversification following these lizards phylogeography. Despite a significant variation, islets/populations were characterized by a remarkably low bacterial uniqueness (2.4% of total OTUs) and a minor differential enrichment of taxa, indicating a negligible impact of local inputs and important host common constraints. Overall, the extant pattern of similarity/dissimilarity among islets is compatible with partial retention of the ancestral mainland microbial pool, with differences among islets potentially explained by a differential loss of bacteria following population fragmentation and bottlenecks (i.e. ecological drift). While more quantitative data are needed to validate this hypothesis, this study unveils the importance of considering both neutral and niche-driven processes in driving contemporary patterns of gut metacommunity diversity. island, drift, dispersal, inheritance, population, allopatry INTRODUCTION A multitude of factors shape the gut microbiota composition and challenge its stability over ecological and evolutionary times. Most critical among them are the host constraints (e.g. gut morphology, physiology and immune system) (Rawls et al.2006; Ley et al.2008; Benson et al.2010; Moeller et al.2016) and environmental inputs (e.g. diet and parasites) (Muegge et al.2011; David et al.2014; Delsuc et al.2014; Aivelo and Norberg 2016; Kohl et al.2017). Although diet and host genotype have been consistently found to be the most influential determinants, large part of animal intraspecific microbial variation remains unexplained by adaptive variables. Hence, the role of neutral processes (drift and dispersal) is gaining more attention (Costello et al.2012; Lankau, Hong and Mackie 2012; Nemergut et al.2013; Zeng et al.2015; Burns et al.2016; Martinson, Douglas and Jaenike 2017; Shoemaker, Locey and Lennon 2017). Population-level studies of gut microbiotas dynamics might represent the right scale to explore the influence of both neutral and selective forces on metacommunity structure. While substantial interindividual microbial variation might exist, populations do typically carry discrete gut metacommunities (Yatsunenko et al.2012; Song et al.2013), resulting from a restricted, although not exclusive circulation of gut bacteria favored by parental transmission (Funkhouser and Bordenstein 2013), preferential dispersal by shared ecological niche (e.g. common diet), reproductive strategy (including parental care) and social behavior (Nemergut et al.2013; Song et al.2013; Ren et al.2016). How gut metacommunities respond to restrictions in ‘free’ circulation among host members, for instance in the event of population split, remains to date poorly understood/investigated in natural systems (Lankau, Hong and Mackie 2012), owing to the concomitant effect of adaptive variables (ecological traits). Allopatric but ecologically similar populations found in islands might provide excellent natural systems to explore the impact of geographic barriers and host demographic history on gut microbiota diversity. In the event of fragmentation of an original insular population following isolation of some specimens into novel islets, distinct factors can act over time to shape the microbiota compositional structure (for simplicity, hereafter use of microbial/microbiota will refer to bacteria only). These include the conservation of taxa present in the original pool by stable transmission through host generations (i.e. microbiota inheritance) (Funkhouser and Bordenstein 2013; Zeng et al.2015), selective and stochastic changes in taxa relative abundances (the latter known as ecological drift) (Chase and Myers 2011; Lankau, Hong and Mackie 2012), acquisition of novel taxa from local pools (i.e. local exposure) and dispersal of taxa among populations (Nemergut et al.2013). Of these, ecological drift during and following population isolation and local exposure should drive microbiota divergence among populations (Lankau, Hong and Mackie 2012), whereas microbiota inheritance and dispersal of common taxa across islands should maintain/increase microbial similarity among populations (Zeng et al.2015). Disentangling the contribution of each of these processes on extant patterns of similarity/dissimilarity of gut microbiota would provide a better understanding of the compositional stability and dynamics of these microbial communities over spatial and temporal changes. We explored this scenario using allopatric populations of the lizard species Podarcis lilfordi inhabiting coastal islets around Menorca, one of the Balearic Islands (Fig. 1). The species originated during the re-flooding of the Mediterranean at the end of the Messinian salinity crisis, about 5.33 Mya (Krijgsman et al.1999), with a first cladogenic event (2.88–2.66 Ma) leading to the separation between P. lilfordi (currently present in coastal islets of Mallorca and Menorca, and Cabrera) and P. pityusensis (Ibiza, Formentera and adjacent coastal islets) (Brown et al.2008; Terrasa et al.2009). To date, endemic populations of P. lilfordi inhabit 16 of the 19 Menorca islets (Pérez-Mellado 1989; Terrasa et al.2009) (Fig. 1). According to mtDNA data, these Menorca populations were largely panmictic during the Pleistocene (Terrasa et al.2009); the following sea level rise during the past 10 000 years (Holocene) (Vecchi et al.2016) led to coastal fragmentation and progressive isolation of islet populations up to date (Pretus, Marques and Pérez-Mellado 2004; Brown et al.2008; Vecchi et al.2016). Archaeological data indicate that the mainland stock became extinct ∼2000 years ago subsequent to the introduction of competitive species and foreign predators (Alcover, Moyá-Solà and Pons-Moyá 1981; Sanders 1984; Pretus, Marques and Pérez-Mellado 2004), while islets provided valuable refugia. Extant populations have been geographically and reproductively isolated for the past 2000 years, although time since separation most likely ranges between 4000 and 9000 years, according to bathymetry with respect to mainland and Holocene sea level curves (Vecchi et al.2016). While morphologically diverse, suggesting incipient speciation (Salvador 2009), mtDNA data show low genetic differentiation among P. lilfordi populations, in accordance with their recent segregation (Terrasa et al.2009). The presence of additional vertebrates (birds, rabbits and humans) is restricted to few islets (Pérez-Mellado 1989; Salvador 2009). Several studies indicate that a dietary shift occurred with variable intensity in the different populations, from an ancestral insectivorous niche to a current omnivorous niche, including the uptake of a large proportion of arthropods and vegetable matter, with seasonal variations (Saez and Traveset 1995; Traveset and Riera 2005; Perez-Cembranos, Leon and Perez-Mellado 2016). Figure 1. View largeDownload slide Map of Menorca and sampled islets (circled). Population size according to Pérez-Mellado et al. (2008) is provided below each islet (no estimates available for Ravl, although expected population size is likely below 1000 individuals, personal observations). Bars below each islet correspond to a 50-mt scale. Islet abbreviations: AdGr (Addaia Gran), Aire (Aire), Colm (Colom), PrCv (Porros de Cavalleria), Ravl (Ravells), Rei (Rei) and Tosq (Tosqueta). Figure 1. View largeDownload slide Map of Menorca and sampled islets (circled). Population size according to Pérez-Mellado et al. (2008) is provided below each islet (no estimates available for Ravl, although expected population size is likely below 1000 individuals, personal observations). Bars below each islet correspond to a 50-mt scale. Islet abbreviations: AdGr (Addaia Gran), Aire (Aire), Colm (Colom), PrCv (Porros de Cavalleria), Ravl (Ravells), Rei (Rei) and Tosq (Tosqueta). Overall, two aspects make these populations a rare system for studying natural processes shaping the vertebrate gut microbiota: their genetic and trophic similarity, which allowed us to detract the differential effect of these two major factors affecting the microbiota, and their distribution into several allopatric populations descending from a single ancestral population. Here we aimed to address some key questions: Was the geographic isolation and subsequent exposure to distinct local pools of bacteria sufficient to drive significant divergence among enteric communities in these populations? Do their gut microbial communities still carry an identifiable signature of their past history, i.e. a core component that has been putatively inherited from the ancestral mainland population? What is the actual impact of ecological drift and dispersal in explaining the current microbiota variation? To explore these questions, we sampled whole intestines from P. lilfordi specimens from seven Menorca islets and performed comparative analyses of their gut microbial community structure. We specifically partitioned their gut microbiota into shared and variable components and investigated forces responsible for shaping contemporary patterns of similarity/dissimilarity. MATERIALS AND METHODS Sample collection and morphometric body measures A total of 72 individuals from 7 of the 16 populated islets surrounding the major Menorca island (3 to 18 individuals per population/islet) were sampled in summer 2015, between 27 July and 1 August (Fig. 1 and Table S1, Supporting Information). Summer represents the season when trophic resources are most limited for these lizard populations. Sampled islets were chosen to vary along a gradient of bathymetry (from deep to shallow channels between islets and Menorca) and distance to the coastline (together providing a proxy for time since population isolation), and for their different exposure to the open sea (including both islets found in open sea and within protected bay areas). Specimens were collected using food traps containing fruits or fruit juices. For each specimen, the following morphometric measurements were taken: number of ventral-row scales (VS), weight (W), snout-vent length (SVL), head width (HW) and head length (HL). Of the 72 specimens, 33 adult males were chosen for microbiota sequencing (five for population, except for Rei, for which three specimens were sampled), while the remaining specimens were immediately released (Table S1, Supporting Information). Adulthood stage was assigned when SVL ≥ 51 mm for females and SVL ≥ 55 mm for males, according to the mean values for the smallest subspecies as described in Salvador (2009). Sex was assigned primarily according to number of VS (<27 for males) and validated by body size and femoral pores (males are larger than females and show pores with visible lipophilic compounds (Salvador 2009) (Table S1, Supporting Information). Statistical correlations were assessed through the cor.test function in the R Stats package version 3.2.2. Specimen sampling and manipulation were carried out in accordance with the ethics guidelines and recommendations of the Species Protection Service (Department of Agriculture, Environment and Territory, Government of the Balearic Islands), under permit No. CEP 31/2015 given to LB and JLP. Microbiota DNA extractions and 16S rRNA Illumina sequencing Sampled specimens for the microbiota study were kept in individual plastic containers for up to 6 h, and then transferred to –20°C until death. The full digestive tracts were dissected with sterile instruments and preserved in RNAlater. Samples were immediately stored at 4°C for 24 h and then kept at –20°C until further processing. Sex was then carefully reassessed, revealing the inclusion of three adult females in the final sample set (Table S1, Supporting Information). Under a sterile hood, the entire digestive tracts were wiped from RNAlater using a Kimwipes filter paper (Sigma-Aldrich) and lengths measured with a Vernier Caliper (Table S1, Supporting Information). Stomachs and intestines were separated according to their different tissue characteristics (stomachs had visibly ticker walls). Stomach contents were removed for examination under a stereoscope. The intact intestines were transferred to tubes with a 800-μL lysis buffer (see Baldo et al.2015), chopped using sterile scissors for facilitating the lysis process, vigorously shacked and incubated at 70°C for 30 min with gentle agitation. After addition of 0.3 g of 0.1-mm zirconia-silica beads, tube contents were homogenized with a Precellys Evolution bead beating instrument (Bertin Technologies) at 5500 rpm for 2 × 45 s, and centrifuged at 4°C for 5 min at 16 000 × g. Supernatant was stored in a separate tube. Additional 300 μL of lysis buffer was added to the pellet, and samples were homogenized again in the Precellys instrument for 2 × 30 s. Stored supernatant was placed back into the tube with beads, 20 μL of proteinase K (to a final of 10 mg/mL) was added and samples incubated at 55°C for 1 h at 450 rpm. Extracted DNA was purified according to Baldo et al. (2015). Briefly, residues were eliminated with ammonium acetate and DNA precipitated with isopropanol and cleaned with standard ethanol washes followed by Qiagen columns (DNeasy Blood and Tissue kit). The region V3-V4 (∼460 bp) of 16S rRNA was amplified using the following non-barcoded primers including TruthSeq adapters: S-D-Bact-0341-b-S-17 (5΄-CTACACGACGCTCTTCCGATCT CCTACGGGNGGCWGCAG-3΄) and S-D-Bact-0785-a-A-21 (5΄-CAGACGTGT GCTCTTCCGATCT GACTACHVGGGTATCTAATCC -3΄). Amplicons from individual specimens were barcoded and pooled according to Baldo et al. (2017). The final library was sequenced on Illumina MiSeq v3 (600-cycle cartridge, 300 paired-end reads) at the Centre for Genomic Regulation in Barcelona (Spain). Sequence analyses After demultiplexing, paired-end reads were merged and quality-filtered using the scripts make.contigs and screen.seqs (maxambig = 0, maxhomop = 8, maxlength = 490) implemented in mothur v.1.25.0 (Schloss et al.2009). Chimeras were identified using uchime_ref in vsearch v1.4.4 (VSEARCH GitHub repository) against the Chimera Slayer reference database (‘gold’) in the Broad Microbiome Utilities and discarded. Sequences were input into the QIIME pipeline (Caporaso et al.2010), using Macqiime 1.9.1., and clustered into operational taxonomic units (OTUs) at 97% similarity through the command pick_open_reference_otus.py. Taxonomy was assigned with the UCLUST method against the Greengenes gg_13_8_otus with RDP classifier (80% confidence). Mitochondrial and chloroplast-derived OTUs were discarded, retaining only bacterial OTUs with ≥10 total sequences. Sequences were aligned with PyNAST default parameters (Caporaso et al.2010) using the reference database ‘core_set_aligned.fasta’ from Greengenes website and default lanemask. Alignments were further trimmed to ensure overlapping ends, and a phylogenetic tree was built with FastTree (Price, Dehal and Arkin 2009). Alpha diversity measures (phylogenetic diversity (PD) whole tree and Shannon Index) were estimated on multiple rarefied biom tables to an even depth of 38 160 reads (10 iterations). Statistical differences in alpha diversity across islets were calculated with a non-parametric two-sample t-test (999 permutations). Beta diversity and PERMANOVA analyses Unweighted Unifrac and binary Jaccard distance matrices were built on the rarefied OTU table (38 160 sequences). Statistical correlations between microbiota distances and categorical variables (‘islet’, ‘age’, presence of ‘rats’ and stable colonies of ‘seagulls’, and ‘exposure’ to open sea, see Table S1, Supporting Information) were assessed by permutational multivariate analysis of variance (PERMANOVA, Anderson 2001; Anderson, Gorley and Clarke 2008) for PRIMER v7 (Clarke and Gorley 2015). All factors were considered fixed save for ‘islet’, which was defined as a random factor and included in all designs as nested within the levels of fixed factors. Categorical ‘age’ provides an approximation of each islet isolation time and was largely assigned according to bathymetry (as depth of the shallowest channel between islets and main island) choosing a threshold of ≥9 m to discriminate between old and more recently separated islets, in accordance with the speed of sea level rising, which slowed down after 7500 yr BP (Vecchi et al.2016) (Table S1, Supporting Information). Bathymetry and distance from mainland were also combined to define different age categories, grouping Aire and PrCv (old) against all remaining islets (recent) (see the Results section). To summarize microbial metacommunity diversity by islet, OTU counts, with and without previous rarefaction, were binned by islet and the data were used to calculate ‘uniqueness’ by islet (i.e. taxa and OTUs occurring only in one islet to the exclusion of all other islets). To test whether microbiota distances correlated with the geographic distances among islets, the original OTU table was filtered for OTUs above 50 total counts, reads were binned by islet and rarefied to the minimum depth per islet. Jaccard and unweighted Unifrac distance matrices were built on this OTU table and tested against geographic (i.e. linear) distances across islets with a Mantel test. Islet distances from mainland were instead approximated by the ‘age’ category and tested with PERMANOVA. Taxa and OTUs enrichment among islets To statistically compare differences in taxa and OTU relative abundances among islets, we used the script group_significance.py in QIIME with the non-parametric Kruskal–Wallis test. In all cases, a false discovery rate (FDR) correction was applied to assess significant differences. Differentially enriched OTUs (i.e. indicator OTUs) were estimated on the original OTU table rarefied to 38 160 reads. For differentially enriched taxa (i.e. indicator taxa), OTU counts were binned to each taxon level. Results were filtered for OTU/taxa with total sum of means across islets greater than 10. Core analyses We used the non-rarefied OTU table to estimate core OTUs and taxa (from phylum to species, estimated by binning OTU counts at each taxon level, as above) at different depths (shared by 100%, 90% and 80% of total specimens), retaining those taxa and OTUs with representation in all islets. Core OTUs were searched with batch BLASTN against the local nt database (cut-off e-value = 30), retrieving the first two best hits. Outputs were parsed using in-house perl scripts. Metadata associated with each BLASTN hit was extracted from Genbank records with Biopython scripts. The 134 core OTUs present in at least 90% of the specimens were used to explore a potential signature of within-population diversification. For this purpose, all reads associated with core OTUs were pooled and reclustered at 98% and 99% of similarity (now referred as to 98% and 99% OTUs or more generally as to sub-OTUs), picking the most abundant as representative and retaining OTUs with at least 50 reads total and presence in at least two specimens. A binary Jaccard distance matrix was built on each sub-OTU abundance table, rarefied to the minimum depth coverage across samples (13 600 and 11 000 for 98% and 99% OTUs, respectively), and significance of clustering was tested with PERMANOVA, as for 97% OTUs (see above). The 23 core OTUs that occurred in all specimens (representing the 100% core) were analyzed individually for their sub-OTU contribution to the islet-based structuring. Specifically, reads from each OTU were reclustered at 99% similarity and significance tested with permutational analysis of variance using adonis (999 permutations) with the script compare_categories.py in QIIME. RESULTS Morphometrics and stomach content The large majority of the captured specimens were adult males, putatively because they were less shy in approaching food sources, although a skewed sex ratio cannot be ruled out (Table S1, Supporting Information). Males were on average significantly larger than females (weight, W was 8.2 g ± 2.17 for males, 5.4 g ± 1.55 for females) and had longer SVL (65.2 mm ± 4.20 for males, 58.6 mm ± 3.68 for females) and larger HW/SVL and HL/SVL values (Table S1, Supporting Information). Considering males only, Addaia Gran (AdGr) hosted the smallest lizards, which showed among the largest HL and HW/SVL ratios, while Aire hosted the biggest and longest specimens. These results are compatible with previous findings (Salvador 2009). Overall, SVL positively correlated with W (Pearson's r = 0.677, P < 0.0001), HW (r = 0.6735581, P < 0.0001) and HL (r = 0.863, P < 0.0001), but not with intestine length (r = 0.132, P = 0.333). Stomach content analyses indicated widespread and abundant presence of arthropods, particularly insects (especially ants) in most specimens and all islets, compatible with a dominant role of insectivory during the dry season. Sequence summary We characterized the gut microbiota of 33 total specimens (30 males and 3 females), with three to five representatives per population (Table S1, Supporting Information). After quality filtering, removal of the three OTUs found in the blank control (water only) and OTUs with less than 10 reads total, we obtained a total of 2 014 276 high-quality reads and an average of 60 081 reads per specimen (ranging between 38 163 and 106 206), corresponding to 6380 OTUs. These were assigned to 17 phyla, 32 classes (6378 annotated OTUs), 58 orders (6369), 99 families (4975) and 147 genera (2321). Rarefaction curves based on the number of observed OTUs indicated that we captured most of the sample diversity, although we did not reach saturation at the chosen rarefied depth of 38 160 sequences (Fig. S1, Supporting Information). The taxonomic composition of the P. lilfordi gut microbiota is highly conserved Podarcis lilfordi showed the typical vertebrate phyla composition, with three phyla representing altogether 97% of the microbiota: Firmicutes (average 65% across islets), Bacteroidetes (18%) and Proteobacteria (14%) (Fig. 2). The proportion of Firmicutes was highly comparable across islets, while the representation of both Bacteroidetes and Proteobacteria varied to a greater extent. Firmicutes and Bacteroidetes were represented by six and five orders, respectively, while being dominated by a single order in all islets: Clostridiales (Firmicutes) and Bacteroidales (Bacteroidetes) (Fig. 2). Two families within both phyla showed the highest relative abundances: Lachnospiraceae and Ruminococcaceae within Firmicutes, and Porphyromonadaceae and Bacteroidaceae within Bacteroidetes (Fig. 2). Patterns of family relative abundance within islet were remarkably conserved, especially for Firmicutes. On the other hand, Proteobacteria were more diverse, including a larger number of represented orders (24) and families (33). Among the most represented taxa within this phylum were four orders (Desulfobrionales, Enterobacteriales, Campylobacteriales and Pseudomonadales) and three families (Helicobacteraceae, Desulfovibrionaceae and Enterobacteriaceae), present in all islets (Fig. 2). Proteobacteria varied widely in their quantitative profiles across islets at all taxonomic levels. Nonetheless, alpha diversity estimates for this phylum were not significantly different among islets (two-sample t-test P > 0.05, both Shannon and PD whole tree for all islet pairwise). Figure 2. View largeDownload slide Order and family-level taxonomic composition of the three major phyla, Firmicutes, Bacteroidetes and Proteobacteria, as proportions of OTUs. Values were estimated on the rarefied dataset and averaged by islet. Legends list only the 10 most abundant taxa. Figure 2. View largeDownload slide Order and family-level taxonomic composition of the three major phyla, Firmicutes, Bacteroidetes and Proteobacteria, as proportions of OTUs. Values were estimated on the rarefied dataset and averaged by islet. Legends list only the 10 most abundant taxa. As a whole, gut communities did not show significantly distinct taxonomic diversity and richness across islets, according to both Shannon and PD-whole tree (Fig. S2, Supplementary; P > 0.05 at both indexes for any islet pairwise), although a visual trend indicates Porros de Cavalleria (PrCv) as the islet hosting the largest diversity. Diversity and richness did not correlate with intestinal length (Pearson's r = –0.092, P = 0.49 for PD whole tree and –0.194, P = 0.15 for Shannon). Overall, the major taxonomic structure of the P. lilfordi gut microbiota was remarkably conserved across islets, with Proteobacteria content representing the most variable quantitative trait. Islet as the major determinant in the P. lilfordi microbiota clustering According to both unweighted Unifrac and Jaccard distances on 97% OTUs, ‘islet’ was the main significant factor accounting for 13.5% of the community variation, as measured by partition into components of variation (PERMANOVA, P = 0.0001, Table 1). There were no significant effects of age (but see below), exposure or presence of rats (P > 0.05). The presence of seagull colonies was barely significant, accounting for only 3.5% of microbiota variation. The islet Rei was grouped with the oldest islets according to bathymetry, but it presents several unique features, including its location within a natural harbor at short distance from the mainland and in a context of a largely anthropized environment (the islet hosted a paleochristian basilica in the 6th century and a hospital in the 19th century). When both bathymetry and distance from the mainland were used as grouping factors, thus excluding Rei as an ‘ancient’ islet and retaining only Aire and PrCv, ‘age’ became a significant factor (P = 0.0001, Table 1). The only interaction that could be tested (age by exposure) was not significant. Microbiota distances across islets did not correlate with geographic distances across islets (Mantel test, r = 0.00588, P = 0.980 for Jaccard and r = –0.02849, P = 0.906 for Unifrac), but perhaps indirectly with their distances from the mainland, here approximated by the ‘age’ category (see the Methods section). The ‘islet’ and ‘age’ effects are also illustrated by a principal coordinate analysis based on unweighted Unifrac distances, depicting the specimens clustering according to the islet origin (color-coded, Fig. 3). The two oldest islets, Aire and PrCv, clearly diverged from the remaining islets along PC2. Figure 3. View largeDownload slide Principal coordinate analysis according to unweighted Unifrac distances. Distances were calculated for OTUs at 97% similarity using the rarefied dataset (38 160 sequences). Dots correspond to single specimens. Specimens within each islet are connected to the median coordinates for each islet. Figure 3. View largeDownload slide Principal coordinate analysis according to unweighted Unifrac distances. Distances were calculated for OTUs at 97% similarity using the rarefied dataset (38 160 sequences). Dots correspond to single specimens. Specimens within each islet are connected to the median coordinates for each islet. Table 1. Results of PERMANOVA analyses, separately for the whole dataset and core OTUs after reclustering at 98 and 99% similarity.     Whole dataset  Core (90% of specimens)      Unweighed Unifrac  Jaccard 97%  Jaccard sub-OTU 98%  Jaccard sub-OTU 99%  Model  Term  P-value  Variation  P-value  Variation  P-value  Variation  P-value  Variation        estimate    estimate    estimate    estimate  Islet  Islet  0.0001  0.028  0.0001  0.035  0.0001  0.029  0.0001  0.031  Age + islet(A)  Age  N.S. [0.0001]  0.004 [0.011]  N.S. [0.0001]  0.006 [0.012]  N.S. [0.0001]  0.005 [0.012]  N.S.  0.003 [0.008]    Islet(A)  0.0001  0.024 [0.022]  0.0001  0.029 [0.028]  0.0001  0.024 [0.023]  0.0001  0.028 [0.027]  Exposure + islet(E)  Exposure  N.S.  0.000  N.S.  0.000  N.S.  0.000  N.S.  0.000    Islet(E)  0.0001  0.028  0.0001  0.034  0.0001  0.029  0.0001  0.031  Age + exposure + islet(A, E)  Age  N.S.  0.002 [0.010]  N.S.  0.005 [0.012]  N.S. [0.0286]  0.002 [0.012]  N.S.  0.001 [0.008]    Exposure  N.S.  0.003 [0.000]  N.S.  0.003 [0.000]  N.S.  0.003 [0.001]  N.S.  0.002 [0.001]    A*E  N.S.  0.003  N.S.  0.002  N.S.  0.003  N.S.  0.003    Islet(A, E)  0.0001  0.022  0.0001  0.027  0.0001 [N.T.]  0.022  0.0001  0.026  Rats + islet(R)  Rats  N.S.  0.001  N.S.  0.001  N.S.  0.002  N.S.  0.002    Islet(R)  0.0001  0.027  0.0001  0.034  0.0001  0.028  0.0001  0.030  Seagulls + islet(S)  Seagulls  0.0465  0.007  0.0388  0.008  0.0299  0.008  0.0340  0.007    Islet(S)  0.0001  0.024  0.0001  0.030  0.0001  0.025  0.0001  0.027    Residual    0.178    0.201    0.161    0.175      Whole dataset  Core (90% of specimens)      Unweighed Unifrac  Jaccard 97%  Jaccard sub-OTU 98%  Jaccard sub-OTU 99%  Model  Term  P-value  Variation  P-value  Variation  P-value  Variation  P-value  Variation        estimate    estimate    estimate    estimate  Islet  Islet  0.0001  0.028  0.0001  0.035  0.0001  0.029  0.0001  0.031  Age + islet(A)  Age  N.S. [0.0001]  0.004 [0.011]  N.S. [0.0001]  0.006 [0.012]  N.S. [0.0001]  0.005 [0.012]  N.S.  0.003 [0.008]    Islet(A)  0.0001  0.024 [0.022]  0.0001  0.029 [0.028]  0.0001  0.024 [0.023]  0.0001  0.028 [0.027]  Exposure + islet(E)  Exposure  N.S.  0.000  N.S.  0.000  N.S.  0.000  N.S.  0.000    Islet(E)  0.0001  0.028  0.0001  0.034  0.0001  0.029  0.0001  0.031  Age + exposure + islet(A, E)  Age  N.S.  0.002 [0.010]  N.S.  0.005 [0.012]  N.S. [0.0286]  0.002 [0.012]  N.S.  0.001 [0.008]    Exposure  N.S.  0.003 [0.000]  N.S.  0.003 [0.000]  N.S.  0.003 [0.001]  N.S.  0.002 [0.001]    A*E  N.S.  0.003  N.S.  0.002  N.S.  0.003  N.S.  0.003    Islet(A, E)  0.0001  0.022  0.0001  0.027  0.0001 [N.T.]  0.022  0.0001  0.026  Rats + islet(R)  Rats  N.S.  0.001  N.S.  0.001  N.S.  0.002  N.S.  0.002    Islet(R)  0.0001  0.027  0.0001  0.034  0.0001  0.028  0.0001  0.030  Seagulls + islet(S)  Seagulls  0.0465  0.007  0.0388  0.008  0.0299  0.008  0.0340  0.007    Islet(S)  0.0001  0.024  0.0001  0.030  0.0001  0.025  0.0001  0.027    Residual    0.178    0.201    0.161    0.175  N.S.: not significant (P > 0.05). N.T.: no test performed. Rarefaction depth used: whole dataset (38 160 sequences); core sub-OTU 98% (13 600), sub-OTU 99% (11 000). In brackets are estimates for ‘age’ when excluding Rei from the ‘old’ category (see Table S1, Supporting Information). Residual are common to all models. View Large Despite a main significant islet effect, a large fraction of the P. lilfordi gut microbiota was common to all populations: according to Jaccard pairwise distances across islets, mean distance between any two islets was only 0.38, with the Aire and PrCv representing the closest pair (mean Jaccard = 0.28). Hereafter, we investigated both shared and variable microbiota components driving the above pattern. Low microbiota uniqueness and differential enrichment among islets To explore the variable component responsible for the observed clustering and distances, we considered both taxa/OTU that are islet-specific in terms of presence/absence (i.e. uniqueness) and those that are differentially enriched. Taxa and OTUs unique to one islet can either represent novel gut acquisitions from the local pool, contamination by allochthonous bacteria, or simply loss of these taxa in all other islets. Total number of OTUs per islet, without rarefaction, ranged between 3013 and 4816 for a comparable sampling effort (five individuals per islet, except for Rei, with three individuals sampled) (Table 2). A total of 1093 OTUs (17% of total OTUs) were shared among all islets (although not all specimens). Uniqueness by islet was remarkably low, accounting for only between 1.1% and 5.2% of an islet's OTU content (average 2.4%) and showing limited distribution (1–2 individuals). This suggests that we either failed to sample them in most specimens or, more likely, that they represent part of the allochthonous component. Of the 666 total unique OTUs detected (Table 2), none corresponded to unique taxa when checked against the unfiltered OTU table including all reads. In all cases, most unique OTUs occurred in Aire, with Tosqueta (Tosq) hosting the lowest uniqueness. Overall, at least 95% of each islet OTU content was shared between two or more islets. The same analyses performed on the rarefied OTU table returned highly comparable results. Table 2. OTU uniqueness and shared component.   Number of OTUs      Islet  Totala  Unique (%)  Shared %b    Tosq  3287  37 (1.13)  98.87    Rei  3013  61 (2.02)  97.98    Ravl  3370  71 (2.11)  97.89    Colm  4334  82 (1.89)  98.11    PrCv  4816  87 (1.81)  98.19    AdGr  4133  110 (2.66)  97.34    Aire  4219  218 (5.17)  94.83    Average  3881.7  95.1 (2.40)  97.60      Total  Unique  Sharedb  Commonc  90% of specimens  All islets  6380  666  5714  1093  134    Number of OTUs      Islet  Totala  Unique (%)  Shared %b    Tosq  3287  37 (1.13)  98.87    Rei  3013  61 (2.02)  97.98    Ravl  3370  71 (2.11)  97.89    Colm  4334  82 (1.89)  98.11    PrCv  4816  87 (1.81)  98.19    AdGr  4133  110 (2.66)  97.34    Aire  4219  218 (5.17)  94.83    Average  3881.7  95.1 (2.40)  97.60      Total  Unique  Sharedb  Commonc  90% of specimens  All islets  6380  666  5714  1093  134  aof OTUs was estimated from a non-rarefied dataset, considering OTUs > 10 total sequences. bby two or more islets. cin all islets, but not all specimens. View Large As for presence/absence, differences in relative abundances of taxa/OTUs among islets were also minor. According to enrichment analyses, only a few taxa (Table S2, Supporting Information) and OTUs (total of 30, Table 3) showed a significant differential representation among islets (kw, FDR P < 0.05). Among indicator taxa (Table S2, Supporting Information), the phylum Synergistetes and its only recognized family Synergistaceae were a signature of Aire and PrCv, with only a few reads detected in AdGr and Colom (Colm). Aire had a particularly high representation of this family in terms of reads (376) and number of specimens (all). The islet of Rei showed a particular enrichment in Oceanospirillales and Halomonadaceae, while Ravells (Ravl) was enriched in the genus Turicibacter. Table 3. Differentially enriched OTUs among islets (Kruskal Wallis, FDR P < 0.05).   Mean of reads    Finest taxonomy  BLASTN Hits  OTU_IDa  Aire  AdGr  Rei  Tosq  Ravl  Colm  PrCv  Specificity  Phylum  Order  Family  Genus/species  Accession number  Identity (%)  Hostb  N.Ref.OTU895  516.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Peptococcaceae    GQ502520  91.8  Termite  721569  26.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    CU926622  98.9    175037  22.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Christensenellaceae    FJ362815  97.7  Human  NC.Ref.OTU185722  20.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      FJ535004  96.8    NC.Ref.OTU27689  20.6  0  0  0  0  0  0  Aire  Bacteroidetes  Bacteroidales  Bacteroidaceae  Bacteroides  HQ745181  96.7  Human  NC.Ref.OTU224538  18.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Peptococcaceae    GQ502520  90.8  Termite  NC.Ref.OTU27099  17.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      LC028698  97.3  Buffalo  NC.Ref.OTU212111  15  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    GQ294881  96.4  Mouse  566434  14  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      FJ535004  97.1    NC.Ref.OTU43839  10.4  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    AB218341  92.3  Dugong  NC.Ref.OTU59527  0  0  0  0  0  0  10.6  PrCv  Firmicutes  Clostridiales      AF371819  96.8  Swine  688934  0  0  0  0  25.2  0  0  Ravl  Proteobacteria  Enterobacteriales  Enterobacteriaceae  Serratia  KY124195  98.3    NC.Ref.OTU106014  0  0  0  0  22.8  0  0  Ravl  Proteobacteria  Enterobacteriales  Enterobacteriaceae    FR848151  98.1  Millipedes  182970  128.6  0  0  0  0  0  664.2  Aire, PrCv  Verrucomicrobia  Verrucomicrobiales  Verrucomicrobiaceae  Akkermansia  LT629973  99.1  Python  N.Ref.OTU217  46.4  0  0  0  0  0  66.2  Aire, PrCv  Firmicutes  Clostridiales      KU506346  94.4  Pig  N.Ref.OTU651  28  0  0  0  0  0  73.4  Aire, PrCv  Bacteroidetes  Bacteroidales  [Barnesiellaceae]    HQ792288  91.5  Human  NC.Ref.OTU124139  108  0  0  0  0  0  6  Aire, PrCv  Firmicutes  Clostridiales      AB606338  95.7  Mouse  NC.Ref.OTU110653  67.6  0  0  0  0  0  0.4  Aire, PrCv  Firmicutes  Clostridiales      JQ248100  97.3  Chicken  NC.Ref.OTU229167  2.2  0  0  0  0  0  16.8  Aire, PrCv  Firmicutes  Clostridiales  Ruminococcaceae    HM123952  95.5  Human  3970651  1041.6  0  0  0  0  23  0  Aire, Colm  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  JX457157  99.3  Cockroach  NC.Ref.OTU12308  15.6  0  0  0  0  0.2  0  Aire, Colm  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  JX457157  96.7  Cockroach  NC.Ref.OTU231957  7.4  0  12  0  0  0  0  Aire, Rei  Firmicutes  Clostridiales  Ruminococcaceae    CU919465  94.4    NC.Ref.OTU25986  0  0  0  0  100.4  0  3.2  Ravl, PrCv  Firmicutes  Clostridiales      AB606272  94.3  Mouse  347529  25  6.2  0  0  329.2  5.2  26.6  Several  Firmicutes  Turicibacterales  Turicibacteraceae  Turicibacter  LT624967  99.1    727945  8.4  2.6  0  0  11  1.4  2.2  Several  Firmicutes  Clostridiales  Eubacteriaceae  Anaerofustis  LT223622  99.3  Human  NC.Ref.OTU117227  1.4  0  0  0  0.2  0  16.4  Several  Firmicutes  Clostridiales  Lachnospiraceae  Coprococcus  EU779145  96.6  Sheep  NC.Ref.OTU63406  7  42.2  0  0  0  0  54.2  Several  Bacteroidetes  Bacteroidales  Bacteroidaceae  Bacteroides  GU105959  92.8  Human  N.Ref.OTU255  10  0  0  0.2  0  14.2  7.2  Several  Firmicutes  Clostridiales  [Mogibacteriaceae]  Anaerovorax  JQ084279  94.8  Mouse  N.Ref.OTU41  274.4  134  3.7  0  0  0  86.2  Several  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  HQ767608  94.3  Human  N.Ref.OTU544  251.6  0.2  46.7  0  0  0  51.6  Several  Bacteroidetes  Bacteroidales  Porphyromonadaceae  P. distasonis  HQ797049  97.4  Human    Mean of reads    Finest taxonomy  BLASTN Hits  OTU_IDa  Aire  AdGr  Rei  Tosq  Ravl  Colm  PrCv  Specificity  Phylum  Order  Family  Genus/species  Accession number  Identity (%)  Hostb  N.Ref.OTU895  516.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Peptococcaceae    GQ502520  91.8  Termite  721569  26.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    CU926622  98.9    175037  22.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Christensenellaceae    FJ362815  97.7  Human  NC.Ref.OTU185722  20.8  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      FJ535004  96.8    NC.Ref.OTU27689  20.6  0  0  0  0  0  0  Aire  Bacteroidetes  Bacteroidales  Bacteroidaceae  Bacteroides  HQ745181  96.7  Human  NC.Ref.OTU224538  18.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Peptococcaceae    GQ502520  90.8  Termite  NC.Ref.OTU27099  17.2  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      LC028698  97.3  Buffalo  NC.Ref.OTU212111  15  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    GQ294881  96.4  Mouse  566434  14  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales      FJ535004  97.1    NC.Ref.OTU43839  10.4  0  0  0  0  0  0  Aire  Firmicutes  Clostridiales  Ruminococcaceae    AB218341  92.3  Dugong  NC.Ref.OTU59527  0  0  0  0  0  0  10.6  PrCv  Firmicutes  Clostridiales      AF371819  96.8  Swine  688934  0  0  0  0  25.2  0  0  Ravl  Proteobacteria  Enterobacteriales  Enterobacteriaceae  Serratia  KY124195  98.3    NC.Ref.OTU106014  0  0  0  0  22.8  0  0  Ravl  Proteobacteria  Enterobacteriales  Enterobacteriaceae    FR848151  98.1  Millipedes  182970  128.6  0  0  0  0  0  664.2  Aire, PrCv  Verrucomicrobia  Verrucomicrobiales  Verrucomicrobiaceae  Akkermansia  LT629973  99.1  Python  N.Ref.OTU217  46.4  0  0  0  0  0  66.2  Aire, PrCv  Firmicutes  Clostridiales      KU506346  94.4  Pig  N.Ref.OTU651  28  0  0  0  0  0  73.4  Aire, PrCv  Bacteroidetes  Bacteroidales  [Barnesiellaceae]    HQ792288  91.5  Human  NC.Ref.OTU124139  108  0  0  0  0  0  6  Aire, PrCv  Firmicutes  Clostridiales      AB606338  95.7  Mouse  NC.Ref.OTU110653  67.6  0  0  0  0  0  0.4  Aire, PrCv  Firmicutes  Clostridiales      JQ248100  97.3  Chicken  NC.Ref.OTU229167  2.2  0  0  0  0  0  16.8  Aire, PrCv  Firmicutes  Clostridiales  Ruminococcaceae    HM123952  95.5  Human  3970651  1041.6  0  0  0  0  23  0  Aire, Colm  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  JX457157  99.3  Cockroach  NC.Ref.OTU12308  15.6  0  0  0  0  0.2  0  Aire, Colm  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  JX457157  96.7  Cockroach  NC.Ref.OTU231957  7.4  0  12  0  0  0  0  Aire, Rei  Firmicutes  Clostridiales  Ruminococcaceae    CU919465  94.4    NC.Ref.OTU25986  0  0  0  0  100.4  0  3.2  Ravl, PrCv  Firmicutes  Clostridiales      AB606272  94.3  Mouse  347529  25  6.2  0  0  329.2  5.2  26.6  Several  Firmicutes  Turicibacterales  Turicibacteraceae  Turicibacter  LT624967  99.1    727945  8.4  2.6  0  0  11  1.4  2.2  Several  Firmicutes  Clostridiales  Eubacteriaceae  Anaerofustis  LT223622  99.3  Human  NC.Ref.OTU117227  1.4  0  0  0  0.2  0  16.4  Several  Firmicutes  Clostridiales  Lachnospiraceae  Coprococcus  EU779145  96.6  Sheep  NC.Ref.OTU63406  7  42.2  0  0  0  0  54.2  Several  Bacteroidetes  Bacteroidales  Bacteroidaceae  Bacteroides  GU105959  92.8  Human  N.Ref.OTU255  10  0  0  0.2  0  14.2  7.2  Several  Firmicutes  Clostridiales  [Mogibacteriaceae]  Anaerovorax  JQ084279  94.8  Mouse  N.Ref.OTU41  274.4  134  3.7  0  0  0  86.2  Several  Bacteroidetes  Bacteroidales  Porphyromonadaceae  Parabacteroides  HQ767608  94.3  Human  N.Ref.OTU544  251.6  0.2  46.7  0  0  0  51.6  Several  Bacteroidetes  Bacteroidales  Porphyromonadaceae  P. distasonis  HQ797049  97.4  Human  Taxonomy in brackets is not officially recognized (under revision). aOTUs were filtered for total sum of means > 10. bPresence in gut or feces samples. Blanks correspond to putative environmental hits. View Large Among indicator OTUs (Table 3), the majority belonged to the order Clostridiales (Firmicutes). Except for six putative environmental OTUs, all others had best hits to gut/feces of animals, typically vertebrates. Of them, 13 were single islet-specific, 10 being uniquely found in Aire. One Aire-specific OTU (N.Ref.OTU895) was present in high abundance (516 mean reads) in most specimens and had best BLASTN hit with low identity (91.8%) to the termite gut. Six OTUs were specific to Aire and PrCv, three of them with large mean read abundance in both islets (Table 3). OTU-182970, corresponding to a member of the genus Akkermansia, was particularly abundant in PrCv and had the best BLASTN hit to the reticulated python gut. This is overall the variable component responsible of the islet-based compositional structure of the P. lilfordi microbiota. Hereafter, we focused on the large shared component, the core microbiota. Core composition of the P. lilfordi gut microbiota We considered core taxa those with presence in at least 90% of specimens and all islets (calculated on the non-rarefied dataset). The core comprised 7 phyla, 12 classes, 13 orders, 25 families, 27 genera and 7 species (for species we considered presence in at least 80% of specimens) (Table 4). Notably, the only represented order within Cyanobacteria was YS2, a recognized member of the novel non-photosynthetic phylum Melainabacteria (Di Rienzi et al.2013) that occurred at very low abundances. Among core genera, most abundant were Bacteroides (average 10% of total reads per specimen) and Parabacteroides (4.6%), followed by Oscillospira, Dorea, Ruminococcus, Desulfovibrio and Coprococcus (all above 1%) (Fig. S3, Supporting Information). Table 4. List of core taxa (90%) by their finest taxonomic resolution, with number of core OTUs per taxon shown in parentheses (total = 134). Phylum  Class  Order  Family  Genus  Species1  Firmicutes (118)  Erysipelotrichi (5)  Erysipelotrichales (5)  Erysipelotrichaceae (5)  Coprobacillus (1)            PSB-M-3            Holdemania      Clostridia (113)  Clostridiales (113)  Lachnospiraceae (37)  Ruminococcus (8)  [Ruminococcus] gnavus          Blautia (2)  Blautia producta          Dorea (10)            Coprococcus            Roseburia          Ruminococcaceae (51)  Oscillospira (24)            Anaerotruncus          Eubacteriaceae (1)  Pseudoramibacter (1)            Anaerofustis          Clostridiaceae  Clostridium          Dehalobacteriaceae  Dehalobacterium          Christensenellaceae (4)  Christensenella (1)          [Mogibacteriaceae] (3)  Anaerovorax          Peptostreptococcaceae (1)            Peptococcaceae (3)            Veillonellaceae            [Tissierellaceae]        Bacilli  Lactobacillales  Lactobacillaceae  Lactobacillus          Leuconostocaceae          Bacillales  Bacillaceae  Bacillus    Actinobacteria (6)  Coriobacteriia (6)  Coriobacteriales (6)  Coriobacteriaceae (6)  Adlercreutzia (1)    Proteobacteria (6)  Deltaproteobacteria (6)  Desulfovibrionales (6)  Desulfovibrionaceae (6)  Desulfovibrio (5)            Bilophila      Gammaproteobacteria  Enterobacteriales  Enterobacteriaceae        Alphaproteobacteria  Campylobacterales  Helicobacteraceae        Epsilonproteobacteria          Bacteroidetes (4)  Bacteroidia (4)  Bacteroidales (4)  Bacteroidaceae (3)  Bacteroides (3)  Bacteroides fragilis            Bacteroides ovatus        Porphyromonadaceae (1)  Parabacteroides (1)  Parabacteroides distasonis            Parabacteroides gordonii        Rikenellaceae  Rikenella          [Odoribacteraceae]  Odoribacter    Verrucomicrobia  Verrucomicrobiae  Verrucomicrobiales  Verrucomicrobiaceae  Akkermansia  Akkermansia muciniphila  Tenericutes  Mollicutes  Anaeroplasmatales  Anaeroplasmataceae  Anaeroplasma        RF39        Cyanobacteria (Melainabacteria)  4C0d-2  YS2        Phylum  Class  Order  Family  Genus  Species1  Firmicutes (118)  Erysipelotrichi (5)  Erysipelotrichales (5)  Erysipelotrichaceae (5)  Coprobacillus (1)            PSB-M-3            Holdemania      Clostridia (113)  Clostridiales (113)  Lachnospiraceae (37)  Ruminococcus (8)  [Ruminococcus] gnavus          Blautia (2)  Blautia producta          Dorea (10)            Coprococcus            Roseburia          Ruminococcaceae (51)  Oscillospira (24)            Anaerotruncus          Eubacteriaceae (1)  Pseudoramibacter (1)            Anaerofustis          Clostridiaceae  Clostridium          Dehalobacteriaceae  Dehalobacterium          Christensenellaceae (4)  Christensenella (1)          [Mogibacteriaceae] (3)  Anaerovorax          Peptostreptococcaceae (1)            Peptococcaceae (3)            Veillonellaceae            [Tissierellaceae]        Bacilli  Lactobacillales  Lactobacillaceae  Lactobacillus          Leuconostocaceae          Bacillales  Bacillaceae  Bacillus    Actinobacteria (6)  Coriobacteriia (6)  Coriobacteriales (6)  Coriobacteriaceae (6)  Adlercreutzia (1)    Proteobacteria (6)  Deltaproteobacteria (6)  Desulfovibrionales (6)  Desulfovibrionaceae (6)  Desulfovibrio (5)            Bilophila      Gammaproteobacteria  Enterobacteriales  Enterobacteriaceae        Alphaproteobacteria  Campylobacterales  Helicobacteraceae        Epsilonproteobacteria          Bacteroidetes (4)  Bacteroidia (4)  Bacteroidales (4)  Bacteroidaceae (3)  Bacteroides (3)  Bacteroides fragilis            Bacteroides ovatus        Porphyromonadaceae (1)  Parabacteroides (1)  Parabacteroides distasonis            Parabacteroides gordonii        Rikenellaceae  Rikenella          [Odoribacteraceae]  Odoribacter    Verrucomicrobia  Verrucomicrobiae  Verrucomicrobiales  Verrucomicrobiaceae  Akkermansia  Akkermansia muciniphila  Tenericutes  Mollicutes  Anaeroplasmatales  Anaeroplasmataceae  Anaeroplasma        RF39        Cyanobacteria (Melainabacteria)  4C0d-2  YS2        Taxonomy in brackets is not officially recognized (under revision). Not all core OTUs were identified below order level. aPresence in 80% of the specimens and all islets. View Large At OTU level, 23 OTUs were present in 100% of the specimens, 134 at 90% and 303 at 80%. Four phyla hosted the 134 core OTUs at 90% prevalence, Firmicutes, Proteobacteria, Actinobacteria and Bacteroidetes (Table 4), with the majority of core OTUs (118) belonging to the phylum Firmicutes and, within this, to the order Clostridiales (113), families Ruminococcaceae (51), Lachnospiraceae (37), and genera Oscillosphira (24), Ruminococcus (8) (both Ruminococcaceae) and Dorea (10, Lachnospiraceae) (Table 4). Despite the fact that Bacteroidetes and Proteobacteria were among the most abundant phyla, they hosted only four and six core OTUs, respectively. All Proteobacteria OTUs belonged to the family Desulfovibronaceae, five of which were classified to the genus Desulfovibrio. All six Actinobacteria core OTUs were Coriobacteriaceae. Overall, these 134 core OTUs accounted for only 2.1% of the total number of OTUs detected, with individual core OTUs typically occurring at low abundances, ranging on average across all specimens between 0.01% and 2.7% of total reads per specimen (Table S3, Supporting Information). The most abundant core OTU was N.Ref.OTU202 (2.7%, classified as Bacteroides), followed by OTU-113 542 (1.6%, Desulfovibrio), N.Ref.OTU980 (1.6%, Erysipelotrichaceae) and OTU-194 286 (1.5%, Parabacteroides gordonii). Altogether, the 134 core OTUs contributed 34.8% of the total reads across all samples and 35.1% average reads per specimen. According to best BLASTN hits and taxonomic classification (Table S3, Supporting Information), core OTUs (i) typically belonged to bacteria known as anaerobic or strictly anaerobic, suggesting a limitation in their habitat range and potential for dispersal; (ii) were not classified as pathogens at present abundance; and (iii) in most cases matched to stool/gut from humans and other vertebrates, rarely invertebrates (representing the main food source). These features suggest large specificity to the vertebrate gut, making these OTUs good candidates to investigate microbiota heritability across host generations and the effect of long-term geographic isolation. Sub-OTU resolution recovers islet signal To address the above scenario, we exploited the sub-OTU diversity of the core (sequence variants within each OTU as inferred by re-clustering of sequences at 98 and 99% similarity). If core OTUs represent a conserved component inherited from the ancestral microbial pool and then ‘confined’ in each islet after population splits, we would expect the sub-OTU variants to carry a signature of their islet-restricted diversification, including local drift and mutation/recombination events. On the contrary, if acquisition of this core occurred by more recent dispersal across islets, the sub-OTU variants should present a stochastic distribution or correlate with geographic distances, but not with islet structure. For this purpose, we used individual and pooled core OTUs, which carry no information for discrimination among specimens/islets, and reclustered all reads at 98 and 99% similarity (the latter approximating the haplotype level). PERMANOVA analyses on all previously tested categorical variables indicated ‘islet’ as the main significant factor (Table 1, P = 0.0001), explaining 15% of the core microbiota sub-OTU variance, a proportion comparable to the one explained by clustering of the whole dataset at 97% similarity. According to Jaccard distance matrices built on de novo 98% and 99% OTUs, distances were significantly shorter within islets than between islets (two-sided Student's two-sample t-test, Bonferroni-corrected P < 0.001). Moreover, mean Jaccard distances across islets mirrored the 97% OTU-based distances: in all cases, the closest islet pair was Aire-PrCv (mean Jaccard = 0.25 and 0.30 at 98% and 99% clustering, respectively). Again, microbiota distances did not correlate with geographic distances (Mantel r = –0.08838, P = 0.752 for 99% clustering and r = –0.10154, P = 0.674 for 98% clustering). To understand the individual contribution of core OTUs to this pattern, the 23 OTUs with presence in 100% of the specimens were individually analyzed (99% OTU clustering, Table S4, Supporting Information). Intra-OTU variation analysis revealed that 18 of the 23 OTUs carried a significant signal by islet (adonis, P < 0.05), explaining between 14% and 23% of the total sub-OTU variation. The five OTUs with no significance had their best BLASTN hits to putative environmental locations, human skin and a novel bacterium from the human gut (i.e. Emergencia timonensis). DISCUSSION Here we provide the first characterization of the P. lilfordi gut microbiota, exploring its compositional variation in adult males from seven isolated populations inhabiting coastal islets around Menorca. Taking advantage of the peculiar phylogeographic history of these populations and the expected minor diversifying role of two crucial variables in microbiota studies (host genotype and diet), we performed a preliminary investigation of the forces driving current patterns of microbiota similarity/dissimilarity among populations, including inheritance, ecological drift, dispersal and local environmental inputs. Geographic separation and age drive the diversity and structuring of the Podarcis microbiota Compatible with a significant spatial and temporal effect, islet was the most significant factor in structuring the Podarcis microbiota, accounting for 15% of the observed microbiota variation (Table 1). Aire in particular, with its geographic remoteness and morphological diversity (it hosts a unique melanic form) (Salvador 2009), presented the most distinct microbial community as indicated by Unifrac distances (Fig. 3) as well as by taxa/OTU uniqueness and enrichment (Tables 2 and 3). Among the fixed factors, within which islet was nested, only ‘age’ was significant when considering Aire and PrCv as the oldest islets. Most of the microbiota variation remained unexplained by all other factors tested, including geographic distances among islets. Consistent with a minor although significant islet effect, microbiota beta-diversity distances across populations were relatively short (average pairwise Jaccard distance = 0.38) and partly associated with time since islet/population isolation. Although the phylogeographic history of the Menorcan populations is far from being fully resolved (and out of the scope of this study), Aire and PrCv represent the two first islets/populations to separate from mainland about 8000–9000 years (Pérez-Mellado 1989), in agreement with both bathymetry and curves of sea level rise in the Mediterranean during Holocene (Vecchi et al.2016). Despite being the most geographically distant islets within our sample set (Fig. 1), their gut microbial composition showed the shortest pairwise Jaccard distances and presented several unique commonalities (Table 3 and Table S2, Supporting Information). Retention of ancestral taxa in these two islets or convergence driven by similar environmental opportunities (open sea and scarce vegetation) could both be invoked. However, none of the ecological variables tested with PERMANOVA were strictly associated with these two islets (Table S1, Supporting Information), largely excluding convergence by niche selection. In terms of gut community similarities across islets, OTU uniqueness per islet/population was remarkably low, with the large majority of OTUs (98% on average, Table 2) occurring in more than one islet. Quantitative differences in common taxa and OTUs representation across islets were also minor (Table 3 and Table S2, Supporting Information) and associated with bacteria belonging to core taxa (i.e. order Clostridiales and family Ruminococcaceae, Table 4), suggesting minor functional community differences among islets. While, to some extent, this large microbial membership and structural similarity might result from common host genotypic constraints and lack of major dietary differences (all populations were largely omnivorous), important differences in ecological niches exist among islets: some show rich vegetation and large anthropic impact (Colm and Rei), while others are characterized by more natural environments, with scarce vegetation and large exposure to the open sea (Aire and PrCv). Overall, this suggests that the striking compositional similarity observed is, at least in part, independent from the environmental inputs/pressures. Ancestral retention versus dispersal effect To explore the origin of the shared taxonomic content among populations, we looked at the nature of core OTUs (defined as presence in 90% of the specimens and all islets). Core OTUs could either represent a conserved inherited component since divergence from the mainland population or an independent acquisition after isolation through taxa dispersal across islets followed by gut selective retention (common niche selection). Excluding recent host immigration events (at least during the past 2000 years), as supported by bathymetric data and the degree of morphological and pigmentation diversity among these populations (Salvador 2009), dispersal should be largely restricted to other biotic or abiotic vectors able to cross the sea barrier. Invertebrates consumed with the diet represent the most likely biotic vectors, particularly flying insects (i.e. all populations are predominantly insectivorous, especially during the summer season) and occasional mollusks and crustaceans captured along the shores. According to best BLASTN hits and taxonomic classification, however, the majority of core OTUs typically matched gut/stool samples of vertebrates, with rare invertebrate hits (Table S3, Supporting Information). This suggests that despite invertebrates (mostly insects) represent the largest component of these lizards’ diet, they are not a major source of gut bacteria. This is in accordance with previous findings showing that invertebrate communities did not significantly contribute to the lizard gut microbiota (Ren et al.2016). Likewise, vertebrates (including rats, rabbits, seagull colonies and humans) are improbable vectors of pervasive dispersal of shared OTUs, as their presence is largely restricted to a few islets (Table S1, Supporting Information; except for Colm and Rei, none of the islets hosts stable human settlements). Aside from animal vectors, recent studies have shown that the phyllosphere microbiome might play a substantial role in shaping the gut microbiota in both lizards (Kohl et al.2017) and desert woodrats (Kohl and Dearing 2012); in this respect, P. lilfordi populations are known to consume vegetable matter (Salvador 2009; Perez-Cembranos, Leon and Perez-Mellado 2016), a finding that is also confirmed by stomach content analysis from this study, revealing plant material and flowers in the majority of samples, and by the presence of several fiber-degrader taxa in the gut (see below). While our current lack of microbial samples from the surrounding environment does not allow us to test the impact of the phyllosphere in these lizards’ microbiota, both best BLASTN matches and large differences in vegetation cover among islets make this scenario improbable. A second important observation is that the 98% and 99% OTU-based structure of the core presented a significant clustering by islet and age (Table 1), which is compatible with an islet-restricted diversification of core OTUs (putatively driven by loss of variants, mutation and recombination events). This hypothesis clearly needs to be corroborated by strain-level analyses. Moreover, these data might suffer from undersampling in some populations/individuals, despite the high sequence depth coverage per specimen obtained (38 160 sequences after rarefaction). Finally, whether any of these sub-OTU variants are under selection or exist in a transitory neutral state remains to be investigated. Indirect evidence of gut microbiome heritability Clearly, the presence of shared OTUs, that are unlikely the result of dispersal, indirectly stands for the inheritance of a microbial component over relatively large temporal and spatial scales. Up to date, levels and mechanisms of gut microbiota inheritance in lizards, and reptiles in general, remain unclear (Colston 2017). Gut microbiotas from Anolis lizards showed a remarkably high intraspecific variation, with only 10% on average of shared OTUs between any two conspecifics in nature (Ren et al.2016), suggesting low microbiota heritability in these lizards. Although their large generalist diet was proposed as a major diversifying factor, this is unlike our findings for P. lilfordi, also a generalist species, where any two conspecifics within an islet share on average 31% of their OTUs and populations are allopatric. Methodological differences, such as the use of full intestines in our study versus fecal samples for Anolis, could be partly responsible for the observed pattern. In contrast, for viviparous lizards of genera Liolaemus and Phymaturus, recent evidence indicated around 34% of microbiome heritability in captivity (Kohl et al.2017), an estimate that nicely matches the percentage of core microbial sequences observed in Podarcis (34.8%). While the Podarcis core summed up to more than one-third of the microbiota, individual core OTUs typically occurred at low abundances (mean = 0.26 ± 0.4% over total reads per specimen, Table S3, Supporting Information), suggesting that, while being more vulnerable to stochastic loss (Chase and Myers 2011), their persistence over time is sustained by important selective forces and putative high inheritability. This is compatible with the general understanding that core gut microbiotas in vertebrates are rather stable in the absence of differential adaptive pressures (such as type of metabolism) (Ley et al.2008; Lozupone et al.2012). In this respect, the overall taxonomic composition of the P. lilfordi gut microbiota mirrored that found in other vertebrates (Fig. 2), with the majority of core families and genera detected (Table 4) representing stable residents of human and mouse guts (Xiao et al.2015), and being consistently found in other lizards (Ren et al.2016; Kohl et al.2017). Core taxa were particularly enriched in bacteria known as active fiber degraders, including the families Ruminococcaceae and Lachnospiraceae (Fig. 2, both Clostridiales) (Biddle et al.2013) and the genera Bacteroides, Blautia, Akkermansia, Oscillospira and Desulfovibrio (Fig. S3, Supporting Informtion). Bacteroides, for instance, which was the most conspicuous genus in our samples, also makes up a very large proportion of the mammalian gut where it is an active degrader of plant material (Ley et al.2008; Xiao et al.2015). Oscillospira, also highly abundant in humans, where it shows high heritability and is associated with leanness, is largely enriched in plant-fed lizards (Kohl et al.2017). Overall, the P. lilfordi core composition is compatible with major conserved functions in vertebrates and the consumption of plant material, suggesting common selective pressures for its maintenance across host generations. Nonetheless, considering that parental care in reptiles is rare and not documented in the genus Podarcis (oviparous), means of bacteria recycling at each new host generation remain unclear. Coprophagy, cannibalism or other types of social behavior have been proposed as mechanisms for bacteria transmission (Colston and Jackson 2016; Perez-Cembranos, Leon and Perez-Mellado 2016; Colston 2017); however, they might suffer from much stochasticity, compromising their ability to preserve a metacommunity structure over time. Covering the eggs with a bacterial layer by means of fecal contamination during oviposition through the single urogenital opening present in lizards (the cloaca) might represent a more reliable mechanism for a stable microbiota transmission. Overall, although our results provided important indirect evidence for gut microbiome heritability in P. lilfordi, confirming previous data for lizards (Kohl et al.2017), this topic definitely demands more empirical investigation. Islet microbiotas as nested subsets of the ancestral mainland microbial pool: a hypothesis To summarize, our main findings indicated that (i) there is a large taxonomic microbiota conservation among all seven allopatric populations studied, down to OTU level, with a reduced impact of local exposure (i.e. uniqueness is negligible); (ii) the minor variation detected is explained by ‘islet’ and ‘age’, supporting a role of the host phylogeographic history in shaping microbiota composition; (iii) dispersal of OTUs among islets appears to be low; and (iv) core variation at 98%–99%-OTU level carries a significant ‘islet’ and ‘age’ signature, compatible with an islet-restricted evolution. Altogether, these findings suggest that persistence of taxa from the mainland stock and ecological drift might have played a substantial role in shaping current pattern of microbiota similarity/dissimilarity among P. lilfordi populations. This intriguing hypothesis clearly requires more extensive quantitative data, targeting additional specimens, islets and microbial markers, as well as ad hoc statistical tools to test drift in host-associated microbiota across populations, a field still in its infancy. Although overlooked in current literature, host population demographic history is likely to represent a crucial factor in explaining current gut metacommunity diversity, as host drift could, to some extent, cause microbiota ecological drift (Orrock and Watling 2010; Chase and Myers 2011). Up to date, only few studies have attempted to address this topic. In the Galápagos archipelago, processes of ecological drift were proposed to explain intraspecific differences in marine iguanas (Lankau, Hong and Mackie 2012). Likewise, in two Anolis lizard species, the complex pattern of similarity/dissimilarity observed, only partially attributed to geography and with no differences among ‘ecomorphs’ (Ren et al.2016), could be potentially explained by chance variation driven by drift and dispersal. In our study system, population densities varied greatly across islets (Fig. 1), with Aire presenting one of the highest population densities of lizards/vertebrates worldwide (Pérez-Mellado et al.2008). Although past population demography for these islets is poorly known (Pérez-Mellado et al.2008), the largest populations (Aire, PrCv and Colm) carried slightly higher diversities (PD whole tree) compared to smaller islets (Tosq, Ravl and Rei), though differences were not significant (Fig. S2, Supporting Information). The small population found in Tosqueta, in particular, which showed the strongest genetic drift among all 16 islets according to unpublished microsatellite data (personal observations), carried the least diverse microbial community and lowest uniqueness, despite the extensive vegetation cover and ample opportunities for local inputs in this islet. These observations, although preliminary, point to an interesting correlation between population demographic history and drift effect in gut metacommunities, which could potentially explain a large part of the extant variation, a topic that should definitely warrant further investigation. CONCLUSIONS Overall, the minor dietary and genotypic differences among our sampled populations (suggesting weak diversifying selection) and the geographically restricted areas to which they are confined, reducing the exposure to microbial pools from other vertebrates, have most likely contributed to the large compositional stability of P. lilfordi gut microbiota, ultimately revealing the strength of drift. We envision that, owing to a progressive within-islet microbiota evolution, in a context of no host immigration events, the core OTU signature of these populations’ past host legacy will gradually disappear, while resemblance at higher taxonomic levels will persist for a longer time. Integrating our data with a fine-scale resolution of P. lilfordi phylogeographic history and a bacterial strain-level characterization that targets core OTUs (under development) might offer a way to disentangle the dynamics of this microbiota component over such short evolutionary time scales. Finally, extending the sampling to all remaining populations from Balearic Islands, including the closely related species P. pityusensis, will provide important comparative data for a most integrative understanding of the gut microbiota diversification in these lizards. Data accessibility Raw sequences are available in the Sequence Read Archive (SRA) database at NCBI under BioProject ID PRJNA429123PRJNA429123. SUPPLEMENTARY DATA Supplementary data are available at FEMSEC online. Acknowledgements We thank Sam Pons, Jordi García-Darocas and Montserrat Farré for aid in sampling, and the Animal Biology department (University of Barcelona) for use of equipment. FUNDING This work was supported by funding from the Spanish Ministerio de Economía y Competitividad to LB (RYC-2012-11443). Conflicts of Interest. None declared. REFERENCES Aivelo T, Norberg A. Parasite-microbiota interactions potentially affect intestinal communities inwildmammals. J Anim Ecol  2017, doi: 10.1111/1365-2656.12708. Alcover JA, Moyá-Solà S, Pons-Moyá J. Les quimeres del passat. Els vertebrats fòssils del Plio-Quaternari de les Balears i Pitiüses . Palma, Espanya: Editorial Moll, 1981. Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecol  2001; 26: 32– 46. Anderson MJ, Gorley RN, Clarke KR. PERMANOVA+ for PRIMER: Guide to Software and Statistical Methods . Plymouth, UK: PRIMER-E, 2008. Baldo L, Pretus JL, Riera JL et al.   Convergence of gut microbiotas in the adaptive radiations of African cichlid fishes. ISME J  2017; 11: 1975– 87. Google Scholar CrossRef Search ADS PubMed  Baldo L, Riera JL, Tooming-Klunderud A et al.   Gut microbiota dynamics during dietary shift in eastern african cichlid fishes. PLoS One  2015; 10: e0127462. Google Scholar CrossRef Search ADS PubMed  Benson AK, Kelly SA, Legge R et al.   Individuality in gut microbiota composition is a complex polygenic trait shaped by multiple environmental and host genetic factors. P Natl Acad Sci USA  2010; 107: 18933– 8. Google Scholar CrossRef Search ADS   Biddle A, Stewart L, Blanchard J et al.   Untangling the genetic basis of fibrolytic specialization by lachnospiraceae and ruminococcaceae in diverse gut communities. Divers Distrib  2013; 5: 627– 40. Google Scholar CrossRef Search ADS   Brown RP, Terrasa B, Perez-Mellado V et al.   Bayesian estimation of post-Messinian divergence times in Balearic Island lizards. Mol Phyl Evol  2008; 48: 350– 8. Google Scholar CrossRef Search ADS   Burns AR, Stephens WZ, Stagaman K et al.   Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. ISME J  2016; 10: 655– 64. Google Scholar CrossRef Search ADS PubMed  Caporaso JG, Kuczynski J, Stombaugh J et al.   QIIME allows analysis of high-throughput community sequencing data. Nat Methods  2010; 7: 335– 6. Google Scholar CrossRef Search ADS PubMed  Chase JM, Myers JA. Disentangling the importance of ecological niches from stochastic processes across scales. Philos T Roy Soc B  2011; 366: 2351– 63. Google Scholar CrossRef Search ADS   Clarke KR, Gorley RN. PRIMER v7: User Manual/Tutorial . Plymouth, UK: PRIMER-E, 2015. Colston TJ. Gut microbiome transmission in lizards. Mol Ecol  2017; 26: 972– 4. Google Scholar CrossRef Search ADS PubMed  Colston TJ, Jackson CR. Microbiome evolution along divergent branches of the vertebrate tree of life: what is known and unknown. Mol Ecol  2016; 25: 3776– 800. Google Scholar CrossRef Search ADS PubMed  Costello EK, Stagaman K, Dethlefsen L et al.   The application of ecological theory toward an understanding of the human microbiome. Science  2012; 336: 1255– 62. Google Scholar CrossRef Search ADS PubMed  David LA, Maurice CF, Carmody RN et al.   Diet rapidly and reproducibly alters the human gut microbiome. Nature  2014; 505: 559– 63. Google Scholar CrossRef Search ADS PubMed  Delsuc F, Metcalf JL, Wegener Parfrey L et al.   Convergence of gut microbiomes in myrmecophagous mammals. Mol Ecol  2014; 23: 1301– 17. Google Scholar CrossRef Search ADS PubMed  Di Rienzi SC, Sharon I, Wrighton KC et al.   The human gut and groundwater harbor non-photosynthetic bacteria belonging to a new candidate phylum sibling to Cyanobacteria. eLife  2013; 2: e01102. Google Scholar CrossRef Search ADS PubMed  Funkhouser LJ, Bordenstein SR. Mom knows best: the universality of maternal microbial transmission. PLoS Biol  2013; 11: e1001631. Google Scholar CrossRef Search ADS PubMed  Kohl KD, Brun A, Magallanes M et al.   Gut microbial ecology of lizards: insights into diversity in the wild, effects of captivity, variation across gut regions and transmission. Mol Ecol  2017; 26: 1175– 89. Google Scholar CrossRef Search ADS PubMed  Kohl KD, Dearing MD. Experience matters: prior exposure to plant toxins enhances diversity of gut microbes in herbivores. Ecol Lett  2012; 15: 1008– 15. Google Scholar CrossRef Search ADS PubMed  Krijgsman W, Hilgen FJ, Raffi Y et al.   Chronology, causes and progression of the Messinian salinity crisis. Nature  1999; 400: 652– 5. Google Scholar CrossRef Search ADS   Lankau EW, Hong PY, Mackie RI. Ecological drift and local exposures drive enteric bacterial community differences within species of Galapagos iguanas. Mol Ecol  2012; 21: 1779– 88. Google Scholar CrossRef Search ADS PubMed  Ley RE, Lozupone CA, Hamady M et al.   Worlds within worlds: evolution of the vertebrate gut microbiota. Nat Rev Microbiol  2008; 6: 776– 88. Google Scholar CrossRef Search ADS PubMed  Lozupone CA, Stombaugh JI, Gordon JI et al.   Diversity, stability and resilience of the human gut microbiota. Nature  2012; 489: 220– 30. Google Scholar CrossRef Search ADS PubMed  Martinson VG, Douglas AE, Jaenike J. Community structure of the gut microbiota in sympatric species of wild Drosophila. Ecol Lett  2017; 20: 629– 39. Google Scholar CrossRef Search ADS PubMed  Moeller AH, Caro-Quintero A, Mjungu D et al.   Cospeciation of gut microbiota with hominids. Science  2016; 353: 380– 2. Google Scholar CrossRef Search ADS PubMed  Muegge BD, Kuczynski J, Knights D et al.   Diet drives convergence in gut microbiome functions across mammalian phylogeny and within humans. Science  2011; 332: 970– 4. Google Scholar CrossRef Search ADS PubMed  Nemergut DR, Schmidt SK, Fukami T et al.   Patterns and processes of microbial community assembly. Microbiol Mol Biol R  2013; 77: 342– 56. Google Scholar CrossRef Search ADS   Orrock JL, Watling JI. Local community size mediates ecological drift and competition in metacommunities. Proc Biol Sci  2010; 277: 2185– 91. Google Scholar CrossRef Search ADS PubMed  Perez-Cembranos A, Leon A, Perez-Mellado V. Omnivory of an insular lizard: sources of variation in the diet of Podarcis lilfordi (Squamata, Lacertidae). PLoS One  2016; 11: e0148947. Google Scholar CrossRef Search ADS PubMed  Pérez-Mellado V. Estudio ecológico de la lagartija balear Podarcis lilfordi (Günther, 1874) en Menorca. Rev de Menorca  1989; 53: 455– 511. Pérez-Mellado V, Hernández-Estévez JA, García-Díez T et al.   Population density in Podarcis lilfordi (Squamata, Lacertidae), a lizard species endemic to small islets in the Balearic Islands (Spain). Amphibia-Reptilia  2008; 29: 49– 60. Google Scholar CrossRef Search ADS   Pretus JL, Marqués R, Pérez-Mellado V. Holocene sea level rise 975 and range of fragmentation of Podarcis lilfordi on Minorcan islets: a vicariance scenario reviewed through a mtDNA tree. In Pérez-Mellado V, Riera N, Perera A, (eds). The Biology of lacertid lizards: Evolutionary and Ecological Perspectives . Institut Menorquí d'Estudis . 2004; 8: 279– 291. Price MN, Dehal PS, Arkin AP. FastTree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol Biol Evol  2009; 26: 1641– 50. Google Scholar CrossRef Search ADS PubMed  Rawls JF, Mahowald MA, Ley RE et al.   Reciprocal gut microbiota transplants from zebrafish and mice to germ-free recipients reveal host habitat selection. Cell  2006; 127: 423– 33. Google Scholar CrossRef Search ADS PubMed  Ren T, Kahrl AF, Wu M et al.   Does adaptive radiation of a host lineage promote ecological diversity of its bacterial communities? A test using gut microbiota of Anolis lizards. Mol Ecol  2016; 25: 4793– 804. Google Scholar CrossRef Search ADS PubMed  Saez ET, Traveset A. Fruit and nectar feeding by Podarcis lilfordi (Lacertidae) on Cabrera archipelago (Balearic islands). Herpetol Rev  1995; 26: 121– 3. Salvador A. Lagartija balear – Podarcis lilfordi (Günther, 1874). In Salvador A, Marco A., (eds). Museo Nacional de Ciencias Naturales, Madrid . Enciclopedia Virtual de los Vertebrados Españoles, 2009. http://www.vertebradosibericos.org. Sanders ED. Evidence Concerning Late Survival and Extinction of Endemic Amphibia and Reptilia from the Bronze and Iron Age settlement of Torralba d’En Salort (Alayor, Menorca) . Palma de Mallorca: Monografies Cientìfiques, Palma de Mallorca. 1984. Schloss PD, Westcott SL, Ryabin T et al.   Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microb  2009; 75: 7537– 41. Google Scholar CrossRef Search ADS   Shoemaker W, Locey KJ, Lennon JT. A macroecological theory of microbial biodiversity. Nat Ecol Evol  2017; 1: 107. Google Scholar CrossRef Search ADS PubMed  Song SJ, Lauber C, Costello EK et al.   Cohabiting family members share microbiota with one another and with their dogs. eLife  2013; 2: e00458. Google Scholar PubMed  Terrasa B, Perez-Mellado V, Brown RP et al.   Foundations for conservation of intraspecific genetic diversity revealed by analysis of phylogeographical structure in the endangered endemic lizard Podarcis lilfordi. Divers Distrib  2009; 15: 207– 21. Google Scholar CrossRef Search ADS   Traveset A, Riera N. Disruption of a Plant-Lizard seed dispersal system and its ecological effects on a threatened endemic plant in the Balearic islands. Conserv Biol  2005; 19: 421– 31. Google Scholar CrossRef Search ADS   Vecchi M, Marriner N, Morhange C et al.   Multiproxy assessment of Holocene relative sea-level changes in the western Mediterranean: Sea-level variability and improvements in the definition of the isostatic signal. Earth-Sci Rev  2016; 155: 172– 97. Google Scholar CrossRef Search ADS   Xiao L, Feng Q, Liang S et al.   A catalog of the mouse gut metagenome. Nat Biotechnol  2015; 33: 1103– 8. Google Scholar CrossRef Search ADS PubMed  Yatsunenko T, Rey FE, Manary MJ et al.   Human gut microbiome viewed across age and geography. Nature  2012; 486: 222– 7. Google Scholar PubMed  Zeng Q, Sukumaran J, Wu S et al.   Neutral models of microbiome evolution. PLoS Comput Biol  2015; 11: e1004365. Google Scholar CrossRef Search ADS PubMed  © FEMS 2017. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com

Journal

FEMS Microbiology EcologyOxford University Press

Published: Feb 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off