Genetic differentiation predicts body size divergence between island and mainland populations of common wall lizards (Podarcis muralis)

Small-bodied vertebrates sometimes evolve gigantism on islands, but there is a lack of consistent association with ecological factors or island characteristics. One possible reason is that, even if the ecological conditions are right, body size might fail to diverge on islands that were isolated recently or if there is gene flow between islands and the mainland. We studied body size, ventral colour polymorphism and genetic structure across nine islands and adjacent mainland populations of common wall lizards ( Podarcis muralis ) off the western coast of France. Population genetic data suggested that island populations might have maintained gene flow after their geographical isolation from the mainland. Island lizards were larger and heavier than mainland lizards on average, but the extent of gigantism varied substantially between islands. Island size and distance from the mainland were poor predictors of body size, but lizards from populations that were highly genetically differentiated from the mainland were larger than lizards from less differentiated populations. Colour morphs that were rare on the mainland tended to be more common on islands. We propose that genetic isolation or bottlenecks promote body size evolution in island lizards, which makes it challenging to identify ecological causes of island gigantism without complementary genetic information.


INTRODUCTION
Animals that live on islands occasionally diverge from mainland populations in body size, colour or other characters (Grant, 1965;Case, 1978;Lomolino, 1985Lomolino, , 2005. Conspicuous instances of body size divergence between island and mainland populations of mammals, birds and reptiles have been explained by a number of ecological factors, including dietary shifts, increased intraspecific competition and reduced predation (e.g. Foster, 1964;Clegg & Owens, 2002;Boback, 2003;Boyer & Jetz, 2010;Aubret, 2012). Nonetheless, the evidence that these factors are sufficient to cause an increase or decrease in body size remains equivocal. Indeed, populations on islands that are ecologically similar often differ in body size, and the best predictor of body size divergence between island and mainland populations for one species is often a poor predictor for another species (Meiri, 2007;Meiri et al., 2011;Itescu et al., 2018;Lokatis & Jeschke, 2018).
One possible reason for this apparent empirical inconsistency is that the evolutionary history of island populations plays an important role in phenotypic divergence but is rarely considered. For example, founder effects may both promote and retard opportunities for phenotypic divergence, and gene flow between mainland and island populations and between island populations can prevent natural selection from having a substantial effect. Although studies of island gigantism commonly consider island size and isolation to be proxies of ecological causes of selection (Lomolino, 2005), such features may also correlate with genetic differentiation from the mainland. Studies that combine phenotypic and genetic data are therefore an important starting point for identifying the conditions that are conducive for island gigantism.
Wall lizards of the genus Podarcis are widespread in southern Europe, and several species are common on, or endemic to, Mediterranean islands. Large body size on small islands has been documented in Podarcis (Pafilis et al., 2009;Runemark et al., 2015), but it does not appear to be a highly repeatable phenomenon (Itescu et al., 2018). One European species, Podarcis muralis, has significantly expanded its range from the Mediterranean basin, and its native range now includes north-western France (Salvi et al., 2013). This species also inhabits several islands of variable sizes in the Atlantic (Fig. 1). It is likely that these populations became established as the islands separated from the mainland by rising sea levels ~5000-7000 years ago (Ters, 1986;Pluet & Pirazzoli, 1991). However, the proximity to the mainland means that more recent isolation, colonization or occasional gene flow between islands and between island and mainland populations cannot be ruled out, and these processes might contribute to patterns of phenotypic divergence observed between populations.
The aim of this study was to establish the relationship between genetic isolation and phenotypic, in particular body size, divergence between island and mainland lizard populations. To this end, we collected genetic and phenotypic data for nine island and adjacent mainland populations in western France to test: (1) whether the genetic structure on islands reflects their colonization history and island characteristics (i.e. area and distance from mainland); (2) whether island populations exhibit phenotypic divergence from their mainland counterparts; and (3) whether there is an association between body size divergence, island area, geographical isolation and genetic differentiation.  Table 1 for details).

