Faecal Microbiota Divergence in Allopatric Populations of Podarcis lilfordi and P. pityusensis, Two Lizard Species Endemic to the Balearic Islands

Gut microbial communities provide essential functions to their hosts and are known to influence both their ecology and evolution. However, our knowledge of these complex associations is still very limited in reptiles. Here we report the 16S rRNA gene faecal microbiota profiles of two lizard species endemic to the Balearic archipelago (Podarcis lilfordi and P. pityusensis), encompassing their allopatric range of distribution through a noninvasive sampling, as an alternative to previous studies that implied killing specimens of these IUCN endangered and near-threatened species, respectively. Both lizard species showed a faecal microbiome composition consistent with their omnivorous trophic ecology, with a high representation of cellulolytic bacteria taxa. We also identified species-specific core microbiota signatures and retrieved lizard species, islet ascription, and seasonality as the main factors in explaining bacterial community composition. The different Balearic Podarcis populations are characterised by harbouring a high proportion of unique bacterial taxa, thus reinforcing their view as unique and divergent evolutionary entities. Supplementary Information The online version contains supplementary material available at 10.1007/s00248-022-02019-3.


Introduction
Recent advances in sequencing technologies and analytical methodologies are improving our understanding of the microbiome in host evolution [1]. Early evolutionary conceptions considering animals and plants as autonomous entities are being challenged by the holobiont point of view, which also considers their numerous microbial symbionts and their genomes [2]. The sum of the genetic information of the host and its associated microbiota has been termed the hologenome [3,4], whose variation can influence phenotypes upon which natural selection and/or genetic drift can operate [2]. Indeed, microbial communities can have a deep impact on host diversification by acting as environmental factors with selective effects [5], influencing many aspects of host evolutionary history such as adaptation to resource utilisation [6], resistance to pathogens [7], control of nutrient inputs [8], tissue and organ development [9], life history strategy [10] and behaviour [11], among many others. The relevance of this interaction is such that the evolutionary history of host species cannot be fully understood without addressing the study of its associated microbiota. However, despite the effort made over the last decade to characterise microbial-host associations, we still know relatively little about the evolutionary and ecological processes shaping them, particularly in non-model and non-captive organisms.
The lizard species endemic to the Balearic Islands, Podarcis lilfordi and P. pityusensis, represent an interesting model to study host-microbiota associations since they are sister taxa with a nonoverlapping distribution, consisting of multiple allopatric populations restricted to coastal islands and islets of Mallorca, Menorca, and Cabrera in the case of P. lilfordi, and the main islands of Ibiza and Formentera in P. pityusensis. Podarcis lilfordi became extinct on both the main Mallorca and Menorca islands due to the anthropic introduction of foreign predators and/or competitors [12,13]. The feeding ecology of both species is well-known [14][15][16][17] and is marked by the scarcity and unpredictability of diet resources in their isolated habitats, which probably determined the adoption of omnivory [17]. Balearic lizards are active foragers exploiting a wide range of animal prey, plant tissues, carrion, and marine subsidies, showing even low levels of cannibalism of juvenile conspecifics [17]. On the other hand, the phylogeographic relationships of the group have also been profusely studied, revealing an evolutionary origin linked to eustatic sea-level changes associated with the reflooding of the Mediterranean at the end of the Messinian Salinity Crisis that occurred 5.33 Ma ago [18][19][20]. The Menorca lineage represents the earliest cladogenetic event (2.6 Ma) within P. lilfordi, followed by the differentiation of the West Mallorca lineage (2.0 Ma), Cabrera (1.2 Ma), and the remaining populations in northern and southern Mallorca islets [20] (see Fig. 1 for island configuration in the Balearic archipelago). Regarding P. pityusensis, the Ibiza and Formentera populations have been reported to be genetically distinct, with a divergence estimated to have occurred ca. 0.111-0.295 Ma ago [21].
Although lizards are a globally distributed and speciesrich group within vertebrates [22,23], they have been largely overlooked in terms of gut microbiota [24]. More than 90% of such studies within vertebrate hosts have been carried out in mammal species, with reptiles being the least investigated group [25]. To date, a single study on the microbiota of the Balearic lizards has been addressed [26]. This work was focused on a very limited portion of the distribution of one of the species (seven Menorcan populations of P. lilfordi) and based on sex-biased sampling (30 males and 3 females) from a single season (summer) using methodologies that imply killing specimens of this endangered (IUCN) species. The study reported a significant effect of the geographical distribution (islet) in explaining the 13.5% of the bacterial community variation and a very low bacterial uniqueness in each lizard population, suggesting the retention of the ancestral mainland microbial pool and the occurrence of stochastic population processes as the main factors shaping the gut microbiota.
Here we aim to further explore the evolutionary history of these two endemic lizard species by using Illumina 16S rRNA sequencing to characterise their faecal microbiota across different seasons of the year and through a sex-balanced and noninvasive sampling. Specifically, we address the following questions: (i) what is the gut bacterial composition of these two endemic species? (ii) Are there species-specific microbiome core signatures? (iii) Do microbial communities mirror the allopatric distribution of their respective lizard populations? (iv) What is the degree of uniqueness of each Podarcis population from a microbiome perspective? (v) Do the microbiome communities reflect the trophic ecology of the host populations?

