Microsatellites were previously hypothesized to be selectively neutral (Charlesworth et al., 1994). Under this hypothesis, we would expect heterozygosity, but not mean allele length, to be positively correlated with the effective population size N_e (Kimura and Crow, 1964).

Empirical evidence (Amos et al., 1996, 2008) suggests that “heterozygous sites show increased mutability due to an extra round of DNA replication when they are recognized and “repaired” in heteroduplex DNA formed during synapsis” (Amos, 2016). This hypothesis suggests that differences in heterozygosity, due for example to population size or demographic processes, will lead to higher relative mutation rates (Amos et al., 2008). Therefore, compared with small populations, in larger populations, where heterozygosity is higher, mean allele length will be higher if loci are expanding (i.e., increasing in number of alleles) and lower if loci are contracting (i.e., decreasing in number of alleles). Because microsatellite length is finite, whereas many loci may have room to expand in length, many others can only contract. Thus, under this hypothesis, we expect a strong correlation between effective population size (or observed heterozygosity) and mean allele length, with the sign of the correlation coefficient being either positive or negative (Amos et al., 2008; Amos, 2010).

Many empirical studies have shown that environmental stress induces a state of hypermutation in bacteria (Harris et al., 1997; Galhardo et al., 2007), as well as in eukaryotes (Heidenreich, 2007; Bristow and Hill, 2008; Sharp and Agrawal, 2012). There are several possible mechanisms to explain the occurrence and benefits of hypermutation under stressful environmental conditions (Forche et al., 2011; Rosenberg et al., 2012; Torres-Barcelo et al., 2013). For example, if maladapted individuals increase their mutation rate, then the population adapts faster and the overall population mean fitness increases (Bjedov et al., 2003; Ram and Hadany, 2012; Ram et al., 2018). Thus, stress-induced mutagenesis implies that more mutations are expected to accumulate in a population during extreme environmental events, regardless of population size and heterozygosity level. Under this hypothesis, we therefore expect allele diversity to positively correlate with environmental stress, while allele length can either increase or decrease with stress.

Microsatellite allele length may be directly selected by biological, environmental, or climatic variables (Vigouroux et al., 2003; Amos and Clarke, 2008; Amos et al., 2008). Allele length may affect fitness by regulating the expression and structure of nearby genes (e.g., Schmidt and Anderson, 2006; Gao et al., 2013). One example of this is that of the decrease in allele lengths of maize with increase in altitude across a wide geographic range (Vigouroux et al., 2003). Under this hypothesis, we expect allele lengths of some loci to correlate with an environmental gradient or with environmental stability, independent of genetic diversity or heterozygosity. The sign of the correlation coefficient may be positive or negative, depending on the locus.

For this study, we collected P. syriacus tadpole samples from 22 vernal pools located along a steep climatic gradient in Israel (Figure 1). We used eight published microsatellite markers to estimate mean allele lengths, allelic richness, Shannon’s diversity, and observed heterozygosity for each sampled pool. To explain the observed variation in allele length along the climatic gradient in Israel and to test for a fit for each of the hypotheses presented above, we used dbRDA and mixed models.

The climatic processes affecting mtDNA variation along a steep gradient in Israel have been studied by Munwes et al. (2010). Here, we used a similar set of vernal pools, with the eastern spadefoot toad (Pelobates syriacus) as our model system. This species is a fossorial toad, spending the hot and dry periods in deep burrows in the soil that it digs with its hind legs. It occupies terrestrial habitats in open uncultivated lands such as sparse forests, steppe, semi-desert, and rocky areas. Spawning sites in this species comprise stagnant temporary water bodies, as well as large permanent pools. In Israel, the species is known to breed only in ephemeral pools and the adults are philopatric, breeding in the same vernal pool in consecutive years (Gafny, 1986).

