# Low neutral and immunogenetic diversity in northern fringe populations of the green toad Bufotes viridis: implications for conservation

Jacob Höglund, Lucas Bolender, Maria Cortazar-Chinarro, Sara Meurling, Anssi Laurila, Adam Hermaniuk, Christophe Dufresnes

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Jacob Höglund, Lucas Bolender, Maria Cortazar-Chinarro, Sara Meurling, Anssi Laurila, et al.. Low neutral and immunogenetic diversity in northern fringe populations of the green toad Bufotes viridis: implications for conservation. Conservation Genetics, 2022, 23 (1), pp.139-149. 10.1007/s10592-021-01407-5. hal-04549764

HAL Id: hal-04549764

https://hal.science/hal-04549764

Submitted on 17 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Low neutral and immunogenetic diversity in northern fringe populations of the green toad Bufotes viridis: implications for conservation

Jacob Höglund, Lucas Bolender, Maria Cortazar-Chinarro, Sara Meurling,

Anssi Laurila, Adam Hermaniuk, Christophe Dufresnes

► To cite this version:

Jacob Höglund, Lucas Bolender, Maria Cortazar-Chinarro, Sara Meurling, Anssi Laurila, et al. Low neutral and immunogenetic diversity in northern fringe populations of the green toad Bufotes viridis: implications for conservation. Conservation Genetics, 2022, 23 (1), pp.139-149. 10.1007/s10592-021-01407-5. hal-04549764

HAL Id: hal-04549764

https://hal.science/hal-04549764v1

Submitted on 17 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Low neutral and immunogenetic diversity in northern fringe populations of the green toad *Bufotes viridis*: implications for conservation

Jacob Höglund<fnref n="1" /> · Lucas Bolender<fnref n="1" />,<fnref n="4" /> · Maria Cortazar-Chinarro<fnref n="1" /> · Sara Meurling<fnref n="1" /> · Anssi Laurila<fnref n="1" /> · Adam Hermaniuk<fnref n="2" /> · Christophe Dufresnes<fnref n="3" />

Received: 5 October 2020 / Accepted: 5 October 2021 / Published online: 23 October 2021
© The Author(s) 2021

## Abstract

Genetic variation is often lower at high latitudes, which may compromise the adaptability and hence survival of organisms. Here we show that genetic variability is negatively correlated with northern latitude in European green toads (*Bufotes viridis*). The result holds true for both putatively neutral microsatellite variation and supposedly adaptive MHC Class IIB variation. In particular, our findings have bearing on the conservation status of this species in Sweden, on the northern limit of its distribution where local populations are small and fragmented. These genetically impoverished populations are closely related to other populations found around the Baltic Sea basin. The low neutral and adaptive variation in these fringe populations compared to population at central ranges confirms a pattern shared across all other amphibians so far studied. In Sweden, the situation of green toads is of concern as the remaining populations may not have the evolutionary potential to cope with present and future environmental challenges.

**Keywords** Adaptive diversity · Neutral diversity · Green toad · Clinical variation · Fringe population

## Introduction

A general trend in latitudinal studies of genetic variation is the loss of genetic variability towards the poles (e.g. Adams and Hadly 2012; Miraldo et al. 2016; Smith et al. 2017). In the northern hemisphere, populations of many species become smaller and more fragmented with increasing latitude (e.g. Wan et al. 2018). This observation is explained by harsher environmental conditions and increasing seasonality

with longer and harder winters towards the north (Vellend 2005; Vellend and Geber 2005). Furthermore, climatic conditions have oscillated during recent evolutionary times. During the Quaternary, large parts of the northern hemisphere were completely covered by glaciers and thus uninhabitable for any form of life (Hewitt 1996, 1999). Since the retreat of the glaciers after the last glacial maximum (LGM) approximately 20.000 years ago, northern species have expanded across the formerly uninhabitable areas. While doing so populations have lost variation during the recolonization process (Hewitt 2000, 2004).

This loss of genetic diversity towards the poles has consequences for the conservation of contemporary populations and species. Since adaptation to changing conditions requires genetic variation, small and fragmented populations are predicted to be more vulnerable to local extinction (McMahon et al. 2014). This prediction is borne out in comparative studies of threatened and non-threatened species (Spielman et al. 2004; Li et al. 2016).

Genetic diversity can be assayed in various ways (Horegger et al. 2006; Rödin-Mörch et al. 2019). On the one hand, molecular surveys of natural populations may consider putatively neutral loci, i.e. loci evolving mainly through the

<figure><img src="image_3.png" /></figure>

<footnote>J Jacob Höglund
jacob.hoglund@ebc.uu.se</footnote>

<footnote n="1">¹ Department of Ecology and Genetics, Uppsala University, Norbyvägen 18D, 75236 Uppsala, Sweden</footnote>

<footnote n="2">² Department of Evolutionary and Physiological Ecology, Faculty of Biology, University of Białystok, Ciołkowskiego 1J, 15-245 Białystok, Poland</footnote>

<footnote n="3">³ LASER, College of Biology and the Environment, Nanjing Forestry University, Nanjing, People's Republic of China</footnote>

<footnote n="4">⁴ Present Address: Advanced Research and Diagnostic Laboratory, Department of Laboratory Medicine and Pathology, University of Minnesota, 1200 Washington Avenue S., Minneapolis, MN 55415, USA</footnote>
# Low neutral and immunogenetic diversity in northern fringe populations of the green toad *Bufotes viridis*: implications for conservation

Jacob Höglund<fnref n="1" /> · Lucas Bolender<fnref n="1" />,<fnref n="4" /> · Maria Cortazar-Chinarro<fnref n="1" /> · Sara Meurling<fnref n="1" /> · Anssi Laurila<fnref n="1" /> · Adam Hermaniuk<fnref n="2" /> · Christophe Dufresnes<fnref n="3" />

Received: 5 October 2020 / Accepted: 5 October 2021 / Published online: 23 October 2021
© The Author(s) 2021

## Abstract

Genetic variation is often lower at high latitudes, which may compromise the adaptability and hence survival of organisms. Here we show that genetic variability is negatively correlated with northern latitude in European green toads (*Bufotes viridis*). The result holds true for both putatively neutral microsatellite variation and supposedly adaptive MHC Class IIB variation. In particular, our findings have bearing on the conservation status of this species in Sweden, on the northern limit of its distribution where local populations are small and fragmented. These genetically impoverished populations are closely related to other populations found around the Baltic Sea basin. The low neutral and adaptive variation in these fringe populations compared to population at central ranges confirms a pattern shared across all other amphibians so far studied. In Sweden, the situation of green toads is of concern as the remaining populations may not have the evolutionary potential to cope with present and future environmental challenges.

**Keywords** Adaptive diversity · Neutral diversity · Green toad · Clinical variation · Fringe population

## Introduction

A general trend in latitudinal studies of genetic variation is the loss of genetic variability towards the poles (e.g. Adams and Hadly 2012; Miraldo et al. 2016; Smith et al. 2017). In the northern hemisphere, populations of many species become smaller and more fragmented with increasing latitude (e.g. Wan et al. 2018). This observation is explained by harsher environmental conditions and increasing seasonality

with longer and harder winters towards the north (Vellend 2005; Vellend and Geber 2005). Furthermore, climatic conditions have oscillated during recent evolutionary times. During the Quaternary, large parts of the northern hemisphere were completely covered by glaciers and thus uninhabitable for any form of life (Hewitt 1996, 1999). Since the retreat of the glaciers after the last glacial maximum (LGM) approximately 20.000 years ago, northern species have expanded across the formerly uninhabitable areas. While doing so populations have lost variation during the recolonization process (Hewitt 2000, 2004).

This loss of genetic diversity towards the poles has consequences for the conservation of contemporary populations and species. Since adaptation to changing conditions requires genetic variation, small and fragmented populations are predicted to be more vulnerable to local extinction (McMahon et al. 2014). This prediction is borne out in comparative studies of threatened and non-threatened species (Spielman et al. 2004; Li et al. 2016).

Genetic diversity can be assayed in various ways (Holderegger et al. 2006; Rödin-Mörch et al. 2019). On the one hand, molecular surveys of natural populations may consider putatively neutral loci, i.e. loci evolving mainly through the

<figure><img src="image_2.png" /></figure>

<footnote>J Jacob Höglund
jacob.hoglund@ebc.uu.se</footnote>

<footnote n="1">¹ Department of Ecology and Genetics, Uppsala University, Norbyvägen 18D, 75236 Uppsala, Sweden</footnote>

<footnote n="2">² Department of Evolutionary and Physiological Ecology, Faculty of Biology, University of Białystok, Ciołkowskiego 1J, 15-245 Białystok, Poland</footnote>

<footnote n="3">³ LASER, College of Biology and the Environment, Nanjing Forestry University, Nanjing, People's Republic of China</footnote>

<footnote n="4">⁴ Present Address: Advanced Research and Diagnostic Laboratory, Department of Laboratory Medicine and Pathology, University of Minnesota, 1200 Washington Avenue S., Minneapolis, MN 55415, USA</footnote>
effects of the stochastic process of genetic drift, such as the mitochondrial control region (d-loop) and nuclear short tandem repeats (referred to as microsatellites). On the other hand, direct assessments of the evolutionary potential of species should rely on loci also subjected to the effects of natural selection, such as the vertebrate major histocompatibility complex (MHC). The MHC encompasses a series of genes encoding cellular proteins involved in immune response, and features highly polymorphic loci, in particular the class II B genes, that serve to assist in the presentation of pathogens. As such, they may be recognized by T cells and trigger the production of antibodies (Sommer 2005; Piertney and Oliver 2006). With higher immunogenetic variation, populations and individuals are consequently better protected against diseases and parasites. For example, resistance to the fungal pathogen *Batrachochytrium dendrobatidis* has been found to be associated with amino acid residue positions that are located within MHC class II (Fu and Waldman 2017).

