Can neutral genetic differentiation explain geographical variation in body size of the natterjack toad, Epidalea calamita?

Acta Herpetologica 18(2): 69-78, 2023 ISSN 1827-9635 (print) © Firenze University Press ISSN 1827-9643 (online) www.fupress.com/ah DOI: 10.36253/a_h-13774 Can neutral genetic differentiation explain geographical variation in body size of the natterjack toad, Epidalea calamita? Federico Marangoni Department of Evolutionary Ecology, Estación Biológica de Doñana, CSIC, Avda. Américo Vespucio s/n, 41092 Sevilla, Spain Present address: Departamento de Biología, Facultad de Ciencias Exactas y Naturales y Agrimensura. Universidad Nacional del Nordeste (FACENA-UNNE) and Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET). Av. Libertad 5470, 3400 Corrientes, Argentina E-mail: fedemarango[email protected]m Submitted on: 2022, 14th September; revised on: 2023, 25th April; accepted on: 2023, 10th July Editor: Raoul Manenti Abstract. Population genetic studies are crucial for evolutionary biologists because the population is the basic substrate on which evolution is forged. However little empirical evidence has been able to demonstrate the role that isolation and gene flow play in maintaining differentiation in populations at short geographic scales. Epidalea calamita exhibits a steep variation in body size and reproductive traits in southwestern Spain, associated with changes in the geological substrate. This implies a decrease of 70.9% of body mass and 28.5% in snout-vent length, on a microgeographic scale of only 60 km. Previous results from both metamorphic and juvenile common garden experiments showed that genetic differentiation may be a causal determinant of geographic variation in adult. This study tested whether neutral genetic differentiation can explain the geographical variation in the body size observed in E. calamita. It was addressed analyzing the level of genetic structuring and gene flow among populations along the cline, comparing the genetic diversity between and within populations, as well as between ecological environments. The study showed that the geographic variation in body size observed in E. calamita has evolved in absence of geographic isolation, with moderate gene flow connecting the populations. Thus, neutral genetic differentiation cannot explain the geographical variation observed. Future studies are needed on the interaction between the genetic component with the environmental factors and will be necessary to analyze the contribution of the maternal effects in the origin and evolution of the geographical variation in the body size observed in E. calamita from southern Spain. Keywords. Epidalea, FST, microsatellite loci, population differentiation, body size. INTRODUCTION Geographic variation in phenotypic and genetic characteristics among species’ populations is a phenomenon that has been very well documented since the 1950s (Stebbins, 1950; Mayr, 1963; Harper, 1977). Adaptive explanations for the evolution and maintenance of geographic variation in body size have been put forward, particularly considering macrogeographical patterns as a response to environmental gradients (Bergmann, 1847; Ray, 1960; Lindsey, 1966; Adams and Church, 2008; Ashton, 2002; Cvetkovic et al., 2009; Sinsch et al., 2010). Nevertheless, few studies have evaluated it at smaller geographical scales with, in many cases, lack of genetic isolation between populations (e.g., Skelly, 2004; Gomez-Mestre and Tejedo, 2004; Lee et al., 2020; Albert and GarcíaNavas, 2022), and thus, what causes and maintains those patterns is still not well understood. Therefore, to infer on the processes causing these patterns, we need to know the mechanisms underlying the observed phenotypic variation, how they are connected to genetic variation, and how they interact with other traits and the environ- 70 Federico Marangoni ment (Stearns, 1989). This will also help us to understand the evolutionary significance of the geographic variation, which ultimately can lead to the formation of new species (Endler, 1977; Foster and Endler, 1999). The agents that change the gene frequencies of populations, that is, the factors of evolution, are mutation, genetic drift, gene flow, and natural selection (Slatkin, 1987). While drift and selection tend to increase population differentiation, gene flow promotes homogenization among connected populations, either increasing or decreasing the genetic diversity of the system (Lenormand, 2002). Thus, gene flow is a major component of population structure because it determines the extent to which each local population of a species is an independent evolutionary unit. If there is a high gene flow between local populations, then all the populations evolve together; whereas in the presence of low gene flow, each population evolves almost independently (Slatkin, 1985). The natterjack toad (Epidalea calamita) populations from southwestern Spain exhibit a steep variation in body size and reproductive traits associated with changes in the geological substrate (Marangoni et al., 2008). This implies a decrease of 70.9% of body mass and 28.5% in snout-vent length, on a micro-geographic scale of only 60 km (Fig. 1). Previous studies suggested that considerable genetic differentiation may be the mechanism underlying the observed geographic variation in metamorphic traits in E. calamita (Marangoni, 2006) and Pelobates cultripes (Marangoni and Tejedo, 2008), which exhibit the same geographic pattern of adult body size variation (Marangoni et al., 2008; Lee et al., 2020). Moreover, the study of age structure and growth pattern across populations of E. calamita suggests that both environmental variations in resources availability associated with the sandy substrate (Marangoni, 2023), but also different growth and maturity pathways, may happen in response to contrasting selective pressures (Marangoni et al., 2021). Two hypotheses could be suggested to explain the evolution and maintenance of the observed cline in body size and reproductive parameters in E. calamita (Marangoni et al., 2008). In the first place, it could be expected that in the presence of gene flow between populations, the alleles favored by a selection pressure of intensity s do not decrease their high frequency, because the rate of immigration m of other alleles is lower than the intensity of selection s, that is, m < s (Slatkin, 1985). The second hypothesis would be that the differentiation between the populations along the cline had occurred in a context of relative population isolation or in the presence of scarce gene flow. Little empirical evidence has been able to demonstrate the role that gene flow plays in maintaining differentiation in populations at short geographical scales. It has been suggested that when gene flow is not homogeneous, evolutionary differentiation can be rapid and can occur on small spatial scales (Kennington et al., 2003; Garant et al., 2005; Postma and van Noordwijk, 2005). In previous studies on E. calamita that exhibited local adaptation to osmotically stressful environments, microsatellite markers revealed little population differentiation, lack of an isolation-by-distance pattern, and moderate gene flow connecting the populations (Gomez-Mestre 2001; Gomez-Mestre and Tejedo, 2004). Present study assess whether neutral genetic differentiation can explain the geographical variation in the body size, age and reproductive parameters observed in natterjack toad, from southern Spain (Marangoni et al., 2008; 2021). This comparison between patterns of population genetic differentiation in neutral markers and quantitative traits can provide valuable insights into the mechanisms driving variation within populations (Leinonen et al., 2007). Thus, this comparison can help us understand evolutionary forces by discerning the relative influences of genetic drift and natural selection on the evolution of quantitative traits (Gomez-Mestre and Tejedo, 2004; Knopp et al., 2007; Páez-Vacas et al., 2001). Secondly, it aids in identifying adaptive traits, for example those traits which might are related to adaptation to different environmental conditions (Gomez-Mestre and Tejedo, 2003; Luque et al., 20015). In addition, it also allows us to understand the distribution of genetic effects within the architecture of quantitative traits (see Leinonen et al. 2007 and references therein). The main goals in the present study were: i) analyze the level of genetic structuring and gene flow in E. calamita populations along the geographic variation in body size observed (Marangoni et al., 2008), ii) compare genetic diversity between and within populations, as well as between ecological environments and iii) test for the existence of an isolation-by-distance pattern of populations differentiation. We expect that the observed geographical variation in the body size of natterjack toad has been evolved in absence of isolation-by-distance, and with a gene flow connecting the populations. MATERIAL AND METHODS Populations Seven populations of Epidalea calamita were selected, representing the cline in body size observed in a previous study (Marangoni et al., 2008), which encompass three areas with different geological substrates. These populations included: two Large-bodied populations, Pedroso (L1) and Navas (L2), from Sierra Morena (old hercinic granites schist soils): four Small-bodied popula- 71 Neutral genetic differentiation in Epidalea calamita tion, Juncosilla (S1), Bodegones (S2), Abalario (S3) and Reserva Biológica de Doñana (S4), from the Doñana area (quaternary sandy eolian deposits): and one population with intermediate body size, Sanlúcar (I), hereafter), geographically located between Sierra and Doñana (mixture of clays and sand) (Fig. 1). A more detailed description of the environments and the life history traits (body size, age, growth patterns and reproductive output) of the populations selected, is available in Marangoni and Tejedo (2008) and in Marangoni et al. (2008, 2021). Sampling and molecular genetic analysis Twenty recently laid clutches of similar larval stage (Gosner stage 10, Gosner, 1960) from each of the seven population of Epidalea calamita were sampled during the breeding season (January 2003). Each clutch (fullsib families) were brought to the laboratory and kept separately in plastic trays filled with dechlorinated tap water until tadpoles reached Gosner’s stage 25 (Gosner, 1960), to be included in a common garden experiFig. 1. Location and geological substrate of the studied Epidalea calamita populations. Abbreviated names of sampling localities, geographic coordinates (Coordinates UTM x/y in meters, Datum European 1950, Spain and Portugal, Zone: 30), elevation, and sample size are as follows: L1, 255170/4190574, 395 m, n = (39); L2, 229255/4187617, 420 m, n = (43); I, 213349/4144548, 34 m, n (22); S1, 203208/4124509, 23 m, n = (44); S2, 175577/4120711, 32 m, n = (45); S3, 174267/4115417, 63 m, n = (35); S4, 188450/4102197, 24 m, n = (42). Circles: Sierra (Paleozoic, Granite + Schist rocks), square: Intermediate (Miocene-Pliocene, Clay + Sandy), triangle: Doñana (Holocene, Sandy soil). Below each abbreviated names of sampling localities are indicated the mean body size (snout-vent length in mm/weight in g) from Marangoni et al. (2008). 72 Federico Marangoni ments (Marangoni, 2006). Tissue samples for the present genetic study were obtained by cutting the tip of the tail of 15-20 tadpoles from each population. They were randomly taken from a sample in which the twenty clutches from each population were previously mixed, to maximize the chances of sampling unrelated tadpoles. In present study were used the same eight microsatellite loci that previously were analyzed by Gomez-Mestre (2001) and Gomez-Mestre and Tejedo (2004), which included Bcal 1, Bcal 2, Bcal 3, Bcal 4, Bcal 5, Bcal 7 (Rowe et al. 1997), Bcal 10, and Bcal 11 (Rowe et al., 2000). The 5’ to 3’-primers were labeled with a color fluorophore, either HEX, TET, or FAM (Gomez-Mestre, 2001). DNA was obtained using an DNA Dneasy Tissue Extraction Kit (QIAGEN). The DNA was extracted by digesting each sample, approximately 25 mg of tissue minced into small pieces, to which 180 µl of ATL buffer and 20 µl of proteinase K were added. The samples incubated for 12 hours at 55 °C. Afterwards, if the tissue was not completely digested, an additional 20 µl of proteinase K was added and incubated for another 3 hours. Once the samples were fully digested, they were shaken for 15 seconds, 200 µl of buffer AL was added and incubated at 70 °C for 10 minutes. Subsequently, after adding 200 µl of ethanol, the mixture was centrifuged for 1 minute at 8000 rpm using the columns provided by the kit. Transferring the columns to other 2 ml tubes, 500 µl of buffer AW1 was added and centrifuged once more. This last step was repeated one more time, but adding 500 µl of buffer AW2, and centrifuged for 3 min at 15,000 rpm. Finally, once the columns were transferred to 1.5 ml Eppendorf tubes, 200 µl of AE buffer was added, incubated at room temperature for 1 min and centrifuged at 8000 rpm. This step was repeated twice, obtaining a final volume of 400 µl of buffer and the DNA extracted from the sample, which was stored at –20 °C until further procedures. Loci amplification was conducted using polymerase chain reactions (PCRs) of 1.5 ml total volume with 5 ml of DNA. The PCR amplifications followed a 66-50 °C touchdown procedure. PCR products of each of the eight loci for each individual sampled were aliquoted and mixed according to their abundance. An aliquot of 1.5 ml of the resulting mix was added to 13 ml of formamide plus 0.3 ml of Tamra 500 (Applied Biosystems, Foster City, CA) standard. Samples were analyzed in a fluorescence-based automatic fragment analyzer (ABI-PRISM 310 Genetic Analyzer, Applied Biosystems). Allele sizes for each locus were resolved by comparison of the peaks obtained to those yielded by the Tamra 500 standard using the software GeneScan version 3.1.2 (Applied Biosystems). All procedures described were performed in the Laboratorio de Ecología Molecular, at the Estación Biológica de Doñana (EBD-CSIC), Seville, Spain. Permits for capture and sampling of E. calamita, including all ethical considerations, were acquired from the regional authorities. Statistics Arlequin software packages (Schneider et al., 2000) was used to perform a nested molecular analysis of variance (nested AMOVA) comparing genetic diversity between and within populations, as well as between ecological environments population groups (Large-bodied populations from Sierra and Small-bodied populations from Doñana). This analysis provides genetic structure parameters in the form of F statistics (Wright, 1951; Excoffier et al., 1992). The significance of these statistics is assessed with the null distributions of the variance components, since both are highly correlated (Excoffier et al., 1992). Using the Arlequin software packages, were analyzed the allele frequencies, mean number of alleles per locus, and observed and expected heterozygosity. The significance of each of the variance components in subsequent analysis was estimated using 10,000 permutations. Genepop (version 3.4 online; Raymond and Rousset, 1995) was used to: i) analyze the population differentiation computing FST estimators (Weir and Cockerham, 1984) and their confidence intervals, ii) tests for linkage disequilibrium between loci using a Fisher exact test using Markov chains, and iii) analyze the existence of isolation-by-distance through Mantel tests carried out between matrices of log-transformed geographic distances and odds-transformed genetic distances (FST/[1-FST]; Rousset, 1997). Mantel Tests were also performed using Isolation By Distance Web Service (IBDWS) Version 2.5 (Jensen et al., 2005). Departures of the allelic frequencies from Hardy-Weinberg expectations were performed with the GENETIX Version 4.04 program (Belkhir et al., 2000) and allelic richness was estimated using the FSTAT program (Goudet, 1995). Finally, was conducted the analysis at a significance level of α = 0.05, and applied the DunnSidak sequential correction of the level of significance for multiple tests (Sokal and Rohlf, 1995) when necessary. RESULTS No large differences between populations were found regarding genetic diversity and allelic richness (Table 1). Moreover, no allele had a frequency greater than 95%, which indicates that all analyzed loci were polymorphic. The mean number of alleles per locus ranged between 8 and 25 (Bcal 2 and Bcal 4, respectively), while in the remaining loci were: Bcal 1 = 23, Bcal 3 = 17, Bcal 4 = 25, 73 Neutral genetic differentiation in Epidalea calamita Bcal 5 = 19 and Bcal 7 = 18. Additionally, the mean number of alleles per locus in each studied population ranged between 9.75 (I) and 12.75 (S4). No private alleles were found, indicating that each allele was found in two or more populations. However, allele frequencies significantly differed across populations for all loci. The analysis of linkage disequilibrium showed that only 19 of all possible comparisons between pairs of loci in each population (28 × 7 = 119) were significant for α = 0.05, which represented 9.6%. However, all of them lost their significance after applying the Dunn-Sidak significance level correction for multiple comparisons. Genepop estimates for pairs of loci taking all populations together did not yield any significant evidence of linkage disequilibrium (for α = 0.05), so the eight loci were treated as independent henceforth. Furthermore, were detected significant departures from Hardy-Weinberg expectations in different loci across all populations (Table 1). Bcal 11 and Bcal 2 were departed in all analyzed populations, while Bcal 4 in six of them. The population with the highest number of loci (7) outside the H-W expectations was S1, while in S2 we found the lowest (4). The rest of the populations presented five loci outside the H-W balance. Observed heterozygosity in the studied loci was generally lower than expected under Hardy-Weinberg equilibrium, indicating a deficiency of heterozygotes. Regarding the relation between geographic distance (km) and, estimated distance genetics (Fst) and gene flow (Nm) of the pairwise comparisons between populations are shown in Table 2. There were significant differences at 13 of the 21 pairwise comparison between populations, before Bonferroni sequential correction (shown in bold in the Table 2). Nevertheless, all of them lost statistical significance after the Bonferroni correction was applied (critical level of significance was P = 0.0023; Sokal and Rohlf, 1995). Neutral genetic distance among populations was not correlated with geographical distance (r = 0.140, P = 0.217), rejecting the hypothesis of isolation by distance (Fig. 2). Two groups were clearly differentiated, one containing the pairwise comparisons between nearby populations (SS or LL), and another group containing the pairwise comparisons of distant populations (SL or LS), however these two groups did not show significant differences in the mean values of genetic distances (Fig. 2). For example, the S4 population had similar genetic distances when compared with both its most distant population (L2, Fst = 0.056/110 km) and its closest one (S3, Fst = 0.042/19 km). I also did not find any significant covariation between geographic distance and gene flow (y = 2.027 - 0,003*x; r = -0.129; P = 0.578; r2 = 0.017). The two AMOVA performed showed that the 95.4% of overall variation was held within populations, whereTable 1. Molecular variation in populations of Epidalea calamita from southern Spain. Populations (mean number of alleles), number of alleles per loci (NA), expected (He) and observed (Ho) heterozygosity, genetic diversity (GD) and allelic richness (AR). Significant deviations from Hardy-Weinberg equilibrium are marked with an asterisk. Locus L1 (11.1) L2 (13) I (9.7) S1 (12) S2 (12.4) S3 (12.5) S4 (15.2) NA HeHoNA HeHoNA HeHoNA HeHoNA HeHoNA HeHoNA HeHo Bcal µ1 17 0.901 0.897 18 0.891 0.837 14 0.888 0.818 15 0.846 0.704* 14 0.884 0.777 13 0.863 0.885 14 0.855 0.785* Bcal µ2 7 0.762 0.575* 5 0.500 0.138* 5 0.555 0.200* 5 0.651 0.232* 6 0.507 0.157* 6 0.527 0.416* 6 0.559 0.114* Bcal µ3 11 0.871 0.281* 14 0.878 0.333* 10 0.879 0.272* 13 0.886 0.405* 15 0.902 0.342* 14 0.875 0.428* 14 0.898 0.323* Bcal µ4 16 0.889 0.550* 17 0.880 0.555* 12 0.885 0.381* 16 0.862 0.780 15 0.894 0.522* 17 0.859 0.512* 18 0.801 0.642* Bcal µ5 8 0.735 0.525* 12 0.752 0.767 10 0.792 0.736 12 0.818 0.697* 11 0.803 0.681 11 0.836 0.666* 14 0.860 0.795 Bcal µ7 11 0.807 0.717 13 0.792 0.613* 6 0.420 0.450 11 0.616 0.613* 14 0.752 0.600 13 0.717 0.717 13 0.773 0.744 Bcal µ10 9 0.854 0.906 11 0.840 0.795 9 0.807 0.727* 9 0.839 0.846* 11 0.832 0.750 10 0.810 0.815 8 0.817 0.825 Bcal µ11 10 0.756 0.240* 14 0.826 0.282* 12 0.883 0.181* 14 0.887 0.620* 13 0.851 0.410* 16 0.897 0.552* 15 0.920 0.700* continuation: GD AR GD AR GD AR GD AR GD AR GD AR GD AR Bcal µ1 0.913 12.948 0.903 12.308 0.911 12.298 0.858 10.063 0.895 10.496 0.876 10.558 0.867 9.970 Bcal µ2 0.777 5.825 0.513 4.019 0.588 5.000 0.664 4.831 0.519 5.255 0.537 4.709 0.574 5.641 Bcal µ3 0.895 9.396 0.901 11.591 0.915 9.538 0.906 10.917 0.922 12.125 0.895 10.935 0.921 11.889 Bcal µ4 0.905 11.969 0.897 11.400 0.920 10.959 0.874 11.265 0.910 12.094 0.875 11.283 0.813 11.268 Bcal µ5 0.748 6.127 0.761 8.076 0.816 9.241 0.830 8.996 0.814 8.535 0.850 8.457 0.872 10.435 Bcal µ7 0.820 8.420 0.804 9.089 0.430 5.236 0.623 6.924 0.763 8.291 0.727 8.643 0.783 8.468 Bcal µ10 0.867 8.088 0.851 8.565 0.829 7.853 0.851 8.097 0.843 8.521 0.821 7.685 0.828 7.464 Bcal µ11 0.783 8.652 0.845 9.845 0.921 11.274 0.908 12.174 0.868 10.248 0.914 11.952 0.935 12.912 74 Federico Marangoni as the lowest variance components were associated to the differences between environments (Sierra Morena vs Doñana), which suggests a lack of population substructuring (Table 3). DISCUSSION This study showed the highest allelic diversity found in Epidalea calamita populations. This was greater than that previously found by Gomez-Mestre and Tejedo (2004) in Spanish populations, and even greater than that found in the northernmost populations of the specie distribution (Rowe et al., 1998; Beebee and Rowe, 2000) (Table 4). The differences in the allelic diversity between this study and that of Gómez-Mestre and Tejedo (2004), particularly in the S1 population (Fresh 1 in GómezMestre and Tejedo, 2004), could be attributed to a significant increase in the sample size in this study. However, I found no difference in the L1 population (Fresh 3 in Gómez-Mestre and Tejedo, 2004) (Table 4). The lower allelic diversity found in the British populations also supports the hypothesis that the Iberian Peninsula constituted a Pleistocene refuge for E. calamita (as it was for other species in other Mediterranean peninsulas; Hewitt, 1996; Taberlet et al., 1998). From this refuge, the species would have expanded rapidly to north and east during the postglacial stage (Beebee and Rowe, 2000), resulting in a pattern of high levels of genetic diversity in populations derived from southern refuge and a progressive loss of diversity in recolonized areas to the north (Avise, 1994). Despite the high variability, there were significant departures from Hardy-Weinberg expectations in different loci across all populations, due to deficiency of heterozygotes. It is possible that there is a certain degree of variability in the primer pairing regions, so that their sequence would not be completely homologous and would fail to amplify some alleles. In this case, by amplifying only one of the two alleles present, the proportion Table 2. Geographical distance (km), distance genetic (Fst) and gene flow (Nm) between populations of E. calamita. Significant values before Bonferroni sequential correction are shown in bold. Pairwise comparison Geographic distance Genetic distance Gene flow S2-S3 5.5 0.006 41.7 S3-S4 19.4 0.042 5.86 S2-4 22.5 0.039 6.35 S1-I 22.5 0.094 2.65 L1-L2 26.1 0.042 5.88 S1-S4 26.8 0.06 4.14 S1-S2 27.9 0.02 12.13 S1-S3 30.3 0.025 9.66 S2-I 44.7 0.064 3.86 L2-I 45.9 0.021 11.39 S3-I 48.7 0.079 3.13 S4-I 49.1 0.024 10.18 L1-I 62.2 0.07 3.56 L2-S1 68.3 0.068 3.63 L1-S1 84.1 0.032 7.76 L2-S2 85.8 0.046 5.41 L2-S3 90.8 0.062 4.03 L2-S4 94.7 0.025 9.64 L1-S2 105.9 0.015 16.31 L1-S3 110.4 0.038 6.49 L1-S4 110.7 0.056 4.41 Figure 2. Pairwise comparisons between populations from the three environments studied, relating the geographic distance and estimated genetic distance. Two groups are clearly differentiated, one (inner square with continues lines) containing the pairwise comparisons between nearby populations (SS or LL), and another group (inner square with broken lines) containing the pairwise comparisons of distant populations (SL or LS). Nevertheless, mean values of genetic distances (middle lines) are not significantly different. L = large-bodied population (Sierra), I = intermediate and S = smallbodied population (Doñana). Table 3. Nested molecular analysis of variance, where df stands for degrees of freedom, SS for sum of squares, Varcomp for variance components, and %Var for proportion of total variance accounted for by each source. Environments are Hercinic (L1 and L2) and Sandy soils (S2-S4). Source of variation df SS Varcomp % Var Between environments 1 19.985 0.023 0.69 Among populations within environments 4 58.991 0.133 3.91 Within populations 508 1660.179 3.268 95.4 Total 513 1739.156 3.425 100 75 Neutral genetic differentiation in Epidalea calamita of heterozygotes may have been underestimated, which would move the observed frequencies away from those expected for a Hardy-Weinberg equilibrium situation (Gómez-Mestre, 2001). However, we cannot rule out the homogenizing effect that gene flow (see below) may have in the departures from Hardy-Weinberg expectations observed. In addition, we could consider that the genes responsible for the expression of body size are potentially under selection and need to be studied. Considering that some estimates of gene flow between populations were quite high, the lack of a relation between geographical distance and the degree of genetic differentiation (Fig. 2) is well exemplified by the L1 population. This population is genetically more similar to the S2 population, located 105 km apart, than to the L2 population, located only 26 km apart. As most of the observed variability corresponds to within-population differences (95.4%), this suggests that the populations are not structured. In addition neutral genetic differentiation cannot explain the geographical variation in body size observed, since only 0.6% of the total variance could be attributed to differences between the two environments. The evolution and geographic variation of the Mediterranean herpetofauna has been influenced by a succession of geographic barriers to faunistic exchange over the last 23 myr (López-Martínez, 1989). The Guadalquivir River basin in the present study area has been suggested as major factor in the speciation processes in amphibians and as a barrier to dispersal of Salamandra salamandra (García-París et al., 1998). This intracontinental barrier has also been suggested as responsible of the geographic variation pattern in water salinity tolerance among Epidalea calamita populations in southern Spain (GomezMestre, 2001). Nevertheless, since all our studied populations are all geographically located on the west bank of the Guadalquivir River, without population in the east, we cannot evaluate the hypothesis that considers the river as a barrier to gene flow. A possible explanation for the lack of isolation-bydistance could be associated with the high dispersive capability that Bufonidae species can potentially present (e.g., 1.3 km/night, in Rhinella marina, Leblois et al., 2000), which could result in high gene flow despite populations being geographically distant. However, although this hypothesis is probable, it does not explain the high difference between closer populations (e.g., L1 and L2, or S3 and S4, Table 2). Alternatively, other less conspicuous physical barriers could exist between closer populations and need to be evaluated in futures studies. So, given the pattern of the gene flow observed in our study, there are some possible explanations for the maintenance of the geographical variation in the body size of Epidalea calamita. Populations subject to selection in two or more environmental patches, completely connected by gene flow, may develop phenotypic plasticity or adaptive reaction norms (Schmalhausen, 1949; Bradshaw, 1965), such that genetically similar individuals express different phenotypes in each environment. Then, it could be expected that the alleles that cause different phenotypes in E. calamita between Sierra and Doñana environments can evolve by natural selection, if the plastic response of the phenotype produces an increase in biological fitness (Via and Lande, 1985), as was demonstrated in Pelobates cultripes using reciprocal transplant experiments (Marangoni, 2006). This and other studies made in E. calamita and newts have shown that the dwarfism could be involved in response to some common environmental factor in Doñana (Díaz-Paniagua et al., 1996; Díaz-Paniagua and Mateo, 1999). It is clear that sandy soil substrates, directly or indirectly impose a strong effect (e.g., by imposing Table 4. Genetic diversity. Mean number of alleles per locus (MAPL), percentage of polymorphic loci (P95), and expected (He) and observed (Ho) heterozygosity of E. calamita throughout its distribution range. PS (in bold): present study, 1: Gómez-Mestre and Tejedo (2004), 2: Beebee and Rowe (2000). * and #, large and smallbodied population respectively, are showing the same population studied in PS and 1, using the same eight microsatellite loci. Population N MAPL P95 He Ho Source L1 (Spain)* 39 11.12 100 0.822 0.586 PS. Fresh 3 (Spain)* 22 11.25 100 0.767 0.677 1 S4 (Spain) 42 12.75 100 0.810 0.616 PS. S2 (Spain) 45 12.37 100 0.803 0.530 PS. S3 (Spain) 35 12.5 100 0.798 0.624 PS. S1 (Spain)# 44 12 100 0.801 0.612 PS. Fresh 1 (Spain)# 22 10.5 100 0.698 0.581 1 L2 (Spain) 43 13 100 0.795 0.540 PS. I (Spain) 22 9.75 100 0.764 0.471 PS. Saline 1 (Spain) 23 9.38 100 0.664 0.626 1 Fresh 2 (Spain) 20 7.5 100 0.629 0.534 1 Velez (Spain) 11 4.88 100 0.689 0.607 2 Brittany (France) 32 4.38 87.5 0.491 0.355 2 Boulogne (France) 15 3.88 75 0.461 0.455 2 Ooy-Polder (The Netherlands) 40 5.13 87.5 0.520 0.466 2 Kerry (Ireland) 40 2.38 62.5 0.344 0.335 2 Merseyside (England) 40 2.63 62.5 0.294 0.295 2 Cumbria (England) 40 3.75 75 0.391 0.344 2 E/SE (England) 40 2.50 75 0.352 0.289 2 Texel (The Netherlands) 40 2.63 62.5 0.367 0.430 2 Sweden 40 1.63 25 0.119 0.144 2 Poland 40 2.00 62.5 0.245 0.283 2 76 Federico Marangoni high energetic costs of maintaining water balance or limiting availability of food resources; Marangoni, 2023), on adult body size, age and growth pattern (Marangoni et al. 2008, 2021). In addition, it could also be possible that gene flow would have prevented differentiation of these E. calamita populations for neutral loci, while intense selection would have maintained differences in adaptive traits of size and reproduction as has been suggested for other processes of adaptive divergence (Bensch et al., 1999; Gómez-Mestre and Tejedo, 2004). Recently, in the first study reporting genetic diversity estimates in the Iberian endemic pygmy newt (Triturus pygmaeus) from Doñana environment, Albert and García-Navas (2022) showed differences in genetic variability between temporary and permanent ponds. The authors suggested that the pond connectivity may constitute a more important factor than hydroperiod length in determining the genetic diversity and viability of pygmy newt populations. Moreover, given the geologically recent formation of the Guadalquivir basin, which may have occurred in mid-Holocene (about 5,000 years ago), it is possible that current genetic diversity patterns still reflect the historical distribution and gene flow among populations and that the effect of current gene flow (or lack thereof) is still not visible. In conclusion, the neutral genetic differentiation cannot explain the geographical variation in the body size of natterjack toad, in accordance with the previous study by Gomez-Mestre and Tejedo (2004). Therefore, it is suggested that future studies are needed on the interaction between the genetic component with the environmental factors, and life history traits (e.g., age and growth pattern, Marangoni et al., 2021; food resources, Marangoni, 2023) at both larval and juvenile stages. Finally, it is considered necessary to investigate other potential source of both withinand between-population components of variance, in addition to purely additive genetic variance, since that preliminary analyses showed that maternal effects may potentially contribute to the origin and evolution of the geographical variation in body size observed in E. calamita (Marangoni, unpublished data), and others amphibians (Bernardo, 1996; Mousseau and Fox, 1998; Räsänen et al., 2003, 2005). ACKNOWLEDGEMENTS I thank Fernando Campos for his invaluable help with the field wok. I thank Ivan Gomez-Mestre for his help in optimizing the molecular assays in the Laboratorio de Genética Molecular (LEM), Estación Biológica de Doñana (EBD-CSIC), and training me in microsatellite analysis. Dan Cogălniceanu, Marco Katzenberger and two anonymous reviewers provided helpful comments on an earlier draft of the manuscript. We thank Helder Duarte, a native speaker, for correcting the English draft of this manuscript. This work was supported by grant PB96– 0861 from Dirección General de Investigación Científica y Técnica conceded to M. Tejedo. Thanks also to the Consejería de Medio Ambiente de la Junta de Andalucía and the Reserva Biológica de Doñana, for providing the corresponding permits and facilities. All animal handling was conducted in accordance with the legal standards of Spain. REFERENCES Adams, D.C., Church, J.O. (2008): Amphibians do not follow Bergmann’s rule. Evolution. 62: 413-420. Albert, E.M., García-Navas, V. (2022): Population structure and genetic diversity of the threatened pygmy newt Triturus pygmaeus in a network of natural and artificial ponds. Conserv. Genet. 23: 575-588. Ashton, K.G. (2002): Do amphibians follow Bergmann’s rule? Can. J. Zool. 80: 708-716. Avise, J.C. (1994): Molecular markers, natural history and evolution. Kluwer Academic Publishers, Boston. Beebee, T.J.C., Rowe, G. (2000): Microsatellite analysis of natterjack toad Bufo calamita Laurenti populations: consequences of dispersal from a Pleistocene refugium. Biol. J. Linn. Soc. 69: 367-381. Belkhir, K., Borsa, P., Chikhi, L., Raufaste, N., Bonhomme, F. (2002): GENETIX 4.04, logiciel sous Windows TM pour la génétique des populations. Laboratoire Génome, Populations, Interactions, CNRS UMR 5000, Université de Montpellier II, Montpellier. Bensch, S., Andersson, T., Akesson, S. (1999): Morphological and molecular variation across a migratory divide in willow warblers, Phylloscopus trochilus. Evolution 53: 1925-1935. Bergmann, C. (1847): Über die Verhältnisse der Wärmeökonomie der Thiere zu ihrer Grösse. Göttinger Studien. 3: 595. Bernardo, J. (1996): Maternal effects in animal ecology. Am. Zool. 36: 83-105. Bradshaw, A.D. (1965): Evolutionary significance of phenotypic plasticity in plants. Adv. Genet. 13: 115-155. Cvetkovic, D., Tomasevic, N., Ficetola, G.F., CrnobrnjaIsailovic, J., Miaud, C. (2009): Bergmann’s rule in amphibians: combining demographic and ecological parameters to explain body size variation among populations in the common toad Bufo bufo. J. Zool. Syst. Evol. Res. 47: 171-180. 77 Neutral genetic differentiation in Epidalea calamita Endler, J.A. (1977): Geographic variation, speciation, and clines. Princeton University Pess, Princeton, New Jersey. Excoffier, J., Smouse, P., Quattro, J. (1992): Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics 136: 343359. Foster, S.A., Endler, J.A. (1999): Thoughts on geographic variation in behavior. In: Geographic variation in behaviour, pp. 287-307. Foster, S.A, Endler, J.A, Eds, Oxford University Press, New York, Oxford. Garant, D., Kruuk, L.E.B, Wilkin, T.A., McCleery, R.H., Sheldon, B.C. (2005): Evolution driven by differential dispersal within a wild bird population. Nature 433: 60-65. García-París, M., Alcobendas, M., Alberch, P. (1998): Influence of the Guadalquivir river basin on Mitocondrial DNA evolution of Salamandra salamandra (Caudata: Salamandridae) from Southern Spain. Copeia 1: 174-176. Gomez-Mestre, I. (2001): Variación geográfica y adaptación local al estrés osmótico en el sapo corredor, Bufo calamita. Unpublished doctoral dissertation. University of Seville, Spain. Gomez-Mestre, I., Tejedo, M. (2003). Local adaptation of an anuran amphibian to osmotically stressful environments. Evolution 57: 1889-1899. Gomez-Mestre I., Tejedo, M. (2004): Contrasting patterns of quantitative and neutral genetic variation in locally adapted populations of the natterjack toad Bufo calamita. Evolution. 58: 2343-2352. Goudet, J. (2001): FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3). Available at: http://www.unil.ch/izea/softwares/fstathtml. Harper, J. (1977): Population biology of plants. Academic Press, New York, USA. Hewitt, G.M. (1996): Some genetic consequences of ice ages, and their role in divergence and speciation. Biol. J. Linn. Soc. 58: 247-276. Jensen, J.L., Bohonak, A.J., Kelley, S.T. (2005): Isolation by distance, web service. BMC Genetics 6: 13. Kennington, W. J., Gockel, J., Partridge, L. (2002). Testing for asymmetrical gene flow in a Drosophila melanogaster body-size cline. Genetics. 165: 667-673. Knopp, T., Cano, J.M., Crochet, P.A., Merilä, J. (2007). Contrasting levels of variation in neutral and quantitative genetic loci on island populations of moor frogs (Rana arvalis). Conserv. Genet. 8: 45-56. Leblois, R., Rousset, F., Tikel, D., Moritz, C. (2000). Absence of evidence for isolation by distance in an expanding cane toad (Bufo marinus) population: an individual-based analysis of microsatellite genotypes. Mol. Ecol. 9: 1905-1909. Lee, H.J., Broggi, J., Sánchez-Montes, G., Díaz-Paniagua, C., Gomez-Mestre, I. (2020): Dwarfism in close continental amphibian populations despite lack of genetic isolation. Oikos. 129: 1243-1256. Leinonen, T., O’Hara, R.B., Cano, J.M,, Merilä, J. (2208). Comparative studies of quantitative trait and neutral marker divergence: a meta-analysis. J. Evol. Biol. 21: 1-17. Lenormand, T. (2002): Gene flow and the limits to natural selection. Trends Ecol. Evol. 17: 183-189. Lindsey, C.C. (1966): Body sizes of poikilotherm vertebrates at different latitudes. Evolution. 20: 456-465. López-Martínez, N. (1989): Tendencias en paleobiogeografia. El futuro de la biogeografía del pasado. In: Paleontologia, pp. 271-296. Aguirre, E., Coord, C.S.I.C., Madrid, Spain. Luquet, E., Léna, J.P., Miaud, C., Plénet, S. (2015). Phenotypic divergence of the common toad (Bufo bufo) along an altitudinal gradient: evidence for local adaptation. Heredity (Edinb) 114: 69-79. Marangoni, F. 2006. Variación clinal en el tamaño del cuerpo a escala microgeográfica en dos especies de anuros (Pelobates cultripes y Bufo calamita). Unpublished doctoral dissertation. University of Seville, Spain. Marangoni, F. (2023). Food Availability’s Role in the reduction in body size associated with sandy substrate in Pelobates cultripes. FACENA 33: In press. Marangoni, F., Tejedo, M. (2008): Variation in body size and metamorphic traits of Iberian spadefoot toads over a short geographic distance. J. Zool. 275: 97-105. Marangoni, F., Tejedo, M., Cogălniceanu, D. (2021): Can age and growth patterns explain the geographical variation in the body size of two toad species? An. Acad. Bras. Cienc. 93: e20190470. Marangoni, F., Tejedo, M., Gomez-Mestre, I. (2008): Extreme reduction in body size and reproductive output associated with sandy substrates in two anuran species. Amphibia-Reptilia 29: 541-553. Mayr, E. (1963): Animal Species and Evolution. Harvard University Press. Cambridge, UK. Mousseau, T.E., Fox, C.W. (1998). Maternal effects as adaptations. Oxford University Press, New York. Páez-Vacas, M.I., Trumbo, D.R, Funk, W.C. (2022) Contrasting environmental drivers of genetic and phenotypic divergence in an Andean poison frog (Epipedobates anthonyi). Heredity (Edinb) 128: 33-44. Postma, E., Noordwijk van, A.J. (2005): Gene flow a large genetic difference in clutch size at a small spatial scale. Nature 433: 65-68.