Phylogenetic, Microbiome, and Diet Characterisation of Wall Lizards in the Columbretes Archipelago (Spain): Clues for Their Conservation

Abstract

The Columbretes archipelago consists of a group of small volcanic islets located in the western Mediterranean near the east of the Iberian Peninsula. Four of its islands are inhabited by the wall lizard Podarcis liolepis, whose populations have been considered vulnerable. The purpose of this study was to assess the level of genetic diversity across the archipelago and the possible evolutionary origin of the Columbretes lizard populations. Additionally, we investigated the evolutionary ecology of these populations using a DNA-based metabarcoding approach to characterise both their microbiota and trophic interactions. The genetic results reported very low genetic diversity and corroborated the conspecificity between insular populations and P. liolepis from the mainland (Peñagolosa region). The results of the metabarcoding analyses based on faecal samples were in accordance with an omnivorous ecology, suggesting that specific microbiota communities in the insular populations might be correlated with differences in host ecology and phylogeny. These results are a valuable contribution to the current understanding of the evolution of Columbretes’ lizards and provide important information for conservation management.

1. Introduction

The Columbretes archipelago consists of a group of small islets occupying a total area of 19 Ha in the western Mediterranean Sea. These islets are relatively close to the eastern Iberian Peninsula (50 km) and to Ibiza Island (Balearic Archipelago) (100 km). Its volcanic origin is dated between 1–0.3 Ma [1,2], which corresponds to the period after the Pleistocene volcanic episodes. Possible land connections could have existed between the archipelago and the mainland during the Quaternary glaciations, when the Mediterranean Sea suffered fluctuations in sea level, since the channel between the Iberian Peninsula and the Columbretes has a depth of approximately 90–100 m [3,4].

The archipelago comprises four different groups of islets (Figure 1): (1) Grossa (13 ha), Mascarat, Senyoreta, and Mancolibre; (2) Ferrera (1.5 ha), Espinosa, Bauzá, Valdes, and Navarrete; (3) Foradada (1.6 ha), Lobo (0.5 ha), and Méndez Nuñez; and (4) Carallot (0.1 ha), Cerquero, Churruca, and Baleato. Until 1982, some islands (Foradada, Ferrera, and Carallot) were used for military exercises with live ammunition, and these activities may have had negative impacts on the flora and fauna of the islands [5].

The terrestrial vertebrate fauna consists of some breeding bird species and the endemic lizard first named as Lacerta muralis atrata by Boscá (1916), even though there was a previous citation under the name Podarcis muralis fusca (Salvator 1895; see more details in [6]). The modern name employed is Podarcis hispanica atrata (Boscá, 1916), and later on was considered a new species, named Podarcis atrata [7,8]. Recently, it has been considered conspecific with Podarcis liolepis (Boulenger, 1905), a common species in the north-eastern Iberian Peninsula [9,10,11,12,13,14,15]. In the eastern area of the Iberian Peninsula, three new lineages have been described within the so-called Podarcis hispanicus (Steindachner, 1870) complex: Podarcis galerai, Podarcis hispanicus sensu stricto (Albacete/Murcia), and P. hispanicus of Valencia lineage [11,15,16,17]. The Columbretes lizard populations are present only on four of the islets in the archipelago: Grossa, Mancolibre, Foradada, and Lobo. It has been postulated that their ancestors may have arrived on the Columbretes Islands from the mainland during the Quaternary glaciations [7,8], when the Mediterranean Sea experienced level changes [3,4].

The fauna of the archipelago suffered a dramatic episode during the construction of the lighthouse of Grossa Island in 1855, with a drastic campaign for the extermination of the Lataste’s Viper, Vipera latastei, that involved a bushfire on Grossa Island. The effects of the intentional fire used to eliminate the viper populations on lizard populations have yet to be studied. Fortunately, the Columbretes Islands were declared a Natural Park in 1988, and a marine reserve in 1990. Although the conservation status of this insular species of Podarcis is still not well defined due to its recent taxonomic change, it is considered vulnerable by national institutions [5,18]. The high fragmentation of the different insular populations, the reduced habitat, and the low number of individuals have been highlighted as characteristics that make this species/population susceptible to extinction [19,20].

There are few studies on the biology of these insular populations that provide information about their phylogenetic origin, their trophic interactions, or the composition of their microbiota. The only available study about the genetic diversity of the lizards from the Columbretes archipelago dates from 1998 [7], and this is only based on a 306 bp fragment of the mitochondrial gene encoding the protein cytochrome b, amplified by means of universal primers [21]. Although this study refers to the insular Podarcis as a separated species, Podarcis atrata, according to genetic differentiation, it was suggested that populations from Columbretes Grossa and Mancolibre are less differentiated than those of Foradada and Lobo and seem to have retained mainland haplotypes whose origins have not been clearly established [7].