The distribution range of P. syriacus extends to southern Romania in the north, Greece in the west, Iran in the east, and Israel in the south (Agasyan et al., 2009). In Israel, it is found from the Golan Heights and Galilee in the north to throughout the central coastal plain, which includes its present southernmost breeding site—Robert’s pool near Ashdod (Munwes et al., 2010). Although Israel is at the southernmost edge of its distribution, there is a steep gradient in environmental conditions within Israel (Figure 1). The northern regions of the country are relatively colder and wetter, with a more stable climate. The vernal pools in the north are filled every year and retain water for long periods (7–10 months) with relatively low variability among years. In contrast, the coastal plain is more arid, with a higher inter-annual variability in precipitation and its vernal pools consequently retain water for shorter periods (up to 5 months) and in some years do not fill at all (Munwes et al., 2010). Spadefoot toad species display extreme variation in larval period duration, due in part to the evolution of thyroid hormone (TH) physiology (Hollar et al., 2011). Eastern spadefoot toad tadpoles exhibit plasticity in the timing of metamorphosis (Székely et al., 2010). Furthermore, environmental factors have been shown to affect developmental traits in the tadpoles, conducing to rapid metamorphosis in drying ephemeral ponds (Storz and Travis, 2007).

The northern (i.e., core-like) and southern (i.e., edge; Munwes et al., 2010) populations of the spadefoot toad in Israel differ in several life-history traits (Gafny, 1986). Breeding season in the coastal plain is from early January to late February, while in the Upper Galilee it is from mid-March to mid-April. Females lay a long string of eggs (mean ± SD: 8,954 ± 5,229 eggs, n = 19; Gafny, 1986), which develop to metamorphosis in about 80 days in the coastal plain and 180 days in the Upper Galilee and Golan Heights. On average, in the coastal plain only about 30% of pool tadpole populations successfully metamorphose, and of those about 45% die later from dehydration. The main reason for the low survival rates in the coastal plain is the early desiccation of the vernal pools, before the tadpoles can achieve metamorphosis or find appropriate shelter. In contrast, at least 70% of the Galilee and Golan Heights tadpole populations successfully metamorphose (Gafny, 1986). Additionally, post-metamorphic individuals are also smaller in the coastal plain (1.3–4.0 g) relative to those at the same developmental stage in the north of Israel (4.5–9.3 g; Gafny, 1986). Body size and stored fat at metamorphosis are thought to be important correlates of fitness in amphibians (Morey and Reznick, 2000). Post-metamorphic individuals in the coastal plain must tolerate the dry and unfavorable conditions for much longer than in the northern regions. The coastal plain, which is the southernmost edge of the species’ range, is thus also ecologically marginal (Munwes et al., 2010).

The study was conducted in 22 vernal pools known as breeding sites for Pelobates syriacus in Israel (Table 1, Supplementary Table 1, and Figure 1). The sampling sites spanned northern and central Israel, including the southernmost known breeding site for this species. Most DNA samples were collected over a period of 4 years (79%; 2012–2015). We also incorporated DNA samples collected between 2006 and 2010 (21%) by Munwes et al. (2010).

We visited each vernal pool 2–3 times annually, during winter and spring (Oct-May), and each time sampled 15–20 tadpoles. The tadpoles were measured (TL—snout to end of tail). Newly hatched tadpoles (≤3 cm) were collected whole and tail-tip samples were taken from larger tadpoles. All samples were stored in 96% alcohol and kept frozen until analysis. To maximize the chance of sampling tadpoles from different spawns (i.e., individuals that are not full siblings), we spatially sampled from around each pool. Most of the pools were small (diameter < 50,m) and shallow (< 1,m), enabling extensive sampling. To further verify that our samples were not composed of highly related individuals, we conducted a relatedness analysis using ML-Relate (Kalinowski et al., 2006). We identified full sibs using the 95% confidence set for R = 0.5 approach. We also mapped out the matrix of relatedness between all tadpoles sampled in each region using PROXCAL multidimensional scaling (Commandeur and Heiser, 1993).