Sampling
The sampling design included 17 localities encompassing the distribution range of P. lilfordi and P. pityusensis in the Balearic archipelago (Fig. 1, Table 1, and Table S1a). A total of 242 faecal samples from P. lilfordi (140) and P. pityusensis (102) were collected between the spring of 2016 and the autumn of 2017. Specimens were captured by noosing and fresh faeces were collected in absolute ethanol vials directly from the animals before releasing them back to their habitats. Individuals were sexed based on the number and size of the femoral pores [27]. The samples were immediately preserved at 4 °C in the field and upon arrival at the laboratory, where they were stored at − 20 °C until DNA extraction.

DNA Extraction, 16S rRNA Library Preparation, and Sequencing
Total DNA was extracted from individual samples using the ISOLATE Fecal DNA kit (Bioline, London, UK) following the manufacturer protocol, and their concentrations were quantified using Qubit fluorometric quantitation (Ther-moFisher, Foster City, CA, USA). For cost-effective reasons and given that our study seeks to describe the bacterial composition of each lizard population as a whole, samples from each island/islet with matching sex, season, and collecting year were pooled in equimolar concentrations, obtaining a final volume of 30 μl at 20 ng/μl per sample (see Table 1 for details on the number of faecal samples pooled per location). A total of 48 samples were submitted to the Roy J. Carver Biotechnology Center (University of Illinois, USA) for amplification of the V4 region of 16S rRNA in a microfluidic high-throughput multiplexed PCR platform (Fluidigm). The primer set 515F (5′-GTG CCA GCMGCC GCG GTAA-3′) and 806R (5′-GGA CTA CHVGGG TWT CTAAT-3′) [28] were used, and CS1 and CS2 Fluidigm universal tags,

Sequence Analyses
QIIME2 version 2020.2 [29] was used for read demultiplexing and subsequent filtering and denoising using the DADA2 pipeline [30]. Sequences were grouped into amplicon sequence variants (ASVs or 100% identity OTUs), and the taxonomic assignment was performed using the q2-featureclassifier plugin with a pre-trained naive Bayes classifier [31] implemented in QIIME2 against the SILVA database release 132 [32]. ASVs identified as chloroplasts or mitochondria, and those with undetermined phylum annotation were excluded for downstream analyses. Sequences were aligned with MAFFT [33] under the default FFT-NS-1 algorithm, and a phylogenetic tree was inferred with FastTree [34] using the QIIME2 plugin align-to-tree-mafft-fasttree.

Datasets
To avoid comparing samples from different seasons and to dissect the potential effect of the host lizard species (P. lilfordi / P. pityusensis) and geography (islet) on the faecal microbiota, we subdivided the data into three datasets by collecting seasons: spring, summer, and autumn. In addition, we further explored the data by analysing those samples with matching collecting seasons: spring + summer dataset [Espardell (Fig. 1f), Na Gorra (Fig. 1e), Aire and Colom ( Fig. 1b)], and summer + autumn dataset [Dragonera ( Fig. 1d), Cabrera, Esclatasang and Na Foradada ( Fig. 1g)].