MATERIAL AND METHODS
The common wall lizard (P. muralis) is a small diurnal lizard. It is widespread in southern and western Europe throughout a range of natural and highly disturbed habitats, including villages or cities. The present study was conducted in western France, where the species is common on the mainland and on several islands in the Atlantic (Fig. 1). Lizards were caught opportunistically by noosing, and snout-to-vent length (SVL) and total length were measured to the closest millimetre with a ruler, and head size to the closest 0.1 mm with a pair of callipers. Body mass was recorded with a digital field balance to the nearest 0.01 g. The tail tip (~5 mm) was removed with surgical scissors and stored in 95% ethanol for genetic analysis. We also scored the sex of each lizard and recorded whether the tail was regenerated or recently lost. The ventral surface of each lizard was photographed and later assigned to one of five colour morphs: three 'pure' morphs (white, yellow or orange) and two mosaic morphs (white-orange or orange-yellow; white-yellow does not appear to exist in this species). These colour morphs have been well studied in P. muralis (e.g. i de Lanuza et al., 2013;Sacchi et al., 2013). The expression of yellow and orange coloration is known to be associated strongly with two genetic loci, whereas the genetic basis (if any) of the mosaic morphs is unknown (Andrade et al., 2019). In total, we collected phenotypic data for 1115 lizards from nine islands and 349 lizards from 11 of the 20 mainland locations (Table 1). For the island Belle Île en Mer, lizards were sampled in four locations, and for the island Île d'Yeu in two locations, but these were pooled into one sample per island in the subsequent analysis. Lizards were released at the site of capture after sampling.