We used eight published microsatellite loci (Cohen et al., 2021; Supplementary Table 2), a set of genetic markers that is particularly suitable for individual-level analysis (Wan et al., 2004). Genomic DNA samples from P. syriacus were sent to the Evolutionary Genetics Core Facility (EGCF) at Cornell University for genomic library construction via Illumina MiSeq sequencing. Five samples from pools across Israel were used to construct the microsite library. The EGCF provided a msatcommander (Faircloth, 2008) output file with thousands of potential microsatellite primer pairs. We identified polymorphic loci even before the PCR testing step by sorting the output data file. We looked for loci with tetramer motifs, size range of 150–400 bp, annealing temperature of 58–60°C, and 7–10 alleles. The primers we finally selected were tested using a 25 μl PCR mix containing 12.5 μl KAPA2G Fast ReadyMix (KAPA Biosystems), 10 μM of forward and reverse primers, 1 μl of DNA and 9 μl of PCR-grad water. The PCR cycling protocol included an initial denaturation step (3 min, 95°C), followed by 35 cycles of denaturation (15 s, 95°C), annealing (15 s, 57°C), and extension (15 s, 72°C), and a final extension step (3 min, 72°C). After preliminary testing, we retained eight optimal primers (Supplementary Table 2). Forward primers were labeled with fluorescent dyes (6-FAM, Tamra, HEX, and ROX) for use in two multiplex reactions. Genotyping was performed using an ABI 3500xl Genetic Analyzer with the GeneScanT 600 LIZ size standard (Applied Biosystems). Peaks were scored manually in GeneMarker (version 1.97, SoftGenetics).

We tested the presence of null allele, stuttering, and large allele dropout using Micro-Checker (version 2.2.3; Van Oosterhout et al., 2004). Deviation from Hardy-Weinberg equilibrium (HWE) and evidence of linkage disequilibrium (LD) were calculated using Genepop (version 4.4; Rousset, 2008).

We evaluated the overall population structure and isolation by distance (IBD) using the Bayesian methods implemented in the program STRUCTURE 2.3.3 (Pritchard et al., 2000). We tested two possible ancestry and allele frequency model combinations (i.e., admixture/correlated model and no admixture/independent model), and performed four independent runs for each. Following an initial test, we set a burn-in period of 20,000 iterations and 20,000 MCMC repeats for each run. We set K (the number of possible populations) to between 4 and 15. The most likely K was inferred from ΔK (Evanno et al., 2005) as implemented in the program Structure Harvester (Earl and von Holdt, 2012). Graphical representation of structure results was produced using Distruct 1.1 (Rosenberg, 2004).

We used Nei’s corrected average of allele length differences between populations (D_A; Nei and Li, 1979), implemented in Arlequin (Excoffier et al., 2005), to evaluate the genetic distance between all pools. Nei’s corrected average of allele length differences between populations (D_A; Nei and Li, 1979) is a Rst-like distance, and for our dataset Nei’s D_A was almost linearly identical to the Goldstein et al. (1995) (δμ)² distance (r² = 0.988, P < 0.0001), which is the squared pairwise difference in mean allele length summed over all loci. Since we were interested in the differences in allele length between pairs of pools, (δμ)² distance would usually be the most appropriate measure to use. However, because (δμ)² distance is highly sensitive to outliers, we used Nei’s D_A instead. Using the D_A matrix, we explored the overall genetic structure with PROXSCAL multidimensional scaling (MDS; Commandeur and Heiser, 1993) and minimum spanning network (MSN; Excoffier and Smouse, 1994).

We examined genetic diversity using three common measures: observed heterozygosity (Ho), number of alleles per locus (allelic richness; Na), and Shannon’s diversity (H; Pielou, 1966). These diversity measures were calculated per pool using GeneAlEx (version 6.5; Peakall and Smouse, 2012). Since the number of samples (N) significantly correlated with Na (r² = 0.37, P = 0.002) and H (r² = 0.18, P = 0.049), but not with Ho (r² < 0.05, P = 0.866) or Ne (r² < 0.05, P ≥ 0.708), we adjusted pool Na and H values for n = 5 using the rarefaction approach of Chao et al. (2014) implemented in the Past program (version 4; Hammer et al., 2001).