Diversity Analyses
Rarefaction curves were explored using the rarecurve function implemented in the R package vegan [35]. The ASV table was rarefied using the rarefy_even_depth option in the R package phyloseq [36] by subsampling the data to the even depth defined by the minimum number of sequences per sample. This value was set to 10,815 in the P. lilfordi dataset and 19,159 in the case of P. pityusensis. When analysing the dataset of both species combined, the lowest value was chosen. Alpha diversity was measured through the estimation of the absolute number of observed ASVs, Shannon and Simpson indexes in phyloseq, and Faith's phylogenetic diversity (PD) in the picante R package [37]. Differences in terms of alpha diversity at both inter-and intra-specific levels were assessed using the Kruskal-Wallis test.
Although the pooled nature of the samples erases the interindividual variability within each island/islet, we carried out exploratory beta diversity analyses as a proxy to understand the extent of change in bacterial community composition across lizard populations. To delve into it, we performed a permutational multivariate analysis of variance (PERMANOVA) using the adonis function implemented in the vegan package and based on both unweighted and weighted UniFrac and Bray-Curtis distance matrices [35]. Specifically, we tentatively explored potential correlations between community microbiota distances and the categorical variables associated to our samples, namely species, sex and islet. The later variable was also explored at a broader scale by assigning each sampling locality (i.e., islet/population) to a main island district in the Balearic archipelago (Mallorca, Menorca, Cabrera, Ibiza, or Formentera). The unweighted UniFrac matrices were also used to perform PCoA ordination analyses. Finally, we tested for correlation between the variable number of faecal pellets pooled per sample and microbiome community composition distances (Unweighted and Weighted UniFrac matrices) by conducting Mantel tests based on the Spearman correlation method and performing 9999 replicates.

Core Microbiota
Shared microbiome taxa were independently investigated for each one of the two lizard species with the core members function of the microbiome R package [38], setting the prevalence threshold at 90%. The intersection of the resulting sets was interpreted as the core of the lineage P. lilfordi + P. pityusensis, and their differences were, respectively, considered the core signature of each one of the two Podarcis species.

Uniqueness and Shared Microbial Taxa
UpSetR package [39] was used to calculate both unique and shared ASVs by population and/or species.

Sequence Summary
We obtained 16S rRNA sequence data from 48 faecal pools representing 140 specimens of P. lilfordi and 102 of P. pityusensis ( Fig. 1 and Table 1). After the quality filtering stage and the removal of features identified as chloroplasts, mitochondria, or Eukaryota and those with undetermined phylum annotation, the dataset was reduced to 3,488,522 high-quality reads with an average value of 72,678 reads per sample (range = 11,107-336,352). The reads could be ascribed to 3542 ASVs representing 20 phyla, 41 classes, 93 orders, 190 families, and 425 genera of microbiome taxa. A total of 369 ASVs were unique to P. lilfordi and 766 to P. pityusensis. The sampling depth values used in the rarefaction analyses allowed us to get most of the diversity of the samples, as shown by the asymptotic trend of the resulting curves (Fig. S1). The taxonomy of rarefied and the nonrarefied ASV tables are provided as Table S1b, c, and d.

