Environmental correlates of the European common toad hybrid zone

The interplay between intrinsic (development, physiology, behavior) and extrinsic (landscape features, climate) factors determines the outcome of admixture processes in hybrid zones, in a continuum from complete genetic merger to full reproductive isolation. Here we assess the role of environmental correlates in shaping admixture patterns in the long hybrid zone formed by the toads Bufo bufo and B. spinosus in western Europe. We used species-specific diagnostic SNP markers to genotype 6584 individuals from 514 localities to describe the contact zone and tested for association with topographic, bioclimatic and land use variables. Variables related to temperature and precipitation contributed to accurately predict the distribution of pure populations of each species, but the models did not perform well in areas where genetically admixed populations occur. A sliding window approach proved useful to identify different sets of variables that are important in different sections of this long and heterogeneous hybrid zone, and offers good potential to predict the fate of moving contact zones in global change scenarios.


Introduction
The study of hybrid zones, or areas where individuals from two genetically differentiated population lineages meet and produce hybrid offspring, is critical to our understanding of how new species form and remain distinct through evolutionary time (Barton & Hewitt, 1985;Harrison & Larson, 2014). The outcome of this admixture process ranges from complete genetic merger to full reproductive isolation with strong selection against hybrid individuals, and is a result of the interaction between intrinsic and extrinsic factors (Mallet, 2005;Abbott et al., 2013). Among the former are developmental, physiological and behavioral traits contributing to reproductive isolation, whereas extrinsic factors include those associated with the spatial and climatic context of the hybrid zone.
Hybrid zones often extend for hundreds of kilometers across heterogeneous landscapes, with features like mountains, rivers, soil types or roads differentially restricting dispersal and, consequentially, the degree of admixture in alternate sections of the contact zone. Similarly, spatially varying climatic conditions along the hybrid zone potentially impose different selective pressures on parental species and their hybrids in different transects. While topography is generally more conserved through time, climatic changes operate at temporal scales that are relevant for the study of hybrid zone dynamics and have been invoked to explain patterns of asymmetrical introgression in moving hybrid zones. The study of the complex interplay of topographic and climatic factors varying over time and space in hybrid zones can thus illuminate the relative role of extrinsic factors in reproductive isolation, explain discordant genetic patterns in different sections of a hybrid zone, and help predict the fate of hybrid zones (i.e., the stability of reproductive isolation through time) in the face of climatic changes (McQuillan & Rice, 2015;Taylor et al., 2015;Hunter et al., 2017;Ryan et al., 2018).
The European common toad hybrid zone is an emerging model system in the study of speciation.
The common (Bufo bufo) and spined (Bufo spinosus) toads meet along a ca. 900 km hybrid zone from the Atlantic coast of northwest France to the Mediterranean coast of southeast France and northwest Italy (Recuero et al., 2012;Arntzen et al., 2018). This hybrid zone has been previously characterized with mitochondrial and nuclear DNA markers, which have revealed different patterns of genetic admixture across the hybrid zone, with transects in northwestern France resembling a classical tension zone model with selection against hybrids and those in the southeast of France and northwest of Italy characterized by strong cyto-nuclear discordance, probably as a result of spatial displacement of the hybrid zone in the past due to climatic changes (Arntzen et al., 2016(Arntzen et al., , 2017van Riemsdijk et al., 2019ab).
While previous studies have discussed the potential role of intrinsic factors in maintaining reproductive isolation across the hybrid zone, the relative role of environmental correlates has not been assessed so far, with the exception of certain topographic features reported to be associated with genetic transitions in different section of the contact zone (Arntzen et al., 2017(Arntzen et al., , 2018. Here we use environmental data and a comprehensive genetic dataset to investigate the role of extrinsic factors in maintaining species borders through space and time.

Materials and methods
A total of 6584 tissue samples was obtained from adult and larval toads from 514 localities across western Europe, with the emphasis on the broad region of species contact across France and Italy. Sampling was somewhat thin at the Atlantic coast and it was unsuccessful over the lower cols of the French Alps (col de Montgenèvre at 1854 m a.s.l., colle della Maddalena at 1996 m, col de Tende at 1870 m). All individuals were studied for the four species-diagnostic SNP nuclear markers bdnf, pomc, rag1 and rpl3 (Arntzen et al., 2016). Individuals with data missing for more than one marker were excluded. Sample size ranged from 1-358 with an average of 12.8. A total of 610 individuals was studied in duplicate. Results were not the same in two cases (0.08%) and involved the score of '22' (i.e. homozygous for the spinosus-allele) versus the heterozygous '12' condition. Both observations were subsequently treated as unavailable. The total of missing data was 2.0%. A subset of 2524 individuals from 185 localities was studied for another 27 diagnostic SNP markers (as in van Riemsdijk et al., 2019b;Arntzen, 2019) to determine the center of the two species contact zone more accurately. These latter localities were arranged in eight transects with hybrid zone centers located in France or Italy. For locality information, sample sizes and molecular species identifications see supplementary information I [not provided]. Mitochondrial DNA and toad morphology were not studied because for these characters species diagnostic performances break down over parts of the species' parapatric range border (Arntzen et al., 2017(Arntzen et al., , 2018. The genetic data were subjected to Bayesian clustering with Structure software (Pritchard et al., 2000) with a predetermined K-value of two, corresponding to the number of species involved. The position of the mutual range border of B. bufo and B. spinosus was estimated by spatial interpolation of the Structure Q-values, to obtain the Q=0.5 isoline. In this procedure data were weighted by sample size and by the number of markers employed. We used Dirichlet cells with ILWIS (ILWIS, 2009) and linear interpolation with MyStat (Systat, 2007), which methods implied a weak and a strong spatial smoothing, respectively.
Environmental data considered include altitude and 19 climate variables from the WorldClim global climate database v2, available at http://www.worldclim.org (Hijmans et al., 2005; see also Fick & Hijmans, 2017). The parameter 'slope' was derived from altitude through a set of filter operations. Soil properties data are from the ESDA European soil Database v2.0, available at https://esdac.jrc.ec.europa.eu (see also Panagos et al., 2012) and were used as far as parameter values could be ordered. Vegetation and land use data were from the CORINE land cover database of the European Environment Agency, available at https://www.eea.europa.eu/publications/COR0-landcover. Data were grouped in the four classes 'crop growing', 'forestation', 'pasture' and 'other' (including e.g. wetlands, urban, industrial and the transport network).
To identify and subsequently reduce colinearity among the environmental variables we constructed the half-matrix of their pairwise Pearson correlation coefficients. This matrix was subjected to clustering with UPGMA. Variables were retained using criteria of partial independence at r<0.7 and selected in alphanumeric order (supplemental information II) [not provided]. The 20 variables considered for the construction of species distribution models are listed in table 1.
We produced two-species distribution models, in which both species' environmental attributes are contrasted under the assumption of a distribution in parapatry, as follows. First, we sampled environmental data for the localities of the genetically investigated toad populations within the area -5 to 12 degrees eastern longitude and 42 to 53 degrees northern latitude. Populations with Q<0.2 that represent B. bufo and Q>0.8 that represent B. spinosus were about equally frequent (B. bufo N=196, B. spinosus N=208). Data for 47 populations with intermediate Q-values were analyzed separately. Second, environmental data were sampled in a string of 17 overlapping and adjoining circular windows, positioned over the species contact zone, as defined by the Q=0.5 isoline (figure 1). A window diameter of 100 km was chosen to widely encompass the ca. 50 km wide B. bufo -B. spinosus hybrid zone (Arntzen et al., 2016(Arntzen et al., , 2017(Arntzen et al., , 2018van Riemsdijk et al., 2019b). In consideration of scale, the string of windows followed the smooth curved species delineation (and not the more angular one; see below). The analyses involved 100 randomly chosen data points per window. The sampling localities were interpreted as pseudo-presences of B. bufo and B. spinosus if at the northeastern or at the southwestern side of the inferred species border, respectively.

Results
The nuclear genetic data confirm that the hybrid zone of the common toad, B. bufo and the spined toad, B.

Discussion
We used two approaches to assess the role of environmental variables in shaping the common-spined toad hybrid zone. While both performed relatively well at predicting where pure populations will be, the global model performed poorly when genetically admixed populations were included. This may reflect the heterogeneity of the hybrid zone, where the spatial context of successive transects varies greatly in association with topography. Our results also highlight the sliding window approach as most useful for dealing with long, complex hybrid zones, because different variables will have different impact on hybrid zone structure in different sections of the zone, and these are highlighted as such.
Among selected environmental variables, most were related to temperature (bio01, bio03, bio06, bio08, bio09) or precipitation (bio12), but soil type also entered the global model. While the former may generally be associated with a differentiation between a northern (cooler, more humid) and a southern (warmer, drier and more extreme) species, the role of soil type is not obvious, and may have an indirect association with the distribution of the two toad species (e.g. through larval or postmetamorphic associations with pH and vegetal communities). However, land use categories did not make it into the models, but perhaps they operate at a smaller spatial scale and have little or no bearing on the wider position of the hybrid zone.
High model fit was also associated with higher altitudes for B. spinosus than for B. bufo. This observation is not compatible with the climate preferences alluded to in the above, with as a further complication that the critical altitude increases with an increasingly continental climate. Low model fit tended to be associated with less populated areas and with rivers. These are areas where dispersal is reduced and where hybrid zones are predicted to settle (Barton & Hewitt, 1985). The highest records for B. bufo are from the Swiss Alps at 2300 m (Malkmus & Grossenbacher, 2013) where the species is expanding its elevation range (Luscher et al., 2016). The highest records for B. spinosus are from the Pyrenees at 2600 m with, however, low population densities above 1500 m (Balcells, 1975;Ortiz Santaliestra, 2014). We were unable to sample from France to Italy and it is currently unclear to what degree the species are (or were) in contact, other than along the Mediterranean coast.
Our study shows that rivers can have a larger impact than climatic/environmental variables, suggesting that the presence of a river can break down environmental associations with toad presence.
The low numerical model fit at rivers may simply reflect the fact that these linear structures were not incorporated in our two-dimensional modelling approach. Historical factors can also overrule the importance of environmental variables, for instance when different regions involve different intraspecific lineages, such as the northern Balkan lineage of B. bufo and its southern, Apennine counterpart (Garcia-Porta et al., 2012;Arntzen et al., 2017). The failure of the model to accurately predict the species present in Normandy and Brittany may be a 'preoccupancy effects', where a species is present in an area that is more environmentally favourable for the other species just because it has been long part of its range.
Alternatively, B. bufo might have occupied these areas but did not make it because the colonization route was inaccessible.
Predicting the fate of moving hybrid zones is challenging because of the complex interaction between topography and environmental variables, but our sliding window approach, coupled with phylogeographic analyses, would seem to offer great potential to make predictions on species advances and retreats in a global change scenario.  Center.

Figure 2
Two-species global distribution model for Bufo toads in western Europe. The colour scale runs from deep red for B. spinosus (P b is zero) to deep blue for B. bufo (P b at unity). The solid black line represents the center of the species' hybrid zone from molecular data, as in figure 1. Outlined circular windows are those for which model fit is less than good (AUC<0.8, windows 1, 6-8 and 14-16). France where the species border coincides with the Rhône and the lower Isère river at window 13. Note the paucity of data for the high Alps at windows 15 and 16 (see also Lescure and de Massary, 2012;Arntzen et al., 2017).

Figure 4
Average values for eight selected environmental variables over 17 windows that follow the Bufo bufo -B.
spinosus hybrid zone from the Atlantic coast (window 1) to the Mediterranean (window 17). Units are as in table 1; see also Hijmans et al., 2005 signals are likely to be absent or void, either from poor sampling (window 1), the presence of rivers (windows 5-9, 13-14), or a thin or absent species' contact (windows 16-17) (see figure 3).