Two of our hypotheses posited that allele length variations are associated with environmental changes and/or stability. We therefore chose to focus on three annual climatic variables of high relevance for spadefoot toad breeding behavior and tadpole survival (Gafny, 1986): the mean daily maximum temperature in January (MJT), the annual rainfall (AR), and the number of days with more than 0.1 mm rain (RD). We obtained climate data from the Israel Meteorological Service (IMS) for the past 20 years (1999–2018), from 127 permanent stations that record temperature and 559 stations that record rainfall across the country. For each year, the annual rainfall and number of rainy days were recorded from August 1 to July 31 of the following year. For each IMS station, we calculated the mean and the coefficient of variation (CV) for MJT, AR, and RD. The mean and CV data were then projected over a map of Israel by ordinary kriging and extraction by mask using ArcGIS (version 10, ESRI). These mean and CV projections (Figure 1) are added to provide the reader with a comprehensive view of the climatic gradient in temperature and rainfall and the climate stability across Israel. For each vernal pool, we calculated the mean and the coefficient of variation (CV) for MJT, AR, and RD from the closest IMS stations (Supplementary Table 1).

Pools in the north of Israel are situated on heavy soil types (terra rosa, rendzina, basaltic protogrumusols, basaltic), while the pools on the coastal plain are usually associated with lighter soil types (sand dunes, hamra, grumusols; Ravikovitch, 1992). These differences in soil types could exert a differential selection pressure on P. syriacus in regard to burrowing into the soil, food type, and nutrient availability, etc. We used the four soil characteristics found in key soil types in Israel (proportion of clay, silt, fine sand, and gross sand; Ravikovitch, 1992) to represent the habitat in the vicinity of the pools (Munwes et al., 2010).

We used Principal Component Analysis (PCA) to reduce the number of predictors and retained the components with an eigenvalue ≥ 1 (Kaiser criterion; Kaiser, 1960). Geographic distance between pools was represented in our analysis by a single predictor, the first component of the PCA on northern and eastern UTM coordinates (thereafter Geography PC1; eigenvalue = 1.9, % of variance = 96.1). Using PCA on the four soil (i.e., habitat) variables, the first two components, soil PC1 (eigenvalue = 2.3, % of variance = 57.1) and soil PC2 (eigenvalue = 1.6, % of variance = 39.8), accounted together for 97.9% of the variance in soil types. Soil PC1 correlated with clay and silt (r≅ 0.90), and soil PC2 with sand particle size (r≅ 0.79). PCA collapsed the mean daily maximum temperature in January (MJT), the annual rainfall (AR), and the number of days with more than 0.1 mm rain (RD) into a single climate predictor (thereafter climate PC1; eigenvalue = 2.6, % of variance = 86.4). The climate CV variables were collapsed into two principal components: climate CV PC1 (eigenvalue = 1.4, % of variance = 45.2), which correlated with the CV in rainfall and number of rainy days (r≅0.81), and climate CV PC2 (eigenvalue = 1.1, % of variance = 36.4), which correlated with the CV in maximum January temperature (r = 0.94). The two climate CV components together accounted for 81.6% of the variance in the climate CV variables.

To address collinearity among the predictors in our model, we used the Variance Inflation Factor (VIF), and any predictor with VIF ≥ 10 was excluded from further use. Collinearity was below the above threshold for all predictors (VIF ranged 1.3–7.9; Table 2).