Alpha Diversity
Faecal samples from P. lilfordi and P. pityusensis showed similar average alpha diversity values ( Table 2) with no significant differences between both lizard species (Table S3). At the population level, most indicators reported Es Vedrà as the location with the highest average microbiome diversity values, followed by Na Gorra, Bosc, Espardell, and Formentera, all of them belonging to P. pityusensis. Oppositely, two P. pityusensis populations received the lowest diversity values (Bleda and Vaixell). Podarcis lilfordi populations were scored with relatively moderate values for all indexes, standing out Na Foradada, Cabrera, and Esclatasang as the most diverse locations, while the Menorcan populations of Colom and Porros islets received the lowest alpha diversity values ( Table 2).
When considering only samples collected in the same season and lizard species, all datasets showed similar average diversity indexes ( Table 2) and no significant  (Table S4). However, inter-dataset comparisons revealed significant differences in terms of alpha diversity indices (observed number of ASVs, Chao1, and phylogenetic diversity) between the spring samples of P. lilfordi and those of P. pityusensis (Table S3 and

Beta Diversity
Mantel tests for correlation between differences in terms of the number of faecal pellets pooled on each sample and microbiome community distances yielded non-significant results regardless of whether the analyses were based on the Unweighted or the Weighted UniFrac matrices from the entire dataset (i.e., all 48 samples) or from each of the subsets by season (spring, summer, and autumn datasets) ( Table S5a-i, Table S6). We therefore tentatively explored if the microbiota community composition from the analysed pooled samples varied across species, island/islet, main island district, season and/or sex. PERMANOVA analyses yielded significant microbiome composition differences between the two lizard species and, also, in relation to the geographical distribution of the samples, regardless the collecting season and of whether the analyses were based on weighted or unweighted UniFraq matrices (Table 3). Among the significant variables, island/islet was reported as the most important source of variance in explaining bacterial community composition, with R 2 values ranging from 0.36 to 0.48, followed by the main island district variable (R 2 range: 0.06 to 0.21). When analysing the datasets conformed by samples with matching collecting seasons (i.e., spring + summer and summer + autumn datasets), PERMANOVA tests consistently retrieved the species, the island/islet, and the main island district as significant variables in explaining the microbial community composition of the samples, and in addition both the season variable and the interaction of population:season were also retrieved as significant in some of the analyses ( Table 3). The ordination analysis of the bacterial community distances by island/islet and species through principal coordinate analysis (PCoA) based on the unweighted UniFrac matrices from the samples collected in the same season of the year resulted in the clustering of the samples from the same location with few exceptions (Fig. 4).

Core Composition
The number of core ASVs present in at least 90% of the samples of P. lilfordi and P. pityusensis was, respectively,    Table 3 Results of PERMANOVA analyses based on both unweighted and weighted UniFrac distance matrices for each season.

ASV Uniqueness
We explored the presence/absence of population-specific ASVs/taxa (Table 4). A total of 2042 ASVs were found in P. lilfordi, while 2384 were detected in P. pityusensis. Within P.

Non-Destructive Sampling as an Alternative to Killing Threatened Podarcis Lizards for Microbiome Analyses
Here we report the first characterisation of the Podarcis pityusensis faecal microbiota and expand the geographic sampling of P. lilfordi (Baldo et   Although the different methodologies implemented in both studies make it difficult to carry out differential abundance analyses, comparing our results to past studies of the same lizard species shows that faecal samples can recover a high proportion of the taxa seen in gut samples. Faecal sampling has been used in microbiome studies on reptiles (e.g., [40][41][42]) and has been demonstrated to provide a comprehensive understanding of their hindgut bacterial communities [24]. Consistently, our findings reinforce this view and suggest that there is no need to kill lizards to characterise their microbiota at least at high taxonomic levels, which is even more relevant when dealing with species included in the IUCN Red List, such as the Balearic lizards P. lilfordi (endangered) and P. pityusensis (near threatened).

The Faecal Microbiota of the Balearic Podarcis
Most lizard species are regarded as primarily feeding on invertebrates (e.g., [43,44]), and less than 2% of the species are known to exploit plants as their sole food source [45]. However, many species do consume plant tissues under conditions of prey scarcity, a behaviour that is more frequently observed in island taxa (see [46] and references therein), including the Podarcis species endemic to the Balearic archipelago [16,47]. Even though herbivory is rare among reptiles, they are known to perform hindgut fermentation either in the cecum or in the large colon/intestine, as also occurs in many lineages of mammals and birds [48,49]. Indeed, there exists evidence of rapid adaption to plant diet through the acquisition of cecal valves, which slow down food flow and act as fermenting chambers in other Podarcis lizard species from the Mediterranean [50]. Like other vertebrates, reptiles lack the endogenous glycoside hydrolases needed to effectively hydrolyse and to ferment the complex plant polymers found in celluloses and hemicelluloses [51], and therefore rely on specialised bacterial communities [52]. Previous studies on the gastrointestinal microbiota of herbivore reptiles have reported a high prevalence of cellulolytic bacteria belonging to the phyla Bacteroidetes and Firmicutes (e.g., [24,52,53]), a pattern that is also matched by our results. Most Firmicutes taxa detected in the faecal microbiota of both P. lilfordi and P. pityusensis belong to class Clostridia order Clostridiales, a lineage of bacteria that includes most Firmicutes in both mammalian and reptilian herbivores [53]. Within Clostridiales, our results are highly represented by families Lachnospiraceae and Ruminococcaceae, both typically found in the gut microbiota of animals and known to decompose complex plant material [54]. Examples of Ruminococcaceae genera reported here with known implications in fibre digestion include Oscillospira and Ruminoccocus [55], represented in our dataset by 4 and 26 distinct ASVs, respectively (Table S1). Regarding Bacteroidetes, both Podarcis species showed a high proportion of microbes from families Bacteroidaceae, Porphyromonadaceae, Rikenellaceae, and Odoribacteraceae, all of them previously reported from the gut microbiota of reptiles with herbivorous habits (e.g., [26,40,56]). Within Bacteroidaceae, our results have reported the presence of 104 different ASVs from genus Bacteroides (Table S1), a lineage of active degraders of plant material [57,58]. The gut microbiotas of P. lilfordi and P. pityusensis also include bacterial lineages that are consistent with their omnivorous ecology. This would be the case of the members of phylum Deferribacteres, which seem to be absent in the microbiota of generalist herbivorous lizards but is found in species that exploits both animal preys and plants [24], or the Clostridiales vadinBB60 group, a Firmicutes lineage with a high prevalence in carnivorous reptiles [59]. The relatively low prevalence of Actinobacteria and Proteobacteria in Podarcis samples is consistent with previous reports from other lizard species [24]. Within the latter our results report the presence of 20 distinct ASVs from Desulfovibrio (Table S1), a genus whose abundance has been correlated with fibre digestion [60] and that may play an important role in herbivorous lizards [24].

Species-Specific Microbiome Signatures
Our exploratory beta diversity analyses pointed to significant differences in bacterial community composition between P. lilfordi and P. pityusensis samples regardless of the dataset (spring, summer, or spring + summer) and the input matrix (weighted or unweighted UniFrac) ( Table 3). Such differences are also evident at the level of alpha diversity, where we have found significant differences between samples of both species collected in spring, being P. pityusensis the species exhibiting the highest diversity values. Although these results should be interpreted with caution due to the pooled nature of the samples, there is evidence that microbiota analyses based on mixed samples are a viable measure to consider in population-level studies, providing estimates of the community-level diversity that are highly correlated with diversity estimates using individually sequenced samples [61][62][63][64][65][66]. The core microbiomes of both Podarcis species intersected in 13 ASVs from several genera commonly found in the gut communities of other lizard species (Alistipes, Bacteroides, Breznakia, Desulfovibrio, Odoribacter, Oscillibacter, Parabacteroides, and Ruminococcaceae UBA1819; e.g., [67][68][69] and also showed a symmetric difference consisting of bacterial taxa with high prevalence levels in either P. lilfordi or P. pityusensis. In the former we detected 11 species-specific core taxa that could be classified to genus (Alistipes, Anaeroplasma, Bacteroides, Bacteroides, Desulfovibrio, Rikenellaceae dgA-11 gut group, Helicobacter, and Parabacteroides). Taxa from these genera were also identified as core members of P. lilfordi gut microbiota by [26] excepting Alistipes, Breznakia, Oscillibacter, Ruminococcaceae UBA1819, and Rikenellaceae dgA-11 gut group. Podarcis pityusensis also showed a distinctive core taxa assemblage consistent with bacteria genera commonly found in the gut microbiome of other lizard species (Bacteroides, Coprococcus 3, Dielma, Erysipelatoclostridium, Eubacterium, Odoribacter, Parabacteroides, Robinsoniella, and Romboutsia) [42,67,70]. All these findings are consistent with the long-term geographic isolation of both P. lilfordi and P. pityusensis lineages that started ca. 5.33 Ma ago and subsequent geographical, ecological and evolutionary divergence [20].

Factors Shaping the Faecal Microbiota of the Balearic Podarcis Lizards
Extant P. lilfordi populations are restricted to coastal islands and islets of Mallorca, Menorca, and the Cabrera archipelago, which acted as refugia after the extinction of mainland populations ca. 2000 years ago due to the anthropic introduction of predators and competitors [12,13]. Such isolation could be even older according to Holocene sea-level data from the western Mediterranean Sea [71]. Therefore, each islet would constitute a particular evolutionary scenario with an independent demographic history and linked to different ecological conditions. Indeed, both morphological and genetic differentiation of P. lilfordi populations have led to consider them as independent evolutionary significant units [72]. In this regard, our results suggest that the allopatric status of the Balearic Podarcis populations could also have shaped their gut microbiota. Our exploratory PER-MANOVA analyses retrieved the islet adscription and main island district as the main factor in explaining bacterial community composition. Consistently, each individual population showed a high proportion of non-shared ASVs (e.g., up to 32.77% of the ASVs found in the Porros population were exclusive from this islet). This result is even more significant if we consider that rarefaction curves of the samples reached the plateau stage, thus ensuring that comparisons are based on a comprehensive view of their respective ASV diversity. All these results reinforce the consideration of every single Podarcis population as a unique evolutionary entity hosting a singular gut microbiota. PERMANOVA analyses based on datasets confirmed by samples from populations with matching collecting seasons: (spring + summer and summer + autumn) allowed us to estimate the effect of the seasons on the gut microbiota of the Balearic Podarcis. In general, samples from the same population (island/islet) and collecting season showed similar composition in terms of microbial communities. In this regard, the seasonality of food has been demonstrated to partially affect the diet composition of P. lilfordi lizards [73]. Resources availability is usually higher in spring when both arthropod prey and plants are abundant, while summer is characterised by the consumption of vegetal tissues and ants, and autumn is marked by food scarcity [73,74]. Previous studies have shown diet composition affects the lizard's gut microbiota [68]. Our results not only point in the same direction but also highlight the importance of the interaction between the specific islets and the seasons of the year, thus suggesting that the different Balearic Podarcis populations inhabit unique ecological scenarios.
The evolutionary history of Podarcis in the Balearic archipelago is marked by the occurrence of allopatric processes leading to the extant isolated populations that exhibit genetic, morphological, ecological, and ethological differences, in many cases unique to one population [72]. Such isolation is likely to have impeded both the gene flow between populations and the dispersal of gut symbionts, thus favoring the parallel divergence of the hosts and their gut microbiota. This divergence process could have been reinforced by a housing effect similar to that reported from other animal groups where genetically divergent individuals/species inhabiting the same place share more bacteria taxa than individuals living apart [75,76]. In this regard, the consumption of conspecifics (cannibalism) reported in P. lilfordi [17] could facilitate horizontal gut microbiome transmission among sympatric individuals as proposed for other organisms [77]. Environmental bacteria acquired from prey and plant material, or even from coprophagy, are also known to constitute an important fraction of the gut microbiome in lizards [24,78], which is an islet system that could contribute to homogenizing the bacterial communities of the individuals. In addition, vertical transmission of gut microbiota during birth has been demonstrated in other lizard species [79]. In turn, microbial communities could drive host evolution by influencing key aspects such as adaptation to resource utilisation and behaviour [5,6,11], thus consolidating the divergence between lizard lineages/populations. All of the above suggest that allopatric divergence of hosts coupled with both the limited dispersal of gut symbionts and the ecological idiosyncrasy of their isolated habitats could have shaped the faecal microbiota of these two species of endangered lizards of the Balearic Islands. Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This study was possible thanks to the project: CGL2015-68139-C2-1-P, "Dinámica de la variación genética y respuesta adaptativa en las Podarcis insulares," financed by of the Ministerio Español de Economia y competitividad and European Regional Development Fund (ERDF). IA supported by FPI CAIB 2017 and research funds from the Conselleria d'Educació, Cultura i Universitats (Govern de les Illes Balears, Spain), co-financed by the ERDF.
Data Availability Raw sequences are available in the Sequence Read Archive (SRA) database at NCBI under BioProject ID PRJNA693423.

Declarations
Ethics Approval Not applicable.

Conflict of Interest
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.