Molecular analyses
To estimate population genetic structure and differentiation, we relied on 16 microsatellite markers and a 655 bp region of the cytochrome b (Cytb) gene in the mitochondrial genome (mtDNA). Genetic data (microsatellites and Cytb sequences) have been published for 14 of the 20 mainland populations (Michaelides et al., 2015b). For the remaining six mainland and 12 island locations (nine islands; no tissue from one location on Île d'Yeu), we generated genetic data following the same protocol (see Michaelides et al., 2015b). In brief, genomic DNA was extracted (N = 853;   and nine by Heathcote et al., (2014). The reactions were carried out in a total volume of 11 μL reaction mix containing 1 μL of genomic DNA, 5 μL of Qiagen MasterMix, 0.2 μL of each primer (forward and reverse in equal concentrations) and 3.8 μL (for multiplexes 1, 2, 3 and 5) or 3.6 μL (for multiplex 4) of PCR-grade H 2 O. The PCR conditions were as follows: 15 min of initialization step at 95 °C, 26 cycles of 30 s at 94 °C, 90 s at 57 °C (for multiplexes 1, 2 and 3) or 55 °C (for multiplexes 4 and 5) and 1 min at 72 °C, and a final extension step of 20 min at 60 °C. The 5′-end of each forward primer was labelled with a fluorescent dye (6-FAM, HEX or NED). The PCR products were run with an internal ladder (red ROX-500), on an ABI3730XL DNA analyser (Thermo Fisher Scientific). We scored alleles in Geneious v.6.1.7, and any ambiguous peaks were repeated to confirm the genotype. For a subset of samples (between five and 14 individuals per population), partial mtDNA Cytb gene was sequenced by PCR using the primer pair LGlulk (5′-AACCGCCTGTTGTCTTCAACTA-3′) and Hpod (5′-GGTGGAATGGGATTTTGTCTG-3′) (Podnar et al., 2007;Schulte et al., 2012;Michaelides et al., 2013). Amplifications were carried out in a total volume of 15 μL consisting of 7.5 μL of MyTaq HS Mix (Bioline), 0.45 μL (8 pmol) of each primer (Eurofins), 4.6 μL PCR-grade H 2 O and 2 μL template DNA. The PCR conditions were as follows: an initial denaturation step at 94 °C for 1 min, followed by 35 cycles at 94 °C for 1 min, 53 °C for 45 s and 72 °C for 1 min, and a final extension step at 72 °C for 10 min. The PCR products were purified using the ExoSAP-IT Cleanup Reagent (Applied Biosystems). Sequencing reactions were carried out with BigDye Terminator v.3.1 Ready Reaction kit (Applied Biosystems) in both directions. Products were precipitated with isopropanol and analysed on an ABI3500 automated capillary sequencer (Applied Biosystems). Sequence reads were corrected by eye and mapped to a reference sequence in Geneious v.6.1.7 (Drummond et al., 2011). Final sequences were aligned and trimmed into a uniform length of 655 bp in BioEdit Sequence alignment editor (Hall, 1999). Newly obtained sequences were submitted to GenBank under the accession numbers MK618780-MK618788.

Genetic data analysis
We used the mitochondrial and microsatellite data to perform a number of population genetic analyses: (1) to characterize the genetic diversity within the sampled populations; (2) to infer how these populations cluster together; (3) to calculate their degree of divergence; (4) to test whether populations experienced a bottleneck event; (5) to estimate the timing of isolation of island populations from their mainland counterparts, and subsequent gene flow and migration rates after colonization; and (6) to estimate effective population size.

Genetic diversity and structure
First, the microsatellite data were inspected for the presence of null alleles, effects of stuttering and large allele dropout using Micro-checker v.2.2.3 (Van Oosterhout et al., 2004). One locus (PmurC275_278) displayed an allele dropout pattern, and we discarded this locus from further analyses. We found no evidence of stuttering, but some loci showed signs of null alleles. However, as these were not present across all populations, we retained the remaining 15 loci for further analyses.
Basic genetic diversity indices, observed and expected heterozygosity (H O and H E ) were calculated with Genalex v.6.5 (Peakall & Smouse, 2006 and allelic richness (A R ) with FSTAT v.2.9.3 (Goudet, 2002). The inbreeding coefficient (F IS ), which estimates the proportion of variance in the population that is contained in an individual, was calculated with Arlequin (Excoffier et al., 2007). To examine the genetic structure of the sampled populations, we used a Bayesian inference clustering method implemented in STRUCTURE v.2.3.4 (Pritchard et al., 2000), with an admixture model (Falush et al., 2003) and correlated allele frequencies. This method finds the number of genetic clusters that best fit the data. Ten independent runs were made, with a burn-in of 100 000 iterations and a run length of 10 6 iterations for a number of genetic clusters (K) = 1-15. The best K was determined according to the method described by Evanno et al. (2005). To characterize the genetic structure further, we used a discriminant analysis of principal components (DAPC) implemented in the R package 'adegenet' (Jombart et al., 2010;Jombart & Ahmed, 2011). This approach, in contrast to STRUCTURE, uses coefficients of the alleles in linear combinations and maximizes between-group variance, while minimizing within-group variance, without the assumption of Hardy-Weinberg equilibrium (Jombart et al., 2010).
Pairwise genetic distances (F ST , G″ ST and Jost's D ST ) between island and mainland populations were calculated in Genalex v.6.5. For these analyses, the mainland populations were grouped according to the highest membership score based on the STRUCTURE analysis (see Results), and each island was then compared separately with the relevant mainland population group. We also calculated pairwise genetic differentiation (F ST ) between all populations (island and mainland) in Arlequin (Excoffier et al., 2007).
We used two approaches to detect whether island and mainland populations had undergone genetic bottlenecks after colonization and range expansion, respectively. First, we calculated the degree of heterozygosity excess, which occurs because of the loss of rare alleles shortly after bottlenecks, using BOTTLENECK (Piry et al., 1999). We used a twophase mutation model, with 95% stepwise and 5% nonstepwise mutations. The significance of heterozygosity excess was then calculated using Wilcoxon tests. Second, we calculated Garza and Williamson's index (M-ratio), by dividing the number of alleles in a population (k) by the range in allele size (r) (Garza & Williamson, 2001) in Arlequin (Excoffier et al., 2007;Excoffier & Lischer, 2010). The M-ratio should be small in populations that have recently been through a bottleneck and larger in stationary populations.

tiMinG of isolation and subsequent Gene flow
To estimate the timing of isolation and subsequent gene flow between the islands and the mainland, each of the islands was compared with the relevant mainland population. As above, the mainland populations were grouped based on the highest membership score of the STRUCTURE analysis (see Results). To test whether the relative timing of isolation of the islands differed from each other, we used an approximate Bayesian computation (ABC) analysis in DIYABC software (Cornuet et al., 2014) on the microsatellite data. For this purpose, we assumed, separately for each island, a scenario where the ancient mainland population (with effective population size Ne A ) splits into island and mainland populations (with effective population sizes of Ne 1 and Ne 2 , respectively) at time t, without any subsequent gene flow. The coalescent-based algorithm simulates datasets for a predefined scenario and compares the resulting summary statistics with those of the observed data. Model priors were kept at default values, with the exception of individual locus mutation rates, where prior limits were set to 1 × 10 -7 to 1 × 10 -2 . The summary statistics used in ABC were the mean number of alleles, mean genetic diversity and mean size variance (one-sample summary statistics), together with pairwise F ST values, classification index (Rannala & Mountain, 1997) and a genetic distance value, (dμ) 2 , developed by Goldstein et al. (1995) (two-sample summary statistics). Pre-evaluation of parameter priors was carried out with preliminary runs before the final 1 × 10 6 datasets were simulated and the posterior distributions examined.
Potential migration from the mainland after the initial colonization of the islands was assessed from the microsatellite data using MIGRATE v.3.6.11 (Beerli & Felsenstein, 2001). MIGRATE estimates effective population sizes and past migration rates between the sampled populations, assuming an asymmetric migration matrix model and asymmetric population sizes. We used Bayesian inference with a Brownian motion mutation model to estimate the parameters jointly. The population size parameter, θ (θ = 4Neµ, where Ne is the effective population size and µ the mutation rate), and the migration rate parameter, θM (estimate of number of immigrants per generation), from the mainland to the island were estimated for each of the islands. We used F ST estimates as starting parameters. Relative mutation rates were estimated from the data, and uniform parameter priors (for θ and M) were set to range from zero to 1000 (in steps of 100). We ran 1 × 10 6 genealogies at a sampling increment of 50 iterations for each locus and discarded the first 10 000 as burn-in. An adaptive heating scheme using four simultaneous Markov chains was applied to increase the search efficiency (Geyer & Thompson, 1995).

PhenotyPic data analyses
Statistical analyses were performed in R v.3.5.1. We first tested for differences in body size between island and mainland populations. To do so, we analysed the variation in body mass (mass) and snout-to-vent length (SVL) with linear mixed effects models (LMMs), with location (island/mainland), sex (female/male) and their interaction as fixed effects and population as a random effect. The incidence of tail loss was analysed using a generalized linear mixed model (GLMM) with the same random and fixed effects, a binomial distribution and a logit link function. All mixed effects models were analysed using the lme4 package (Bates et al., 2015). The significance of fixed effects was evaluated using F-tests, with degrees of freedom approximated using the Kenward-Roger method (for LMMs), or likelihood ratio tests (for GLMMs). We used DHARMa for model diagnostics (Hartig, 2018).
Second, we tested whether there were any associations between the mean body size (mass and SVL) and island area, distance to the mainland and genetic differentiation (F ST ) from the mainland. Given that island size and distance to the mainland are predicted to affect population divergence in body size both through their association with putative causes of selection and through independent effects on genetic differentiation (e.g. drift, migration), it would be desirable to use structural equation modelling to calculate direct and indirect causal effects. However, this is not feasible with only nine islands, and we therefore used a multiple regression approach with both backward and forward selection (stepAIC in the MASS package; Venables & Ripley, 2002). Island area and F ST were log 10 -transformed, and all input variables were standardized (z-transformed) before analysis (Schielzeth, 2010). Differences in the frequency of ventral colour polymorphism between populations were tested using χ 2 tests (MASS package). For this analysis, we combined the orange-white and orange morphs and eliminated the very rare orange-yellow morph to avoid the many zeros that would otherwise occur. Morph frequencies were visualized in a ternary diagram using the R software package 'ggtern' (Hamilton, 2018), an extension to ggplot2 (Wickham, 2009).

Genetic analyses
We found nine mtDNA haplotypes across the island populations. All the sequenced individuals in the Glénan islands had a Cytb haplotype identical to that found almost uniformly in all the mainland populations ( Supporting Information, Fig. S1). In contrast, Île d'Yeu, and especially Belle Île en Mer, exhibited more structure in the mtDNA ( Supporting Information, Fig.  S1). Excluding the one very common haplotype present all over the mainland, there were no haplotypes shared exclusively between individual islands and mainland populations (Supporting Information, Table S2).
Genetic diversity in terms of heterozygosity varied between 0.12 (population GGU) and 0.73 (PU and NE) and was generally lower within the island populations (average H E = 0.48 for all island populations) than in the mainland populations (average H E = 0.61 in the northwestern populations, H E = 0.70 in the southern populations, Mann-Whitney U-test mainland vs. island; P = 0.006; Table 1). Allelic richness varied between 1.45 (population GGU) and 5.33 (NE) and was lower within the island populations (average A R = 3.41 for all island populations) than in the mainland populations (average A R = 4.13 in the north-western populations, A R = 4.97 in the southern populations, Mann-Whitney U-test mainland vs. island; P = 0.021; Table 1). The estimates of genetic diversity were higher for large islands than for small islands (Spearman rank order correlation: H E , r = 0.80, P = 0.014; A R = 0.88, P < 0.01).
The Bayesian clustering analysis STRUCTURE indicated that two genetic clusters best captured the genetic structure in our data (Supporting Information, Fig. S2A), separating the southern France populations from the north-western populations. However, the likelihood was almost as high for three clusters, further dividing the north-western populations into two. Given that we were mainly interested in the genetic structure in this region (Fig. 1), we ran the analysis separately for the northern populations (i.e. excluding BA, RQ, LI, SG, OU, FC, NE and LU). For this dataset, two clusters were inferred (Supporting Information, Fig. S2B). Given that there was considerable structure among the north-western populations, we considered the number of genetic clusters to be three in the whole dataset, consisting of two north-western clusters (I and II) and a southern cluster (III; see Supporting Information, Table  S3; Fig. 2  north-western populations. All southern populations formed one very uniform cluster (Q = 0.97-0.99; Fig.  2). The DAPC clustering analysis, which maximizes the between-population variation, showed similar genetic structuring, distinguishing the southern and the north-western populations, with population DN as a mixture population between the three main clusters, but also some further subdivision within north-western mainland and island populations (Fig. 3).
To calculate average genetic distances of the islands to the mainland, we divided the mainland populations into relevant genetic clusters based on the STRUCTURE results (Table 2). For this purpose, we considered populations as part of the genetic cluster to which they showed the highest membership ( Fig. 2; see Supporting Information, Table S3 for population average Q-values). Genetic differentiation from the mainland varied substantially between islands, and increased with decreasing island size (Pearson's product-moment correlation on F ST and island area, both log 10 -transformed; r = −0.73, P = 0.02). There was low genetic differentiation between island populations (F ST < 0.07, with GGU being an exception; Supporting Information, Fig. S3).
The population divergence times estimated by DIYABC were very similar between all the islands, with substantial overlap of the 95% confidence intervals (Table 3). Both the lowest and the highest estimate were found within the Glénan islands [lowest, Guéotec (GGU) 362-6800 years ago, mode 747 years; and highest, Glénan Le Loc'h (GLL) 1120-8730 years ago, mode 3220 years; Table 3]. The estimates of gene flow (number of immigrants per generation from the mainland) were also similar among the islands (Table  4). The estimate was lowest for Glénan Guéotec (GGU; mode = 1.67; confidence interval CI 97.5% = 0-18) and highest for Belle Île en Mer (BEL; mode = 9.67; CI 97.5% = 0-26), but all estimates included zero migrants within their 97.5% confidence interval. The mode values for effective population size estimates (Ne) were lower in all the island populations when compared with their respective mainland populations and reached statistical significance for YSS, GBV and GGU (Table 3).
We used two tests to detect genetic bottlenecks. There was no significant heterozygosity excess (P HE ) across all populations (Supporting Information, Table  S4), but the M-ratio was significantly lower in island than in mainland populations (Mann-Whitney U-test mainland vs. island, P = 0.0002).
Body mass increased significantly with genetic differentiation from the mainland, whereas there was no significant effect of island area or distance from the mainland (Table 5). Snout-to-vent length showed a very similar effect of genetic differentiation after model selection (Table 5).
The frequency of ventral colour polymorphism was highly variable between populations (χ 2 = 96.9, P < 0.001, d.f. = 2). Morphs that expressed yellow and/ or orange were more common on most of the islands than on the mainland, with one instance of apparent fixation of the orange-white mosaic morph on the island of Guéotec (Fig. 5).

DISCUSSION
There has been considerable controversy over the evolution of island gigantism in small-bodied vertebrates (Lomolino, 2005;Meiri et al., 2005Meiri et al., , 2006Meiri et al., , 2008Meiri, 2007;Lokatis & Jeschke, 2018). Although island gigantism was once considered sufficiently common to be described as a rule (e.g. van Valen, 1973;Whittaker, 1998), more recent analyses across a broader range of taxa suggest that gigantism is occasional rather than the norm, and that no single ecological factor is consistently associated with increases in body size (Lomolino, 2005;Itescu et al., 2018;Lokatis & Jeschke, 2018). Our results demonstrate that common wall lizards in island populations in the Atlantic tend to be larger and heavier than lizards in mainland populations in western France. However, there was considerable variation between populations in the extent of gigantism. Some populations showed no increase or a weak increase in body size, whereas lizards on other islands were substantially larger than the average mainland population, with one island population (Guéotec) reaching an average 20% increase in SVL and twice the body mass compared with the mainland average. We show that this variation in gigantism is positively associated with the extent of neutral genetic differentiation between islands and the mainland. Although the ecological conditions that cause this variation in body size among island populations were not identified in the present study, our results suggest that genetic isolation or bottlenecks might promote opportunities for body size divergence between islands and the mainland.

Genetic differentiation and body size variation
in wall lizards in western france The islands off the coast of western France diverged from the mainland about 5000-7000 years ago as a result of rising sea levels. The estimates of divergence from the mainland based on microsatellite data are imprecise but tended towards a more recent isolation. Although the estimate of immigration, bottlenecks and genetic differentiation appears to rule out recent colonization, there could have been significant gene flow from the mainland even after the geographical separation from the mainland, drawing the divergence estimate closer to the present. Indeed, the unique cytb haplotypes in the Belle Île en Mer, despite the close genetic relationship to the mainland, might reflect early separation from an unsampled or extinct population, followed by more recent gene flow. Although islands that are further apart are more divergent genetically, the pattern of population differentiation also suggests that there might have been contemporary gene flow among the Glénan islands. Given that mitochondrial identities did not reveal any exclusive island-mainland connections in our relatively sparsely sampled dataset, they provide limited insight into the history of the populations. In contrast, wall lizard populations on the Channel Islands show a more consistent signal of isolation from the mainland, despite that the islands were isolated geographically at about the same time as the islands in the present study (Michaelides et al., 2015a). To some extent, this difference might be a consequence of the fact that the Channel Islands are at the northern extreme of the species distribution, which might historically have had lower lizard densities and a patchy population distribution (as is the case today; Lescure & Massary, 2012), hence reduced genetic diversity compared with the focal region in the present study (see Michaelides et al., 2015aMichaelides et al., , 2016. Genetic analysis using more informative markers would be required to identify the demographic history of the island populations and the extent to which island and mainland populations have exchanged genes over the past 5000 years. Genetic diversity and genetic differentiation from the mainland differed substantially between island populations, and both were negatively correlated with island size. The finding that populations on small islands lose genetic diversity and differentiate is expected, but it might also have consequences for the extent of phenotypic divergence. Indeed, genetic differentiation was a good predictor of body size divergence, even when controlling for island size and distance from the mainland. This relationship should be interpreted with caution given the small sample size, the correlation between island size and genetic differentiation, and the fact that several processes could contribute to high genetic differentiation between islands and the mainland. Nevertheless, one possible Table 3. Posterior distributions for the parameters from the approximate Bayesian analysis (DIYABC; Cornuet et al., 2014) interpretation is that genetic isolation promotes the capacity to respond to selection for increased body size. Indeed, the most exaggerated body size is found on an island [Guéotec (GGU)] in the centre of the Glénan archipelago, with very low estimates of immigration rate and high genetic differentiation from nearby islands. This population was also unusual in that one of the, usually rare, ventral colour morphs has gone to fixation. Given that Guéotec is not substantially smaller than some of the other islands, it is possible that this population has not only experienced limited gene flow, but also that it has been subjected to a more severe bottleneck than other populations. If so, the lack of heterozygosity excess suggests that this bottleneck was not very recent. A temporarily low population size can contribute to phenotypic divergence, partly owing to drift but also because it can increase the additive genetic variance and thus facilitate a response to selection (e.g. Whitlock et al., 1993). Both processes might be important for explaining why lizards on this island are so large bodied. However, it appears unlikely that the very large body size can be explained by founder effects and drift alone, nor can this explain the positive relationship between genetic differentiation and body size across the nine island populations. Further inference is, however, difficult to make on the basis of these data, not least because multiple processes can give rise to the same estimate of genetic differentiation.

unknown causes of selection on body size
Although the large body size of wall lizards on islands might be explained, in part, by reduced mortality and increased food availability, plasticity is unlikely to account completely for the phenotypic divergence observed in this study (e.g. the body sizes observed on some of the islands are never reached in captive animals of mainland origin; T. Uller, pers. obs.). There was no increase in sexual size dimorphism on islands, suggesting that natural selection is likely to be more important than sexual selection. However, more data on sexually selected characters would be needed to rule out sexual selection completely. Many studies of island body size evolution use island size and distance from the mainland as proxies of ecological factors that cause selection on body size (reviewed by Lomolino, 2005;Itescu et al., 2018), but for these populations island size and isolation from the mainland were poor predictors of body size once genetic differentiation was taken into account. Furthermore, although we did not quantify how similar the islands are in terms of their ecology, the composition of possible predators is known to show limited variation between islands (F. Aubret, pers. obs.). Tail loss is often used as an index of predation (e.g. Bateman & Fleming, 2011), but also tends to increase at high population densities as a result of intraspecific aggression (Itescu et al., 2017). Given that tail loss was not more common on islands than on the mainland, neither predation rate nor intraspecific competition stands out as a strong candidate for explaining the large body size on islands.
Thus, the ecological factors that explain the increase in body size in these island populations remain to be identified. Although island populations of reptiles have previously been shown to vary in the extent of gigantism relative to mainland populations, there does not appear to have been any explicit test of whether or not instances of gigantism are associated with neutral genetic differentiation. It is possible that a failure to identify consistent ecological drivers of body size evolution in island reptiles (Itescu et al., 2018) is explained, in part, by a lack of attention to genetic isolation. It would be worth revisiting some of the more well-studied species to see whether genetic analysis can help to establish why gigantism is occasional even across islands that are ecologically very similar. In the case of wall lizards, the recent highquality genome assembly for P. muralis will enable a far more extensive population genetic analysis than the one conducted here, and the search for regions of the genome that are associated with particular phenotypes (Yang et al., 2018;Andrade et al., 2019). In this context, the colour morphs are interesting, because recent work has identified two main loci that determine red or yellow coloration, but failed to identify genetic associations with mosaic morphs (Andrade et al., 2019). Given that both mosaics and the pure colour morphs appear in unexpectedly high  Table S3) and in the island populations. The island populations are ordered from left to right by their average genetic distances from the mainland populations (see Table 2). Error bars represent the standard error.
frequencies on islands, studies of the genetic basis of the morphs within and across islands could help us to understand evolutionary lability in colour morph determination.
suMMary G e n e t i c d i f f e r e n t i a t i o n f r o m t h e m a i n l a n d predicted instances of gigantism across nine island populations of common wall lizards off the coast of western France. We suggest that this reflects that selection for large body size is more effective in island populations that experience very low levels of gene flow, although there might also be a role for bottlenecks. Future studies that attempt to explain why some island populations evolve gigantism should combine studies of ecology with reliable estimates of the timing of divergence and the extent of gene flow between populations.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Figure S1. A 99% maximum parsimony haplotype network for partial Cytb sequences (655 bp), constructed with TCS inference in PopART. Figure S2. Estimation of the numbers of clusters in the genetic data inferred from the STRUCTURE analysis, for all populations (A) and for northwestern populations (B). Figure S3. Population differentiation (F ST ) among island populations (A) and all sampled population pairs (B). Table S1. Microsatellite marker multiplexes used to genotype all individuals in this study. Table S2. Cytochrome b nucleotide diversity values, in terms of mean observed pairwise distance within population (π). Table S3. Population membership coefficient matrix (Q-matrix; assuming K = 3) inferred with the program STRUCTURE. The proportion with highest membership value is presented in bold for each population.