Morphological and genetic diversity of the Balearic lizard, Podarcis lilfordi (Günther, 1874): Is it relevant to its conservation?

To characterize the genetic and morphological diversification of the endangered Balearic lizard Podarcis lilfordi and to assess the relevance of this diversity, and how it is described, to conservation measures.


| INTRODUC TI ON
For a variety of taxa, islands make a very important contribution to biodiversity, out of proportion to their land area in comparison with continents. Island contribution to global diversity is mainly in terms of endemic forms, instead of species diversity (Whittaker & Fernández-Palacios, 2007). This is the case of Balearic Islands, with only three species of autochthonous terrestrial vertebrates, one midwife toad and two lacertid lizards, all of them endemic to Balearic Islands.
According to IUCN assessment, the Balearic lizard, Podarcis lilfordi (Günther, 1874), is an endangered species (Pérez-Mellado & Martínez-Solano, 2009) of Mediterranean wall lizard endemic from the Balearic Islands (Western Mediterranean, Spain). P. lilfordi is a member of the first post-Messinian faunal assemblage reported from Menorca Island (Bover, Quintana, & Alcover, 2008;Bover et al., 2014). During the Holocene, around 2,000 years ago, the spe- Phylogeographical studies indicate that the Balearic lizard separated from its sister species, Podarcis pityusensis (Boscá, 1883), after the refilling of the Mediterranean basin, at the end of the Messinian Salinity Crisis (Brown et al., 2008;Terrasa, Rodríguez, et al., 2009). Menorcan and Mallorcan populations appear to have diverged at the beginning of the Quaternary period, 2.6 my ago (Brown et al., 2008). Sea levels were higher during the Late Pliocene than during the present day (Emig & Geistdoerfer, 2004), precluding a land connection between Mallorca and Menorca. The genetic separation of lizards from Mallorca, including Cabrera archipelago, and Menorca was maintained during later periods of the Pleistocene, despite potential for secondary contact that could have taken place during eustatic Pleistocene sea level changes (Emig & Geistdoerfer, 2004). Consequently, during more than 2.5 my, the reciprocally monophyletic clades of Mallorca, Cabrera and Menorca evolved with little or no introgression. Another major divergence occurred around 2 my in Western Mallorca and the rest of populations around Mallorca Island (Brown et al., 2008).
In addition to these major divergences, within the groups of Mallorca and Menorca there has been a subsequent differentiation of isolated populations forming one of the most surprising arrays of allopatric populations. Thus, the Balearic lizard was recognized as a polytypic species (Huxley, 1940in Mayr, 1963. Taxonomic studies to provide formal recognition of this variation began with Von Bedriaga who defined three infraspecific taxa (Bedriaga, 1879).
However, the validity of this subspecific arrangement is largely debatable, especially if we take into account the genetics of these insular populations ). The construction of taxonomies based on phylogenetic species concept could remove the need for subspecific descriptions (Haig et al., 2006;). An alternative solution would be the use of evolutionary significant units (ESUs) that initially were defined as units that should be reciprocally monophyletic for mtDNA alleles and that show a significant divergence of allele frequencies at nuclear loci (Moritz, 1994). We will explore the arrangement of these ESUs in the case of the Balearic lizard and its consistence with the traditional separation of subspecies.
The temporal process of separation during lineage divergence can accumulate genetic, ecological and morphological changes, resulting in a better adaptation to local environmental conditions. The occupation and use of different habitats by lizards would lead to a divergent selection on traits that define several morphological characteristics as body size, body shape, coloration patterns or scalation characters (Hu, Wu, Ma, Chen, & Ji, 2019;Muñoz et al., 2013;Wollenberg, Wang, Glor, & Losos, 2013). We analyse, within the frame of genetic variability, these morphological characteristics in all extant populations of the Balearic lizard. Nuclear DNA markers are analysed for the first time, and we describe the morphometry, scalation and colour patterns of populations. Morphological traits are also related to different climatic conditions at different geographical locations (Mayr, 1963). In Balearic Islands, due to their latitudinal situation, climatic conditions can be different for the two main clades of Mallorca and Menorca populations. Thus, we tested the potential influence of main climatic conditions of Mallorca and Menorca, that is, rainfall and environmental temperatures, on morphological characters of lizards.
Our goal is to answer the following questions: 1. How valid is the current taxonomic arrangement of populations if we consider their genetic and morphological variability? and 2. What is the relationship between the observed morphological diversity and the ESUs are particularly suitable to describe and recognize such diversity and, especially, to ensure the continuity of the evolutionary process.