Dietary studies are essential for understanding ecosystem functionality [22], the contribution of species to food webs, and the relevance of food resources for their demographic viability [23]. Feeding habits are expected to vary among species depending on their evolutionary history, the characteristics of their microhabitats, and prey availability [24], as well as the foraging strategies they adopt [25]. Lizards have been considered to be generalists in relation to their diet [26]; however, the foraging strategy in island ecosystems is different from those in continental areas mainly due to food-resource scarcity, which promotes the exploitation of marine trophic resources [27,28]. Such use of marine subsidies is very evident in the case of lizard species inhabiting extremely dry coastal areas [29,30]. Little is known about the feeding behaviour of P. liolepis on the Columbretes Islands, but previous studies based on the direct observations of feeding activity or the microscopic analysis of faecal contents [18,31,32] have considered them omnivorous [31], including cannibalism of both eggs and juveniles [33], which is very frequent in small Mediterranean islets [34]. However, these methodological approaches are known to have limitations, because they are difficult to conduct, require a high level of training, and usually contribute a biased and incomplete view of the trophic spectrum [34,35,36,37].

Molecular approaches, such as DNA barcoding, have emerged as a promising tool to investigate complex trophic interactions [38,39]. DNA barcoding also represents an important advance in the molecular identification of microbial taxa, allowing the characterization of microbial communities in complex environmental samples and providing new insights into their responses to environmental factors [40]. It is well established that the complex microbial communities living in the gastrointestinal tract of animals can have a significant impact on host diversification [41], resource adaptation [42], and behaviour [43], as well as other factors. Moreover, there is increasing evidence of the effects of dietary behaviour and habitat on shaping gut microbiota [44,45,46].

Our purpose was to assess the level of intra-archipelago genetic diversity, the origin of the founders of the Columbretes lizard populations, and the genetic differentiation between the lizard species from Columbretes and their sister taxa from the mainland. Furthermore, we investigated some aspects of the ecology of these insular lizard populations through a DNA-based characterization of both their microbiota and trophic interactions from faecal samples. Our findings will improve our knowledge on the Columbretes’ Podarcis populations and could be useful for understanding the relationship between phylogenetic and ecological differences with distinct microbial composition. The results and the genetic resources presented here will provide important information for official conservation bodies.

2. Materials and Methods

2.1. Phylogenetic Analyses

2.1.1. Sampling

Lizards were carefully noosed in their habitats. A total of 32 individuals from Columbretes Grossa, 4 from Foradada, 2 from Lobo and 2 from Mancolibre were sampled in July of 2009. Tail tips were removed and immediately stored in 100% ethanol. All lizards were released at the point of capture. Samples and their locations are described in Table S1 and Figure 1. The field protocols for the capture, handling, and release of lizards were approved by Conselleria de Medi Ambient of the Government of Valencia (REF 2009/143-jvee). This study was approved by the ethical guidelines of the University of the Balearic Islands and Salamanca.

2.1.2. DNA Amplification and Sequencing

DNA extraction was performed by the phenol-chloroform method following the standard protocol described in Sambrook et al. [47]. A mtDNA fragment providing an alignment of 2311 bp length was obtained for 23 specimens from the Columbretes Islands, including partial 12S rRNA, partial cytochrome b (CYTB), partial control region (CR), and two partial subunits of the NADH dehydrogenase gene and associated tRNAs (referred to as ND1, ND2, tRNA_Ile, tRNA_Gln, and tRNA_Met). Primers and amplification conditions were the same as those used in our previous studies on Podarcis [48,49,50,51]. In addition, mtDNA sequences of 56 individuals were retrieved from GenBank [17] (Table S1): 11 from Columbretes; 19 P. liolepis (from Peñagolosa, South of France and La Rioja); 20 belonging to the so-called P. hispanicus complex from Valencia, Albacete/Murcia, and Galera lineages; 6 Podarcis muralis (Laurenti, 1768) (from Asturias, France, Castilla y León, Cantabria and Castellón).

Additionally, a cytochrome oxidase I (COI) fragment (657 bp) was amplified by means of primers LCO1490: GGTCAACAAATCATAAAGATATTGG and HCO2198: TAAACTTCAGGGTGACCAAAAAATCA [52]. The amplification reaction consisted of 35 cycles after initial denaturalization at 94 °C for 4 min. The cycles comprised 30 s at 94 °C, 30 s at 48 °C, and 1 min at 72 °C, followed by a final extension at 72 °C for ten minutes. A total of 24 individuals were studied in this case (Table S1): 10 of them from Columbretes (including 1 from Lobo, Mancolibre and Foradada); 10 P. liolepis from the Peñagolosa area; one from Rioja; three from the south of France (Pêch de Foix and Vaychis). Two COI sequences were also used from GenBank, one P. liolepis from Catalonia, Palamós (accession number MN015111.1) and one P. muralis from Panticosa, Spain (accession number MN015063.1).

