Research article | Open | Published:
Where does diversity come from? Linking geographical patterns of morphological, genetic, and environmental variation in wall lizards
BMC Evolutionary Biologyvolume 18, Article number: 124 (2018)
Understanding how phenotypic variation scales from individuals, through populations, up to species, and how it relates to genetic and environmental factors, is essential for deciphering the evolutionary mechanisms that drive biodiversity. We used two species of Podarcis wall lizards to test whether phenotypic diversity within and divergence across populations follow concordant patterns, and to examine how phenotypic variation responds to genetic and environmental variability across different hierarchical levels of biological organization, in an explicit geographic framework.
We found a general concordance of phenotypic variation across hierarchical levels (i.e. individuals and populations). However, we also found that within-population diversity does not exhibit a coherent geographic structure for most traits, while among-population divergence does, suggesting that different mechanisms may underlie the generation of diversity at these two levels. Furthermore, the association of phenotypic variation with genetic and environmental factors varied extensively between hierarchical levels and across traits, hampering the identification of simple rules to explain what yields diversity.
Our results in some cases comply with general ecological and evolutionary predictions, but in others they are difficult to explain in the geographic framework used, suggesting that habitat characteristics and other regulatory mechanisms may have a more substantial contribution in shaping phenotypic diversity.
Variation is the cornerstone of evolution. It is the raw material for the action of selection and a central component for identifying how phenotypic diversification occurs across all levels of biological organization . The unit of phenotypic variation is the individual, through which fitness is optimized via effects on survival and reproduction. Selective influences and stochastic processes like drift accumulate across individuals, yielding variation at the population level . When intraspecific differentiation is further reinforced through mechanisms of population divergence, reproductive isolation may gradually develop, which produces phenotypic variation among different species . Through this pipeline across levels of biological organization, we may conceptually trace variation from individuals, through populations of the same species, to macroevolutionary diversification across species. Due to its relevance for understanding evolution, phenotypic diversity has been explored at all these hierarchical levels. Comparative studies of phenotypic evolution combine species phenotypic data and phylogenetic trees, to explore, for instance, how shared ancestry and biological traits determine macroevolutionary diversification . Similar questions have been addressed, using a different set of analytical tools, to examine phenotypic variation at the intraspecific level from ecological, evolutionary, and conservation-focused perspectives, as it could determine the potential of a species to adapt to local conditions and variations of these across time (e.g. [5, 6]).
Common to all these approaches is the recognition of the importance of both genetic and environmental factors in shaping phenotypic variation. The connection between genotype and phenotype is a central topic in evolutionary biology, and in recent years our understanding of the genetic architecture underlying many phenotypic traits has been highly enhanced through the application of modern sequencing techniques and quantitative genetic tools . Most usually, however, indirect associations are established between genetic and phenotypic patterns, by using neutral markers to characterize genetic variation. In this type of studies, the hypothesis tested is usually that of a positive correlation between genetic and phenotypic variation . At the intraspecific level, greater degrees of (neutral) genetic divergence between populations are typically associated to lower connectivity between them, and the more differentiated these populations are predicted to be phenotypically, under the null hypothesis of random accumulation of variation. Of course, evolution, and particularly that of the phenotype, rarely occurs under stochastic processes alone. Instead, phenotypic traits are usually also correlated with the environment where organisms live, through effects on physiological, functional and ecological performance . Across populations, local adaptation and phenotypic plasticity are major frameworks for studying associations between the phenotype and the environment, and for elucidating the mechanisms underlying intraspecific phenotypic diversity [10, 11]. As the phenotypic optima are expected to be similar for individuals of the same species that inhabit similar environments, we frequently observe population-level responses to environmental variation, and therefore it is not unusual for patterns of phenotypic divergence to correlate to environmental differences across populations [5, 12].
In all these cases, additional perspectives may be gained by considering the geographic structure of variation. Because both genetic and environmental variation are frequently geographically structured, phenotypic traits also tend to exhibit some geographic pattern (e.g. ). Genetic variation in particular may exhibit geographic structuring due to the spatial dynamics of migration and population differentiation [13, 14], or as a result of direct environmental influences, as described by the isolation-by-distance family of models . In a similar fashion, environmental characteristics exhibit strong spatial autocorrelation at certain geographical scales, following diffusion mathematical models (e.g. ). Similar modelling techniques have been frequently used to examine the spatial structure of phenotypic variation by focusing on differentiation across populations [17,18,19,20]. A less appreciated aspect is that the degree of phenotypic variation within populations could also exhibit a geographic structure, determined by underlying links with genetic and environmental traits. This geographic structuring of both within-population diversity and among-population divergence needs to be statistically taken into account when trying to establish links between phenotypic and genetic or environmental factors. On the other hand, it also provides the opportunity to examine patterns of variation in a geographically explicit framework , which can aid in pinpointing patterns of variation that may remain unaccounted for by the genetic and environmental factors explicitly considered in the examined models.
Here, we investigated how phenotypic variation emerging from individuals is expressed as diversity within and divergence across populations of the same species, and at linking morphological, genetic, and environmental patterns in a geographically explicit analytical framework. For this purpose, we sampled morphological, genetic and environmental variation across populations of two species of wall lizards of the genus Podarcis. These lizards are a particularly intriguing system for exploring the processes involved in phenotypic differentiation, as they exhibit high intraspecific morphological variation, but a relatively conserved morphology across species, which results in high levels of cryptic genetic diversity [22, 23]. Using exactly the same set of individual specimens, we quantified a series of morphological traits and used genetic sampling to characterize recent population connectivity and diversity, and evolutionary history. Based on these individually-matched samples, we addressed the following questions: 1) Are patterns of morphological variation concordant at the within- and among-population levels? This is expected if phenotypic diversity responds to the same evolutionary influences – either stochastic or adaptive – at the individual and population levels. 2) Is diversity within populations and divergence among them geographically structured? Here, geographic structure, particularly if common across scales and traits, would point to direct, causal links between phenotypic, genetic, and environmental variation. 3) Which is the relative contribution of genetic and environmental factors in shaping morphological variation within and across populations? If variation in phenotypic traits is more strongly associated to population dynamics and evolutionary history, this will point to drift as a major factor shaping diversity patterns; alternatively, if it is strongly influenced by environmental factors, this will reflect plastic and adaptive mechanisms. By addressing these questions, we aimed at providing new insights into how phenotypic variation scales up from individuals to populations of the same species, and to examine whether it responds to the same genetic and environmental influences across different hierarchical levels of biological organization.
Patterns of morphological variation
The examination of the principal components (PC) of morphological variation across individuals and across populations indicated that, although the amount of variance captured differed between the two analyses, variation patterns were transversal to both hierarchical levels for all traits, except for head shape, in both species. Indeed, the structure of phenotypic variation as expressed by PC eigenvectors of the individual-based and population-based analyses did not differ significantly for neither linear biometric traits (P. bocagei: θ = 13.8°, p = 0.918; P. vaucheri: θ = 15.0°, p = 0.125), nor scalation (P. bocagei: θ = 32.2°, p = 0.155; P. vaucheri: θ = 14.6°, p = 0.286). The main direction of biometric variation was related to body size, both at individual and population analyses, as indicated by high, equal-sign loadings of all biometric traits on PC1 in both species (Fig. 1). This also resulted in a high eigenvalue for the first PC in both analyses and species, with more than 60% of variance explained. For scalation traits, the first PC generally correlated to higher scale counts for all examined traits, with a lower relative contribution of VSN, again under a common pattern across levels (Fig. 1). Instead, for head shape as quantified using geometric morphometrics, the direction of principal component eigenvectors differed significantly across individuals vs. across populations for both species (P. bocagei: θ = 50.7°, p = 0.001; P. vaucheri: θ = 50.4°, p = 0.001): while in both cases the first PC was related to a relative shortening of the snout, individual variation was also associated to a widening of the posterior region of the head, whereas population differentiation was rather related to a reduction of this area (Fig. 1).
Although interpolations were mostly accurate (Additional file 1: Table A8.1), within-population diversity values for studied traits did not exhibit a geographic structure, with the exception of body size in P. bocagei (Additional file 1: Figure A8). Spatial variation in interpolated traits was inexistent or extremely reduced and frequently presented spatial artefacts (e.g. body size in P. vaucheri; Additional file 1: Figure A8).
By contrast, phenotypic differentiation among populations exhibited a coherent geographic structure for most studied traits. We globally obtained a reasonable number of accurate rasters to illustrate general patterns of geographic variability, except for the case of scalation traits in P. bocagei (three rasters only; Additional file 1: Table A8.2). In both species, some populations were identified as more dissimilar compared to the rest (i.e. they exhibited higher distances in interpolations; Fig. 2), but population divergence patterns varied extensively across traits. For instance, some coastal and inland mountain populations of P. bocagei were more dissimilar than the rest in body size, body shape and head shape; while for P. vaucheri, Rif and Western High Atlas populations differed from the rest in scalation, body size (High Atlas populations only) and shape (Rif populations only) (Fig. 2).
Factors shaping morphological variation
Generalized least-squares (GLS) analyses indicated a significant association between levels of within-population phenotypic and genetic diversity, as well as a significant potential contribution of environmental factors (Additional file 2). Levels of population genetic diversity were relevant for explaining phenotypic diversity for all phenotypic traits – except for head shape – in both species. Interestingly, the association between phenotypic and genetic diversity within populations was not always positive. Indeed, while a positive correlation was observed for scalation in P. bocagei and body size in P. vaucheri, the remaining significant associations (i.e. for body size and shape in P. bocagei; and for body shape and scalation in P. vaucheri) were negative (Fig. 3). Population phenotypic diversity was also influenced by environmental factors, though such associations varied between species (Fig. 3). Variables related to temperature always had a positive effect on phenotypic diversity, while the effect of precipitation was limited to a negative influence on head shape diversity in P. bocagei.
Generalized Dissimilarity Modelling (GDM) indicated that both genetic and environmental factors influence phenotypic differentiation across populations, but the relative contribution of different explanatory variables varied across traits and between species (Fig. 3). Population genetic differentiation at the mtDNA level was associated with the degree of population divergence in scalation in both species, but recent population connectivity, as captured using microsatellite data, did not exhibit a significant contribution. Population divergence in scalation was also influenced by environmental variation across sites and specifically by precipitation in the driest month (Pdry12) and maximum temperature in the warmest month (maxTwarm12), in P. vaucheri, but not in P. bocagei. Variation across populations in body size and shape was associated to environmental, but not genetic, variation. For body size, slope, mean temperature of the wettest quarter, temperature range and productivity (NDVI) were relevant for explaining population divergence in P. vaucheri but not in P. bocagei, for which a significant model could not be obtained based on the considered explanatory variables. By contrast, variation across populations of P. bocagei in body shape was associated mainly to slope, with an additional significant contribution of geographic distance, whereas precipitation in the driest month was the only variable contributing to explaining variation in body shape in P. vaucheri. Finally, divergence across populations in head shape could not be associated to any of the examined genetic or environmental factors in neither species.
Understanding how phenotypic variation is structured across hierarchical levels of biological organization, and determining how it is associated to underlying genetic and environmental factors, is essential for our comprehension of the evolutionary mechanisms that drive biodiversity patterns. Both at the intra- and at the inter-specific levels, neutral evolution through drift-related processes has been predicted to bind genetic and phenotypic variation together, which may translate into common patterns of within-population diversity, among-population divergence and, ultimately, cross-species diversification . At the other extreme of this spectrum, selective mechanisms may cause plastic and adaptive responses that are identifiable as correlations between the environmental niche and phenotypic traits.
The structure of variation from individuals to populations
The accumulation of phenotypic diversity across hierarchical levels may be predicted to cause concordant patterns of variation across individuals, populations and even species, if the mechanisms driving diversification are uniform across these levels. Our results partially support this idea, as general patterns of morphological variation were found to be uniform across individuals and among populations (Fig. 1). For biometric traits, body size appears as a strong integrating factor, which dominates variation both among individuals and across populations. This is a common and well-known pattern in intraspecific morphological studies, as ontogenetic growth , but also possible sampling effects , easily yield high levels of size variation, particularly across individuals . Interestingly, the same was the case in our system for scalation, which also exhibited similar structures of variation at the individual and population level, but not for head shape, where principal vectors of variation differed significantly for individual and population-level analyses (Fig. 1). These results add to our understanding of how phenotypic diversity scales up from individuals, through populations, to among species. Previous studies have compared the structure of between intra- and interspecific differentiation, suggesting that concordance at these levels is prone to occur as a response to selection . Our results suggest that this scaling up of morphological variation can also be traced to the individual level, and individual phenotypic variation then translates into divergence among populations, which has been suggested to have potential effects on ecological success .
The uniformity in the general direction of morphological variation across hierarchical levels contrasts with its geographic structure, where visibly dissimilar patterns were found when examining divergence across populations and diversity within them. Indeed, while coherent geographic patterns where observed for divergence in most morphological traits (Fig. 2), intra-population diversity levels only exhibited geographic structuring in one trait (i.e. body size; Additional file 1: Table A8.2, Figure A8). This result may be partially related to the covariance structure among individual measured traits, which is generally stronger at the population than at the individual level, as expressed by high eigenvalues concentrated across the first PC axes of population-level analyses. This means that while variation across individuals is more prone to vary in direction across traits, population differentiation in some traits is also accompanied by differentiation in others. Therefore, while principal components of variation are concordant at the individual and population levels, there is no agreement between them in geographic structure. This result suggests that, whatever the mechanisms that underlie the links observed across hierarchical levels of phenotypic diversity, these are not associated to geography. The presence of geographic structure in population divergence, particularly in scalation traits and body shape, is typically evoked as providing evidence for local adaptation or phenotypic plasticity as a response to environmental conditions (e.g. [30, 31], see also below), but the geographic structure of within-population diversity has not been frequently explored – to our knowledge. Based on population dynamics, a dissociation between the amount of within-population diversity and among-population divergence may be expected as a result of reduced effective population sizes: if intra-population diversity is low, divergence among populations may increase due to drift and random phenotypic oscillation, as typically observed in e.g. island populations [32, 33], or during the establishment of invasive populations . Nevertheless, if drift alone drives phenotypic variation, divergence among populations is not expected to exhibit a geographic structure, as that observed in our data. Associations of within-population diversity and among-population divergence with genetic and environmental factors shed more light into the possible causes underlying phenotypic diversification.
Diversity within populations: Recent phenotypic dynamics
Our results suggest that variation across individuals within populations results from an amalgamation of recent population dynamics and environmental factors. The detailed mechanisms that drive phenotypic variation at the individual level seem to act at a very local scale, as within-population diversity does not exhibit a coherent geographic pattern in most studied traits (Additional file 1: Table A8.2). Also, such mechanisms appear quite labile, as patterns vary extensively across traits and between species (Fig. 3). Population genetic diversity estimated based on microsatellite analyses had a significant effect on phenotypic diversity for all examined traits except for head shape. Such an association between genetic and morphological diversity is expected under several evolutionary scenarios. First, higher genetic diversity may be intuitively expected to be linked to higher phenotypic diversity, as the first would provide the genetic background for the expression of different phenotypes. On the contrary, a negative association can also be observed if a depletion of genetic diversity interferes with basic developmental mechanisms involved in phenotypic regulation .
Our results support the idea that developmental mechanisms may be important in shaping levels of phenotypic diversity. In addition to a negative association between genetic and phenotypic diversity, we also found a positive contribution of temperature-related variables. Such an effect may also be mediated by development: as oviparous ectotherms, wall lizards are quite exposed to environmental effects, particularly after oviposition and throughout the post-natal period. Studies of variance within and across individuals in lizards exposed to different environmental conditions have shown that some morphological traits are particularly prone to temperature and other environmental effects during development [36,37,38,39,40], with potential relevance for shaping intrapopulation diversity, as well as population divergence and patterns of geographic variation .
Divergence across populations: Different traits, different stories
Patterns of divergence across populations varied extensively across traits, and were associated to different genetic and environmental factors. For scalation, we found that phylogeographic structure and evolutionary history partially explained divergence among populations in both species. This reinforces the long-established perception that these traits are useful for investigating past and emergent differentiation  and for making systematic inferences (e.g. ). Still, variation in scalation traits does not seem to be completely neutral, as it has been related to environmental and ecogeographic variation [17, 19]. Indeed, the number of scales has been frequently associated to thermoregulation and water balance in lizards [43, 44], and it is known to co-vary with the thermal and hydric characteristics of the environment, both within [31, 32] and across species [28, 30] in several lizard groups. Accordingly, we found that divergence among populations of P. vaucheri in scalation is associated to both temperature and precipitation (Fig. 3), and specifically to conditions that may be limiting for these species, as the variables that significantly contributed to explaining morphological variation represent extremes of environmental factors.
This contribution of restrictive environmental variables, together with the lack of significant effects of environmental factors on scalation divergence in P. bocagei, suggest that the magnitude of variation examined here might also be contributing to the observed results. For instance, the distribution range of P. vaucheri, and especially the populations examined here, encompass a wider range of environmental conditions , which may facilitate the detection of statistically significant effects. This does not seem to be limited to scalation, as significant environmental effects on body size and shape divergence were also identified mostly for populations of P. vaucheri (Table 1). However, coherent patterns of geographic variation were identified for all these traits in both species, where some coastal and high-mountain populations exhibited a higher distinctiveness, depending on the trait under examination (Fig. 2). This is in accordance with previous observations in other species of the P. hispanica complex (i.e. P. carbonelli ), and it suggests that possibly other environmental factors, not examined here, or the spatial structuring of genetic variation, may contribute in shaping geographic morphological variation.
On the other hand, head shape divergence does not exhibit an easily explicable geographic pattern, and it could not be associated to either genetic or environmental factors in any of the two species. Local variation observed in our data, as well as geographic gradients not related to environmental differences between sampling sites, suggest that other factors, not investigated here and working at the population level might be important. First, the ecophysical characteristics of the sites sampled may have a contribution. Indeed, head shape is known to exhibit fast evolutionary responses to structural habitat in these lizards, and it varies both across species with different ecological habits [46, 47] and between populations that occupy different microhabitats . Second, developmental processes such as allometry and developmental stability are known to promote fast evolutionary change in head shape throughout ontogeny (e.g. [31, 37]) with important influences on phenotypic diversity .
Overall, our results from two species of wall lizards suggest that establishing simple, global rules to link the phenotype and either genetic or environmental variation is not straightforward. Despite concordant general patterns across hierarchical levels (i.e. individuals and populations), within-population diversity and among-population divergence vary extensively across phenotypic traits, and they exhibit varying levels of geographic structuring and of association to genetic and environmental factors. While some of the traced links comply with ecological predictions and evolutionary intuition, others are difficult to explain in the geographic framework used, suggesting that habitat characteristics and other regulatory mechanisms – not examined here –, may be more influential in shaping patterns of phenotypic diversity.
Podarcis wall lizards are an intriguing model for exploring how variation scales from individuals to populations, eventually determining diversity across species. Ranging around the Mediterranean Basin, species of this genus occupy a wide variety of ecosystems . Their morphology has been shown to respond to different characteristics of the environment including habitat structure [46,47,48], insularity [51, 52], and human-mediated disturbance [36, 37]. Further, geographically-structured morphological variation has been reported within some species (e.g. ), but the causes underlying such structuring have not been explored.
Here, we investigated two species of wall lizards, Podarcis bocagei from Northwestern Iberia, and P. vaucheri from Morocco and Algeria (sensu ). For each of these species, we sampled a total of nine populations (Fig. 4, Additional file 3), in the years between 2005 and 2009. In each population, we captured around 30 adult individuals of both sexes (Additional file 3). For each individual, we recorded several biometric traits (see below); we took high-resolution photographs to record head shape using geometric morphometrics, and to quantify scalation traits; and we removed a small tip of the tail, which was stored in 96% ethanol to be used for genetic analyses. All individuals were processed in the field and released immediately after.
To characterize morphological diversity within and divergence across populations, we considered four sets of traits: body size, body shape, head shape using geometric morphometrics, and scalation. We chose these traits due to their relevance for our hypotheses, as they may be associated to both stochastic (genetic) and adaptive (environmental) factors. Body size, body shape and head shape are known to be of relevance for a series of physiological, ecological and social functions in lizards. On the other hand, scale counts are frequently related to environmental variability, showing geographic variation [17, 19], but they have also been extensively used for systematics (e.g. ), which suggests at least some association to genetic divergence.
To quantify body size and shape, we recorded seven linear biometric traits to the closest 0.01 mm using electronic calipers: snout-vent length, trunk length, head dimensions (i.e. length, width and height), and fore and hind limb length (Fig. 4). These traits capture the general biometric properties of each individual and they are associated to whole-organism performance [53, 54], thus being potentially associated to environmental variation. Because Podarcis wall lizards exhibit marked sexual dimorphism in body size and shape [55, 56], we first centered all observations to the corresponding sex-mean. The first principal component of the matrix of sex-mean-centered biometric traits was then treated as a measure of body size, while the remaining PC axes were treated as a multivariate, size-corrected measure of body shape (see also results).
Because head shape has been shown to be an important component of ecologically-related variation in these lizards (i.e. [46, 48]), we also used geometric morphometrics (GM) to quantify variation in the position of 24 two-dimensional landmarks on the dorsal head surface (Fig. 4). A detailed account of all GM data processing can be found in Additional file 4.
To capture variation in scalation, in each individual we counted the number of gular (GSN) and ventral (VSN) scales, femoral pores (FPN), and supraciliary granules (SCGN) (Fig. 4). These traits were chosen because they are all continuously valued, and they exhibit enough variation across individuals of Podarcis to be informative (i.e. they are not fixed ). All scalation traits were sex-mean-centered prior to statistical analyses.
Based on each of the above sets of traits, we quantified the level of diversity within each population using the mean distance to the population centroid . In addition, we represented the degree of morphological divergence among populations of each species by calculating multivariate (or univariate, in the case of body size), pairwise Euclidean distances between population mean values. All raw data used to calculate diversity within and divergence between populations are available in Additional file 5 (P. bocagei) and Additional file 6 (P. vaucheri).
Population genetics and phylogeographic structure
To characterize recent population dynamics, we genotyped all individuals examined for morphological analyses for nine microsatellite loci, initially developed for P. bocagei . Based on the data obtained, we characterized within-population genetic diversity using the number of alleles per population (Na) and expected heterozygosity (Hexp). We also characterized differentiation between populations by calculating the Cavalli-Sforza and Edwards (1967)  genetic distance (genD, with a correction for null alleles, see Additional file 7) and Fst values according to .
In addition, we examined a 623 bp fragment of the NADH dehydrogenase subunit 4 (ND4) mtDNA gene to describe the phylogeographic structure and deeper evolutionary history of the populations under study. Based on these data, we described the genetic distance among populations using Dxy, the average number of nucleotide substitutions per site between populations ; and estimated differentiation using the Hudson et al. (1992) estimator of Fst , and calculated diversity measures per locality (haplotype diversity, Hd, and nucleotide diversity, π).
A detailed account of laboratory and analytical procedures followed to obtain genetic data, and the corresponding descriptors of population diversity and differentiation, are provided in Additional file 7. Raw individual multilocus genotypes for microsatellite data are available in Additional file 8.
To characterize the average climatic environment of the studied populations, and fluctuations to it, we clipped the 19 bioclimatic variables (at 30 arc-seconds of resolution) included in the Worldclim database  to our study area. To avoid variable collinearity, we reduced this set to six variables that exhibited low correlations (R < 0.7), and which are known to be biologically informative for the species under study [45, 63], including: isothermality (isoT), maximum temperature of the warmest month (maxTwarm12), annual range of temperature (Trange), mean temperature of the wettest quarter (meanTwet4), precipitation of the wettest and driest month (Pwet12 and Pdry12, correspondingly). In addition to these climatic variables, and in order to characterize the local productivity of habitats occupied by wall lizards, we extracted the values of a Normalized Difference Vegetation Index (NDVI) for each locality (see ). We obtained these values by averaging the satellite-derived temporal raster series of the MODIS NDVI imagery (from GIMMS Global Agricultural Monitoring System, http://glam1.gsfc.nasa.gov/; at 250 m of resolution) for the presumed activity period of wall lizards during our sampling periods (late February to early November, in the years between 2005 and 2009). Furthermore, to characterize topographic conditions of each locality we derived the slope variable from altitude (also downloaded from the Worldclim database), using the “slope” function of ArcGIS . Finally, we quantified environmental suitability for each of the studied populations extracting values from ecological niche-based models (ENM) calculated for the whole range of each species using the Maxent software . Further details on ENMs can be consulted in Additional file 9.
The structure of morphological variation
To investigate whether the structure of morphological variation was concordant across individuals and across populations of each species, we performed principal components analyses of each multivariate trait block separately (i.e. biometry, head shape, and scalation), considering a) the matrix of individual observations and b) population means for each trait. We then permuted population-centered residuals of individual observations across populations to examine whether the first principal component of individual morphological variation aligned with that of population divergence, by testing whether the angle θ between these two PC vectors was larger than that expected at random.
Modelling diversity within populations
To visualize the spatial structure of morphological diversity within populations, we spatially interpolated within-population mean distance to centroid for each of the four explored trait sets for each species, using the kriging interpolation method in ArcGIS . A detailed description of the methods used for interpolations is available in Additional file 1.
To evaluate which genetic and environmental factors may contribute in shaping levels of intrapopulation phenotypic diversity, we used Generalized Least Squares (GLS) to fit single-predictor linear models while taking the spatial structure of our sampling into account. Specifically, we estimated GLS regressions with population phenotypic diversity (i.e. mean distance to the centroid, see above) for each of the morphological traits examined as the response; and each of the genetic diversity and environmental factors described above as predictors. We preferred this simple, pairwise approach, instead of fitting multiple-predictor GLS models, to obtain an overall view of which factors may influence intraspecific phenotypic diversity without confounding statistical interactions among predictors.
The most adequate autocorrelation structure for each response-predictor combination was determined by fitting models using a linear, exponential, Gaussian, or rational quadratic structure, and then comparing their fit to the data using Akaike’s Information Criterion (AIC). The model with the lowest AIC for each combination was considered to optimize the expected covariance among observations due to spatial proximity. When GLS models suggested a significant contribution of a genetic or environmental predictor on the level of population diversity for a trait, we examined the strength of such an association by inspecting the Pearson correlation between response and predictor variables.
Modelling divergence across populations
Similar to what we did for intra-population diversity, we used kriging of inter-population phenotypic distances to explore the geographic structure of phenotypic divergence among populations (Additional file 1). In continuation, we used Generalized Dissimilarity Modelling (GDM)  to investigate whether recent population connectivity, evolutionary history, and/or environmental factors contribute in explaining morphological divergence across populations, while taking spatial structure into account. This technique allows to model dissimilarities in morphology, represented as the matrix of pairwise distances between population means, with pairwise matrices representing geographic proximity, genetic differentiation, and environmental dissimilarity across populations. It is particularly suitable for our purposes, because (1) it facilitates the application of regression-type analyses on data that are naturally expressed as distances (e.g. genetic differentiation); (2) it allows the inclusion of non-linear relationships between response and explanatory dissimilarity matrices ; and (3) by working with multivariate distance matrices, one can evaluate the sources of divergence in multivariate phenotypes examining the response of several phenotypic traits simultaneously, similarly to approaches implemented to genomic data .
We fit GDMs to inter-population Euclidean distance matrices for each of the investigated phenotypic traits and for each species separately using the gdm package for R . Predictors included the matrices of Fsts and genetic distances, based on microsatellite (recent population connectivity) and mtDNA data (evolutionary history); the nine environmental variables described above; and the matrix of pairwise geographic distances between sampling sites. We performed a backward-stepwise selection of variables using Monte Carlo permutations to evaluate predictor significance [67, 71]. To quantify the relative contribution of each predictor in explaining variance in phenotypic differentiation for each of the response traits we summed the coefficients of the I-splines, which corresponds to the maximum height attained by the response curve for each predictor [67, 69].
To account for the possibility that environmental resistance, instead of geographic distance, may best represent the expected covariance among populations for phenotypic traits, we repeated all GDM analyses considering pairwise matrices of environmental connectivity, calculated in Circuitscape  based on the inverse of habitat suitability calculated for each species using ENMs in Maxent (see above, and Additional file 9). However, because the results obtained did not vary substantially when using connectivity instead of geographic distance, we only present the results using geographic distances.
Akaike’s Information Criterion
Average number of nucleotide substitutions per site between populations
Ecological niche-based models
Number of femoral pores
Generalized Dissimilarity Modelling
Number of gular scales
Maximum temperature in the warmest month
Mean temperature of the wettest quarter
Number of alleles
Normalized Difference Vegetation Index
Precipitation in the driest month
precipitation of the wettest month
Number of supraciliary granules
Annual range of temperature
Number of ventral scales
Hallgrímsson B, Hall BK. Variation and variability: central concepts in biology. In: Hallgrímsson B, Hall BK, editors. Variation: a central concept in biology. Burlington: Elsevier Academic Press; 2005. p. 1–7.
Lande R. Natural selection and random genetic drift in phenotypic evolution. Evolution. 1976;30:314–34.
Turelli M, Barton NH, Coyne JA. Theory and speciation discuss recent examples of progress in each of these areas. Trends Ecol Evol. 2001;16:330–43.
Harvey PH, Pagel MD. The comparative method in evolutionary biology. Oxford: Oxford University Press; 1991.
Thomassen HA, Buermann W, Milá B, Graham CH, Cameron SE, Schneider CJ, et al. Modeling environmentally associated morphological and genetic variation in a rainforest bird, and its application to conservation prioritization. Evol Appl. 2010;3:1–16.
González-Suárez M, Bacher S, Jeschke JM. Intraspecific trait variation is correlated with establishment success of alien mammals. Am Nat. 2015;185:737–46.
Wray GA. Genomics and the evolution of phenotypic traits. Annu Rev Ecol Evol Syst. 2013;44:51–72.
Orgogozo V, Morizot B, Martin A. The differential view of genotype-phenotype relationships. Front Genet. 2015;6:179.
Arnold SJ. Morphology, performance and fitness. Am Zool. 1983;23:347–61.
Ghalambor CK, McKay JK, Carroll SP, Reznick DN. Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments. Funct Ecol. 2007;21:394–407.
Mullen LM, Vignieri SN, Gore JA, Hoekstra HE. Adaptive basis of geographic variation: genetic, phenotypic and environmental differences among beach mouse populations. Proc R Soc Lond B Biol Sci. 2009;276:3809–18.
Campbell P, Pasch B, Pino JL, Crino OL, Phillips M, Phelps SM. Geographic variation in the songs of neotropical singing mice: testing the relative importance of drift and local adaptation. Evolution 2010; 64(7): 1955–1972.
Pinho C, Kaliontzopoulou A, Harris DJ, Ferrand N. Recent evolutionary history of the Iberian endemic lizards Podarcis bocagei (Seoane, 1884) and Podarcis carbonelli Pérez-Mellado, 1981 (Squamata: Lacertidae) revealed by allozyme and microsatellite markers. Zool J Linnean Soc. 2011;162:184–200.
Martínez-Freiría F, Velo-Antón G, Brito JC. Trapped by climate: interglacial refuge and recent population expansion in the endemic Iberian adder Vipera seoanei. Divers Distrib. 2015;21:331–44.
Wang IJ, Bradburd GS. Isolation by environment. Mol Ecol. 2014;23:5649–62.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.
Martinez-Freiria F, Santos X, Pleguezuelos JM, Lizana M, Brito JC. Geographical patterns of morphological variation and environmental correlates in contact zones: a multi-scale approach using two Mediterranean vipers (Serpentes). J Zool Syst Evol Res. 2009;47:357–67.
Kaliontzopoulou A, Carretero MA, Sillero N. Geographic patterns of morphological variation in the lizard Podarcis carbonelli, a species with fragmented distribution. Herpetol J. 2010;20:41–50.
Martínez-Freiría F, Brito JC. Integrating classical and spatial multivariate analyses for assessing morphological variability in the endemic Iberian viper Vipera seoanei. J Zool Syst Evol Res. 2013;51:122–31.
Santos X, Vidal-García M, Brito JC, Fahd S, Llorente GA, Martínez-Freiría F, Parellada X, Pleguezuelos JM, Sillero N. Phylogeographic and environmental correlates support the cryptic function of the zigzag pattern in a European viper. Evol Ecol. 2014;28:611–26.
Legendre P. Spatial autocorrelation: trouble or new paradigm? Ecology. 1993;74:1659–73.
Kaliontzopoulou A, Pinho C, Harris DJ, Carretero MA. When cryptic diversity blurs the picture: a cautionary tale from Iberian and north African Podarcis wall lizards. Biol J Linn Soc. 2011;103:779–800.
Kaliontzopoulou A, Carretero MA, Llorente GA. Morphology of the Podarcis wall lizards (Squamata: Lacertidae) from the Iberian Peninsula and North Africa: patterns of variation in a putative cryptic species complex. Zool J Linnean Soc. 2012;164:173–93.
Hansen TF, Martins EP. Translating between microevolutionary process and macroevolutionary patterns: the correlation structure of interspecific data. Evolution. 1996;50:1404–17.
Gould SJ. Allometry and size in ontogeny and phylogeny. Biol Rev. 1966;41:587–640.
Prothero J. Methodological aspects of scaling in biology. J Theor Biol. 1986;118:259–86.
Dreyer AP, Ziabari OS, Swanson EM, Chawla A, Frankino WA, Shingleton AW. Cryptic individual scaling relationships and the evolution of morphological scaling. Evolution. 2016;70:1703–16.
Calsbeek R, Knouf JH, Smith TB. Variation in scale numbers is consistent with ecologically based natural selection acting within and between lizard species. Evol Ecol. 2006;20:377–94.
Forsman A, Wennersten L. Inter-individual variation promotes ecological success of populations and species: evidence from experimental and comparative studies. Ecography. 2016;39:630–48.
Oufiero CE, Gartner GEA, Adoplh SC, Garland T. Latitudinal and climatic variation in body size and dorsal scale counts in Sceloporus lizards: a phylogenetic perspective. Evolution. 2011;65:3590–607.
Jaffe AL, Campbell-Staton SC, Losos JB. Geographical variation in morphology and its environmental correlates in a widespread north American lizard, Anolis carolinensis (Squamata: Dactyloidae). Biol J Linn Soc. 2016;117:760–74.
Brown RP, Thorpe RS, Báez M. Parallel within-island microevolution of lizards on neighbouring islands. Nature. 1991;352:60–2.
Jordan MA, Snell HL. Historical fragmentation of islands and genetic drift in populations of Galápagos lava lizards (Microlophus albemarlensis complex). Mol Ecol. 2008;17:1224–37.
Kolbe JJ, Larson A, Losos JB. Differential admixture shapes morphological variation among invasive populations of the lizard Anolis sagrei. Mol Ecol. 2007;16:1579–91.
Klingenberg CP, Nijhout HF. Genetics of fluctuating asymmetry: a developmental model of developmental instability. Evolution. 1999;53:358–75.
Lazić MM, Kaliontzopoulou A, Carretero MA, Crnobrnja-Isailović J. Lizards from urban areas are more asymmetric: using fluctuating asymmetry to evaluate environmental disturbance. PLoS One. 2013;8:e84190.
Lazić MM, Carretero MA, Crnobrnja-Isailović J, Kaliontzopoulou A. Effects of environmental disturbance on phenotypic variation: an integrated assessment of canalization, developmental stability, modularity, and allometry in lizard head shape. Am Nat. 2015;185:44–58.
Braña F, Ji X. Influence of incubation temperature on morphology, locomotor performance, and early growth of hatchling wall lizards (Podarcis muralis). J Exp Zool. 2000;286:422–33.
Ji X, Lin LH, Luo LG, Lu HL, Gao JF, Han J. Gestation temperature affects sexual phenotype, morphology, locomotor performance, and growth of neonatal brown forest skinks, Sphenomorphus indicus. Biol J Linn Soc. 2006;88:453–63.
Lazić MM, Carretero MA, Crnobrnja-Isailović J, Kaliontzopoulou A. Postnatal dynamics of developmental stability and canalization of lizard head shape under different environmental conditions. Evol Biol. 2016;43:368–79.
Qualls CP, Shine R. Geographic variation in lizard phenotypes: importance of the incubation environment. Biol J Linn Soc. 1998;64:477–91.
Arnold EN, Arribas O, Carranza S. Systematics of the Palaearctic and oriental lizard tribe Lacertini (Squamata: Lacertidae: Lacertinae), with descriptions of eight new genera. Zootaxa. 2007;1430:1–86.
Alibardi L. Adaptation to the land: the skin of reptiles in comparison to that of amphibians and endotherm amniotes. J Exp Zool Part B - Mol Dev Evol. 2003;298B:12–41.
Soulé M, Kerfoot WC. On the climatic determination of scale size in a lizard. Syst Zool. 1972;21:97–105.
Caeiro-Dias G, Luís C, Pinho C, Crochet P-A, Sillero N, Kaliontzopoulou A. Lack of congruence of genetic and niche divergence in Podarcis hispanicus complex: allopatric speciation as the main driving factor. J Zool Syst Evol Res. In press
Kaliontzopoulou A, Carretero MA, Adams DC. Ecomorphological variation in male and female wall lizards and the macroevolution of sexual dimorphism in relation to habitat use. J Evol Biol. 2015;28:80–94.
Urošević A, Ljubisavljević K, Jelić D, Ivanović A. Variation in the cranium shape of wall lizards (Podarcis spp.): effects of phylogenetic constraints, allometric constraints and ecology. Zoology. 2012;115:207–16.
Kaliontzopoulou A, Carretero MA, Llorente GA. Intraspecific ecomorphological variation: linear and geometric morphometrics reveal habitat-related patterns within Podarcis bocagei wall lizards. J Evol Biol. 2010;23:1234–44.
Piras P, Salvi D, Ferrara G, Maiorino L, Delfino M, Pedde L, Kotsakis T. The role of post-natal ontogeny in the evolution of phenotypic diversity in Podarcis lizards. J Evol Biol. 2011;24:2705–20.
Harris DJ, Arnold EN. Relationships of wall lizards, Podarcis (Reptilia: Lacertidae) based on mithocondrial DNA sequences. Copeia. 1999;1999:749–54.
Raia P, Guarino FM, Turano M, Polese G, Rippa D, Carotenuto F, et al. The blue lizard spandrel and the island syndrome. BMC Evol Biol. 2010;10:289.
Runemark A, Hansson B, Pafilis P, Valakos ED, Svensson EI. Island biology and morphological divergence of the Skyros wall lizard Podarcis gaigeae: a combined role for local selection and genetic drift on color morph frequency divergence? BMC Evol Biol. 2010;10:269.
Kaliontzopoulou A, Adams DC, Meijden A, Perera A, Carretero MA. Relationships between head morphology, bite performance and ecology in two species of Podarcis wall lizards. Evol Ecol. 2012b;26:825–45.
Kaliontzopoulou A, Bandeira V, Carretero MA. Sexual dimorphism in locomotor performance and its relation to morphology in wall lizards (Podarcis bocagei). J Zool. 2013;289:294–302.
Kaliontzopoulou A, Carretero MA, Llorente GA. Head shape allometry and proximate causes of head sexual dimorphism in Podarcis lizards: joining linear and geometric morphometrics. Biol J Linn Soc. 2008;93:111–24.
Kaliontzopoulou A, Carretero MA, Llorente GA. Sexual dimorphism in traits related to locomotion: ontogenetic patterns of variation in Podarcis wall lizards. Biol J Linn Soc. 2010;99:530–43.
Foote M. Contributions of individual taxa to overall morphological disparity. Paleobiology. 1993;19:403–19.
Pinho C, Sequeira F, Godinho R, Harris DJ, Ferrand N. Isolation and characterization of nine microsatellite loci in Podarcis bocagei (Squamata: Lacertidae). Mol Ecol Notes. 2004;4:286–8.
Cavalli-Sforza LL, Edwards AWF. Phylogenetic analysis: models and estimation procedures. Evolution. 1967;32:550–70.
Weir BS. Genetic data analysis II. Sunderland, Mass: Sinauer Associates; 1996.
Nei M. Molecular evolutionary genetics. New York: Columbia University Press; 1987.
Hudson RR, Slatkin M, Maddison WP. Estimation of levels of gene flow from DNA sequence data. Genetics. 1992;132:583–9.
Kaliontzopoulou A, Brito JC, Carretero MA, Larbes S, Harris DJ. Modelling the partially unknown distribution of wall lizards (Podarcis) in North Africa: ecological affinities, potential areas of occurrence, and methodological constraints. Can J Zool. 2008;86:992–1001.
Pettorelli N, Vik JO, Mysterud A, Gaillard JM, Tucker CJ, Stenseth NC. Using the satellite-derived NDVI to assess ecological responses to environmental change. Trends Ecol Evol. 2005;20:503–10.
ESRI. ArcGIS Desktop: Release 10.3. Redlands CA; 2014.
Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006;190:231–59.
Ferrier S, Manion G, Elith J, Richardson K. Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment. Divers Distrib. 2007;13:252–64.
Thomassen HA, Cheviron ZA, Freedman AH, Harrigan RJ, Wayne RK, Smith TB. Spatial modelling and landscape-level approaches for visualizing intra-specific variation. Mol Ecol. 2010;19:3532–48.
Fitzpatrick MC, Keller SR. Ecological genomics meets community-level modelling of biodiversity: mapping the genomic landscape of current and future environmental adaptation. Ecol Lett. 2014;18:1–16.
Manion G, Lisk M, Ferrier S, Nieto-Lugilde D, Fitzpatrick MC. gdm: Functions for Generalized Dissimilarity Modeling. R package version 1.3.2. https://CRAN.R-project.org/package=gdm; 2017.
Fitzpatrick MC, Sanders NJ, Ferrier S, Longino JT, Weiser MD, Dunn RR. Forecasting the future of biodiversity: a test of single- and multi-species models for ants in North America. Ecography. 2011;34:836–47.
McRae BH, Dickson BG, Keitt TH, Shah VB. Using circuit theory to model connectivity in ecology and conservation. Ecology. 2008;10:2712–24.
We thank I. Freitas, V. Gomes, A. Perera, C. Rato, and Y. Simone for useful discussion on a previous version of the manuscript; S. Carranza and an anonymous reviewer for constructive comments during the review process; R. Godinho and S. Lopes for help with microsatellite genotyping; and colleagues from CIBIO for help during fieldwork.
AK and CP were supported by IF contracts (IF/00641/2014/CP1256/CT0008 and IF/01597/2014/CP1256/CT0009 respectively), and FMF by a post-doctoral grant (SFRH/BPD/109119/2015), all under the Programa Operacional Potencial Humano – Quadro de Referência Estratégico Nacional funds from the European Social Fund and Portuguese Ministério da Educação e Ciência. The funding body had no involvement in the design of the study, the collection, analysis, and interpretation of data, or in writing the manuscript.
Availability of data and materials
All data generated or analysed during this study are included in this published article [and its additional files].
Consent for publication.
Data collection for this study was conducted following national laws and the guidelines of the Ethical Commission of the University of Porto (Comissão de Ética da Universidade do Porto, CEUP). All specimens used were captured under scientific capture and handling permits from Instituto para a Conservação da Natureza e Biodiversidade (ICNB, Portugal) and from Le Haut Commissaire aux Eaux et Forets et à la Lutte Contre la Désertification (Morocco) to AK.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Interpolation of morphological diversity within and differentiation across populations, provides details on the methods used and the results obtained to study the geographic structure of morphological variation. (DOCX 169 kb)
Results of Generalized Least Squares analysis performed to identify environmental and genetic variables that contribute significantly in explaining differences across populations in levels of intra-populational morphological diversity while taking geographic distance into account. (DOCX 73 kb)
Geographic coordinates in the WGS1984 system, number of individuals sampled per sex (Nm: males, Nf: females) and in total (Ntot), and abbreviation (CODE) for the examined populations. (DOCX 16 kb)
Geometric morphometric procedures for quantifying variation in dorsal head shape, provides a detailed description of all steps involved in obtaining head shape data, including landmark definitions and digitizing, superimposition of landmark coordinates and other pre-processing applied to GM data. (DOCX 13 kb)
Sex-corrected individual values for all morphological traits for populations of P. bocagei, provides the data used for all analyses for this species, after mean-centering by sex, including biometric and scalation variables as described in Methods, and superimposed, sex-mean-centered landmark coordinates (X1 – Y24), as well as the population of origin in accordance with Additional file 3. (TXT 218 kb)
Sex-corrected individual values for all morphological traits for populations of P. vaucheri, provides the data used for all analyses for this species, after mean-centering by sex, including biometric and scalation variables as described in Methods, and superimposed, sex-mean-centered landmark coordinates (X1 – Y24), as well as the population of origin in accordance with Additional file 3. (TXT 203 kb)
Genetic analyses, provides a detailed account of of laboratory and analytical procedures followed to obtain genetic data, and the corresponding descriptors of population diversity and differentiation. (DOCX 81 kb)
Individual multilocus genotypes for the nine microsatellite loci. Each allele is coded as its length, given in bp. Missing data is represented by “0”. (XLSX 66 kb)
Maxent models for quantifying environmental suitability, provides a description of methods and results related to ENM used to obtain habitat suitability maps for each species. (DOCX 727 kb)