Variation in body size and sexual size dimorphism in the most widely ranging lizard: testing the effects of reproductive mode and climate

Abstract Reproductive mode, ancestry, and climate are hypothesized to determine body size variation in reptiles but their effects have rarely been estimated simultaneously, especially at the intraspecific level. The common lizard (Zootoca vivipara) occupies almost the entire Northern Eurasia and includes viviparous and oviparous lineages, thus representing an excellent model for such studies. Using body length data for >10,000 individuals from 72 geographically distinct populations over the species' range, we analyzed how sex‐specific adult body size and sexual size dimorphism (SSD) is associated with reproductive mode, lineage identity, and several climatic variables. Variation in male size was low and poorly explained by our predictors. In contrast, female size and SSD varied considerably, demonstrating significant effects of reproductive mode and particularly seasonality. Populations of the western oviparous lineage (northern Spain, south‐western France) exhibited a smaller female size and less female‐biased SSD than those of the western viviparous (France to Eastern Europe) and the eastern viviparous (Eastern Europe to Far East) lineages; this pattern persisted even after controlling for climatic effects. The phenotypic response to seasonality was complex: across the lineages, as well as within the eastern viviparous lineage, female size and SSD increase with increasing seasonality, whereas the western viviparous lineage followed the opposing trends. Altogether, viviparous populations seem to follow a saw‐tooth geographic cline, which might reflect the nonmonotonic relationship of body size at maturity in females with the length of activity season. This relationship is predicted to arise in perennial ectotherms as a response to environmental constraints caused by seasonality of growth and reproduction. The SSD allometry followed the converse of Rensch's rule, a rare pattern for amniotes. Our results provide the first evidence of opposing body size—climate relationships in intraspecific units.