Differentiating among the working hypotheses required us to associate allele length variation with both environmental variation and stability (proxy for environmental and habitat selection) and geographical distance (proxy for reproductive isolation and founder effects). We used distance-based redundancy analysis (dbRDA; Legendre and Anderson, 1999; McArdle and Anderson, 2001) to explore the association between Nei’s corrected average allele length differences between populations (D_A) and a set of climatic and environmental variables. The dbRDA is an ordination procedure, which is constrained to find linear combinations of the predictor variables that explain the greatest variation in the D_A matrix. We used marginal tests to explore the amount of variance explained by each predictor alone while ignoring all others, and sequential tests to explore the amount of variance explained by each individual predictor added in a specified order. In the sequential tests, we added geographic distance first and the climatic variables last. P-values were evaluated using randomizations. The dbRDA models were calculated using PERMANOVA + (version 7, PRIMER-E Ltd.).

We also explored the effect of geographic distance, soil type, and climate on the observed heterozygosity (Ho), rarefaction-adjusted Shannon’s diversity (H_R), and rarefaction-adjusted allelic richness (N_R) using the general linear mixed-model (GLMM). In these models, locus, geography PC1, soil type PC1 and PC2, climate PC1, and climate CV PC1 and PC2 were the predictors, and region (Table 1) was set as a random effect. We used a normal distribution and the identity link function. P-values were evaluated using randomizations. We used the total effect as a measure of variable importance (i.e., effect size) following the Saltelli (2002) method. The total effect reflects the relative contribution of each predictor, both alone and in combination with other factors, in a way that is independent of the model type and fitting method. The GLMM models were fitted using JMP Pro (version 16, SAS Inc.).

We estimated contemporary effective population size (Ne) from the microsatellite data using two approaches: The molecular co-ancestry method of Nomura (2008), and the bias-corrected linkage disequilibrium (LD) method (Waples and Do, 2010). Both methods are implemented in the program NeEstimator (version 2; Do et al., 2014).

We also estimated recent asymmetric migration rates among pools using the software BayesAss (version 3.0, Wilson and Rannala, 2003). BayesAss implements an MCMC algorithm for Bayesian inference of asymmetric rates of recent immigration among populations. We ran the program with different mixing parameters until the acceptance rates were between 20 and 60%, as recommended by the manual. The final mixing parameters were dM = 0.3, dA = 0.5, and dF = 0.6 for the coastal plain pools and dM = 0.1, dA = 0.5, and dF = 0.6 for the northern pools. Each run was set at 10⁶ MCMC iterations, 50,000 burn-in period, and a sampling interval of 100.

We sampled 22 vernal pools spanning 200 km along a sharp climatic gradient (Table 1 and Figure 1). The tadpoles we sampled were not composed of clusters of full sibs. First, each geographical region showed a wide distribution of allele frequencies, which in most cases were not strongly positively or negatively skewed (Supplementary Figure 1). Second, the frequency of full sib relationships among tadpoles (95% confidence set for relatedness = 0.5) was very low (1.37% in the coastal plain, n = 13,781 pairwise comparisons; 2.6% in the Lower Galilee, n = 3,152; 2.0% in the Upper Galilee, n = 2,831; 0.6% in the Golan Heights, n = 93,123). Last, mapping out the relatedness among tadpoles using PROXCAL MDS revealed a wide scatter within most rain pools (Supplementary Figure 2), suggesting that most of the individuals were unrelated.

The frequency of null allele, deviation from Hardy-Weinberg equilibrium (HWE) and evidence of linkage disequilibrium (LD) are presented in Supplementary Tables 3–5. We did not detect evidence for stuttering and large allele dropout. Locus Psy 2 showed a significant heterozygote deficit and was the only one suspected of having null alleles.

The polymorphic microsatellite loci we used had an allele length range of 134–384 bp, and the difference between the longest and shortest allele within each locus ranged from 68 to 154 bp. A distinct decrease was evident in mean (± SD) standardized allele length from south (coastal plain: 257.9 ± 47.7) to north (Golan Heights: 240.2 ± 55.9): that is, from the edge of the species distribution to the core (Figure 2 and Supplementary Figure 1). While a comparison of mean allele length between regions showed a significant difference [F_{(3, 31)} = 15.6, P ≤ 0.0001], the significant interaction in this model [region*locus: F_{(21, 144)} = 3.1, P ≤ 0.0001] implied that the trend of change in the regional means varies between loci. For example, the difference between northern and southern regions is not as strong in some loci (i.e., Psy 6 and 17), and Lower Galilee has a low mean in Psy 3 but a high mean in all other loci (Figure 2).