K E Y W O R D S
Balearic Islands, conservation, Lacertidae, morphology, phylogeny, Podarcis lilfordi, scalation evolutionary history of populations (phylogeny) or their current climatic conditions?

| Specimens
We studied all known populations ( Figure 1) of P. lilfordi, although we did not obtain genetic, morphometric or coloration information from all of them. Samples from Cabrera Island were obtained from four sites: Cabrera Port (the area around the Cabrera Bay), Cabrera far (the Ansiola Peninsula lighthouse, at the south-western corner of the Island), Morro den Tià (another Western Peninsula of the Island) and Miranda (central part of Cabrera Island, Figure 1). Note that the population studied from the Bay of Palma (Mallorca), at Porrassa Islet, is known to have been introduced, as well as the small population of Colonia de Sant Jordi in Southern Mallorca.

| Climatic characteristics
The aridity of the Balearic Islands was estimated with the I q index (Sahin, 2012). In this index, the ratio of annual precipitation to annual mean specific humidity (Sh) is employed. The mean specific humidity can be easily computed with mean temperature, relative humidity and local pressure which are the most commonly measured meteorological data (Sahin, 2012). Due to the limited geographical distribution of P. lilfordi, we were only able to compare climatic data from three weather stations: the Port of Palma de Mallorca (Palma F I G U R E 1 The distribution of major clades of Podarcis lilfordi in Balearic Islands (Ibiza Island is not represented in the map). We showed all known populations of the Balearic lizard in Menorca, Mallorca and Cabrera Archipelago (see Table 1 for a list of described subspecies and their distribution) TA B L E 1 Known populations of the Balearic lizard and described subspecies (see Figure 1 for locations of populations under study)
The two weather stations on Mallorca Island are very close. We employed data from these weather stations freely provided by AEMET Open Data (Agencia Estatal de Meteorología, Spain). Aridity indices and rainfall were calculated with annual data from 1981 to 2018.

| DNA extraction, amplification and sequencing
Total genomic DNA was extracted from tail tips following standard protocols (González et al., 1996). A total of 104 (Table 1) individuals of P. lilfordi have been studied for four non-overlapping mtDNA fragments: (a) partial 12S rRNA, (b) two partial fragments of cytochrome b (CYTB), (c) partial control region (CR) and (d) 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 are the same as those used in our previous studies of Podarcis . Additionally, 92 (Table 1) individuals were sequenced   for COI fragment using primers LCO-1490: GGT CAA CAA ATC ATA   AAG ATA TTG G and HCO-2198: TAA ACT TCA GGG TGA CCA AAA AAT CA (Folmer, 1994

| Divergence and phylogenetic analyses
Distance-based method (p-distance) based on mitochondrial data (mtDNA and COI fragment) between subspecies was calculated in mega 7 (Kumar, Stecher, & Tamura, 2016) to stablish the level of sequence divergence. We used a dataset including 46 P. lilfordi populations for genetic analyses ( 0.01. Results were analysed in tracer v1.6 (Rambaut, Suchard, Xie, & Drummond, 2014) to assess convergence and effective sample sizes (ESS) for all parameters. A burnin of 25% was applied, and the phylogenetic trees were visualized and edited using figtree v1.4.2 (Rambaut, 2012).
All body measurements were done by the same researcher (VPM).
Morphometric and scalation characters were compared among the three groups of populations (Mallorca, Cabrera and Menorca) using ANOVA type, Lawley-Hotelling type, Bartlett-Nanda-Pillai type and Wilks' lambda type test statistics, as implemented in the R-package "npmv" (Burchett, Ellis, Harrar, & Bathke, 2017). These statistics were then used as the basis for permutation or randomization tests. Nonparametric relative effects of each morphometric and scalation character were also tested. These effects give an indication of stochastic superiority; that is, they measure the probability that a value obtained from one group is larger than a value randomly selected from the whole dataset. In agreement with default settings of npmv, we employed the results of Wilks' lambda in pairwise comparisons, and because N > 30, we used F approximation (Burchett et al., 2017).
In addition, we employed a non-metric multidimensional scaling (NMDS) to establish morphological divergence, separately for adult females and males, among 43 populations of P. lilfordi.
The method aims to depict the inherent pattern of a dissimilarity matrix in a geometric picture with a minimum number of dimensions (Clover, 1979). We used the metaMDS function from the "vegan" R-package (Oksanen et al., 2018). This function runs NMDS several times from random starting configurations, compares results and stops after detecting two similar minimum stress solutions (Oksanen et al., 2018). The goodness-of-fit of the ordination was assessed by the coefficients of determination (R 2 ) for the linear and nonlinear regressions of the NMDS distances on the original ones (Borcard, Gillet, & Legendre, 2011). Finally, we recorded the stress values of NMDS (Zuur, Ieno, & Smith, 2007). All calculations were done within R environment (R Core Team, 2019).
Dorsal colouration of lizards was analysed with the "colordistance" R-package (Weller, 2019;Weller & Westneat, 2019). This package is an objective comparative tool to colour profiling and comparison of digital images. The standard RGB image analysis cannot provide the wavelength resolution of reflectance spectrophotometry. However, the use of digital images can reflect the visual sensitivities of several species (Losey et al., 2003;Weller & Westneat, 2019).
Our colour analysis was restricted to 30 populations of adult males and 27 populations of adult females of the Balearic lizard from which we had good JPEG dorsal colour images. We employed 2-15 individual images per population. Non-background pixels were binned to read images into the R environment as 3D arrays. We used colour histograms as the binning method in a 5 × 2 × 3 hue-saturation-value (HSV) colour space. HSV is the most suitable method when different digital cameras or variable light conditions were employed to obtain images (Weller & Westneat, 2019). For calculating the distance between one binned image and another, we employed the earth mover's distance or Wasserstein metric (EMD) method. In this way, colour histograms from each image or from a group of images (i.e., from all individuals of a given population) are compared with a symmetrical distance matrix that is plotted as its corresponding heatmap. In the heatmap, we used default colours ranging from yellow (least similar) to blue (most similar).

| Comparative analysis
We did comparative analyses of morphological traits using the mtDNA tree ( Figure 2). Due to the low levels of divergence and consequently the high number of polytomies, we were unable to obtain an ultrametric and fully dichotomous tree for most comparative analyses. Thus, we restricted our analysis to the comparison of morphological continuous variables with the "caper" package (Orme et al., 2018), using the method of Pagel (1992) to calculate contrasts at polytomies. We compared in this way all continuous variables of morphometry and scalation using SVL as covariate (Orme et al., 2018). In addition, to study the presence/absence of melanic coloration, we employed the D statistics (Fritz & Purvis, 2010) for binary variables of the "caper" package. In order to standardize the effects of phylogeny size and prevalence, "phylo.d" uses two simulated null models: a model of phylogenetic randomness, where trait values are randomly shuffled relative to the tips of the phylogeny, and the Brownian threshold model, where a continuous trait is evolved along the phylogeny under a Brownian process and then converted to a binary trait using a threshold that reproduces the relative prevalence of the observed trait. D typically varies between 0 and 1. A D value close to zero indicates that the binary trait evolves on a tree following the Brownian model (i.e., strong phylogenetic signal).

| Climatic characteristics in the Balearic Islands
We found significantly lower aridity (I q index) in Menorca

| Genetic analysis
We have sequenced four fragments plus COI of mtDNA genome. The In Porrassa Island inhabits an introduced population, it is included in phylogenetic trees and exhibits a high proximity to Cabrera
Thus, all analyses were done separately for adult males and females. Morphology shows a transition from Menorca to Mallorca and Cabrera characterized by increased body size (see Figures 5 and 6 for SVL, and Tables S1-S6), and a decrease in adjusted head dimensions (Tables S1-S6) There is an overall significant difference in body dimensions and scalation between Mallorca and Menorca lizards (p < .05) and Cabrera and Menorca lizards (p < .05), but not between Mallorca and Cabrera archipelagos (p > .05). In adult females, we also reject the hypothesis of equality among Cabrera, Mallorca and Menorca (Table 6, p < .001). We find significant differences in body dimensions and scalation between Mallorca and Menorca lizards (p < .05) and Cabrera and Menorca lizards (p < .05), but again, we do not find significant differences between Mallorca and Cabrera (p > .05).
The analysis of relative effects shows higher values in terms of probabilities for all morphometric and scalation characters of male lizards from Cabrera, with lowest values for all characters from male lizards from Menorca (Table 7). Similar tendencies are observed in females (Table 8). In males, higher probabilities correspond to the length of intact tails, body mass, head height and pileus length of lizards from Cabrera. In females, all probabilities of effects of each variable are lower than in males, indicating a closer distance of their morphologies among the three groups of islands. Only tail lengths of Cabrera females showed a value above 0.80 (Table 8).

| Patterns of lizard coloration
The

| The evolutionary framework of the Balearic lizard
The Balearic lizard was a member of the first post-Messinian faunal assemblage reported from Menorca Island (Bover et al., 2008(Bover et al., , 2014. This fauna spread during the Early and Middle Pliocene, including a lizard (Podarcis sp., Bailón, 2004) which could be the ancestor of TA B L E 2 Evolutionary divergence over sequences using p-distances among defined subspecies of P. lilfordi based on the four mtDNA fragments

| Influence of climate in morphological features
Bergmann's rule (Bergmann, 1847) proposes that in birds and mammals, body size is greater in high latitudes for adaptive reasons: the lower surface-volume ratio in endotherms of greater body size confers thermoregulatory advantages, as heat losses are minimized in colder environments. The issue is more complex in ectotherms.
In Squamata, the existence of an inverse Bergmann's rule has been proposed (Oufiero, Gartner, Adolph, & Garland, 2011 and references therein), with larger body sizes in lower latitude individuals, that is, in warmer environments (Ashton & Feldman, 2003). But there are also data that demonstrate an opposite trend in some groups of reptiles (Sears & Angilletta, 2004 Numerous studies have addressed the adaptive value of the size and number of scales in relation to climatic variables (Arnold, 1973;Bogert, 1949;Calsbeek, Knouft, & Smith, 2006;Hellmich, 1951;Soulé, 1966;Soulé & Kerfoot, 1972). In some cases, lizards from larger islands showed fewer and larger scales, a trend interpreted as an adaptive response to the higher average temperatures on larger islands. Thus, in localities where climate is warmer, we can expect larger scales. In addition, the experimental work of Regal (1975) supported the hypothesis that enlarged scales of lizards inhabiting warmer habitats may function as heat shields. In this way, scales can be considered a protection against abrasion and water loss (Alibardi, 2003;Walker & Liem, 1994). Thus, populations from arid habitats tend to have fewer and larger scales to reduce skin exposure and the amount of evaporative water loss (see, e.g., Hu et al., 2019 and references therein). A fewer number of scales reduce skin exposure, while a larger number of smaller scales increase the exposed skin surface and, consequently, the rate of dehydration (Calsbeek et al., 2006;Neilson, 2002).

| Taxonomic arrangement of P. lilfordi according to genetic and morphological traits
As expected, if we compare the arrangement of described subspecies of the Balearic lizard with genetic, morphological and colouration results, it is clearly a lack of congruence of present-day accepted subspecies (Table 1). In this study, we only employed mtDNA and one nuclear gene and we were unable to find a congruent arrange-  Mones and Mel islets (Table 1). Both are morphologically well separated (Figures 7 and 8), but with a strong genetic similarity to populations apparently unrelated, as in the case of Addaia Gran and Mel (Figures 2 and 3). Some of the described subspecies of P. lilfordi are morphologically and genetically differentiated, as those from Dragonera, Malgrat islets and Toro in Mallorca coast (Figures 2 and 3). However, it is not the case of several populations of the Balearic lizard described as different subspecies, as P. lilfordi lilfordi from Aire Island and P. lilfordi balearica from Rei Island  and this study).
Insular populations of terrestrial vertebrates encounter a physical impediment to gene flow between populations, and it is therefore expected that such populations may diverge in isolation (Mayr, 1963).
Being frequently smaller in population size, the fixation of neutral genes is likely to take place in a faster way than in continents. However,  Rei Colom conservation utility, as proxies for the sub-structure found within species (Phillimore & Owens, 2006). However, from a cladistic viewpoint, the usefulness of subspecies is a controversial issue (Vinarski, 2015;Wilson & Brown, 1953;Winker, 2010

| The application of evolutionary significant units to Balearic lizards
Solving this dilemma can only be achieved abandoning the very concept of subspecies in the case of the Balearic lizard. As an alternative, we propose the use of evolutionary significant units (ESUs, Moritz, 1994Moritz, , 2002. In addition to its original genetic definition (Moritz, 1994), an ESU was also defined as a group of organisms substantially reproductively isolated that represent an important component in the evolutionary legacy of the species (Waples, 1995). This is the case of every microinsular population of P. lilfordi, physically isolated even from closest islets. An ESU can often correspond to species or subspecies boundaries in classical taxonomy, but in some cases, it can be extended to isolated populations (Karl & Bowen, 1999).
The ESU has been proposed as a unit of conservation (Moritz, 1994;Vogler & DeSalle, 1992;Waples, 1995). We have to evaluate each population along two axes of diversity which might be described as molecular genetics and adaptive morphology. Adaptive diversity is the diversity of morphological and ecological traits we observe today in different populations of P. lilfordi. This diversity represents the raw material for future evolution. Consequently, we have to recognize every isolated population of the Balearic lizard as an ESU (Pérez-Mellado, 2008).
In this study, we have only addressed the description and possible origin of some of the morphological features of the Balearic lizard. However, it remains to elucidate the origin and causality of a myriad of ethological and ecological features that, in one way or another, characterize each population (Pérez-Mellado, 2009 and references therein). This is the task that implies the complete preservation of the microevolutionary process in the islets where the Balearic lizard survives for thousands of years. The potential taxonomic identity of each population is not its most relevant aspect, and it is paradoxical that it is the main aspect that has attracted the attention of experts and amateurs.   In short, in each population of P. lilfordi we observe an evolutionarily independent history linked to different ecological conditions, extremely variable population sizes and availability of trophic resources, or the presence of particular competitors and predators.
This situation has led to adaptive responses, in many cases unique to each population, that are resolved in morphological, ethological and ecological features so far only observed in that populations. It is impossible to foresee what the future knowledge of numerous populations holds for which there is no reliable data on ecology and behaviour of their lizards.
Consequently, the recognition of all the populations as ESU implies that each and every one of the populations of the Balearic lizard is fully relevant to its conservation. We do not want to protect only the product of evolution, that is, the extraordinary variability of the polytypic Balearic lizard today observed. We want to preserve the future process of evolution of the Balearic lizard (Moritz, 2002).
The Balearic lizard is listed by IUCN as endangered (criteria: B1ab ( Paradoxically, the regional conservation managers of the Balearic Islands consider that the species is only vulnerable (Viada, 2006).
Probably, this unrealistic assessment derives from the fact that some specific populations, such as those of the Dragonera or Cabrera Islands, have a large number of lizards. The complete preservation of all the extant populations will be the only guarantee of conservation of the unique evolutionary process of this species (Pérez-Mellado, 2008).

ACK N OWLED G EM ENTS
This study was possible thanks to a grant for the project entitled: "Evolución morfológica y genética de la lagartija balear Podarcis lilfordi (Günther, 1874)

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available on request from the corresponding author.