Latitudinal and altitudinal clines in body size are the most widely observed ecogeographic patterns (e.g., Ashton & Feldman, 2003;Blanckenhorn & Demont, 2004). Temperature is often assumed to be the principal determinant of ecogeographic body size clines, because temperature covaries consistently with latitude and altitude, and it strongly affects vital processes in the organisms (Angilletta, 2009). Yet, a number of recent studies have found that water availability (precipitation, humidity), and particularly seasonality (within-year variation in temperature or precipitation), often explain a higher proportion of body size variation than does mean temperature (e.g., Ashton, 2001;Çağlar, Karacaoğlu, Kuyucu, & Sağlam, 2014;Stillwell, Morse, & Fox, 2007). For each of these factors, multiple mechanistic hypotheses have been proposed (see below). However, rigorous testing of such hypotheses is often impeded by collinearity between climatic variables (Millien et al., 2006) which is particularly common within limited geographic areas. Specifically, colder environments are often associated with a greater seasonality (Aragón & Fitze, 2014;Chown & Klok, 2003;Körner, 2000). Furthermore, the pattern of the relationship between a phenotypic trait and a climatic covariate can be nonmonotonic, such as inverted U clines (Hjernquist et al., 2012) or saw-tooth patterns (Masaki, 1967;Mousseau, 1997). Yet, such more complex patterns are unlikely to be revealed within limited spatial and environmental ranges (Ashton & Feldman, 2003;Gaston, Chown, & Evans, 2008).
Wide-ranging species present promising models for studying phenotype-climate relationships, because the variation of target traits can be documented for numerous geographically distinct populations exhibiting a wide range of climates and diverse combinations of putative predictors (Roitberg et al., 2013(Roitberg et al., , 2015 cf. Meiri, Yom-Tov, & Geffen, 2007;Jetz, Ashton, & La Sorte, 2009).
The problem is that wide-ranging species often consist of several phylogeographic lineages. Pooling samples from different lineages may lead to a spurious trait-climate correlation (Romano & Ficetola, 2010) or obscure true relationships (Ashton, 2001). Therefore, even though body size is both phenotypically plastic and evolutionary malleable (Falconer, 1989;Green & Middleton, 2013;Jetz et al., 2009;Millien et al., 2006), the effects of current environment should be examined jointly with those of ancestry (Ashton, 2004;Diniz-Filho, 2008;Gaston et al., 2008). Yet, comprehensive range-wide studies of this kind have rarely been conducted on widespread species, even for fundamentally important traits such as body size (Angilletta, Niewiarowski, et al., 2004;Horváthová et al., 2013;Roitberg et al., 2013). Furthermore, although studies of intraspecific body size variation usually consider both sexes, they seldom explore how body size differences between males and females (i.e., sexual size dimorphism, SSD) vary along geographic or climatic gradients (Laiolo, Illera, & Obeso, F I G U R E 1 Geographic ranges of different clades of Zootoca vivipara, their phylogenetic relationships (after Surget-Groba et al., 2006), and our study sites. Details for study samples  are given in Appendix A7 (Table A1). Two relic lineages (central viviparous I and II), sporadically distributed in the south of Central Europe, are not shown. Note that the distribution border of the eastern oviparous clade (Z. v. carniolica) is actually jagged so that samples 25 and 26 may represent mixtures of different clades (W. Mayer, unpublished data; see also Cornetti et al., 2015;Lindtke et al., 2010) and are excluded from our analyses of climate on adult body size and sexual size dimorphism in this species. Our specific hypotheses and their predictions are summarized in Table 1 and presented in detail below. Other things being equal, we consider explanations based on plasticity as more parsimonious than hypotheses implying genetic adaptation (Chown & Klok, 2003;Madsen & Shine, 1993;Roitberg et al., 2013).
An alternative but not mutually exclusive hypothesis for Bergmann's clines is the temperature-size rule, a predominant pattern of developmental plasticity in ectotherms. This rule postulates that individuals growing at lower temperature mature later but at a larger size than conspecifics growing at higher temperature (Angilletta, 2009;Atkinson, 1994). (1) Dehydration resistance hypothesis a (2) Immediate negative effect of rainfall on lizard activity, food intake, and hence on body growth b

Body size increases with P2 (Prediction 4)
Delayed positive effect of rainfall on habitat quality, including food availability b Length of activity season Seasonality (here, Mean winter temperature, T1) Body size increases with T1 (converse pseudo-Bergmann's cline, Prediction 5) Adolph and Porter (1993) "null model" (age at maturity is constant) b Body size decreases with T1 (pseudo-Bergmann's cline, Prediction 6) (1) Adolph and Porter (1996) "main model" (modal age at maturity shifts abruptly as season length reaches a threshold) b (2) Starvation resistance hypothesis a Viviparity is associated with: (1) lower reproduction frequency; (2) higher gestation costs; (3) stronger maternal body-volume constraints on reproductive output Note: See text for details and references.
Abbreviation: SSD, sexual size dimorphism. a Hypotheses based on genetic adaptation.
b Hypotheses based on plasticity.
c Adolph & Porter's "main model" actually predicts a marked decrease of body size with T1 around the threshold values resulting in a saw-tooth cline whose overall linear trend is decreasing body size with T1. d Female size varies more than male size among populations.
e Male size varies more than female size among populations.

| Hypotheses related to water availability
Considering that large individuals have a reduced surface-to-volume ratio and overall higher absolute water content compared to small individuals, an adaptive dehydration resistance hypothesis predicts a negative correlation of body size with precipitation or humidity (Prediction 3) (see Stillwell et al., 2007 for references). Summer precipitation may also directly affect growth and body size in ectotherms; specifically in lizards and insects, there can be immediate, negative effects on insolation and consequently on animal's activity, food intake and thus body growth and positive, delayed effects on habitat productivity, which enhance foraging opportunity of lizards in later times (reviewed by Çağlar et al., 2014;Le Galliard et al., 2010). The negative effects correspond to our Prediction 3, while the positive ones to Prediction 4.

| Hypotheses related to seasonality
Among multiple hypotheses relating body size variation to the length of annual activity (or inactivity) the models of Porter (1993, 1996) seem particularly relevant for our study since they were developed specifically for lizards as perennial terrestrial ectotherms with advanced behavioral thermoregulation. The basic point of their reasoning is that annual growth increment in such organisms is mainly determined by the length of activity season rather than environmental temperature. Their "null physiological model", that is, the basic model assuming the absence of other factors, predicts smaller body size in more seasonal climates (Prediction 5; Adolph & Porter, 1993) which reduce energy acquisition opportunities. This pattern can be termed "converse pseudo-Bergmann's cline" to distinguish from the true Bergmann's and true converse Bergmann's clines which relate to environmental temperature.
The "main" model by Adolph and Porter (1996) additionally considers a discontinuous variation in the age at maturity (the age at the first reproduction), which is inherent to perennial ectotherms in seasonal climates. It predicts that this variation may reverse the body size cline expected by the null model. As the length of activity season decreases below some threshold, the modal group of subadult individuals cannot attain an appropriate body size within the reproductive season of the same year of life as their conspecifics in an environment allowing longer activity season. Under such constraints, the subadults invest available energy into further growth and start reproduction in the following season. This disproportional prolongation of juvenile growth may overcompensate the shortening of the annual activity period and enhance the typical size at maturity, at least within some range of climates. The predicted pattern is increasing adult body size in more seasonal climates ("pseudo-Bergmann's" cline; Prediction 6).
The main model predicts a saw-tooth cline (Adolph & Porter, 1996: Figure 4) whose overall linear trend is likely a pseudo-Bergmann cline as well. As reproduction is expected to more strongly inhibit body growth in females than in males, earlier maturation might be responsible for smaller female relative to male size in warmer or less seasonal climates (see Roitberg & Smirina, 2006 for indirect evidence in another lacertid lizard, Lacerta agilis). Thus, specifically the effect predicted by the main model may be female-biased.
Note that the underlying mechanism of Adolph & Porter's models is a direct response to environmental constraints which is not necessarily accompanied by genetic divergence.
Prediction 6 is also made by the adaptive fasting endurance, or starvation resistance hypothesis (e.g., Aragón & Fitze, 2014;Ficetola et al., 2010). Its version that applies to temperate zone reptiles explains pseudo-Bergmann clines via ability of larger-sized animals to acquire and carry larger fat reserves relative to metabolic needs than smaller-sized animals, this advantage being more important in climates with longer winters (Ashton, 2001;Litzgus et al., 2004).

| Hypotheses related to sex-specific selection or plasticity
Two distinct hypotheses related to sex-specific selection predict more female-biased SSD in colder and more seasonal climates. The extended fecundity-advantage hypothesis (reviewed by Cox, Skelly, & John-Alder, 2003; see also Angilletta, Steury, & Sears, 2004;Litzgus & Smith, 2010;Roitberg et al., 2015) suggests that reduced reproduction frequency should select for higher fecundity and thus for larger females (Prediction 7a).
The small male advantage hypothesis (reviewed by Blanckenhorn, 2000Blanckenhorn, , 2005Zamudio, 1998;Cox et al., 2003) argues that at low population densities the importance of male-male agonistic interactions, which select for larger body size, should decrease, while the disadvantage of lower mobility (which is often associated with large size) should increase. This hypothesis can be extended for a wider range of conditions. For instance in ectotherms, cold or highly seasonal climates, which reduce energy acquisition opportunities (Congdon, 1989), should also select against energetically costly aggressive behavior (and decrease the benefits of large male size) while increasing the small size-associated advantage of lower resource demands. Thus, we predict smaller male size in cold or highly seasonal climates (Prediction 7b).

| Study species
Zootoca vivipara is a small (adult snout-vent length 40-80 mm), ground-dwelling, insectivorous, heliothermic lizard. Compared to most other lizards Z. vivipara shows a high resistance to low temperatures and a low resistance to desiccation (Reichling, 1957). It prefers humid habitats, mostly in the forest vegetation zone (see Figure 1 and Roitberg et al., 2013 for further details and references).

| Body size data
We used the snout-vent length (SVL), the primary measure of body size in lizards and snakes (Roitberg et al., 2011 and references therein), as a proxy for overall structural size. We summarized original and published SVL data for 19,935 common lizards from 240 localities combined in 72 geographically distinct study samples (populations); these cover a major part of the species range ( Figure 1; Appendix A7: Table A1). The original data come from museum samples or from previous studies mostly performed for parasitological monitoring . Additional data were extracted from published histograms (e.g., Pilorge, 1987); in few cases, we also considered summary statistics for sex-specific adult SVL. When data for the same site were found in several studies, all unique samples for each sex were included to increase sample size. In total, 5,055 males and 6,474 females were considered as adults and constituted our study samples (see below).

| Data analysis
Within localities, samples from different years were pooled to increase sample sizes and to apply a standard approach across all data.
Whenever reasonable sample sizes were available we used strictly local samples, both for original and published data. When local sample sizes were too small, however, we pooled them into compound samples for larger geographic areas ( Figure 1) and used in our analyses weighted means for the study traits and nonweighted means for climatic variables. See Roitberg et al. (2013) for our criteria for pooling samples.
To test robustness of our analyses of the variation in adult body size and SSD to potential confounding factors (see Section 2.4 below) all main analyses were performed for two sets of data. Data set 1 included the animals assigned to adults by the primary researcher (Appendix A7: Table A1, A). Data set 2 was based on single inclusion criteria across samples: body length equal to or exceeding 45 mm for males and 48 mm for females (Table A1, B). These thresholds are close to typical minimum SVL of mature common lizards reported for most viviparous (Avery, 1975;Cavin, 1993;Orlova, 1975;Pilorge & Xavier, 1981) and some oviparous (Sinervo et al., 2007) populations studied. Furthermore, each main analysis was run for mean values and additionally for the 80th percentiles of the size distributions (see Roitberg, 2007 for details and references on higher percentiles as useful estimators of population's typical adult body size in indeterminate growers; see also Case, 1976). Thus, we used four metrics for each of our target traits, that is, male size, female size, and SSD.
Means and percentiles of sex-specific SVL were ln-transformed for all analyses except SSD.
Sexual size dimorphism was quantified with the index: SDI = (size of larger sex/size of smaller sex) − 1, conventionally expressed as positive if females are larger and negative if males are larger (Lovich & Gibbons, 1992). This index shows several favorable properties (Lovich & Gibbons, 1992;Smith, 1999) and has become a standard metric for studies on sexual dimorphism (Fairbairn, Blanckenhorn, & Szekely, 2007). SSD allometry was quantified with the slope of major axis regression (model II) of log(male SVL) on log(female SVL) (Fairbairn, 1997). The slopes (β) and their 95% confidence intervals were computed with the lmodel2 package (Legendre, 2013) in R (R Core Team, 2018). They were tested against the null hypothesis of β = 1 (isometry). The pattern with β > 1 is referred to as Rensch's rule, and that with β < 1 as converse Rensch's allometry (Table 1).
The following bioclimatic indices were used as explanatory vari- To simultaneously analyze categorical (clade identity, hereafter Clade) and continuous (climatic variables) effects, as well as their interactions, on the variation among population means (or percentiles) of a target trait, we used general linear models (GLMs). We determined the best combination of predictors of each target trait using an information-theoretic approach (Burnham & Anderson, 2002) based on the Akaike's Information Criterion corrected for small sample sizes (AICc). Prior to model fitting, all continuous input variables were standardized to mean = 0 and SD = 1 to improve the interpretability of main effects in the presence of significant interactions (Schielzeth, 2010). We then fitted models encompassing all possible combinations of input variables and their first-order interactions, including an intercept-only model, calculating for each combination the AICc score. The interactions were included for explorative purposes.
Models with ΔAICc ≤ 2 were considered candidate models and used for further analysis. As the Akaike's criterion may select overly complex models, we considered a complex model as a candidate model only if its AICc was lower than AICc of all simpler models nested in the complex model (Richards, Whittingham, & Stephens, 2011). Model selection was performed in R version 3.4.3 using the "MuMIn" package (Bartoń, 2017). Considering collinearity between Clade and T1 (Appendix A6) we refrained from model averaging, as recommended by Freckleton (2011). Instead, we summarized our candidate models verbally. We also performed some additional analyses exploring the effect of clade identity on "climate-corrected" body size.
To evaluate the effect of reproductive mode the western oviparous clade (the only oviparous clade in our data set) was used as the reference level of the factor Clade. Furthermore, all GLM analyses were repeated for two viviparous clades only. Comparing the effects of Clade in the two data sets (three clades vs. two viviparous clades) provided additional evaluations of the effect of reproductive mode.
Values of all response and explanatory variables for 72 study samples are provided in Appendix A7 (Tables A1 and A2).

| Methodological caveats
In species with continuing growth after maturity, numerous factors unrelated to geographic variation, such as local and temporal fluctuations in the abiotic (e.g., temperature and humidity) and/ or biotic (e.g., food resources) environment, can affect patterns of growth, maturation, and survival of different cohorts and thus body size distribution in a particular study sample (Kratochvíl & Frynta, 2002;Roitberg, 2007;Shine, 1994Shine, , 2005Stamps, 1993;Watkins, 1996). Further biases can come from compiling data of several independent researchers. They may differ in measuring routine, type of material (living vs. freshly euthanized vs. preserved specimens), and in their criteria of separating adults from immature animals, that is, inclusion criteria (these can be based e.g., on body size and color pattern vs. the state of gonads). The biases from the first two factors are expected to be within a few percents (Case, 1976;Roitberg et al., 2011; E. S. Roitberg, unpublished data; Vervust, Van Dongen, & Van Damme, 2009), and this is much lower than the observed variation within and among our study samples.
Indeed, the factor Type of material was never significant when included in our best models as additional predictor. We also examined potential bias from temporal trends in adult body size (e.g., Chamaillé-Jammes et al., 2006;Green & Middleton, 2013) by adding the factor Time (1950-1990 vs. 1991-2000 vs. 2001-2015); this addition did not improve our best models. The effects of inclusion criteria, and those of temporal variation in the proportion of newly matured animals (Watkins, 1996), were accounted for by using four different metrics for adult body size (see Section 2.3).

| Candidate models for male size
Geographic variation in male SVL was weak (range of sample means 7 mm, min-max 48-55 mm; Appendix A7: Table A1) and poorly explained with our predictors. Even the AICc best-fit models (hereafter top models) explain only 2%-6% of the total variance (Table 2), and in most analyses they do not perform considerably better than the intercept-only model (ΔAICc ≤ 2). Only metric 4 explains 17% (Table 2, model 29) and only in one data set.

| Candidate models for female size
Compared to males, body size variation in females is clearly larger (range of sample means is circa 16 mm, min-max 51-67 mm, Table   A1), with much greater part of this variation being explained by our predictors (26%-49%, Table 3). In the three-clade analyses, all five candidate models include Clade and P2 (summer precipitation); three models include T1 (winter temperature), two of them also including the T1 × P2 interaction (Table 3). In the two-clade analyses, the top models explain a smaller proportion of the total variance than in the three-clade analyses (26%-36% vs. 40%-49%, Table 3); they are less consistent among different metrics, and the total number of candidate models is larger (13 vs. 5). Only around half of the candidate models include Clade and P2 (Table 3). T1, which is present in 11 of 13 candidate models in the two-clade analyses, "replaces" Clade (its collinear counterpart, see Section

2.3) as the most consistent predictor. Model coefficients of T1
were always negative, that is, female size overall increases with decreasing T1. As in the three-clade analyses, a considerable number of candidate models (5 of 13, three of them being the top models, Table 3) include the P2 × T1 interaction. Note that the latter two predictors are also present in the only significant model for male size (Table 3). Model coefficients of P2 were consistently negative (see also Figure 2c), thus supporting our Prediction 3 (Table 1). The P2 × T1 interaction indicates that the major relationship, female size-T1 (Figure 2a), is modulated by covariate P2: in the eastern viviparous clade, this relationship is stronger at higher than at lower values of P2 (Appendix A8: Figure A2).

| Candidate models for SSD and the shared patterns of female size and SSD variation
Sexual size dimorphism was consistently female-biased: SSD index for means varied from 0.04 to 0.27 (Appendix A7: Table A1); that is, females were on average 4%-27% longer in SVL than males. In both the three-clade and the two-clade data sets, the major axis regression slope of log(male SVL) on log(female SVL) was significantly lower than 1 (Table 4). This pattern corresponds to a converse of Rensch's rule; that is, the SSD variation is primarily due to variation in female rather than male size.
The variation in SSD is even better explained by our predictors (up to 58%, Table 5) than that of absolute female size. Model selection shows a high consistency in both the three-and the twoclade data sets. None of the top models for SSD includes P2, an important predictor of the female size models (see above). Another discordance with the female size models is a consistent presence of the T1 × T2 interaction, a predictor infrequently occurring in the models for female size. This interaction indicates that in the eastern viviparous clade, as well as for the whole data set, the negative SSD-T1 relationship is stronger at lower than at higher values of T2 (Appendix A8: Figure A3).
The following patterns are common to the female size and the SSD variation. First, candidate models for both traits frequently include the Clade × T1 interaction (Table 3, Table 5) Table A3). The negative female size-T1 and SSD-T1 correlations for the whole data set are highly significant too (Table A3). The negative correlation corresponds to a pseudo-Bergmann's cline (Table 1: Prediction 6), while the positive one corresponds to its converse (Table 1: Prediction 5). Second, the main effect of T1 is much stronger than that of T2: (a) while ten models include T1 but not T2, no models show the opposite pattern (Table 3, Table 5); (b) in no model, T2 exhibits a significant main effect. Thus, no support for our Predictions 1 or 2 (Table 1) was found. Third, as with female size, clade identity is a consistent predictor of SSD in the three-clade analyses (it occurs in all six candidate models, Table 5) and becomes inconsistent in the two-clade analyses, where four of the six candidate models include climatic predictors only (Table 5); the top models explain a smaller proportion of the total variance in the two-clade analyses than in the three-clade analyses (26%-47% vs. 41%-58%, Table 5).
To further explore this issue, we visualized between clade differences in female size and SSD as they appear with and without controlling for climatic covariates (Figure 3). In line with Prediction 8 (Table 1), the western oviparous clade exhibits a smaller female size and less female-biased SSD than the western and the eastern viviparous clades, whereas differences between the two viviparous clades become insignificant when controlling for climatic variables.

| D ISCUSS I ON
We investigated intraspecific divergence of adult body length in a wide-ranging lizard across the temperate Eurasia. Using general linear models, we tested the effects of mean summer temperature, summer precipitation, seasonality, and reproductive mode/lineage identity on female size, male size, and SSD. We found a moderate effect of reproductive mode and precipitation and a strong but complex effect of seasonality. The latter differed drastically between the lineages, being also modulated by precipitation and especially by temperature. Female size and SSD varied stronger than male size, and virtually all the effects were strongly female-biased. Below, we relate the revealed body size patterns to several evolutionary and ecological hypotheses (Table 1).

TA B L E 3
AICc-selected models (ΔAICc ≤ 2) for female size (SVL) F I G U R E 2 (a-c) Female body size (snout-vent length, SVL) and (d-f) sexual size dimorphism (SSD) in Zootoca vivipara plotted against three climatic predictors. Regression lines are shown whenever the slopes differ from zero at p < .1

| Allometry of SSD
Amniotes, including reptiles, tend to exhibit standard Rensch's allometry meaning that male size varies more than female size (Fairbairn, 1997). This trend is mainly expressed among species (reviewed in Fairbairn et al., 2007) but was also reported for variation among conspecific populations (e.g., Saino & De Bernardi, 1994;Garel, Solberg, Saether, Herfindal, & Høgda, 2006;Aglar & López-Darias, 2016), even in species with overall female-biased SSD (Pearson et al., 2002). In contrast, SSD variation in Z. vivipara follows a converse of Rensch's rule. In Z. vivipara, male size varies not only lower but also qualitatively less regular in relation to our predictors, as compared to female size. This pattern is not solely due to a divergence between oviparous and viviparous populations (see below), as it persists in our two-clade analyses (the western + eastern viviparous clades). Obviously, the revealed SSD allometry is also shaped by a steeper slope of the body size-climate relationship in females versus males, a pattern which is opposite to a prevailing trend found in a meta-analysis of 98 animal species (Blanckenhorn et al., 2006).

| Temperature and water availability during the active season
A strong sexual difference in the extent of body size variation, and in the percentage of this variation explained by our predictors, reduces the relevance of the hypotheses that apply equally to both sexes. This concerns the adaptive hypotheses related to heat acquisition, heat conservation, dehydration resistance, and fasting endurance (Table 1). In accordance with this reasoning, Predictions 1 and 2, which are associated with heat acquisition and heat conservation mechanisms (Table 1), received no support in our analyses: The main effect of the corresponding predictor, mean summer temperature (T2), was consistently weak. The temperature-size rule, which also makes Prediction 2 ( Note: SSD metrics D1, D2, D3, and D4, used as response variables, are based on the corresponding metrics for sex-specific SVL defined in Tables 2  and 3. Other designations as in Table 2.  (Table 1). Previous studies on lizards, including Z. vivipara, made on a small geographic scale in warmer regions (Díaz, Iraeta, Verdú-Ricoy, Siliceo, & Salvador, 2012;Dunham, 1978;Lorenzon, Clobert, & Massot, 2001;Taylor, 2003) reported a positive effect of precipitation or humidity on body size (Prediction 4). The opposite pattern revealed at a large geographic scale (Figure 2c) could mean that in cooler climates, which occur on a major part of the study species range, the negative, immediate effect of precipitation on animal's activity (Table 1) exceeds the positive, delayed effect on habitat productivity (Table 1).

| Seasonality
The strongest and most interesting pattern revealed in this study concerns the relationship of female size and SSD with mean winter temperature (T1) which is tightly correlated to seasonality and used as proxy for the length of activity season. The western viviparous clade exhibits a converse pseudo-Bergmann's cline corresponding to the Adolph and Porter (1993) "null physiological model" (Table 1: Prediction 5). In contrast, a standard pseudo-Bergmann's cline found in the eastern viviparous clade, as well as across the clades, corresponds to their main model (Table 1: Prediction 6) that considers shifts in the age at maturity (Adolph & Porter, 1996). Phenotypic responses predicted by Porter (1993, 1996) models can be female-biased due to sex-differential plasticity of growth and maturation (Fairbairn, 2005;Cox & John-Alder, 2007;Cox & Calsbeek, 2010; see also Table 1).
We hypothesize that geographic variation in age at maturity, specifically in females, is minor in the western viviparous clade, while pronounced in the eastern viviparous clade, at least in our data set. Available data on typical age at maturity in different populations of Z. vivipara appear to be in line with this hypothesis. In virtually all western viviparous populations, which were studied demographically, the modal age at maturity is 2 years (Bauwens & Verheyen, 1987;Heulin, 1985;Pilorge, 1987;Pilorge & Xavier, 1981;Strijbosch & Creemers, 1988;S. Hofmann, unpublished data). Demographic data are scarce for the eastern viviparous clade, but the modal age at maturity is expected to be generally higher and apparently more variable than in the western viviparous clade, because its geographic distribution (Figure 1) results in a strongly higher range of experienced seasonality ( Figure 2a).
Indeed, in the south of West Siberia the typical age at maturity was found to range from 2 to at least 3 years (Bulakhova, Kuranova, & Savelyev, 2007;Epova, Kuranova, Yartsev, & Absalyamova, 2016), being probably even higher in still more severe climates of the Middle and East Siberia. The above considerations allow us to suggest that the opposite correlations of female size (and SSD) with mean winter temperature found in the two major clades of Z. vivipara may reflect a saw-tooth cline along a seasonality gradient predicted by the Adolph and Porter (1996) model (Figure 4).
A stronger effect of seasonality (as estimated by winter temperature) on female body size and particularly SSD in cooler versus warmer summer climates (T1 × T2 interaction; Tables 3 and 5; see also Appendix A8: Figure A3) is beyond our predictions. Perhaps at lower ambient temperatures, lizards use more of the potential activity period predicted by the Adolph and Porter (1993) "null" model to compensate for process limitations (sensu Congdon, 1989). In contrast, in warmer summer climates, lizards might often reduce their activity to avoid predation risk (Werner & Anholt, 1993; see also Sears, 2005), as well as the risk of overheating and dehydration, which are generally higher in such environments.
The fasting endurance hypothesis (larger individuals possess larger fat reserves and survive hibernation better than smaller individuals - Table 1), which also makes Prediction 6, cannot explain the opposing cline in the western viviparous clade and can hardly integrate a strongly female-biased phenotypic response.
The extended "small male advantage hypothesis", which implies a Rensch's allometry of SSD (Prediction 7b), is strongly rejected in the present study that found the opposite, converse Rensch pattern (Table 4). This converse Rensch's allometry (Prediction 7a), as well as the pseudo-Bergmann cline in female size and SSD, corresponds well with the extended fecundity-advantage hypothesis viewing this cline as an adaptive compensation of reduced reproduction frequency

F I G U R E 4
The saw-tooth relationship between population's typical adult female size and seasonality in the lizard Zootoca vivipara, as hypothesized from Porter's (1993, 1996) models. Thin line segments correspond to constant ages at the first reproduction, where the body size-seasonality relationship follows Adolph and Porter's (1993) null model. Thick segments indicate thresholds at which the age at the first reproduction changes abruptly resulting in a reversed body size-seasonality relationship (Adolph & Porter, 1996). See text for explanations (Table 1). However, females of our study species produce a single litter per year in a wide range of climates: repeated clutches occur in considerable frequencies only in lowland oviparous populations, which constitute a tiny portion of our study samples (2 from 72); only exceptional cases of multiple clutches per season are known for viviparous populations (Horváthová et al., 2013). At the same time, no evidence for biennial or intermittent breeding so far exists for Z. vivipara. Furthermore, previous work (Roitberg et al., 2013) found no significant relationship of clutch size or mass with seasonality, the relationship of these reproductive traits with summer temperature being negative (Roitberg et al., 2013 applied climatic variables tightly correlated to those used in the present study; see Section 2).
The above points impair a straightforward application of extended fecundity-advantage hypothesis to the body size patterns presented here (this issue will be addressed in our future work). Also, as most other presented hypotheses, the fecundity-advantage hypothesis does not explain the converse pseudo-Bergmann's cline in the western viviparous clade.
The Porter (1993, 1996) (Adolph & Porter, 1993). A further warming can reverse this trend when a major part of yearling females would reproduce (thereby impeding their further growth) because they reach the threshold size within the "reproduction window" of their second calendar year (Adolph & Porter, 1996).
Another important point can be inferred from the revealed patterns of climatic variation. Ashton and Feldman (2003) based their comprehensive meta-analysis of ecogeographic body size clines in reptiles on mean annual temperature. In our study, mean annual temperature is rather strongly correlated with seasonality and mean winter temperature and only weakly correlated with mean summer temperature (Appendix A3), which might mean that many of the clines reported by Ashton and Feldman (2003)  were consistently positive in warm-water species, while consistently negative in cool-/cold-water species (Rypel, 2014).

| Reproductive mode and lineage identity
In be it the western viviparous clade (Lindtke, Mayer, & Böhme, 2010) or a relic clade the central viviparous II (Recknagel et al., 2018). Yet the differences between oviparous and viviparous populations are not clear-cut and explain a moderate part of the intraspecific variation in female size and SSD in Z. vivipara (Figure 2). This is a likely reason why previous range-wide studies on Z. vivipara (Horváthová et al., 2013;Roitberg et al., 2013), which had less representative data sets, found no significant effect of reproductive mode on female body size.
A problematic point of this study is evaluating the main effect of clade identity (Clade) on the body size variation among viviparous populations. Clade is collinear with T1 (Appendix A6), our proxy for seasonality, which is the strongest predictor. This collinearity arises due to a very low overlap between the western and eastern viviparous clades along the seasonality gradient (Figure 2a,d) resulting from the west-east separation of their geographic distributions (Figure 1).
Further complications come from a significant Clade × T1 interaction (Tables 3 and 5). Under such conditions, the absence of clade identity in some candidate models (Tables 3 and 5

| CON CLUS ION
The present study confirms a major role of seasonality (as compared to mean summer temperature and precipitation) in shaping the geographic body size variation as suggested by previous work in Z. vivipara (Horváthová et al., 2013;Roitberg et al., 2013). We show that the body size response to the seasonality gradient is strongly female-biased, more complex than previously thought, and is parsimoniously interpretable as a saw-tooth cline along a gradient of activity season lengths predicted by Adolph & Porter's models (Table 1; Figure 4). Within the western viviparous clade occurring from France to Eastern Europe, female size and SSD decrease as seasonality increases. Such a response is predicted under constant age at maturity (Adolph & Porter, 1993), an overall likely condition for this region, judging from its low seasonality and from available demographic data. Within the eastern viviparous clade (Eastern Europe to Far East), as well as across the clades, female size and SSD increase with increasing seasonality.
This response is predicted under varying age at maturity (Adolph & Porter, 1996), and such variation, being confirmed by scarce empirical evidence, seems very likely, given the high and strongly Further, our study demonstrates that not only males (e.g., Aglar & López-Darias, 2016) but also females can be a driver of intraspecific divergence in amniotes.

ACK N OWLED G EM ENTS
This paper is dedicated to Werner Mayer (1943-2015) whose studies profoundly contributed to the knowledge of genetic differentiation within our study species and to our deceased coauthor Olga Leontyeva . We are grateful to Werner Mayer for advice on lineage composition of Central European samples, Andrey Reshetnikov for assistance with creating Figure 1, and Julia Kavalerchik for developing numerous scripts in perl aiding at processing of our voluminous data. Vladimir Shitikov kindly helped us with AIC-IT analyses in R; these analyses have benefited from comments by Holger Schielzeth.
We thank Astrid Clasen, Michael Fokt, Regina Shamgunova, Igor Tarasov, Vladimir Yakovlev, and Igor Doronin for sharing their SVL data. We also thank the researchers whose published data were com-

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
ER designed the study. All authors provided morphometric and other data from their fieldwork or museum samples; SH provided climatic data. ER conducted analyses and wrote the manuscript, to which all coauthors contributed and gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
Original individual-based data on body size (SVL), as well as the geographic and climatic data, are archived into the public repository

A PPE N D I X A1 E X TR AC TI O N A N D CO R R EC TI O N O F CLI M ATE DATA
Climatic data (monthly mean minimum and maximum tempera-  (Sears, 2005;Angilletta, Oufiero, & Leaché, 2006;Stillwell et al., 2007;Roitberg et al., 2013) would adjust the temperatures but not the precipitation. To more effectively reduce the resulting deviation in climatic values we used an original approach briefly described below. We generated an array of coordinate values (subsidiary points) around the problematic site. If the discrepancy in altitude is merely due to a small random deviation in coordinate values, as is apparently true for the vast majority of cases, then some of the subsidiary points should be closer to the "true" site than the site with the reported coordinates. Choice of the most relevant point was based on the least difference in the altitude values, the geographic proximity being also accounted for.
For virtually all problematic sites of this study, at least one subsidiary point in close proximity provided a sufficiently close altitude value.

A PPE N D I X A 2 A N A LYS I S O F CLI M ATE DATA
Seasonality was estimated with the coefficient of variation (CV) among monthly values of temperature or precipitations (Dobson & Wigginton, 1996;Ashton, 2001;Çağlar et al., 2014). Temperature data were transformed to Kelvin to allow calculation of the variation coefficient (Ashton, 2001).
The principal component analysis (PCA) was used to summarize the geographic variation of 36 primary climatic variables (monthly mean minimum and maximum temperatures, and monthly mean precipitations): these were reduced to a smaller set of orthogonal vectors, which include a major portion of the total variation. In our previous studies (Roitberg et al., 2013(Roitberg et al., , 2015 the first two principal components, PC1-clim and PC2-clim, were directly used as predictors for body size and SSD in statistical models. In this study, we took advantage of the fact that PC1-clim is tightly correlated to

VA R I A B LE S ACROSS TH E S TU DY S ITE S
The first axis of the principal component analysis of the climatic variables (PC1-clim) explained 55.2% of the total variance among localities ( Figure A1). PC1-clim is strongly and positively correlated with all monthly temperature and precipitation parameters outside the warmest quarter ( Figure A1); PC1-clim is highly correlated with T1 (Spearman rank correlation coefficient, r s = .975, p < .001; N = 72) but not T2 (r s = −.016, p = .893). PC1-clim is also tightly linked to temperature seasonality (r s = −.959, p < .001) and less strongly to precipitation seasonality (r s = −.819, p < .001). As expected, T1 exhibits similar correlations to these seasonality indices (−0.962 and −0.807, respectively). In contrast, the second principal component (PC2-clim, 28.3% of the total variance, Figure A1) is strongly correlated with the monthly values of the warmer season (April-September), the loadings being consistently positive for temperatures and negative for precipitation ( Figure A1). PC2-clim is tightly related to T2 (r s = .959, p < .001) but not T1 (r s = .122, p = .307). Both PC2-clim and T2 show no significant relations to temperature and precipitation seasonality (all r s < .21, all p > .09).

CO N S I D E R I N G TH E E FFEC T S O F E VO LUTI O N A RY LI N E AG E
The phylogeographic study by Surget-Groba et al. (2001 provides reasonably dense covering of border areas between the major clades. Thereby, virtually all our study samples could be

A PPE N D I X A 5 CO N S I D E R I N G TH E E FFEC T S O F S PATI A L AUTO CO R R EL ATI O N S
To test whether the results of our GLMs are biased by spatial autocorrelation we computed spline correlograms based on Moran's I statistic; for this we used the 'ncf' package (Bjornstad, 2013) in R (R Core Team, 2018). Correlograms were computed for the residuals of the two AICc best-fit models of each metric of each study trait. No significant autocorrelation was detected (95% confidence intervals always included zero), suggesting that spatial autocorrelation did not bias our models (Dormann et al., 2007).

Multicollinearity between our input variables was tested by
computing Generalized Variance Inflation Factor, GVIF (Fox & Monette, 1992) with the car package (Fox & Weisberg, 2011) in R (R Core Team, 2018). In the three-clade analyses, GVIF amounted 3.98 for clade identity (df = 2), 3.39 for T1, 1.39 for T2, and 1.08 for P2. Despite their higher GVIF values, both T1 and clade identity were retained in our models as these variables have independent biological relevance. Yet, this collinearity was appreciated when processing the results of model selection (see Section 2: Data analysis).

TA B L E
.175 .304 .071 .089 All three clades .580 .637 .029 .026 .019 .021 .019 .059 .017 .010 For correlations of female size and SSD with T1, significance of the differences between the two viviparous clades is also indicated. All designations of our response and explanatory variables as in Tables  1-3.

S E A SO N A LIT Y R EL ATI O N S H I P S I N CO M M O N LIZ A R DS
Body size at maturity is usually tightly related to overall adult size (Shine, 1990;Adolph & Porter, 1996;Stamps, Mangel, & Phillips, 1998). Yet in some perennial ectotherms, geographic patterns or sexual differences in body size at maturity can be strongly modified at later stages via postmaturity growth and survival (Howard, 1981;Rutherford, 2004;Tracy, 1999). In colder and more seasonal climates, perennial ectotherms tend to grow more slowly but live longer than in milder climates. The form of the adult size-seasonality relationship would apparently be determined by the relative strength of these two opposite trends. Obviously, a pseudo-Bergmann cline more likely arises when the increase in adult survival is "disproportional", i.e. when the life-time activity actually increases in seasonal climates instead of being just distributed over more but shorter seasons. Such an increase may occur if the survival increment caused by reduced lizards' potential annual activity (Adolph & Porter, 1993) is added with an extra increment caused by reduced predation and parasite pressure. A corresponding increase in SSD with increasing longevity is expectable if the post-maturity growth is less decelerated in the larger sex, especially when the larger sex also exhibits a higher adult survival (Watkins, 1996). Both patterns were reported for Zootoca vivipara (growth: Epova et al., 2016;survival: Strijbosch & Creemers, 1988;Uller, Massot, Richard, Lecomte, & Clobert, 2004;Le Galliard et al., 2010). The presented mechanism may strengthen the overall pseudo-Bergmann trend predicted by the Adolph & Porter (1996) model and explain the increase in female size and SSD with increasing seasonality revealed within the eastern viviparous clade and across the studied clades of Z. vivipara (Figure 2a,d). However, this mechanism cannot explain the opposing cline in the western viviparous clade.
The intriguingly disparate body size-climate relationships in the western and the eastern viviparous clades are difficult to explain outside the Adolph & Porter (1996) model. One alternative scenario would be that the western and eastern viviparous clades differ in their inherited thermal reaction norms (Angilletta, 2009) or potential for response to selection, thus being genetically predisposed to develop opposing body size clines along similar climatic gradients.
Currently, such a scenario, implying a substantive role of intrinsic factors, seems unlikely considering that the two clades are closely related  and show no appreciable differences in life-history (Roitberg et al., 2013) and external morphology .

F I G U R E A 2
The relationship between female body size and winter temperature (T1, our proxy of seasonality) in Zootoca vivipara is modulated by summer precipitation (P2) resulting in a significant T1×P2 interaction (Table 3). Fit line on the right panel indicates a significant relationship (p < .01), the corresponding line on the left panel is given for comparison

F I G U R E A 3
The relationship between sexual size dimorphism (SSD) and winter temperature (T1, our proxy of seasonality) in Zootoca vivipara is modulated by summer temperature (T2) resulting in a significant T1 × T2 interaction (Table 5). Fit lines for the western oviparous clade are not shown because of small sample sizes. Note that the opposite pattern of the SSD-Seasonality relationship in the western and eastern viviparous clades persists in both subsets of populations differing in summer temperature