Results of the STRUCTURE analysis using two different models of ancestry and allele frequency were very similar (i.e., similar number of pool clusters). We selected to use the results of the more biologically realistic model, the admixture-correlated model. The most likely clustering according to ΔK included eight clusters (Supplementary Table 6 and Supplementary Figure 3). Pools in the coastal plain were clustered together, except for the most southern pool (CQ). The pools in the Upper Galilee were also grouped into one cluster. The pools in the Golan Heights, however, were divided in this analysis into four clusters.

The Nei’s D_A distance also showed a clear division of the pools according to geography (Figures 1A–C). The PROXCAL MDS plot (Figure 1B) shows three clear clusters (Golan Heights, Upper Galilee, and coastal plain/Lower Galilee), which accounted for 98.5% of the genetic variation (stress = 0.03). Further, the minimum spanning network connected most pools according to their geographical position (Figure 1C). The Lower Galilee (LB, LD) and Samaria (SF) pools deviated from this trend and were clustered with the coastal plain pools.

Analysis of asymmetric migration rate showed that in all pools most breeding individuals were originated from their natal pool (84.6% ± 10.6), and the rate of migration between pools was generally low (1.2% ± 1.9 and 1.7% ± 3.1 of individual were migrants in a given pool for northern and southern pools, respectively; Supplementary Table 7).

In Figure 1, we provided a spatial projection of the three climatic variables (annual maximum January temperature, annual rainfall, and annual rainy days) that most influenced amphibian breeding in the study area. We also present the spatial projection of annual CV, a measure of climatic stability, for these climatic variables. Temperature, rainfall, and number of rainy days all range widely on a gradient from north to south. Mean maximum temperature in January is about 13°C at the northernmost pools and 17.5°C at the southernmost pools (Figure 1D). Rainfall averages about 870 mm at the northern pools but only 430 mm at the southern pools (Figure 1E). The number of rainy days at the northern pools averages about 63 days, while at the southern pools it averages about 43 days (Figure 1F). Maximum temperature during January becomes less stable with increasing elevation, with coastal locations experiencing more stable temperatures among years (Figure 1G). Rainfall (Figure 1H) and the number of rainy days (Figure 1I) become less consistent among years toward the south.

To determine whether geography, soil type (i.e., habitat), and climate characteristics can explain the variation in pool allele length, we tested for association between Nei’s D_A (i.e., a distance matrix between pools based on allele length) and a set of predictors, using dbRDA (Table 2 and Figure 3). The marginal tests showed that all predictors, except climate CV PC1, were significantly correlated with the Nei’s D_A distance, and that geography, climate, and climate CV (PC2) each accounted for at least 50% of the variance in the D_A distance matrix. In the sequential test, we ordered the predictors such that the geography predictor (i.e., isolation by distance and founder effect hypotheses) was entered first into the model, and the climatic predictors (i.e., natural selection and stress-induced mutagenesis hypotheses) were entered last. Geography significantly accounted for 59% and climate accounted for an additional 14% of the variance in the D_A distance. All the other predictors did not account for any unique variation that was not already explained by geography and climate (Table 2). The dbRDA ordination diagram demonstrates that geography and the climatic factors account for separate trajectories in the variation of the D_A projected data (Figure 3). Using the same dbRDA model with the (δμ)² distance revealed both geography (56%, P = 0.001) and climate (13%, P = 0.003) also to be significant effects.