Furthermore, two partial nuclear genes were amplified and sequenced: the apolipoprotein B gene (APOBE28) (25 samples), and the transcription factor gene KIAA2018 (24 samples). The primers and conditions are described in our previous works [51,53] (Table S1). The sequences of two other nuclear markers, melanocortin 1 receptor (MC1R) (665 bp) and recombination activating gene 1 (RAG1) (939 bp), were obtained from GenBank [17] and were also included.

Both strands of all PCR products were sequenced using the BigDye^® Terminator v.3.1 Cycle sequencing kit (Applied Biosystems, Foster City, CA, USA) and then genotyped in an automated ABI 3130 sequencer (Applied Biosystems). The sequences were edited using BioEdit v.7.0.5.2. (Bioedit Company, Manchester, UK) [54] and aligned with the MAFFT online server [55] using the iterative refinement method (FFT-NS-i). For the protein-coding genes, the alignments were verified by translating nucleotide sequences into amino acids. All GenBank accession numbers are indicated in Table S1.

2.1.3. Divergence and Phylogenetic Analyses

The basic genetic-diversity parameters were calculated with DnaSP v.6 software (University of Barcelona, Spain) [56] for the concatenated mitochondrial alignment, the COI fragment alone, and each phased nuclear gene based on Columbretes Islands samples only. DnaSP was also used to calculate Tajima’s D [57] neutrality statistic, which contrasts estimates of θ based on segregating sites (S) and pairwise differences (k) to determine deviation from selective neutrality. In addition, pairwise mismatch distribution to test for population expansion [58] was also carried out using DnaSP.

Best-fit nucleotide-substitution models and partitioning schemes were simultaneously chosen using PartitionFinder v.2.1.1 [59] under the corrected Akaike Information Criterion (cAIC) for the mtDNA dataset (Table S1). The set of models was restricted to those available in BEAST, and the partitioning schemes were manually defined through the “user” option available in the software. Branch lengths of alternative partitions were estimated as “unlinked” to search for the best-fit scheme, which consisted of three partitions: non-coding fragments [GTR + I], 1st and 2nd codon positions of coding regions [HKY + I + G], and 3rd codon positions of coding regions [GTR + I].

We performed phylogenetic analyses using maximum-likelihood (ML) and Bayesian inference (BI) methods based on the mtDNA dataset along with DNA sequences from seven Podarcis from the Balearic Islands, four Podarcis lilfordi (Günther, 1874), and three Podarcis pityusensis (Boscá, 1883) as an outgroup.

Maximum-likelihood analyses were performed using IQ-TREE v.2.0.3 [60]. We applied the partitions and the best-fit substitution model and performed 10⁶ bootstrap replicates based on the ultrafast bootstrap approximation (UFBoot) [61,62] for statistical nodal support.