The European green toad (*Bufotes viridis*) is a widespread anuran amphibian distributed in central and Eastern Europe, reaching southern Sweden and Estonia in its northern fringe (Sillero et al. 2014). In Sweden the species has been listed as one of the most vulnerable amphibians in the country (https://artfakta.se/naturvard/taxon/100022, accessed 1/9 2020). The Swedish distribution is confined to the southernmost provinces, where only a few fragmented populations remain, geographically isolated from the European mainland. The potentially reduced gene pool of these northern isolates creates a significant problem for the conservation of

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Occurrence of the *viridis* (green dots) and *sitibundus* (blue dots) mitochondrial lineages in Central and Eastern Europe, and surrounding regions. The actual range of our target species *B. viridis* is shown in light green background (data from Dufresnes et al. 2019)</figcaption></figure>

*B. viridis* in Sweden, as lack of adaptive potential and lost parasite resistance could prevent survival and reproduction (Höglund 2009).

In the framework of contemporary conservation efforts, the aim of this study was to understand the geographic distribution of putatively neutral vs adaptive variation in *B. viridis*, notably to test whether the fragmented, isolated and small Sweden populations are genetically impoverished compared to the more abundant mainland populations, as expected under the biogeographic paradigms. To this aim, we specifically targeted putatively neutral markers (mitochondrial control region and nuclear microsatellites) and loci that show signals of being subjected to local environmental selection (MHC Class IIB loci).

# Methods

## Mitochondrial barcoding and microsatellite genotyping

One mitochondrial and seven microsatellite markers were genotyped from samples taken all over the range.

For mtDNA, we amplified ~900 bp of the hypervariable mitochondrial control region (CR, also named d-loop in previous work) in 340 new individuals (69 localities, Table S1). PCRs were carried out in 25 µL reactions containing 3 µL of template DNA, 12.5 µL of nanopure water, 7.5 µL of multiplex master mix (Qiagen, containing buffer, dNTPs and hot-start polymerase), 1 µL of primer *CytbA-L* (5'-GAA TYGGRGGWCAACCAGTAGAAGACCC-3'; 10 µM), 1 µL of primer *ControlB-H* (5'-GTCCATTGGAGGTTAAGATC TACCA-3'; 10 µM). The following thermocycling conditions were applied: 95 °C for 15', 35 cycles of 94 °C for 30", 53 °C for 45" and 72 °C for 1'; 72 °C for 5'. Amplicons were Sanger-sequenced in a single direction (*CytbA-L*) and sequences were assigned to one of the two main mitochondrial lineages found in *B. viridis*: the *viridis* lineage and the *sitibundus* lineage, the latter having spread in Europe through past mtDNA introgression by the Anatolian taxon *B. sitibundus* (Dufresnes et al. 2019, 2021). To visualize lineage distributions, our new data was combined to the mtDNA dataset recently compiled for all major *Bufotes* lineages by Dufresnes et al. (2019), and which encompasses 162 localities where the *viridis* and/or *sitibundus* mtDNA occur.

The seven microsatellite markers were D106, Bcalµ10, C203, C218, C205, C223, D105, all previously analyzed in southeastern Europe (Dufresnes et al. 2018), and genotyped here in a similar fashion. Briefly, this included amplification in multiplexed PCRs, genotyping on an ABI3130 genetic analyzer and allele scoring with Genemapper 4.0 (Applied Biosystems); details are available in the supplementary

<figure><img src="image_5.png" /></figure>
effects of the stochastic process of genetic drift, such as the mitochondrial control region (d-loop) and nuclear short tandem repeats (referred to as microsatellites). On the other hand, direct assessments of the evolutionary potential of species should rely on loci also subjected to the effects of natural selection, such as the vertebrate major histocompatibility complex (MHC). The MHC encompasses a series of genes encoding cellular proteins involved in immune response, and features highly polymorphic loci, in particular the class II B genes, that serve to assist in the presentation of pathogens. As such, they may be recognized by T cells and trigger the production of antibodies (Sommer 2005; Piertney and Oliver 2006). With higher immunogenetic variation, populations and individuals are consequently better protected against diseases and parasites. For example, resistance to the fungal pathogen *Batrachochytrium dendrobatidis* has been found to be associated with amino acid residue positions that are located within MHC class II (Fu and Waldman 2017).

The European green toad (*Bufotes viridis*) is a widespread anuran amphibian distributed in central and Eastern Europe, reaching southern Sweden and Estonia in its northern fringe (Sillero et al. 2014). In Sweden the species has been listed as one of the most vulnerable amphibians in the country (https://artfakta.se/naturvard/taxon/100022, accessed 1/9 2020). The Swedish distribution is confined to the southernmost provinces, where only a few fragmented populations remain, geographically isolated from the European mainland. The potentially reduced gene pool of these northern isolates creates a significant problem for the conservation of

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Occurrence of the *viridis* (green dots) and *sitibundus* (blue dots) mitochondrial lineages in Central and Eastern Europe, and surrounding regions. The actual range of our target species *B. viridis* is shown in light green background (data from Dufresnes et al. 2019)</figcaption></figure>

*B. viridis* in Sweden, as lack of adaptive potential and lost parasite resistance could prevent survival and reproduction (Höglund 2009).

In the framework of contemporary conservation efforts, the aim of this study was to understand the geographic distribution of putatively neutral vs adaptive variation in *B. viridis*, notably to test whether the fragmented, isolated and small Sweden populations are genetically impoverished compared to the more abundant mainland populations, as expected under the biogeographic paradigms. To this aim, we specifically targeted putatively neutral markers (mitochondrial control region and nuclear microsatellites) and loci that show signals of being subjected to local environmental selection (MHC Class IIB loci).

# Methods

## Mitochondrial barcoding and microsatellite genotyping

One mitochondrial and seven microsatellite markers were genotyped from samples taken all over the range.

For mtDNA, we amplified ~900 bp of the hypervariable mitochondrial control region (CR, also named d-loop in previous work) in 340 new individuals (69 localities, Table S1). PCRs were carried out in 25 µL reactions containing 3 µL of template DNA, 12.5 µL of nanopure water, 7.5 µL of multiplex master mix (Qiagen, containing buffer, dNTPs and hot-start polymerase), 1 µL of primer CytbA-L (5'-GAA TYGGRGGWCAACCAGTAGAAGACCC-3'; 10 µM), 1 µL of primer ControlB-H (5'-GTCCATTGGAGGTTAAGATC TACCA-3'; 10 µM). The following thermocycling conditions were applied: 95 °C for 15', 35 cycles of 94 °C for 30", 53 °C for 45" and 72 °C for 1'; 72 °C for 5'. Amplicons were Sanger-sequenced in a single direction (CytbA-L) and sequences were assigned to one of the two main mitochondrial lineages found in *B. viridis*: the *viridis* lineage and the *sitibundus* lineage, the latter having spread in Europe through past mtDNA introgression by the Anatolian taxon *B. sitibundus* (Dufresnes et al. 2019, 2021). To visualize lineage distributions, our new data was combined to the mtDNA dataset recently compiled for all major *Bufotes* lineages by Dufresnes et al. (2019), and which encompasses 162 localities where the *viridis* and/or *sitibundus* mtDNA occur.

The seven microsatellite markers were D106, Bcalµ10, C203, C218, C205, C223, D105, all previously analyzed in southeastern Europe (Dufresnes et al. 2018), and genotyped here in a similar fashion. Briefly, this included amplification in multiplexed PCRs, genotyping on an ABI3130 genetic analyzer and allele scoring with Genemapper 4.0 (Applied Biosystems); details are available in the supplementary

<figure><img src="image_4.png" /></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 2** **a** PCA on individual allele frequencies for 776 microsatellite genotypes of *B. viridis* sampled across Europe. Colors distinguish the main regions associated with genetic structure: Balkans/E-Europe (green: Albania, Greece, N-Macedonia, Romania, Russia, Serbia, Ukraine, Belarus, Latvia), northern ranges spanning around the Baltic Sea (purple Denmark, N-Germany, Poland, Sweden) and the isolated population from W-Germany (orange). **b** Bayesian clustering analyses of the microsatellite genotypes into 2 to 5 groups, with the population average ancestries for 5 groups displayed on the map</figcaption></figure>

information of Dufresnes et al. (2018). A total of 336 individuals (65 localities, Table S1) were genotyped with the same panels, bins and operator (CD) as in Dufresnes et al. (2018). For downstream analyses, we could thus merge our new data with their genotype matrix of the SE-European populations of *B. viridis* (their loc. 20–106), for a total of 776 individuals (134 localities).

**Population genetic analyses of mitochondrial and microsatellite data**

Departure from neutrality, possible linkage among loci, genetic structure and diversity were inferred from the combined microsatellite data (7 loci, 776 individuals) by individual and population-based analyses. Neutrality was tested for by examining departures from Hardy–Weinberg assumptions and linkage was examined by the probability of random associations among alleles at the different loci using

the software Fstat 2.9.3 (Goudet 1995). For genetic structure, we first performed a Principal Component Analysis (PCA) on individual allele frequencies with the R package *adegenet* (Jombart 2008). Second, we applied the Bayesian clustering algorithm of STRUCTURE (Pritchard et al. 2000) to assign individuals to genetic groups. The admixture model with correlated allele frequencies was used, with runs of 100,000 MCMC iterations after a burn in of 10,000. We tested K from 1 to 10 groups, each with ten replicates. Likelihood estimates between runs and between K models (Ln(PlK) and ΔK statistics) were summarized and visualized in STRUCTURE Harvester (Earl and von Holdt 2012), runs were combined with CLUMPP (Jakobsson and Rosenberg 2007), and graphical results were produced with DISTRUCT (Rosenberg 2004). Third, for localities with *n* ≥ 5, we computed population differentiation (pairwise Fst) and within-population diversity indices—observed (H_o) and expected
<figure><img src="image_5.png" /><figcaption>**Fig. 2** **a** PCA on individual allele frequencies for 776 microsatellite genotypes of *B. viridis* sampled across Europe. Colors distinguish the main regions associated with genetic structure: Balkans/E-Europe (green: Albania, Greece, N-Macedonia, Romania, Russia, Serbia, Ukraine, Belarus, Latvia), northern ranges spanning around the Baltic Sea (purple Denmark, N-Germany, Poland, Sweden) and the isolated population from W-Germany (orange). **b** Bayesian clustering analyses of the microsatellite genotypes into 2 to 5 groups, with the population average ancestries for 5 groups displayed on the map</figcaption></figure>

information of Dufresnes et al. (2018). A total of 336 individuals (65 localities, Table S1) were genotyped with the same panels, bins and operator (CD) as in Dufresnes et al. (2018). For downstream analyses, we could thus merge our new data with their genotype matrix of the SE-European populations of *B. viridis* (their loc. 20–106), for a total of 776 individuals (134 localities).

## Population genetic analyses of mitochondrial and microsatellite data

Departure from neutrality, possible linkage among loci, genetic structure and diversity were inferred from the combined microsatellite data (7 loci, 776 individuals) by individual and population-based analyses. Neutrality was tested for by examining departures from Hardy–Weinberg assumptions and linkage was examined by the probability of random associations among alleles at the different loci using

the software Fstat 2.9.3 (Goudet 1995). For genetic structure, we first performed a Principal Component Analysis (PCA) on individual allele frequencies with the R package *adegenet* (Jombart 2008). Second, we applied the Bayesian clustering algorithm of STRUCTURE (Pritchard et al. 2000) to assign individuals to genetic groups. The admixture model with correlated allele frequencies was used, with runs of 100,000 MCMC iterations after a burn in of 10,000. We tested K from 1 to 10 groups, each with ten replicates. Likelihood estimates between runs and between K models (Ln(PlK) and ΔK statistics) were summarized and visualized in STRUCTURE Harvester (Earl and von Holdt 2012), runs were combined with CLUMPP (Jakobsson and Rosenberg 2007), and graphical results were produced with DISTRUCT (Rosenberg 2004). Third, for localities with *n* ≥ 5, we computed population differentiation (pairwise Fst) and within-population diversity indices—observed (H_o) and expected
<figure><img src="image_7.png" /><figcaption>**Fig. 3** *a* Within-population genetic diversity, as estimated from the microsatellite dataset for localities of at least 5 genotyped toads. The variation in $H_o$ and $H_e$ estimates are illustrated by color gradients. *b* Relationships between within-population genetic diversity and latitude for these populations. Both are significant (see "Results")</figcaption></figure>

heterozygosity ($H_e$)—using the R package *hierfstat* (Goudet 2005).

## MHC Class II genotyping

A total of 430 samples were independently collected from Belarus, Bulgaria, Denmark, Germany, Greece, Latvia, Moldova, Poland, Romania, Russia, Sweden, and Ukraine

(Table S2). Within Sweden, we obtained 272 samples. DNA was extracted with a modified high salt extraction process (Paxton et al. 1996). DNA quality and concentration were then measured for each sample using electrophoresis with 1.5% agarose and a Thermo Scientific Nanodrop 2000 spectrophotometer.

We used the forward primer 2F347F (GTGACCCTC TGCTCTCCATT) and reverse primer 2R307R (ATAATT

<figure><img src="image_8.png" /></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 3** **a** Within-population genetic diversity, as estimated from the microsatellite dataset for localities of at least 5 genotyped toads. The variation in $H_o$ and $H_e$ estimates are illustrated by color gradients. **b**</figcaption></figure>

Relationships between within-population genetic diversity and latitude for these populations. Both are significant (see “Results”)

heterozygosity ($H_e$)—using the R package *hierfstat* (Goudet 2005).

## MHC Class II genotyping

A total of 430 samples were independently collected from Belarus, Bulgaria, Denmark, Germany, Greece, Latvia, Moldova, Poland, Romania, Russia, Sweden, and Ukraine

(Table S2). Within Sweden, we obtained 272 samples. DNA was extracted with a modified high salt extraction process (Paxton et al. 1996). DNA quality and concentration were then measured for each sample using electrophoresis with 1.5% agarose and a Thermo Scientific Nanodrop 2000 spectrophotometer.

We used the forward primer 2F347F (GTGACCCTC TGCTCTCCATT) and reverse primer 2R307R (ATAATT

<figure><img src="image_7.png" /></figure>
CAGTATATACAGGGTCTCACC) to target and amplify
sequences of 279-282 base pairs of the MHC Class II B
gene as described in Zeisset and Beebee (2014). Samples
were randomly placed into pool groups that determined the
forward and reverse primer tag combinations which attached
to the primers (e.g. NNN-ACTAGATC) for each sample, so
that individual samples could be identified during analysis.
The PCR mixtures were made using 13.9 µl double distilled
water, 4 µl Phusion High Fidelity Buffer, 0.4 µl dNTPS,
0.5 µl each of the required forward and reverse primers,
0.2 µl Phusion High Fidelity DNA Taq polymerase, and
0.5 µl of the extracted DNA. The PCRs were run, and the
results of each sample were measured on a 1.5% agarose gel
after electrophoresis. We included 20% of replicated sam-
ples. Samples were grouped by the strength of the resulting
bands, separating the pools into strong and weak clusters.
Bands within the pre-pools were purified by using a Qia-
gen® MinElute gel extraction procedure (Qiagen, Sollen-
tuna), following the manufacturer's protocol. The pre-pools
were measured for DNA concentration with the Qbit and
consecutively pooled in equimolar amounts. Genetic librar-
ies from every pre-pool (A, B, C and D) were prepared using
ThruPLEX DNA-seq Kit (Rubicon Genomics) following the
manufactures protocol. Sequencing was completed using
Illumina MiSeq at the National Genomic Infrastructure,
SciLifeLab in Uppsala.

MHC data analyses

For MHC variation, final number of individuals used for
data analysis was 330, sourced from 12 countries where
toads were sampled. Sequencing data were extracted from
the raw data and combined into single forward reads using
FLASH (Magoč and Salzberg 2011). Output files from
FLASH were converted into fasta file format with the soft-
ware NextGen Sequence Workbench (Heracle BioSoft).
The pools were dimultiplexed and sequences were paired
by the attached primer tags to the corresponding individu-
als through the jMHC software (Stuglik et al. 2011). The
output files were aligned with the MUSCLE alignment
package and primer tags were trimmed from the sequences
for easier analysis (Edgar 2004). An Excel macro (Lighten
et al. 2014) was used to estimate allele variants per indi-
vidual based on Degree of Change (DOC) calculations.
The DOC criterion calculates cumulative amplicon
sequencing depth among the variants in the amplicon.
Amplicons with total sequence reads ≤300 and with DOC
ai counts ≥ 5 and higher were excluded for consecutively
analyses for quality reasons. All valid allele sequences
selected for further confirmation using similar criteria to
Cortázar-Chinarro et al. (2017). The allele sequences were
then aligned and analysed by eye in MEGA 7 to ensure no
stop codons were present to note insertions and deletions

(Kumar et al. 2018). Abnormal sequences were examined
with a BLAST search and removed to filter out contamina-
tion (Altschul et al. 1990). The final valid allele aligned
sequences were named using structure proposed by Klein
(1975), e.g. BuvaDAB01. The sequences were then con-
structed into a phylogenetic tree with the neighbor-joining
method (Saitou and Nei 1987). An unrooted phylogenetic
network tree of the sequences was also created with the
software SplitsTree4 (Huson and Bryant 2006). The domi-
nant sequences for each sample were assigned based on
locus and the DOC ai values. For each country, nucleotide
diversity (π), observed (H_o) and expected (H_c) heterozygo-
sity were computed with DnaSP v6 (Rozas et al. 2017) and
Arlequin 3.5 (Excoffier and Lischer 2010), respectively.
We calculated heterozygosities for MHC separately and for
the putative loci combined. When comparing populations
in different countries we used the combined estimate (see
Rozsa et al. 2016).

To test for signals of balancing (positive) selection we
calculated dN/dS for the entire amplified exon 2 sequences
and also for the peptide binding region (putative antigen
binding sites according to Bondinas et al. 2007, PBR) and
non-PBR sites using the software MEGA X (Kumar et al.
2018).

Results

Distribution of mitochondrial lineages

Combining our new data (340 new individuals from 69
European localities) with the dataset of Dufresnes et al.
(2019), we mapped the occurrence of the original (*viridis*)
and foreign (*sitibundus*) mtDNA lineages across the range
of the European taxon *B. viridis* (Fig. 1). These show some
latitudinal associations, with *viridis* mtDNA spread on the
western side of the distribution (from Albania – Hungary
– Slovakia – Czechia – Germany), as well as in far-eastern
Europe (Russia – E-Ukraine), while *sitibundus* mitochondria
extend in between, from the Balkan Peninsula (Greece, Ser-
bia, Romania) to the Baltic coast (Poland, Latvia, Denmark,
Sweden).

Genetic structure and diversity in the microsatellite
dataset

The microsatellite data (776 individuals from 134 locali-
ties) showed some departures from Hardy–Weinberg expec-
tations and linkage when analysed by country (Table S3).
However, most of the significant departures and associa-
tions disappeared when analysed by population (Table S3).
Moreover, inheritance analyses never mapped these loci on
the same linkage groups (e.g. Betto-Colliard et al. 2015;
CAGTATATACAGGGTCTCACC) to target and amplify
sequences of 279-282 base pairs of the MHC Class II B
gene as described in Zeisset and Beebee (2014). Samples
were randomly placed into pool groups that determined the
forward and reverse primer tag combinations which attached
to the primers (e.g. NNN-ACTAGATC) for each sample, so
that individual samples could be identified during analysis.
The PCR mixtures were made using 13.9 µl double distilled
water, 4 µl Phusion High Fidelity Buffer, 0.4 µl dNTPS,
0.5 µl each of the required forward and reverse primers,
0.2 µl Phusion High Fidelity DNA Taq polymerase, and
0.5 µl of the extracted DNA. The PCRs were run, and the
results of each sample were measured on a 1.5% agarose gel
after electrophoresis. We included 20% of replicated sam-
ples. Samples were grouped by the strength of the resulting
bands, separating the pools into strong and weak clusters.
Bands within the pre-pools were purified by using a Qia-
gen® MinElute gel extraction procedure (Qiagen, Sollen-
tuna), following the manufacturer's protocol. The pre-pools
were measured for DNA concentration with the Qbit and
consecutively pooled in equimolar amounts. Genetic librar-
ies from every pre-pool (A, B, C and D) were prepared using
ThruPLEX DNA-seq Kit (Rubicon Genomics) following the
manufactures protocol. Sequencing was completed using
Illumina MiSeq at the National Genomic Infrastructure,
SciLifeLab in Uppsala.

MHC data analyses

For MHC variation, final number of individuals used for
data analysis was 330, sourced from 12 countries where
toads were sampled. Sequencing data were extracted from
the raw data and combined into single forward reads using
FLASH (Magoč and Salzberg 2011). Output files from
FLASH were converted into fasta file format with the soft-
ware NextGen Sequence Workbench (Heracle BioSoft).
The pools were dimultiplexed and sequences were paired
by the attached primer tags to the corresponding individu-
als through the jMHC software (Stuglik et al. 2011). The
output files were aligned with the MUSCLE alignment
package and primer tags were trimmed from the sequences
for easier analysis (Edgar 2004). An Excel macro (Lighten
et al. 2014) was used to estimate allele variants per indi-
vidual based on Degree of Change (DOC) calculations.
The DOC criterion calculates cumulative amplicon
sequencing depth among the variants in the amplicon.
Amplicons with total sequence reads ≤300 and with DOC
ai counts ≥ 5 and higher were excluded for consecutively
analyses for quality reasons. All valid allele sequences
selected for further confirmation using similar criteria to
Cortázar-Chinarro et al. (2017). The allele sequences were
then aligned and analysed by eye in MEGA 7 to ensure no
stop codons were present to note insertions and deletions

(Kumar et al. 2018). Abnormal sequences were examined
with a BLAST search and removed to filter out contamina-
tion (Altschul et al. 1990). The final valid allele aligned
sequences were named using structure proposed by Klein
(1975), e.g. BuvaDAB01. The sequences were then con-
structed into a phylogenetic tree with the neighbor-joining
method (Saitou and Nei 1987). An unrooted phylogenetic
network tree of the sequences was also created with the
software SplitsTree4 (Huson and Bryant 2006). The domi-
nant sequences for each sample were assigned based on
locus and the DOC ai values. For each country, nucleotide
diversity (π), observed (H_o) and expected (H_c) heterozygo-
sity were computed with DnaSP v6 (Rozas et al. 2017) and
Arlequin 3.5 (Excoffier and Lischer 2010), respectively.
We calculated heterozygosities for MHC separately and for
the putative loci combined. When comparing populations
in different countries we used the combined estimate (see
Rozsa et al. 2016).

To test for signals of balancing (positive) selection we
calculated dN/dS for the entire amplified exon 2 sequences
and also for the peptide binding region (putative antigen
binding sites according to Bondinas et al. 2007, PBR) and
non-PBR sites using the software MEGA X (Kumar et al.
2018).

Results

Distribution of mitochondrial lineages

Combining our new data (340 new individuals from 69
European localities) with the dataset of Dufresnes et al.
(2019), we mapped the occurrence of the original (viridis)
and foreign (sitibundus) mtDNA lineages across the range
of the European taxon B. viridis (Fig. 1). These show some
latitudinal associations, with viridis mtDNA spread on the
western side of the distribution (from Albania – Hungary
– Slovakia – Czechia – Germany), as well as in far-eastern
Europe (Russia – E-Ukraine), while sitibundus mitochondria
extend in between, from the Balkan Peninsula (Greece, Ser-
bia, Romania) to the Baltic coast (Poland, Latvia, Denmark,
Sweden).

Genetic structure and diversity in the microsatellite
dataset

The microsatellite data (776 individuals from 134 locali-
ties) showed some departures from Hardy–Weinberg expec-
tations and linkage when analysed by country (Table S3).
However, most of the significant departures and associa-
tions disappeared when analysed by population (Table S3).
Moreover, inheritance analyses never mapped these loci on
the same linkage groups (e.g. Betto-Colliard et al. 2015;
a

<figure><img src="image_9.png" /></figure>

b

<figure><img src="image_10.png" /></figure>

c

<figure><img src="image_11.png" /></figure>
a

<figure><img src="image_8.png" /></figure>

b

<figure><img src="image_9.png" /></figure>

c

<figure><img src="image_10.png" /></figure>
►Fig. 4 a neighbor-joining tree (numbers above branches indicate bootstrap support values, when above 50%) and b unrooted phylogenetic network of MHC Class II B sequences; sequences found in Swedish samples are indicated in bold, with asterisks; c geographic occurrence of these sequences (pie size is proportional to the number of samples)

Gerchen et al. 2018). We thus assume the microsatellite loci to be neutral and unlinked, and that instances of departures from neutrality and linkage disequilibrium can be ascribed to population structure within some countries (e.g. Poland).

Accordingly, in the PCA, the first component highlighted differentiation between southern and eastern European ranges compared to the northern ranges (Fig. 2a). The geographic pattern was confirmed by STRUCTURE analyses (Fig. 2b), where these populations are sorted under *K*=2 (the preferred solution of the ΔK statistic, Fig. S1). Increasing the number of STRUCTURE groups highlighted some additional yet widely admixing clusters, in the northern ranges (*K*=3) and the Balkans (*K*=4). Interestingly, the isolated populations from W-Germany stand out as a separate cluster by *K*=5 (which also shows a ΔK peak, Fig. S1). Further increase of *K* did not improve the likelihood nor the ΔK (Fig. S1), and did not reveal additional structure with a meaningful geographic association. These instances of genetic structure were confirmed by inspections of pairwise $F_{st}$ (Fig. S2).

Both within-population diversity indices, namely the observed and expected heterozygosity $H_o$ and $H_e$, recovered similar patterns (Fig. 3a), i.e. higher estimates in southern populations compared to the rest of range, with the lowest values found around the Baltic Sea. Accordingly, there was a significant negative relationship between genetic variation and latitude ($P<0.001$, df=41; linear regressions, with $R^2$ of 0.40 and 0.45 for $H_o$ and $H_e$, respectively, Fig. 3b).

## MHC Class II variability

Of the obtained sequences, 32 were 279 bp in length while 5 (BuvaDAB22–26) were 282 bp due to insertions/deletions identified after alignment (Fig. S3). Across the sequences, 82 sites were found to be polymorphic. All sequences were present in individuals collected from mainland Europe, while 14 of the sequences were not found in Swedish toads (Fig. 4a).

We interpret the two main basal branches of the neighbor-joining tree and the corresponding sequences on the left of the network as belonging to two different loci: locus 1 containing sequences 1 through 21 and 27 through 32, and locus 2 containing sequences 22 through 26 (Fig. 4b). No further divisions of loci could be determined with certainty. Out of all individuals, 29 contain sequences from both loci while the remaining 311 individuals only expressed sequences

contained within locus 1, indicating copy number variation. Ratios of dN/dS showed that the putative PBR harboured more non-synonymous substitutions than the non-PBR both when analysed overall and as separate loci (Table S4). This indicates that the amplified loci have been under positive adaptive selection.

Comparing diversity estimates among the Swedish and other European populations reveal low diversity in the former, especially with respect to $H_o$, which was lower than expected (Table 1; Fig. 5). Furthermore, there was a trend that within populations, MHC- allele sequence diversity tended to be dominated by a single allele sequence towards the northern limits of the distribution (Fig. 4), thus explaining the lower $H_o$.

## Discussion

The main purpose of this study was to investigate the genetic diversity of *Bufotes viridis* green toad populations throughout its vast range, with a particular focus on the isolated and red-listed Swedish populations. The species is one of only 11 anuran taxa found in the country, and like most of them, its occurrence is confined to the southernmost provinces (here Blekinge, Skåne and Öland). Considerable conservation effort has been directed towards protecting *B. viridis* during the last decades, including captive breeding and support releases (Wirén 2006). In 2015, it was listed as vulnerable on the national red list, after previously being classified as critically endangered (https://artfakta.se/naturvard/taxon/100022). Historically, the species has been reported from more than 50 localities, but at present, breeding is confirmed in only about ten of these (Kauri 1948; Andrén and Nilson 1993).

Moreover, it is worth mentioning that Swedish green toads have been affected by taxonomic and even nomenclatural changes in recent years (Dufresnes et al. 2019). *Bufotes systematics* are complicated by far reaching introgression of *B. sitibundus* (previously "*B. variabilis*") mt-DNA into the range of *B. viridis* (Stöck et al. 2006; Dufresnes et al. 2019, 2021, this study). As a consequence, Swedish populations used to be called "*B. variabilis*" (as still listed on the Swedish red list, following Stöck 2005), even though they actually correspond to *B. viridis*, the only species present in continental Europe – expect in Ukraine-Russia (*B. sitibundus*), Italy and Spain (*B. balearicus*) (Dufresnes et al. 2019). In addition, the interspecific status of these taxa remains controversial, and some authors rather consider them as subspecies (Speybroeck et al. 2020), in which case Swedish green toads corresponds to *B. viridis viridis*. Regardless of how the peripheral Swedish populations are called or classified, they are of great interest from a conservation genetic perspective,
◄Fig. 4 a neighbor-joining tree (numbers above branches indicate bootstrap support values, when above 50%) and b unrooted phylogenetic network of MHC Class II B sequences; sequences found in Swedish samples are indicated in bold, with asterisks; c geographic occurrence of these sequences (pie size is proportional to the number of samples)

Gerchen et al. 2018). We thus assume the microsatellite loci to be neutral and unlinked, and that instances of departures from neutrality and linkage disequilibrium can be ascribed to population structure within some countries (e.g. Poland).

Accordingly, in the PCA, the first component highlighted differentiation between southern and eastern European ranges compared to the northern ranges (Fig. 2a). The geographic pattern was confirmed by STRUCTURE analyses (Fig. 2b), where these populations are sorted under *K*=2 (the preferred solution of the ΔK statistic, Fig. S1). Increasing the number of STRUCTURE groups highlighted some additional yet widely admixing clusters, in the northern ranges (*K*=3) and the Balkans (*K*=4). Interestingly, the isolated populations from W-Germany stand out as a separate cluster by *K*=5 (which also shows a ΔK peak, Fig. S1). Further increase of *K* did not improve the likelihood nor the ΔK (Fig. S1), and did not reveal additional structure with a meaningful geographic association. These instances of genetic structure were confirmed by inspections of pairwise $F_{st}$ (Fig. S2).

Both within-population diversity indices, namely the observed and expected heterozygosity $H_o$ and $H_e$, recovered similar patterns (Fig. 3a), i.e. higher estimates in southern populations compared to the rest of range, with the lowest values found around the Baltic Sea. Accordingly, there was a significant negative relationship between genetic variation and latitude ($P<0.001$, df=41; linear regressions, with $R^2$ of 0.40 and 0.45 for $H_o$ and $H_e$, respectively, Fig. 3b).

## MHC Class II variability

Of the obtained sequences, 32 were 279 bp in length while 5 (BuvaDAB22–26) were 282 bp due to insertions/deletions identified after alignment (Fig. S3). Across the sequences, 82 sites were found to be polymorphic. All sequences were present in individuals collected from mainland Europe, while 14 of the sequences were not found in Swedish toads (Fig. 4a).

We interpret the two main basal branches of the neighbor-joining tree and the corresponding sequences on the left of the network as belonging to two different loci: locus 1 containing sequences 1 through 21 and 27 through 32, and locus 2 containing sequences 22 through 26 (Fig. 4b). No further divisions of loci could be determined with certainty. Out of all individuals, 29 contain sequences from both loci while the remaining 311 individuals only expressed sequences

contained within locus 1, indicating copy number variation. Ratios of dN/dS showed that the putative PBR harboured more non-synonymous substitutions than the non-PBR both when analysed overall and as separate loci (Table S4). This indicates that the amplified loci have been under positive adaptive selection.

Comparing diversity estimates among the Swedish and other European populations reveal low diversity in the former, especially with respect to $H_o$, which was lower than expected (Table 1; Fig. 5). Furthermore, there was a trend that within populations, MHC-allele sequence diversity tended to be dominated by a single allele sequence towards the northern limits of the distribution (Fig. 4), thus explaining the lower $H_o$.

## Discussion

The main purpose of this study was to investigate the genetic diversity of *Bufotes viridis* green toad populations throughout its vast range, with a particular focus on the isolated and red-listed Swedish populations. The species is one of only 11 anuran taxa found in the country, and like most of them, its occurrence is confined to the southernmost provinces (here Blekinge, Skåne and Öland). Considerable conservation effort has been directed towards protecting *B. viridis* during the last decades, including captive breeding and support releases (Wirén 2006). In 2015, it was listed as vulnerable on the national red list, after previously being classified as critically endangered (https://artfakta.se/naturvard/taxon/100022). Historically, the species has been reported from more than 50 localities, but at present, breeding is confirmed in only about ten of these (Kauri 1948; Andrén and Nilson 1993).

Moreover, it is worth mentioning that Swedish green toads have been affected by taxonomic and even nomenclatural changes in recent years (Dufresnes et al. 2019). *Bufotes* systematics are complicated by far reaching introgression of *B. sitibundus* (previously “*B. variabilis*”) mt-DNA into the range of *B. viridis* (Stöck et al. 2006; Dufresnes et al. 2019, 2021, this study). As a consequence, Swedish populations used to be called “*B. variabilis*” (as still listed on the Swedish red list, following Stöck 2005), even though they actually correspond to *B. viridis*, the only species present in continental Europe – expect in Ukraine-Russia (*B. sitibundus*), Italy and Spain (*B. balearicus*) (Dufresnes et al. 2019). In addition, the interspecific status of these taxa remains controversial, and some authors rather consider them as subspecies (Speybroeck et al. 2020), in which case Swedish green toads corresponds to *B. viridis viridis*. Regardless of how the peripheral Swedish populations are called or classified, they are of great interest from a conservation genetic perspective,
<table><caption>**Table 1** MHC sequence variation among all populations</caption><thead><tr><th>Population</th><th>n</th><th>nS</th><th>π</th><th>H<sub>o</sub></th><th>H<sub>e</sub></th></tr></thead><tbody><tr><td>Belarus</td><td>2</td><td>3</td><td>0.06691</td><td>1.00</td><td>0.750</td></tr><tr><td>Bulgaria</td><td>2</td><td>3</td><td>0.07348</td><td>0.500*</td><td>1.00</td></tr><tr><td>Denmark</td><td>18</td><td>5</td><td>0.06221</td><td>0.667</td><td>0.477</td></tr><tr><td>Germany</td><td>15</td><td>5</td><td>0.06691</td><td>0.143</td><td>0.00</td></tr><tr><td>Locus 1</td><td>14</td><td>4</td><td>–</td><td>0.286</td><td>0.267</td></tr><tr><td>Locus 2</td><td>1</td><td>1</td><td>–</td><td>0.00</td><td>0.00</td></tr><tr><td>Greece</td><td>5</td><td>6</td><td>0.08058</td><td>0.542</td><td>0.333</td></tr><tr><td>Locus 1</td><td>3</td><td>3</td><td>–</td><td>0.333</td><td>0.667</td></tr><tr><td>Locus 2</td><td>4</td><td>3</td><td>–</td><td>0.750</td><td>0.583</td></tr><tr><td>Latvia</td><td>2</td><td>1</td><td>–</td><td>0.00</td><td>0.00</td></tr><tr><td>Moldova</td><td>1</td><td>2</td><td>0.08273</td><td>0.00</td><td>0.00</td></tr><tr><td>Locus 1</td><td>1</td><td>1</td><td>–</td><td>0.00</td><td>0.00</td></tr><tr><td>Locus 2</td><td>1</td><td>1</td><td>–</td><td>0.00</td><td>0.00</td></tr><tr><td>Poland</td><td>13</td><td>7</td><td>0.07719</td><td>0.308*</td><td>0.756</td></tr><tr><td>Romania</td><td>17</td><td>13</td><td>0.08115</td><td>0.377*</td><td>0.767</td></tr><tr><td>Locus 1</td><td>14</td><td>10</td><td>–</td><td>0.643</td><td>0.868</td></tr><tr><td>Locus 2</td><td>9</td><td>3</td><td>–</td><td>0.111</td><td>0.667</td></tr><tr><td>Russia</td><td>40</td><td>15</td><td>0.07874</td><td>0.375*</td><td>0.733</td></tr><tr><td>Locus 1</td><td>36</td><td>10</td><td>–</td><td>0.750</td><td>0.607</td></tr><tr><td>Locus 2</td><td>7</td><td>5</td><td>–</td><td>0.00</td><td>0.857</td></tr><tr><td>Sweden</td><td>218</td><td>11</td><td>0.07439</td><td>0.103*</td><td>0.614</td></tr><tr><td>Locus 1</td><td>218</td><td>7</td><td>–</td><td>0.206</td><td>0.228</td></tr><tr><td>Locus 2</td><td>4</td><td>4</td><td>–</td><td>0.00</td><td>1.00</td></tr><tr><td>Ukraine</td><td>7</td><td>10</td><td>0.07058</td><td>0.071*</td><td>0.833</td></tr><tr><td>Locus 1</td><td>7</td><td>8</td><td>–</td><td>0.143</td><td>1.00000</td></tr><tr><td>Locus 2</td><td>3</td><td>2</td><td>–</td><td>0.00</td><td>0.667</td></tr></tbody></table>

Populations without specified locus rows only possessed single locus-1 individuals

$H_o$ values marked with a * were significantly lower than the expected n number of samples, nS number of sequences, π nucleotide diversity, $H_o$ observed heterozygosity, $H_e$ expected heterozygosity

notably to compare their diversity with central populations of the range.

In this respect, we show that Swedish green toads harbor among the lowest levels of genetic variation documented in the species across Europe, both at putatively neutral microsatellites and MHC loci. This is also true for mtDNA, where a single haplotype of the *B. sitibundus* lineage (SIT12) has become almost fixed from Anatolia to northern ranges (Dufresnes et al. 2019). Note that the independent origin of the samples analyzed for each marker set prohibits direct statistical comparisons (see Methods). For MHC, our analyses of dN/dS ratios at PBR and non-PBR sites suggest that the sequences analyzed are subjected to selection (see also e.g. Cortazar-Chinarro et al. 2017, for evidence of selection on MHC Class II in Swedish *Rana arvalis*). The same gradients of diversity seem to apply to several amphibian species, at least in anurans, where northern populations show

<figure><img src="image_12.png" /><figcaption>**Fig. 5** Comparison of a observed ($H_o$) and b expected ($H_e$) heterozygosity indices from microsatellite (grey) and MHC Class II B sequences (white) between countries where both data were available for at least 5 toads</figcaption></figure>

lesser immunogenetic variation than southern populations, e.g. natterjack toad *Epidalea calamita* (Zeisset and Beebee 2014; Höglund et al. 2015), common frog *Rana temporaria* (Palo et al. 2004), moor frog *R. arvalis* (Cortazar-Chinarro et al. 2017) and common toad *Bufo bufo* (Meurling 2019). In fact, lower variability in northern ranges was also retrieved from neutral markers e.g. microsatellites in *R. arvalis* (Cortazar-Chinarro et al. 2017), and ddRAD-derived SNPs in
<table><caption>**Table 1** MHC sequence variation among all populations</caption><thead><tr><th>Population</th><th>n</th><th>nS</th><th>π</th><th>H<sub>o</sub></th><th>H<sub>e</sub></th></tr></thead><tbody><tr><td>Belarus</td><td>2</td><td>3</td><td>0.06691</td><td>1.00</td><td>0.750</td></tr><tr><td>Bulgaria</td><td>2</td><td>3</td><td>0.07348</td><td>0.500*</td><td>1.00</td></tr><tr><td>Denmark</td><td>18</td><td>5</td><td>0.06221</td><td>0.667</td><td>0.477</td></tr><tr><td>Germany</td><td>15</td><td>5</td><td>0.06691</td><td>0.143</td><td>0.00</td></tr><tr><td>Locus 1</td><td>14</td><td>4</td><td>–</td><td>0.286</td><td>0.267</td></tr><tr><td>Locus 2</td><td>1</td><td>1</td><td>–</td><td>0.00</td><td>0.00</td></tr><tr><td>Greece</td><td>5</td><td>6</td><td>0.08058</td><td>0.542</td><td>0.333</td></tr><tr><td>Locus 1</td><td>3</td><td>3</td><td>–</td><td>0.333</td><td>0.667</td></tr><tr><td>Locus 2</td><td>4</td><td>3</td><td>–</td><td>0.750</td><td>0.583</td></tr><tr><td>Latvia</td><td>2</td><td>1</td><td>–</td><td>0.00</td><td>0.00</td></tr><tr><td>Moldova</td><td>1</td><td>2</td><td>0.08273</td><td>0.00</td><td>0.00</td></tr><tr><td>Locus 1</td><td>1</td><td>1</td><td>–</td><td>0.00</td><td>0.00</td></tr><tr><td>Locus 2</td><td>1</td><td>1</td><td>–</td><td>0.00</td><td>0.00</td></tr><tr><td>Poland</td><td>13</td><td>7</td><td>0.07719</td><td>0.308*</td><td>0.756</td></tr><tr><td>Romania</td><td>17</td><td>13</td><td>0.08115</td><td>0.377*</td><td>0.767</td></tr><tr><td>Locus 1</td><td>14</td><td>10</td><td>–</td><td>0.643</td><td>0.868</td></tr><tr><td>Locus 2</td><td>9</td><td>3</td><td>–</td><td>0.111</td><td>0.667</td></tr><tr><td>Russia</td><td>40</td><td>15</td><td>0.07874</td><td>0.375*</td><td>0.733</td></tr><tr><td>Locus 1</td><td>36</td><td>10</td><td>–</td><td>0.750</td><td>0.607</td></tr><tr><td>Locus 2</td><td>7</td><td>5</td><td>–</td><td>0.00</td><td>0.857</td></tr><tr><td>Sweden</td><td>218</td><td>11</td><td>0.07439</td><td>0.103*</td><td>0.614</td></tr><tr><td>Locus 1</td><td>218</td><td>7</td><td>–</td><td>0.206</td><td>0.228</td></tr><tr><td>Locus 2</td><td>4</td><td>4</td><td>–</td><td>0.00</td><td>1.00</td></tr><tr><td>Ukraine</td><td>7</td><td>10</td><td>0.07058</td><td>0.071*</td><td>0.833</td></tr><tr><td>Locus 1</td><td>7</td><td>8</td><td>–</td><td>0.143</td><td>1.00000</td></tr><tr><td>Locus 2</td><td>3</td><td>2</td><td>–</td><td>0.00</td><td>0.667</td></tr></tbody></table>

Populations without specified locus rows only possessed single locus-1 individuals

$H_o$ values marked with a * were significantly lower than the expected n number of samples, nS number of sequences, π nucleotide diversity, $H_o$ observed heterozygosity, $H_e$ expected heterozygosity

notably to compare their diversity with central populations of the range.

In this respect, we show that Swedish green toads harbor among the lowest levels of genetic variation documented in the species across Europe, both at putatively neutral microsatellites and MHC loci. This is also true for mtDNA, where a single haplotype of the *B. sitibundus* lineage (SIT12) has become almost fixed from Anatolia to northern ranges (Dufresnes et al. 2019). Note that the independent origin of the samples analyzed for each marker set prohibits direct statistical comparisons (see Methods). For MHC, our analyses of dN/dS ratios at PBR and non-PBR sites suggest that the sequences analyzed are subjected to selection (see also e.g. Cortazar-Chinarro et al. 2017, for evidence of selection on MHC Class II in Swedish *Rana arvalis*). The same gradients of diversity seem to apply to several amphibian species, at least in anurans, where northern populations show

<figure><img src="image_11.png" /><figcaption>**Fig. 5** Comparison of a observed ($H_o$) and b expected ($H_e$) heterozygosity indices from microsatellite (grey) and MHC Class II B sequences (white) between countries where both data were available for at least 5 toads</figcaption></figure>

lesser immunogenetic variation than southern populations, e.g. natterjack toad *Epidalea calamita* (Zeisset and Beebee 2014; Höglund et al. 2015), common frog *Rana temporaria* (Palo et al. 2004), moor frog *R. arvalis* (Cortazar-Chinarro et al. 2017) and common toad *Bufo bufo* (Meurling 2019). In fact, lower variability in northern ranges was also retrieved from neutral markers e.g. microsatellites in *R. arvalis* (Cortazar-Chinarro et al. 2017), and ddRAD-derived SNPs in
R. arvalis (Rödin-Mörch et al. 2019) and B. bufo (Thörn et al. 2021).

These general latitudinal trends can be explained by two non-mutually exclusive hypotheses. First, the central-marginal hypothesis posits that peripheral populations harbor less variation than central core populations (Eckert et al. 2008). Second, populations at northern latitudes are predicted to be less variable because of the recolonization history since the last glaciation (Hewitt 1996, 1999). In essence, both hypotheses rest on the idea that populations at the range edge have been subjected to demographic histories signified by smaller and more fragmented populations, which would result into losses of genetic diversity through the increased effects of genetic drift (Kimura 1983). The patterns observed here for green toads are in line with the recolonization hypothesis, where the southern populations bear the highest levels of genetic diversity. Accordingly, previous phylogeographic studies inferred post-glacial expansions from southern refugia putatively located in the Balkan Peninsula (Stöck et al. 2006; Dufresnes et al. 2019). In addition, the presence of discrete genetic clusters in the northern ranges, corresponding to areas with the most impoverished diversity, could be another effect of drift, as genetic variants distinct from the core populations become fixed (e.g. a few MHC alleles are in much higher proportions in Sweden and Germany than elsewhere). However, the comparatively high genetic diversity in southern Europe may have been inflated by the diffusion of alleles from the closely related *B. sitibundus*. This taxon forms a wide hybrid zone with *B. viridis* across Western Anatolia, which could have been previously located in the Balkans (Dufresnes et al. 2021). The occurrence of foreign alleles through past introgression, as seen at least from mtDNA, would obviously obscure the comparisons.

Northern green toad populations are predicted to be more threatened and less adaptable than southern populations, as shown for most widespread European amphibians (Dufresnes and Perrin 2015). Given the fragmented distribution and low population sizes of Swedish green toads, it is questionable whether it is at all possible for these populations to adapt to local conditions, since the strength of genetic drift is inversely related to effective population size (Kimura 1983). The current captive breeding program has mainly selected individuals originating from a single locality in Sweden, and there is some evidence of adaptive variation in a few larval life history traits (Rogell et al. 2011). Should the breeding program extend to additional populations, it may be advisable to consider the origin of the sources. This is especially important since MHC-variation may be locally adapted to the parasite fauna (Savage and Zamudio 2011, 2016). Furthermore, the unusual pattern of mtDNA diversity, namely mitochondrial capture of many *B. viridis*

populations, is worthy of consideration. Studies on insects have shown clear evidence of cytoplasmic-nuclear incompatibilities (e.g. Dowling et al. 2006, 2007). The far-reaching spread of mtDNA of *B. sitibundus* origin could thus affect the ability of *B. viridis* to persist under local conditions in northern Europe.

To conclude, our study illustrates the strong heterogeneity in the level of genetic diversity across the range of a widespread species that was strongly affected by past climate changes. Importantly, the loss of genetic variants associated with the higher drift experienced by populations that have experienced strong founder effects during post-glacial expansion from refugia and that remain small and fragmented because they occupy more challenging environments on the range periphery, seem to have lost both neutral (microsatellites) and adaptive (MHC) genetic diversity.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10592-021-01407-5.

**Acknowledgements** We thank numerous field workers and curators of herpetological collections for their assistance in obtaining the samples. Gunilla Engström provided assistance with the MHC lab work. The study was supported by a grant from the Swedish Research Council Formas (215-2014-294) to JH. Sampling was done with permission from local authorities.

**Funding** Open access funding provided by Uppsala University.

## Declarations

**Conflict of interest** The authors declare no conflict of interests (MHC sequences will be uploaded on GenBank on acceptance; microsatellite and mitochondrial barcoding data are available in Table S1).

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Altschul SF, Gish W, Miller W (1990) Basic alignment search tool. J Mol Biol 205:403–410

Andrén C, Nilson G (1993) Grönfläckiga Paddan Och Framtiden. Skånes Nat 80:27–30

Adams R, Hadly EA (2012) Genetic diversity within vertebrate species is greater at lower latitudes. Evol Ecol 27:133–143

Bondinas GP, Moustakas AK, Papadoupoulos GK (2007) The spectrum of HLA-DQ and HLA-DR alleles, 2006: a listing
R. arvalis (Rödin-Mörch et al. 2019) and B. bufo (Thörn et al. 2021).

These general latitudinal trends can be explained by two non-mutually exclusive hypotheses. First, the central-marginal hypothesis posits that peripheral populations harbor less variation than central core populations (Eckert et al. 2008). Second, populations at northern latitudes are predicted to be less variable because of the recolonization history since the last glaciation (Hewitt 1996, 1999). In essence, both hypotheses rest on the idea that populations at the range edge have been subjected to demographic histories signified by smaller and more fragmented populations, which would result into losses of genetic diversity through the increased effects of genetic drift (Kimura 1983). The patterns observed here for green toads are in line with the recolonization hypothesis, where the southern populations bear the highest levels of genetic diversity. Accordingly, previous phylogeographic studies inferred post-glacial expansions from southern refugia putatively located in the Balkan Peninsula (Stöck et al. 2006; Dufresnes et al. 2019). In addition, the presence of discrete genetic clusters in the northern ranges, corresponding to areas with the most impoverished diversity, could be another effect of drift, as genetic variants distinct from the core populations become fixed (e.g. a few MHC alleles are in much higher proportions in Sweden and Germany than elsewhere). However, the comparatively high genetic diversity in southern Europe may have been inflated by the diffusion of alleles from the closely related *B. sitibundus*. This taxon forms a wide hybrid zone with *B. viridis* across Western Anatolia, which could have been previously located in the Balkans (Dufresnes et al. 2021). The occurrence of foreign alleles through past introgression, as seen at least from mtDNA, would obviously obscure the comparisons.

Northern green toad populations are predicted to be more threatened and less adaptable than southern populations, as shown for most widespread European amphibians (Dufresnes and Perrin 2015). Given the fragmented distribution and low population sizes of Swedish green toads, it is questionable whether it is at all possible for these populations to adapt to local conditions, since the strength of genetic drift is inversely related to effective population size (Kimura 1983). The current captive breeding program has mainly selected individuals originating from a single locality in Sweden, and there is some evidence of adaptive variation in a few larval life history traits (Rogell et al. 2011). Should the breeding program extend to additional populations, it may be advisable to consider the origin of the sources. This is especially important since MHC-variation may be locally adapted to the parasite fauna (Savage and Zamudio 2011, 2016). Furthermore, the unusual pattern of mtDNA diversity, namely mitochondrial capture of many *B. viridis*

populations, is worthy of consideration. Studies on insects have shown clear evidence of cytoplasmic-nuclear incompatibilities (e.g. Dowling et al. 2006, 2007). The far-reaching spread of mtDNA of *B. sitibundus* origin could thus affect the ability of *B. viridis* to persist under local conditions in northern Europe.

To conclude, our study illustrates the strong heterogeneity in the level of genetic diversity across the range of a widespread species that was strongly affected by past climate changes. Importantly, the loss of genetic variants associated with the higher drift experienced by populations that have experienced strong founder effects during post-glacial expansion from refugia and that remain small and fragmented because they occupy more challenging environments on the range periphery, seem to have lost both neutral (microsatellites) and adaptive (MHC) genetic diversity.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10592-021-01407-5.

**Acknowledgements** We thank numerous field workers and curators of herpetological collections for their assistance in obtaining the samples. Gunilla Engström provided assistance with the MHC lab work. The study was supported by a grant from the Swedish Research Council Formas (215-2014-294) to JH. Sampling was done with permission from local authorities.

**Funding** Open access funding provided by Uppsala University.

## Declarations

**Conflict of interest** The authors declare no conflict of interests (MHC sequences will be uploaded on GenBank on acceptance; microsatellite and mitochondrial barcoding data are available in Table S1).

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Altschul SF, Gish W, Miller W (1990) Basic alignment search tool. J Mol Biol 205:403–410

Andrén C, Nilson G (1993) Grönfläckiga Paddan Och Framtiden. Skånes Nat 80:27–30

Adams R, Hadly EA (2012) Genetic diversity within vertebrate species is greater at lower latitudes. Evol Ecol 27:133–143

Bondinas GP, Moustakas AK, Papadoupoulos GK (2007) The spectrum of HLA-DQ and HLA-DR alleles, 2006: a listing
correlating sequence and structure with function. Immunogen-
etics 59:539–553

Betto-Colliard C, Sermier R, Litvinchuk S, Perrin N, Stöck M (2015) Origin and genome evolution of polyploid green toads in Central Asia: evidence from microsatellite markers. Heredity 114:300–308

Cortazar-Chinarro M, Lattenkamp EZ, Meyer-Lucht Y, Luquet E, Laurila A, HöglundJ, (2017) Drift, selection or migration? Processes affecting genetic differentiation along a latitudinal gradient in an amphibian. BMC Evol Biol 17:189

Dowling DK, Friberg U, Arnqvist G (2006) A comparison of nuclear and cytoplasmatic genetic effects on sperm competitiveness and female remating in a seed beetle. J Evol Biol 20:2113–2125

Dowling DK, Friberg U, Hailer F, Arnqvist G (2007) Intergenomic epistasis for fitness: within-population interactions between cytoplasmic and nuclear genes in *Drosophila melanogaster*. *Genetics* 175:235–244

Dufresnes C, Perrin N (2015) Effect of biogeographic history on population vulnerability in European amphibians. Cons Biol 29:1235–1241

Dufresnes C, Lymberakis P, Kornilios P, Savary R, Perrin N, Stöck M (2018) Phylogeography of Aegean green toads (Bufo viridis subgroup): continental hybrid swarm vs. insular diversification with discovery of a new island endemic. BMC Evol Biol 18:67

Dufresnes C, Mazepa G, Jablonski D, Caliari Oliveira R, Wenseleers T, Shabanov DA, Auer M, Ernst R, Ramírez-Chaves HE, Mulder KP, Simonov E; Tiutenko A, Kryvokhyzha D, Wennekes PL, Zinenko OI, Korshunov OV, Al-Johany AM, Peregentsev EA, Betto-Colliard, Denoël M, Borkin LJ, Skorinov DV, Pasyn-kova RA, Mazanaeva LF, Rosanov JM, Dubey S, Litvinchuk SN (2019) Fifteen shades of green: the evolution of *Bufotes* toads revisited. Mol Phylogen Evol 141:106615

Dufresnes C, Probonas NM, Strachinis I (2021) A reassessment of the diversity of green toads (*Bufotes*) in the circum-Aegean region. Int Zool 16:420–428

Earl DA, von Holdt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Res 4:359–361

Eckert C, Samis K, Lougheed S (2008) Genetic variation across species' geographical ranges: the central–marginal hypothesis and beyond. Mol Ecol 17:1170–1188

Edgar RC (2004) MUSCLE: a multiple sequence alignment with high accuracy and high throughput. Nucl Acids Res 32:1792–1797

Excoffier L, Lischer HEL (2010) Arlequin suite ver. 3.5: a new series of programs to perform population genetic analyses under Linux and Windows. Mol Ecol Res 10:564–567

Fu M, Waldman B (2017) Major histocompatibility complex variation and the evolution of resistance to amphibian chytridiomycosis. Immunogenetics 69:529–536

Gerchen J, Dufresnes C, Stöck M (2018) Introgression across hybrid zones is not mediated by Large X-effects in green toads with undifferentiated sex chromosomes. Am Nat 192:178–188

Goudet J (1995) FSTAT (version 1.2): a computer program to calculate F-statistics. J Hered 86:485–486

Goudet J (2005) *hierfstat*, a package for R to compute and test hierarchical F-statistics. Mol Ecol Notes 5:184–186

Hewitt GM (1996) Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc 58:247–276

Hewitt GM (1999) Post-glacial re-colonization of European biota. Biol J Linn Soc 68:87–112

Hewitt GM (2000) The genetic legacy of the Quaternary ice ages. Nature 405:907–913

Hewitt GM (2004) Genetic consequences of climatic oscillations in the Quaternary. Philos Trans Roy Soc B 359:183–195

Höglund J (2009) Evolutionary conservation genetics. Oxford University Press, Oxford

Höglund J, Wengström Å, Rogell B, Meyer-Lucht Y (2015) Low MHC variation in isolated island populations of the Natterjack toad (*Bufo calamita*). Conserv Genet 16:1007–1010

Holderegger R, Kamm U, Gugerli F (2006) Adaptive vs. neutral genetic diversity: implications for landscape genetics. Landsc Ecol 21:797–807

Huson DH, Bryant D (2006) Application of phylogenetic networks in evolutionary studies. Mol Biol Evol 23:242–252

Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23:1801–1806

Jombart T (2008) *adegenet*: A R package for the multivariate analysis of genetic markers. Bioinformatics 24:1403–1405

Kauri H (1948) Über die Ausbreitung und die Ausbreitungsumstände der Wechselkröte (*Bufo viridis*) im Ostseegebiet. Kungliga Fysio-grafiska Sällskapets Handlingar 59:1–30

Kimura M (1983) The neutral theory of molecular evolution. Cambridge University Press, Cambridge

Klein J (1975) Biology of the mouse histocompatibility-2 complex. Springer, Berlin

Kumar S, Stecher G, Li M, Knyaz C, Tamura K (2018) MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol Biol Evol 35:1547–1549. https://doi.org/10.1093/molbev/msy096.

Li H, Xiang-Ju J, Dai G, Gu Z, Ming C, Yang Z, Ryder OA, Li W-H, Fu Y-X, Zhang Y-P (2016) Large numbers of vertebrates began rapid population decline in the late 19th century. PNAS 113:14079–14085

Lighten J, van Oosterhout C, Paterson IG, McCullan M, Bentzen P (2014) Ultra-deep Illumina sequencing accurately identifies MHC Class IIB alleles and provides evidence for copy number variation in the guppy. Mol Ecol Res 14:753–767

Magoc T, Salzberg SL (2011) FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27:2957–2963

McMahon BJ, Teeling EC, Höglund J (2014) How and why should we implement genomics in conservation? Evol Appl 7:999–1007

Meurling S (2019) *The response in native wildlife to an invading pathogen: Swedish amphibians and Batrachochytrium dendrobatidis*. PhD thesis, Uppsala University. http://uu.divaportal.org/smash/record.jsf?pid=diva2%3A1369696&dswid=3881

Miraldo A, Li S, Borregard MK, Florez-Rodriguez A, Gopalakrishnan S, Nogues-Bravo RMD (2016) An anthropocene map of genetic diversity. Science 353:1532–1535

Palo JU, Schmeller DS, Laurila A, Primmer CR, Kuzmin SL, Merilä J (2004) High degree of population subdivision in a widespread amphibian. Mol Ecol 13:2631–2644

Paxton RJ, Thorén PA, Tengö J, Estoup A. Pamilo P (1996) Mating structure and nestmate relatedness in a communal bee, *Andrena jacobi* (Hymenoptera, Andrenidae), using microsatellites. Mol Ecol 5:511–519

Piertney SB, Oliver MK (2006) The evolutionary ecology of the major histocompatibility complex. Heredity 96:7–21

Pritchard JK, Stephens M, Donelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945–959

Rödin-Mörch P, Luquet E, Meyer-Lucht Y, Richter-Boix A, Höglund J, Laurila A (2019) Latitudinal divergence in a wide-spread amphibian: contrasting patterns of neutral and adaptive genomic variation. Mol Ecol 28:2996–3011

Rogell B, Berglund A, Laurila A, Höglund J (2011) Population divergence of life history traits in the endangered green toad: implications for a support release programme. J Zool 285:46–55
correlating sequence and structure with function. Immunogen-
etics 59:539–553

Betto-Colliard C, Sermier R, Litvinchuk S, Perrin N, Stöck M (2015) Origin and genome evolution of polyploid green toads in Central Asia: evidence from microsatellite markers. Heredity 114:300–308

Cortazar-Chinarro M, Lattenkamp EZ, Meyer-Lucht Y, Luquet E, Laurila A, HöglundJ, (2017) Drift, selection or migration? Processes affecting genetic differentiation along a latitudinal gradient in an amphibian. BMC Evol Biol 17:189

Dowling DK, Friberg U, Arnqvist G (2006) A comparison of nuclear and cytoplasmatic genetic effects on sperm competitiveness and female remating in a seed beetle. J Evol Biol 20:2113–2125

Dowling DK, Friberg U, Hailer F, Arnqvist G (2007) Intergenomic epistasis for fitness: within-population interactions between cytoplasmic and nuclear genes in *Drosophila melanogaster*. *Genetics* 175:235–244

Dufresnes C, Perrin N (2015) Effect of biogeographic history on population vulnerability in European amphibians. Cons Biol 29:1235–1241

Dufresnes C, Lymberakis P, Kornilios P, Savary R, Perrin N, Stöck M (2018) Phylogeography of Aegean green toads (Bufo viridis subgroup): continental hybrid swarm vs. insular diversification with discovery of a new island endemic. BMC Evol Biol 18:67

Dufresnes C, Mazepa G, Jablonski D, Caliari Oliveira R, Wenseleers T, Shabanov DA, Auer M, Ernst R, Ramírez-Chaves HE, Mulder KP, Simonov E; Tiutenko A, Kryvokhyzha D, Wennekes PL, Zinenko OI, Korshunov OV, Al-Johany AM, Peregontsev EA, Betto-Colliard, Denoël M, Borkin LJ, Skorinov DV, Pasyn-kova RA, Mazanaeva LF, Rosanov JM, Dubey S, Litvinchuk SN (2019) Fifteen shades of green: the evolution of *Bufotes* toads revisited. Mol Phylogen Evol 141:106615

Dufresnes C, Probonas NM, Strachinis I (2021) A reassessment of the diversity of green toads (*Bufotes*) in the circum-Aegean region. Int Zool 16:420–428

Earl DA, von Holdt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Res 4:359–361

Eckert C, Samis K, Lougheed S (2008) Genetic variation across species' geographical ranges: the central–marginal hypothesis and beyond. Mol Ecol 17:1170–1188

Edgar RC (2004) MUSCLE: a multiple sequence alignment with high accuracy and high throughput. Nucl Acids Res 32:1792–1797

Excoffier L, Lischer HEL (2010) Arlequin suite ver. 3.5: a new series of programs to perform population genetic analyses under Linux and Windows. Mol Ecol Res 10:564–567

Fu M, Waldman B (2017) Major histocompatibility complex variation and the evolution of resistance to amphibian chytridiomycosis. Immunogenetics 69:529–536

Gerchen J, Dufresnes C, Stöck M (2018) Introgression across hybrid zones is not mediated by Large X-effects in green toads with undifferentiated sex chromosomes. Am Nat 192:178–188

Goudet J (1995) FSTAT (version 1.2): a computer program to calculate F-statistics. J Hered 86:485–486

Goudet J (2005) *hierfstat*, a package for R to compute and test hierarchical F-statistics. Mol Ecol Notes 5:184–186

Hewitt GM (1996) Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc 58:247–276

Hewitt GM (1999) Post-glacial re-colonization of European biota. Biol J Linn Soc 68:87–112

Hewitt GM (2000) The genetic legacy of the Quaternary ice ages. Nature 405:907–913

Hewitt GM (2004) Genetic consequences of climatic oscillations in the Quaternary. Philos Trans Roy Soc B 359:183–195

Höglund J (2009) Evolutionary conservation genetics. Oxford University Press, Oxford

Höglund J, Wengström Å, Rogell B, Meyer-Lucht Y (2015) Low MHC variation in isolated island populations of the Natterjack toad (*Bufo calamita*). Conserv Genet 16:1007–1010

Holderegger R, Kamm U, Gugerli F (2006) Adaptive vs. neutral genetic diversity: implications for landscape genetics. Landsc Ecol 21:797–807

Huson DH, Bryant D (2006) Application of phylogenetic networks in evolutionary studies. Mol Biol Evol 23:242–252

Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23:1801–1806

Jombart T (2008) *adegenet*: A R package for the multivariate analysis of genetic markers. Bioinformatics 24:1403–1405

Kauri H (1948) Über die Ausbreitung und die Ausbreitungsumstände der Wechselkröte (*Bufo viridis*) im Ostseegebiet. Kungliga Fysio-grafiska Sällskapets Handlingar 59:1–30

Kimura M (1983) The neutral theory of molecular evolution. Cambridge University Press, Cambridge

Klein J (1975) Biology of the mouse histocompatibility-2 complex. Springer, Berlin

Kumar S, Stecher G, Li M, Knyaz C, Tamura K (2018) MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol Biol Evol 35:1547–1549. https://doi.org/10.1093/molbev/msy096.

Li H, Xiang-Ju J, Dai G, Gu Z, Ming C, Yang Z, Ryder OA, Li W-H, Fu Y-X, Zhang Y-P (2016) Large numbers of vertebrates began rapid population decline in the late 19th century. PNAS 113:14079–14085

Lighten J, van Oosterhout C, Paterson IG, McCullan M, Bentzen P (2014) Ultra-deep Illumina sequencing accurately identifies MHC Class IIB alleles and provides evidence for copy number variation in the guppy. Mol Ecol Res 14:753–767

Magoc T, Salzberg SL (2011) FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27:2957–2963

McMahon BJ, Teeling EC, Höglund J (2014) How and why should we implement genomics in conservation? Evol Appl 7:999–1007

Meurling S (2019) *The response in native wildlife to an invading pathogen: Swedish amphibians and Batrachochytrium dendrobatidis*. PhD thesis, Uppsala University. http://uu.divaportal.org/smash/record.jsf?pid=diva2%3A1369696&dswid=3881

Miraldo A, Li S, Borregard MK, Florez-Rodriguez A, Gopalakrishnan S, Nogues-Bravo RMD (2016) An anthropocene map of genetic diversity. Science 353:1532–1535

Palo JU, Schmeller DS, Laurila A, Primmer CR, Kuzmin SL, Merilä J (2004) High degree of population subdivision in a widespread amphibian. Mol Ecol 13:2631–2644

Paxton RJ, Thorén PA, Tengö J, Estoup A. Pamilo P (1996) Mating structure and nestmate relatedness in a communal bee, *Andrena jacobi* (Hymenoptera, Andrenidae), using microsatellites. Mol Ecol 5:511–519

Piertney SB, Oliver MK (2006) The evolutionary ecology of the major histocompatibility complex. Heredity 96:7–21

Pritchard JK, Stephens M, Donelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945–959

Rödin-Mörch P, Luquet E, Meyer-Lucht Y, Richter-Boix A, Höglund J, Laurila A (2019) Latitudinal divergence in a wide-spread amphibian: contrasting patterns of neutral and adaptive genomic variation. Mol Ecol 28:2996–3011

Rogell B, Berglund A, Laurila A, Höglund J (2011) Population divergence of life history traits in the endangered green toad: implications for a support release programme. J Zool 285:46–55

<figure><img src="image_12.png" /></figure>
Rosenberg NA (2004) DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes 4:137–138

Rozas J, Ferrer-Mata A, Sanchez-DelBarrio JC, Guiraco-Rico S, Librado P, Ramos-Onsins SE, Sanchez-Gracia A (2017) DNAsp6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol 34:3299–3302

Rozsa J, Strand TM, Montadert M, Kozma R, Höglund J (2016) Effects of a range expansion on adaptive and neutral genetic diversity in dispersal limited Hazel grouse (Bonasa bonasia) in the French Alps. Conserv Genet 17:401–412

Saitou N, Nei M (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol 4:406–425

Savage AE, Zamudio KR (2011) MHC genotypes associate with resistance to a frog killer. PNAS 108:16705–16710

Savage AE, Zamudio KR (2016) Adaptive tolerance to a pathogenic fungus drives major histocompatibility complex evolution in natural amphibian populations. Proc Roy Soc B 283:20153115

Sillero N, Campos J, Bonardi A, Corti C, Creemers R, Crochet PA, Isailovic JC, Denoel M, Ficetola GF, Goncalves J, Kusmin S, Lymberakis P, de Pous P, Rodriguez A, Sindaco R, Speybroeck J, Toxopeus B, Vieites DR, Vences M (2014) Uptaded distribution and biogeography of amphibians and reptilians in Europe. Amph-Rept 35:1–31

Smith BT, Seeholzer GF, Harvey MG, Cuervo AM, Brumfield RT (2017) A latitudinal phylogeographic diversity gradient in birds. PLoS Biol 15:e2001073

Sommer S (2005) The importance of immune gene variability (MHC) in evolutionary ecology and conservation. Front Zool 2(1):1–18

Speybroeck J, Beukema W, Dufresnes C, Fritz U, Jablonski D, Lymberakis P, Martinez-Solano I, Razzeti E, Vamberger M, Vences M, Vörös J, Corchet P-A (2020) Species list of the European herpetofauna - 2020 update by the taxonomic committee of the Societas Europea Herpetologica. Amphibia-Reptilia 41:139–189

Spielman D, Brook BW, Frankham R (2004) Most species are not driven to extinction before genetic factors impact them. PNAS 101:15261–15264

Stöck M (2005) *Preliminary report on the phylogeographic analysis of green toads (Bufo viridis) subgroup from Sweden*. Länsstyrelsen i Skåne län.

Stöck M, Moritz C, Hickerson M, Frynta D, Dusjebayeva T, Eremchenko V, Macey JR, Papenfuss TJ, Wake DB (2006) Evolution of mitochondrial relationships and biogeography of Palearctic green toads (Bufo viridis subgroup) with insights in their plasticity. Mol Phylogenet Evol 41:663–689

Stuglik MT, Radwan J, Babik W (2011) jMHC: software assistant for multilocus genotyping of gene families using next-generation amplicon sequencing. Mol Ecol Resour 11:739–742

Thörn F, Rödin-Mörch P, Cortazar-Chinarro M, Richter-Boix A, Laurial A, Höglund J (2021) The effects of drift and selection on latitudinal genetic variation in Scandinavian common toads (Bufo bufo) following post-glacial colonization. Heredity 126:656–667

Vellend M (2005) Species diversity and genetic diversity: Parallel processes and correlated patterns. Am Nat 166:199–215

Vellend M, Geber MA (2005) Connections between species diversity and genetic diversity. Ecol Lett 8:767–781

Wan HY, Cushman SA, Ganey JL (2018) Habitat fragmentation reduces genetic diversity and connectivity of the Mexican spotted owl: a simulation study using empirical resistance models. Genes 9:403

Wirén M (2006) *Grönfläckig padda (Bufo viridis) i Sverige – Utvärdering av utförda bevarandeåtärder (1994–2005), förslag om framtida åtgärder samt artens tidigare och nuvarande förekomst.* Länsstyrelsen i Skåne och Naturvårdsverket.

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rosenberg NA (2004) DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes 4:137–138

Rozas J, Ferrer-Mata A, Sanchez-DelBarrio JC, Guiraco-Rico S, Librado P, Ramos-Onsins SE, Sanchez-Gracia A (2017) DNAsp6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol 34:3299–3302

Rozsa J, Strand TM, Montadert M, Kozma R, Höglund J (2016) Effects of a range expansion on adaptive and neutral genetic diversity in dispersal limited Hazel grouse (Bonasa bonasia) in the French Alps. Conserv Genet 17:401–412

Saitou N, Nei M (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol 4:406–425

Savage AE, Zamudio KR (2011) MHC genotypes associate with resistance to a frog killer. PNAS 108:16705–16710

Savage AE, Zamudio KR (2016) Adaptive tolerance to a pathogenic fungus drives major histocompatibility complex evolution in natural amphibian populations. Proc Roy Soc B 283:20153115

Sillero N, Campos J, Bonardi A, Corti C, Creemers R, Crochet PA, Isailovic JC, Denoel M, Ficetola GF, Goncalves J, Kusmin S, Lymberakis P, de Pous P, Rodriguez A, Sindaco R, Speybroeck J, Toxopeus B, Vieites DR, Vences M (2014) Uptaded distribution and biogeography of amphibians and reptilians in Europe. Amph-Rept 35:1–31

Smith BT, Seeholzer GF, Harvey MG, Cuervo AM, Brumfield RT (2017) A latitudinal phylogeographic diversity gradient in birds. PLoS Biol 15:e2001073

Sommer S (2005) The importance of immune gene variability (MHC) in evolutionary ecology and conservation. Front Zool 2(1):1–18

Speybroeck J, Beukema W, Dufresnes C, Fritz U, Jablonski D, Lymberakis P, Martinez-Solano I, Razzeti E, Vamberger M, Vences M, Vörös J, Corchet P-A (2020) Species list of the European herpetofauna - 2020 update by the taxonomic committee of the Societas Europea Herpetologica. Amphibia-Reptilia 41:139–189

Spielman D, Brook BW, Frankham R (2004) Most species are not driven to extinction before genetic factors impact them. PNAS 101:15261–15264

Stöck M (2005) *Preliminary report on the phylogeographic analysis of green toads (Bufo viridis) subgroup from Sweden*. Länsstyrelsen i Skåne län.

Stöck M, Moritz C, Hickerson M, Frynta D, Dusjebayeva T, Eremchenko V, Macey JR, Papenfuss TJ, Wake DB (2006) Evolution of mitochondrial relationships and biogeography of Palearctic green toads (Bufo viridis subgroup) with insights in their plasticity. Mol Phylogenet Evol 41:663–689

Stuglik MT, Radwan J, Babik W (2011) jMHC: software assistant for multilocus genotyping of gene families using next-generation amplicon sequencing. Mol Ecol Resour 11:739–742

Thörn F, Rödin-Mörch P, Cortazar-Chinarro M, Richter-Boix A, Laurial A, Höglund J (2021) The effects of drift and selection on latitudinal genetic variation in Scandinavian common toads (Bufo bufo) following post-glacial colonization. Heredity 126:656–667

Vellend M (2005) Species diversity and genetic diversity: Parallel processes and correlated patterns. Am Nat 166:199–215

Vellend M, Geber MA (2005) Connections between species diversity and genetic diversity. Ecol Lett 8:767–781

Wan HY, Cushman SA, Ganey JL (2018) Habitat fragmentation reduces genetic diversity and connectivity of the Mexican spotted owl: a simulation study using empirical resistance models. Genes 9:403

Wirén M (2006) *Grönfläckig padda (Bufo viridis) i Sverige – Utvärdering av utförda bevarandeåtärder (1994–2005), förslag om framtida åtgärder samt artens tidigare och nuvarande förekomst.* Länsstyrelsen i Skåne och Naturvårdsverket.

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.