Additional dbRDA locus-by-locus models revealed that for specific loci the Nei’s D_A distance was significantly influenced by the geography of pools [Psy 1: F₍₂₀₎ = 50.5, P = 0.001; Psy 2: F₍₂₀₎ = 18.1, P = 0.007; Psy 4: F₍₂₀₎ = 16.0, P = 0.007; Psy 5: F₍₂₀₎ = 22.8, P = 0.002], soil composition [Psy 1: PC2, F₍₁₈₎ = 18.8, P = 0.041; Psy 2: PC1, F₍₁₉₎ = 18.9, P = 0.024; Psy 8: PC1, F₍₁₉₎ = 10.9, P = 0.043 and PC2, F₍₁₈₎ = 47.3, P = 0.019], and climate [Psy 3: F₍₁₇₎ = 21.4, P = 0.024; Psy 5: F₍₁₇₎ = 54.2, P = 0.001]. Psy 6 and Psy 17 were not significantly influenced by any of the predictors.

Since the heterozygote instability and the genetic drift hypotheses are both associated with population demography, we also explored the association between the above predictors and Ho, H_R and N_R; all of which could reflect on population size or the presence of a bottleneck (Table 1). Ho was not significantly associated with the Nei’s D_A matrix [dbRDA marginal test: F₍₂₀₎ = 0.0, P = 0.999]. H_R and N_R were not significantly associated with the geography of pools (F = 0.5, P ≥ 0.36 in both models; Table 3). For Ho, only loci Psy 2 and 5 were significantly associated with geography (Supplementary Figure 4). However, both N_R and H_R were positively associated with climate PC1, soil PC1, and climate CV PC1, and negatively associated with soil PC2 (Table 3 and Supplementary Figure 4). Psy 5 was the only locus that did not associate with climate for H_R and N_R, and the only locus negatively associated with climate for Ho (Supplementary Figure 4). In other words, genetic diversity (i.e., H_R and N_R) in seven out of eight loci was associated with climate and soil variation. Most importantly, the total effect in these models showed that climate is by far (0.56–0.57) the most influential effect on H_R and N_R (Table 3). Ho was most affected by the difference between loci (total effect = 0.725), and secondly by the climate variation (total effect = 0.391).

Finally, we also tested for regional differences in effective population size (i.e., N_e) in order to directly address the predictions of the heterozygote instability and the genetic drift hypotheses. N_e estimates (i.e., log-transformed contemporary effective population size; Supplementary Table 8) from the co-ancestry model [GLMM; F_{(3, 15)} = 0.51, P = 0.684] and the bias-corrected linkage disequilibrium (LD) model [F_{(3, 17)} = 0.06, P = 0.978] did not significantly differ between regions. Further, neither of the N_e estimates significantly correlated with Ho [co-ancestry model: r = 0.21, F_{(1, 17)} = 1.02, P = 0.326; linkage disequilibrium model: r = 0.09, F_{(1, 6)} = 0.78, P = 0.409], H_R [co-ancestry model: r = 0.36, F_{(1, 16)} = 02.55, P = 0.129; linkage disequilibrium model: r = 0.33, F_{(1, 17)} = = 2.03, P = 0.172], and N_R [co-ancestry model: r = 0.36, F_{(1, 17)} = 2.48, P = 0.134; linkage disequilibrium model: r = 0.36, F_{(1, 17)} = 2.45, P = 0.136].

Under these hypotheses we expected to observe a correlation between Ne and diversity, and in the case of heterozygote instability, a correlation with mean allele length (the sign of this correlation may be either negative or positive). Our data analyses revealed an increase in allele length from northern to southern populations, but this change in allele length was not significantly correlated with Ho levels nor with N_e. N_e was similar among the four geographic regions using both estimates (i.e., from the co-ancestry and the LD models) and did not significantly correlate with any of the diversity measures. Furthermore, pool tadpole population densities in the south (i.e., edge of distribution) were not significantly lower than those of the pools in the north of Israel (Munwes et al., 2010). The lack of difference in Ne between regions and the lack of correlation between Ne and allele diversity do not fit the key prediction of the genetic drift hypothesis: lower genetic variation in pools of small Ne. Moreover, the lack of correlation between Ne and allele size variation does not fit the prediction of the heterozygote instability hypothesis. Finally, the Ho and diversity measures did not significantly correlate with geographic distance, the effect size for geographic distance was small, and its model estimates did not differ from zero.