Bayesian analyses were performed with MrBayes 3.2.1 [63]. The analyses were run for 10⁷ generations with sampling frequency every 10³ generations. Numbers of independent runs and chains were set to default, two, and four, respectively. A sufficient number of generations was confirmed by examining the stationarity of the log-likelihood (lnL) values of the sampled trees and the value of the average standard deviations of the split frequencies being lower than 0.01. The results were analysed in Tracer v.1.7 [64] to assess the convergence and effective sample sizes (ESS) for all the parameters. A burn-in of 25% was applied, and the phylogenetic trees were visualised and edited using Figtree v.1.4.2 (available from https://github.com/rambaut/figtree/releases), accessed 17 July 2019 [65].

A TCS statistical parsimony network approach [66] with 95% connection limit implemented in the program PopART 1.7 (http://popart.otago.ac.nz, accessed 17 July 2019 [67,68]) was used to infer (i) genealogical relationships between P. liolepis from the Columbretes Islands and from the mainland based on the haplotypes of phased nuclear genes (RAG1, MC1R, KIAA2018 and APOBE28) and mtDNA alignment, using P. muralis as an outgroup and (ii) intra-island diversity based on the RAG1 nuclear gene and mtDNA alignment.

BEAST v.2.5 [69] was used to simultaneously infer the phylogenetic relationships and divergence times between Columbretes’ Podarcis and different lineages of the Iberian Peninsula based on mtDNA-sequence datasets (72 individuals) and seven Podarcis from Balearic Island, which were used for calibration. The calibration was specified from a normal distribution (5.32, 0.01). This calibration was based on knowledge of the timing of the end of the Messinian salinity crisis (5.33 Ma) and the very rapid refilling of the Mediterranean basin that would have separated the two Balearic Island Podarcis (i.e., P. lilfordi and P. pityusensis; see [70]). Partitions and evolutionary models were the same as those used for MrBayes. The *BEAST MCMC sampler was run twice for 5 × 10⁸ generations, sampling every 5000 steps. A relaxed log-normal clock model was specified, and a Yule model was used as the tree prior.

2.2. Metabarcoding Faecal Microbiota and Diet Associations

2.2.1. Sampling

Eight fresh faecal samples from four females and four males of P. liolepis lizards were collected in July of 2018 from the Columbretes Grossa Island for both the metabarcoding analysis of diet content and the molecular characterization of their microbiota. Faecal sampling has been used to analyse the microbiome of P. liolepis, as previous studies in reptiles [71,72,73] have shown that it provides a complete understanding of their hindgut bacterial communities [74], without the need to sacrifice the animal.

Specimens were noosed and fresh faeces were obtained through a gentle abdominal massage and stored directly from lizards’ cloaca in a sterile tube with absolute ethanol. The lizards were immediately released at the capture site. The samples were preserved on ice in the field and then stored at −20 °C upon arrival at the laboratory until DNA extraction. These samples were not obtained from the same specimens used in the phylogeographic analysis.

2.2.2. DNA Extraction, DNA-Library Preparation, and Sequencing

DNA extractions from individual samples were obtained using the Isolate Faecal DNA kit (Bioline, London, UK). Samples were submitted to the Roy J. Carver Biotechnology Center (University of Illinois, Urbana-Champaign, USA) for amplification in a microfluidic high-throughput multiplexed PCR platform (Fluidigm©, San Francisco, CA, USA).

For the amplification of bacterial 16S rDNA, we used the primer set 515F/806R (5′-GTGCCAGCMGCCGCGGTAA-3′/5′-GGACTACHVGGGTWTCTAAT-3′) [75], targeting the V4 hypervariable region. For animal-prey detection, we used two different primer pairs targeting the mitochondrial cytochrome oxidase I gene (COI): the universal pair mlCOlintF/jgHCO2198 (5′-GGWACWGGWTGAACWGTWTAYCCYCC-3′/5′-TANACYTCNGGRTGNCCRAARAAYCA-3′) [76,77] and the arthropod-specific pair ArtF11/ArtR17 (5′-GGNKYNGGNACWGGATGAACWGTNTAYCCNCC-3′/5′-GGRTCAAAAAATGAWGTATTHARATTTCGRTCWGTTA-3′) [78]. For plant-diet detection, we selected two primer sets targeting chloroplast regions: the large subunit of the RuBisCO gene (rbcL), rbcLa_F/rbcLa_R (5′-ATGTCACCACAAACAGAGACTAAAGC-3′/5′-GTAAAATCAAGTCCACCRCG-3′) [79,80], and the psbA-trnH intergenic spacer, psbA3_f/trnHf_05 (5′-GTTATGCATGAACGTAATGCTC-3′/5′-CGCGCATGGTGGATTCACAATCC-3′) [81]. The universal CS1 and CS2 Fluidigm tags along with barcode labels specific to each sample and Illumina adapters i5 and i7 were used, and the resulting amplicons were validated on a Fragment Analyzer (Agilent, Santa Clara, CA, USA) using the HS NGS kit (DNF-474-33). Sequencing was conducted on an Illumina MiSeq v2 platform (Illumina, San Diego, CA, USA), yielding 2 × 250 paired-end reads.

2.2.3. Bioinformatic Analyses

The dietary reads were merged, trimmed, and filtered with Micca v.1.7.2 [82]. For taxonomic assignment, the resulting DNA sequences (OTUs) were combined into 98% similarity clusters with Usearch v.10 [83]. Centroid (representative) sequences were selected for each cluster and subsequently used as BLAST queries to retrieve the top 1000 hits from GenBank. Independent DNA matrices were built for each of the OTU clusters by combining their dietary sequences with the top 1000 hits from their respective BLAST searches. The resulting matrices were aligned with MAFFT [55] and used for maximum-likelihood phylogenetic inference in IQ-TREE [60]. The trees were explored in FigTree [65] to determine the systematic position of the diet sequence according to the highest taxonomic rank, supported by a nodal bootstrap value ≥ 70% [84].

OTUs from multiple markers targeting the same group of organisms (i.e., animals or plants) were combined into a unique list using the python script described in da Silva et al. [85].

The microbial sequences were analysed with QIIME2 v.2020.2 (available from http://qiime.org/), accessed 17 July 2019 [86]. The sequences were demultiplexed and subsequently filtered and denoised using the DADA2 pipeline [87]. Amplicon sequence variants (ASVs = 100% identity OTUs) were taxonomically classified against the SILVA132 database [88]. OTUs were identified as chloroplasts or mitochondria, and those with undetermined phylum annotations were excluded from the dataset. The remaining sequences were aligned with MAFFT [55] under a default FFT-NS-1 algorithm and a phylogenetic tree was inferred with FastTree [89] using the QIIME2 plugin align-to-tree-mafft-fasttree.

A microbial OTU table was rarefied using the rarefy_even_depth option in the R package phyloseq [90] based on the sequencing depth of the sample with the minimum number of reads (sample KK310; 24,981 reads).

3. Results

3.1. Genetic Diversity of Columbretes Islands’ Podarcis Based on Mitochondrial and Nuclear Genes

Only eleven polymorphic sites out of the 2311 bp of the mtDNA fragment sequenced were detected in the 28 genotyped specimens from the Columbretes archipelago, corresponding to an extremely low diversity (pi = 0.0009). We detected 10 haplotypes, eight in Columbretes Grossa, one present in Mancolibre, and one haplotype shared between Foradada and Lobo (

Table 1. Genetic-diversity parameters and neutrality tests for each nuclear gene and for mtDNA alignment and COI genes from the Columbretes archipelago samples (S, polymorphic positions; h, number of haplotypes; Hd, haplotype diversity; π, nucleotide diversity; k, average nucleotide differences).

	N (N Phased)	bp	S	h	Hd	π	k	D Tajima (1989)
Nuclear genes
RAG1	26 (52)	939	12	12	0.704 ± 0.045	0.00192 ± 0.00015	1.802	−0.572 n.s.
MC1R	12 (24)	672	12	9	0.768 ± 0.078	0.00586 ± 0.00067	3.069	−0.154 n.s.
KIAA2018	10 (20)	489	3	3	0.563 ± 0.063	0.00230 ± 0.00026	1.126	0.886 n.s.
APOBE28	26 (52)	665	10	15	0.882 ± 0.027	0.00540 ± 0.00020	3.593	1.770 n.s.
mtDNA
CYTB, NADH, CR, 12S	28	2311	11	10	0.770 ± 0.076	0.00091 ± 0.00017	2.106	−0.833 n.s.
COI	10	657	1	2	0.200 ± 0.015	0.00030 ± 0.00023	0.200	−1.111 n.s.

n.s.: not significant.

). In relation to the COI dataset alone, the ten individuals sequenced (including seven from Columbretes Grossa and one each from Mancolibre, Foradada and Lobo) shared the same haplotype apart from one nucleotide position in an individual from Columbretes Grossa. This low variability was also reflected in the low values of haplotype (Hd = 0.200) and nucleotide (pi = 0.0003) diversity.

Four partial nuclear genes (RAG1, KIAA2018, APOBE28 and MC1R) were sequenced, comprising a total of 2765 bp. They also showed low values of genetic diversity (Table 1) in terms of both haplotype diversity (Hd: 0.563–0.882) and the average of nucleotide differences (k: 1.126–3.593). The APOBE28 gene was the most variable gene with 15 haplotypes (phased), followed by RAG1 with 12 haplotypes, MC1R with 9 haplotypes, and KIAA2018, the least diverse, which only accounted for three haplotypes. Tajima’s D test indicated negative but non-significant deviations from neutrality for mtDNA fragments and for two nuclear genes (RAG1 and MC1R), but positive values for the other nuclear genes. The pairwise mismatch distribution analysis showed a pattern that suggests a growth–decline model rather than a constant-size model (Figure S1).

3.2. Phylogenetic Relationships and Possible Origin of Columbretes Islands Wall Lizard

The phylogenetic relationships between the different islands of the Columbretes archipelago are shown in Figure 2 based on the mtDNA database (Figure 2a) and the nuclear gene RAG1 (Figure 2b) as an example. In both, it is interesting to highlight that Foradada and Lobo share the same haplotype, while the nuclear and mtDNA haplotype of Mancolibre is closer to Columbretes Grossa. These results can be explained by the geographical position of these islands, given that Foradada and Lobo, on the one hand, and Grossa and Mancolibre on the other hand belong to different and distant groups of islands. A maximum of three mutational steps separate both groups of islands. Similar patterns were observed in the three other nuclear genes sequenced.

In our study, samples from the Columbretes Islands form a separate lineage that in turn is nested in a clade that also includes P. liolepis from the mainland, specifically from Peñagolosa (Castellon, Spain) (Figure 3a). The same pattern was revealed in the TCS analysis based on the COI fragment (Figure 3b), where P. liolepis samples from France and from La Rioja and Catalonia in Spain also form a separate cluster with respect to the Columbretes and Peñagolosa samples.

This pattern was more evident in the analyses based on mtDNA, both for phylogenetic trees (ML and BI) (Figure 4) and TCS (Figure 3), than the patterns shown by nuclear genes (Figure 5), indicating that the P. liolepis population from the Columbretes Islands shares some haplotypes with the conspecific mainland population, confirming the genetic proximity between Columbretes and Peñagolosa lizards.

Time-calibrated phylogeny based on mtDNA (Figure 6) suggested a time of divergence between Foradada + Lobo and Columbretes Grossa + Mancolibre at 1.59 Ma (95% highest posterior density, HPD: 2.95–0.52 Ma), coinciding with glacial and interglacial periods during the Pleistocene. The split between the Columbretes Islands’ population and the mainland (Peñagolosa) populations would have also taken place during this period of glaciation (1.77 Ma, 95% HPD: 3.10–0.69 Ma).

3.3. Molecular Diet Composition

Integrating data from multiple taxonomically overlapping markers resulted in a total of 22 diet taxa (12 animals and 10 plants) that could be identified to 6/4 orders and 9/7 families of animals and plants, respectively (

Table 2. List of identified diet taxa to the maximum resolution obtained in the diet composition of Podarcis liolepis from Columbretes Grossa through molecular analysis.

Phylum	Class	Order	Family	Genus	Species	Number of Faecal Samples
Arthropoda
	Arachnida	Araneae	Araneidae			1
	Insecta	Coleoptera				1
			Curculionidae	Trachyphloeus		1
			Tenebrionidae			1
		Hemiptera	Rhyparochromidae	Lamprodema	maurum	1
		Hymenoptera	Formicidae	Pheidole		1
		Lepidoptera				1
			Blastobasidae	Blastobasis		1
			Geometridae			2
			Psychidae			2
	Malacostraca	Isopoda				3
			Armadillidae	Armadillo	officinalis	2
Streptophyta
	Magnoliopsida	Apiales	Apiaceae	Daucus		1
		Caryophyllales	Chenopodiaceae			2
				Beta	vulgaris	1
				Patellifolia		1
				Suaeda		5
			Plumbaginaceae	Limonium		6
			Portulacaceae	Portulaca		2
		Rosales	Rosaceae	Prunus		3
			Ulmaceae	Ulmus		1
		Solanales	Solanaceae	Lycium		1

). Lepidoptera, Isopoda, and Coleoptera were the most common orders identified in terms of prey (Figure 7 and Table 2). At the family level, Armadillidae (isopods), Geometridae, and Psychidae (Lepidoptera) were the most abundant prey (Figure 7). Regarding vegetal food, our study revealed the presence of plant DNA from the phylum Streptophyta in faecal samples (Table 2). Caryophyllales (Figure 7) stands out as the most abundant order, dominated by the families Plumbaginaceae, Chenopodiaceae, and Portulacaceae (Figure 7).

We also identified some diet OTUs at the genus and even the species level. This was the case with Trachyphloeus sp. (Coleoptera), Lamprodema maurum (Hemiptera), Pheidole sp. (Hymenoptera), and Blastobasis sp. (Lepidoptera) in terms of prey, and Daucus sp. (Apiales), Beta vulgaris, Patellifolia sp., Suaeda sp., Limonium sp., Portulaca sp. (Caryophyllales), Prunus sp., Ulmus sp., (Rosales), and Lycium sp. (Solanales) in terms of plants (Table 2).

3.4. Microbiota Analysis

After quality filtering, the dataset was reduced to 482,404 high-quality reads with an average of 60,300.5 reads per sample (range = 24,981–116,862). The reads could be ascribed to 968 OTUs, representing 17 phyla, 31 classes, 71 orders, 116 families, and 209 genera of microbiome taxa.

Two phyla represented altogether almost 82% of the microbiota: Bacteroidetes (average = 44.9; range 32.8–71.4%) and Firmicutes (average = 37.0; range 14.0–45.2%) (Figure 8). Faecal microbiota was dominated by two orders (Figure S2): Bacteroidales (average = 44.9%; range = 32.9–71.3%; Bacteroidetes) and Clostridiales (average = 29.2%; range = 13.0–42.0%; Firmicutes). Within Bacteroidetes, the Bacteroidaceae family showed the highest relative abundance (average = 20.8%; range = 18.1–24.8%), followed by Rikenellaceae (average = 8.4%; range = 3.7–16.2%), Marinifilaceae (average = 6.9%; range = 2.9–19.2%), and Tannerellaceae (average = 6.6%; range = 3.5–14.6%). Regarding Firmicutes, the families Lachnospiraceae (average = 7.9%; range = 2.1–15.5%), Erysipelotrichaceae (average = 6.3%; range = 0.2–16.8%), Peptostreptococcaceae (average = 4.1%; range = 0.0–22.6%), and Ruminococcaceae (average = 3.2%; range = 1.5–4.9%) were found to be the most abundant.

4. Discussion

The very low genetic diversity observed using all datasets in lizard populations from the Columbretes Islands drastically contrasts with the genetic-diversity values obtained by Castilla et al. [7], which observed 95 polymorphic sites, 85 of which were parsimony informative, in a fragment of 306 bp of the CYTB gene. We compared our CYTB sequences with those of 14 haplotypes from Columbretes Islands deposited in GenBank by Castilla et al. [7] (accession numbers AJ004987, AJ004994–AJ004996, AJ004910–AJ004911, AJ004990–AJ004992, AJ224407–AJ224409, and AF052636–AF052637), resulting in a pattern of variation with numerous changes alternating every three nucleotides. This discrepancy between genetic diversity was reported in [7] and our results are not explained by differences in sample size, but could be in accordance with a non-specific amplification of the mtDNA CYTB gene, probably because Castilla et al. [7] utilised less specific primers [21,91] and we used primers that were specific for amphibians and reptiles, as described by Palumbi [92].

Thus, we did not confirm the mtDNA variability previously described in Podarcis from the Columbretes archipelago [7] but rather its extreme homogeneity, indicating a possible recent founder effect at the maternal level or the passage through a recent bottleneck of the population (Table 1). The results obtained in the pairwise mismatch distribution analysis correspond more to a growth–decline population model than to a constant-size population model [93] (Figure S1). This findings together with the star-like haplotype networks, the higher values of haplotype diversity than nucleotide diversity, as well as the tendency indicated by the negative values of Tajima’s D neutrality statistics suggested a population expansion [93] (Table 1 and Figure 2).

The P. liolepis population from the Columbretes Islands shares many haplotypes with the conspecific mainland population, confirming the genetic proximity between Columbretes and Peñagolosa (Castellón) lizards. These results are not in agreement with those obtained by Castilla et al. [7,8], which elevated Columbretes lizards to the species rank based on the high genetic divergence found between them and the samples from the mainland (Valencia). Considering the current deeper knowledge of lizards from the southeast region of the Iberian Peninsula [17], we suggest that Castilla et al. [7] compared Columbretes samples with populations from one of the other lineages of the P. hispanicus complex that inhabit this area with the Valencia lineage and not with specimens of P. liolepis. Thus, we can confirm the inclusion of Columbretes’ lizards within the P. liolepis species and great similarity with individuals from Peñagolosa, especially at the mitochondrial level. The distribution of P. liolepis in the Iberian Peninsula is not yet well defined, so it would be necessary to obtain more samples from more locations to obtain a better picture of P. liolepis lizard populations.

Time-calibrated phylogeny based on mtDNA (Figure 6) provided a time of divergence between Foradada + Lobo and Columbretes Grossa + Mancolibre at 1.59 Ma (95% HPD: 2.95–0.52 Ma), coinciding with an interval of time of several sea-level fluctuations between glacial and interglacial periods [94]. A time-calibrated tree also shows that the Columbretes Islands population, once separated from the mainland (Peñagolosa) populations (~1.77 Ma, 95% HPD: 3.10–0.69 Ma), had no further contact with them and therefore diverged in isolation. Several events leading to a decrease in diversity (bottleneck) followed by a possible population expansion could have occurred, as is corroborated by the growth–decline analysis performed with DnaSP (Figure S1), in the demographic history of Columbretes’ lizards. One of these demographic events may have been caused by the provoked fire that occurred in the 19th century. Genome-wide approaches would help to confirm the phylogenetic relationships and address the results of the expansion analysis.

In this study, we present the first approach to the molecular characterisation of the eating habits of P. liolepis on Grossa Island while avoiding sacrificing specimens and reinforcing that it is not necessary to kill lizards to characterise their diet. This is a key point to consider when dealing with populations subjected to habitat fragmentation and accounting for low numbers of individuals, as may be the case for the wall lizards from the Columbretes archipelago. Our sampling sizes were directly related to lizard abundances on different islets of the Columbretes Archipelago. Only on Grossa Island were we able to estimate a lizard density with line transects (Pérez-Cembranos et al., unpublished data). In Foradada, Mancolibre, and Lobo, lizards were present but with very reduced population sizes (see also, [6]). Thus, even from a conservation viewpoint we were obliged to obtain very small sample sizes.

Most lizard species are considered to feed mainly on invertebrates [95,96]. However, numerous species will consume plant matter, marine resources or even practice cannibalism under conditions of prey scarcity or in extremely dry habitats [30,31], a behaviour that is most frequently observed in island species (see [96] and the references therein). In studies on the diet of the P. hispanica complex from the Iberian Peninsula, there are no mentions of the use of vegetal matter in the diet [97]. Although more studies are needed to characterise the specific microbiota of P. liolepis populations from the mainland, neither of the two current studies on this species mention plant consumption [98,99]. Available dietary studies based on morphological observations of faecal pellets and stomach contents, and on direct observations of the feeding behaviour of P. liolepis from Columbretes Islands, evidenced the consumption of Coleoptera, Lepidoptera, Araneae, Gastropoda, scorpions (Buthus occitanus), as well as plant material and fruits (Opuntia). Some studies even found that some faecal pellets exclusively contained vegetable matter [5,18,31,100].

Our molecular results support the omnivorous behaviour of P. liolepis from the Columbretes Islands, yielding a similar diet composition that shows a preference for Lepidoptera, Coleoptera, and Isopoda. However, we were unable to detect cannibalism (as we have not been able to differentiate between host sequences and cannibalised conspecifics), nor the consumption of scorpions or specific marine prey. Although the marine isopod Ligia italica Fabricius 1798 is very abundant in the Columbretes, it is very difficult to observe these lizards feeding on the vertical, eroded and inaccessible cliffs where these crustaceans live [100]. Castilla et al. [100] reported the opportunistic consumption of marine Isopoda, L. italica, in the Columbretes lizards’ populations under experimental conditions. Our results did not identify specific marine prey, but the order Isopoda is one of the most abundant and some counts have not been identified at the species level.

According to our results, we can say that metabarcoding allows for higher taxonomic precision in the identification of prey than macroscopic methods [101], identifying 16 families, 14 genera, and even 3 species. However, further samples from different seasons, islands/islets, and locations within island/islets should be studied to obtain a comprehensive record of the diet of P. liolepis in the Columbretes archipelago over temporal and spatial scales.

The first characterization of P. liolepis’ faecal microbiota was reported in this study. Although the sample size was reduced, the results were consistent with previous reports in other omnivorous lizards [74], since around 80% of the microbiota found in P. liolepis samples belong to phyla Firmicutes and Bacteroidetes. As mentioned above, P. liolepis has omnivorous feeding habits and in order to process plant polymers found in celluloses and hemicelluloses, these organisms need endogenous glycoside hydrolases [102]. Podarcis liolepis shows a high prevalence within the Firmicutes of the families Lachnospiraceae and Ruminococcaceae, both of which are related to decomposed complex plant material [103]. Furthermore, within Bacteroidetes phyla, P. liolepis showed a high proportion of Bacteroidaceae followed by Rikenellaceae families, both previously reported from the gut microbiota of herbivorous reptiles [73,104,105]. All of these findings indicate that the microbiome composition of P. liolepis lizards is consistent with their omnivorous trophic ecology.

5. Conclusions

Lizard populations from the Columbretes Islands are characterised by their extraordinarily low genetic diversity based on mtDNA and nuclear markers. This finding could be the result of the effect of a bottleneck or founder effect, which are characterised by a drastic reduction in the population size due to environmental variables or evolutionary history. Phylogenetic analysis corroborated that Podarcis populations that inhabit the Columbretes archipelago are conspecific with P. liolepis, specifically those in the Peñagolosa region in Castellón. The divergence time between the insular and mainland forms was dated to 3.10–0.69 Ma, coinciding with a period of several sea-level fluctuations during Pleistocene glaciations, when land connections possibly existed. This work presents, for the first time, the microbial diversity and diet composition of these insular lizard populations based on faecal samples by means of a DNA-metabarcoding approach. The obtained results are in accordance with an omnivorous trophic ecology, with Lepidoptera, Isopoda, and Coleoptera being the most common orders of prey. The consumption of plant matter in their diet was also observed, supporting the behaviour that is most frequently observed in island species. The combination of the results provided by both metabarcoding and phylogeographic analysis suggested that the specific microbiota found in the insular populations of P. liolepis might be correlated with differences in host ecology and phylogeny. The microbiota communities found in faecal samples are in accordance with an omnivorous diet including plant matter, a characteristic feeding strategy that could be the result of adaptation to resource scarcity, which is a situation commonly observed in islands. However, future research is required to determine microbiota and trophic interactions of their sister taxa from the mainland to better understand the relationship between the evolutionary history and the microbiota diversity of these lizard populations. Finally, all the genetic results reported in this study contribute to increasing our knowledge of the Columbretes Islands lizard populations and provide important information for conservation management.