Under this hypothesis we expected allelic diversity to positively correlate with environmental stress, regardless of population size and heterozygosity level, while mean allele length could either increase or decrease with stress. In this study, stress was considered to be correlated with an environmental gradient (from mild to harsh) or with the annual stability of environmental conditions (e.g., temperature, rainfall), measured by coefficient of variation over years. Our results revealed that allelic diversity significantly varied with the climate but not with the climatic stability predictors (Table 2 and Figure 3). The harsher living conditions in the south (i.e., at the distribution edge; Munwes et al., 2010) are expected to result in a higher mutation rate, which is expressed by longer alleles (e.g., Whittaker et al., 2003). The similar range of allele sizes found in the different regions along the gradient was not predicted by this hypothesis but could be explained by the more intense selection in the southern populations that restricted a wider range of alleles. Further, microsatellite diversity in the southern populations, where climatic conditions for P. syriacus survival and breeding are marginal (Munwes et al., 2010), were similar to those of the northern populations, contrasting with the prediction of the stress-induced mutagenesis hypothesis of higher microsatellite diversity in the southern populations.

Under this hypothesis (i.e., direct selection, where alleles are either longer or shorter with gradient), we expected allele length to negatively or positively correlate with a climatic or environmental gradient or with environmental stability, independent of genetic diversity or heterozygosity. Overall, our results revealed that climate had a strong significant effect on allele length. The amount of rainfall is a predictor for the hydroperiod (i.e., the duration of water availability in the vernal pool), and the number of rainy days is a key signal for the onset of breeding in P. syriacus (Gafny, 1986; Munwes et al., 2010). Our results are in line with those of other recent other studies and reviews that have provided evidence of an association between genetic structure and the environment (e.g., Schmidt and Anderson, 2006; McCairns and Bernatchez, 2008; Ortego et al., 2012; Stojak et al., 2019). Further support for an environmental effect on microsatellite variation comes from the study by Cohen et al. (2021), which revealed a significant correlation between climatic variation among successive years, primarily the amount of rainfall and number of rainy days, and both microsatellite allelic composition and diversity (15% of the total explained genetic variance). Cohen et al. (2021) sampled seven rain pools, which we also sampled for the present study, making their earlier study’s results highly relevant. However, explaining our results according to the natural selection hypothesis is currently controversial, as most microsatellite loci in the genome are generally thought to be neutral, with allele length changing due to drift.

Another interesting result was that geographical distance (i.e., IBD) showed a strong correlation with variation in allele length among pools but not with allele diversity. In our system, climate explained much of the variation in pool allele diversity. Notably, P. syriacus is a highly philopatric species (Gafny, 1986); a fact supported by the exceptionally low rates of migration we calculated among vernal pools (Supplementary Table 7). Furthermore, in the recent past (i.e., late nineteenth century) the temporary wetland habitat during winter was vast and continuous, enabling individual toads unrestricted migration to nearby pools (Levin et al., 2009). Over the past 150 years, the wetland habitat, mostly in the coastal plain, has dramatically declined (> 90%), resulting in the physical isolation of temporary pools. Given the toads’ philopatric behavior and the loss of wetlands, we expected a strong signature of isolation by distance, which was indeed observed in the allele length variation. However, we also found that the effect of climate was meaningful (i.e., 14% of the variation in allele length). These findings concur with those of Munwes et al. (2010), who showed that soil type and hydroperiod were both predictors of mtDNA substructure for a similar set of pools, and together accounted for 25% of the variance in genetic distances between pools. Such variation in pool diversity measures is often associated with population size and founder effect, which are also related to the location of pools relative to the core of distribution. In the present study system, we did not detect differences in contemporary Ne, and the variation in allele diversity between pools was largely explained by the climate regime.