Pronounced genetic structure associated with differences in a reproductive trait and climatic barriers in Canadian populations of the western toad ( Anaxyrus boreas )

Pronounced genetic structure associated with differences in a reproductive trait and climatic barriers in Canadian populations of the western toad (Anaxyrus boreas) | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Pronounced genetic structure associated with differences in a reproductive trait and climatic barriers in Canadian populations of the western toad ( Anaxyrus boreas ) View ORCID Profile Jayna C. Bergman , Juan Enciso-Romero , View ORCID Profile Gregory B. Pauly , View ORCID Profile Roseanna Gamlen-Greene , Melissa Todd , Julie A. Lee-Yaw doi: https://doi.org/10.1101/2025.06.17.660222 Jayna C. Bergman 1 Department of Biology, University of Ottawa , Ottawa K1N 9B4, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Jayna C. Bergman For correspondence: jaynabergman{at}gmail.com Juan Enciso-Romero 1 Department of Biology, University of Ottawa , Ottawa K1N 9B4, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site Gregory B. Pauly 2 Department of Herpetology, Natural History Museum of Los Angeles County , USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Gregory B. Pauly Roseanna Gamlen-Greene 3 Department of Biochemistry, University of Otago , Aotearoa New Zealand Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Roseanna Gamlen-Greene Melissa Todd 4 Coast Area Research Section, British Columbia Ministry of Forests , Nanaimo, British Columbia, V9T 6E9, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site Julie A. Lee-Yaw 1 Department of Biology, University of Ottawa , Ottawa K1N 9B4, Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract Identifying genetic groups within species is critical for inventorying biodiversity, understanding species’ distributions, and delineating conservation units. The western toad ( Anaxyrus boreas ) is one of the most widespread amphibians in western Canada and one of just two amphibian species to maintain an extensive range on both sides of the Canadian Rocky Mountains. Curiously, toads on the eastern side of the mountains have a vocal sac and produce a distinct advertisement call whereas toads on the western side of the mountains lack these traits. However, the extent to which these “Calling” and “Non-Calling” populations are genetically and ecologically distinct remains unclear. We used reduced representation sequencing (ddRAD-seq) to genotype individuals from across the Canadian portion of the species’ range. Combining population genomic analysis with the original phenotypic data used to delineate the Calling and Non-Calling populations and with predictions from ecological niche models, we show that Calling and Non-Calling populations of western toads represent distinct genetic groups separated by both differences in a major reproductive trait and by climate barriers. We additionally find evidence for a third, undescribed group in southwestern Alberta and southeastern British Columbia. Our results suggest Calling and Non-Calling populations should be managed as distinct groups and highlight the potential for there to be cryptic and strong genetic structure among northern populations of widespread species. Introduction Characterizing intraspecific genetic structure is fundamental to understanding species’ geographic distributions. In particular, genetic structure provides insight into the history of species’ ranges and the factors influencing colonization (Hewit, 2004; Cairns et al., 2021 ), as well as the extent to which populations are currently connected. Genetic breaks in the range inform our understanding of geographic and ecological barriers that influence dispersal ( Waters et al., 2020 ; Suarez et al., 2022) and thus our understanding of the ability of species to track suitable conditions as the climate changes. Finally, assessments of genetic structure within species’ ranges are important for inventorying biodiversity, allowing for the identification of cryptic taxa ( Bickford et al., 2007 ; Fiser et al., 2017) and the delineation of conservation units for wildlife management ( Funk et al., 2012 ; Coates et al., 2018 ; Forester et al., 2022 ). Genetic substructure is expected to be particularly pronounced in the Northern Hemisphere, where species experienced repeated range expansions and contractions during the Pleistocene glaciations ( Hewitt, 2004 ). Indeed, numerous studies have used molecular markers and methods to assess substructure within species’ ranges in areas impacted by the Pleistocene glaciations, revealing distinct genetic groups within what are otherwise widespread species (see review by: Soltis et al., 1997 ; Hewitt, 2000 ; Brunsfelt et al., 2001; Sommer and Zachos, 2009 ; Shafer et al., 2010 ; Lyman and Edwards, 2022 ; Fu and Wen, 2023 ; Garcia-Rodriguez et al., 2024). However, gaps remain in our understanding of structure in some regions. For example, for some taxonomic groups in North America, genetic studies are heavily biased towards populations in the United States, limiting our understanding of genetic structure in the previously glaciated, northern portion of many species’ ranges ( Lesbarrères et al., 2014 ; Bergman et al. unpublished). In other cases, what is known about genetic structure across Canada comes from earlier studies using mitochondrial markers or a limited number of nuclear markers (e.g. Graham and Burg, 2012 ; Benefer et al., 2013 ; Nobarinezhad et al., 2020 ; but see Lehnert et al., 2023 ; Miller et al., 2024 ; Stacy et al., 2025 ). It is now widely recognized that issues with incomplete lineage sorting, resolution, and/or introgression may limit the conclusions that can be drawn from studies using sparse sampling of the genome to assess genetic structure (e.g. Ballard and Whitlock, 2004 ; Funk et al., 2019 ). Thus, for many northern species, our understanding of intraspecific genetic diversity remains incomplete. Next generation sequencing has provided an opportunity to revisit questions about genetic structure within species’ ranges using information from across the genome and has clarified the number and distribution of genetic groups within the ranges of several species (e.g. Barbosa et al., 2018 ; Dinca et al., 2019; Cairns et al., 2021 ; Gallego-García et al., 2021 ). However, care is needed when using genomic data to delineate meaningful groups. For instance, when there are gaps in the sampling of the geographic range, what is continuous genetic variation associated with isolation-by-distance (IBD) may be misinterpreted as distinct genetic substructure ( Chambers and Hillis, 2020 ; Turbek et al., 2023 ). Ensuring samples are representative of the transition across suspected contact zones is thus important for identifying true genetic breaks within species’ ranges ( Chambers et al., 2022 ; Turbek et al., 2023 ). At the same time, the high resolution of genomic data is likely to reveal fine-scale structure, which may lead to issues with the over splitting of groups ( Sukumaran and Knowles, 2017 ; Coates et al., 2018 ; Chambers et al., 2025 ). Thus, it is ideal to pair genomic assessments of structure with phenotypic (e.g. Gordon et al., 2020 ; Waples and Lindley, 2018 ) and/or ecological (e.g. Campbell et al., 2022 ) data that can provide insight into the biological significance of any observed groups. Amphibians are expected to exhibit strong genetic structure within their ranges because they have low dispersal ability, are often philopatric, and show evidence of limited gene flow over small spatial scales ( Beebee, 2005 ; Zeisset and Beebee, 2008 ). Indeed, amphibians in North America have been shown to consist of distinct genetic groups, with continentally-distributed species often consisting of eastern and western clades (e.g. Lee-Yaw et al., 2008 ; Everson et al., 2021 ). Within the eastern United States, intraspecific genetic structure is often aligned with major geographic features such as the Atlantic Coast of the Florida peninsula, Apalachicola and Mississippi Rivers, and the Appalachian Mountains ( Soltis et al., 2006 ; Pauly, et al., 2007 ; Lyman and Edwards, 2022 ). In the west, sharp genetic breaks also tend to line up with major geographic features, including the Cascade/Coast mountains range, the southern Rocky Mountains, and the Columbia River Basin ( Zeisset and Beebee, 2008 ; Shafer et al., 2010 ). However, most phylogeographic studies of North American amphibians are based on populations in the United States ( Lesbarrères et al., 2014 ; but see Lee-Yaw et al., 2008 ; Wilson et al., 2008 ; Lee-Yaw and Irwin, 2012 ; Cairns et al., 2021 ), limiting our understanding of the features that shape genetic structure across much of North America, as well as efforts to protect genetic diversity in northern areas. Western toads ( Anaxyrus boreas ; also known as Bufo boreas ) are one of the most widely distributed amphibians in western North America and are one of just two amphibian species in Canada to maintain an extensive range on either side of the Canadian Rocky Mountains (CRMs), a barrier that shapes genetic diversity in other taxonomic groups ( Shafer et al., 2010 ; Jensen et al., 2024 ). Furthermore, western toads on either side of the CRMs exhibit differences in key reproductive traits: individuals east of the mountains in Alberta have a vocal sac and pronounced advertisement call, whereas those to the west in British Columbia (BC) lack both traits ( Fig. 1 ; Pauly, 2008 ). To the best of our knowledge, this is the only documented case in anurans where there is a polymorphism in terms of the presence of and ability to produce male advertisement calls within the range of the same species (see Elias-Costa & Faivovich 2025 for more general discussion of the gain and loss of calls across Anurans). Based on these differences, “Calling” and “Non-Calling” populations are recognized as distinct conservation units in Canada (i.e. Designatable Units) by the Committee on the Status of Endangered Wildlife in Canada ( COSEWIC 2012 ). Yet, despite attention to genetic structure elsewhere in the species’ range ( Goebel et al., 2009 ; Switzer et al., 2009 ; Addis et al., 2015 ; Gordon, 2017 ; Gordon et al., 2020 ; US Fish and Wildlife, 2017; Gamlen-Greene, 2022 ; Trumbo et al., 2023 ; Fig. 1A ), the extent to which Calling and Non-Calling populations of western toads are genetically distinct remains unclear. Given conservation concern for western toads across much of their range (e.g. Colorado Parks and Wildlife, New Mexico Department of Game and Fish, Wyoming Department of Game and Fish), including in Canada ( COSEWIC, 2012 ), clarifying genetic structure in this part of the species’ range not only adds to our understanding of genetic diversity in the north but has the potential to guide management plans for the species. Download figure Open in new tab Figure 1. (A) The geographic range of western toads ( Anaxyrus boreas ) in western North America (green shading). Open points represent approximate locations of the sampling sites used to identify mitochondrial groups (purple polygons) from previous work (recreated from Goebel et al., 2009 ). Genetic structure across the Canadian range, including across the Canadian Rocky Mountains (CRMs) is unknown. (B) Locations of the genetic sampling sites for the present study. The putative boundary between the Calling and Non-Calling populations is shown as a red line (redrawn from: Environment and Climate Change Canada, 2020 ). In this study, we evaluate genetic structure across the Canadian range of western toads. We specifically ask whether genomic data support the recognition of the Calling and Non-Calling populations as distinct groups and whether observed genetic boundaries line up with previously described variation in call morphology ( Pauly 2008 ). We additionally use ecological niche models to explore niche overlap between Calling and Non-Calling toads and the potential for gene flow between groups. We discuss the implications of our results for the conservation status of the Calling and Non-Calling populations in Canada and for our broader understanding of the biogeographical factors that shape genetic diversity in western North America. Materials and Methods Tissue sampling Tissue samples were collected between June and August of 2018 to 2023 from coastal BC to the eastern edge of the species’ range in Alberta ( Fig. 1B ; Table S1). Sites across the putative boundary between the Calling and Non-Calling populations were selected to be as close as possible to sites from which individuals were characterized for the presence or absence of vocal sacs (as per Pauly 2008 ). Between one and five individuals were collected from each site, resulting in a total of 46 tissue samples from 39 sites. Whole larvae (or young of the year when larvae were not available) were captured by hand or by dipnet from each site. Individuals were humanely euthanized on site using a 3% buffered MS-222 solution. Samples were preserved in 95% ethanol before being brought to the lab for long-term storage at -80°C. Fieldwork and sample collection was conducted under permits granted by the Government of Alberta (Alberta No. 22-175 and 23-185), Government of British Columbia Ministry of Forests (MRNA16-236743, MRCB23-787858), Parks Canada (WTNP-2023-45325), and the Solutions Table of the Council of the Haida Nation and the Haida Gwaii Natural Resource District (Gamlen-Greene Haida Gwaii Amphibian Research 03-Nov-16). Samples were collected under approved Canadian Council for Animal Care protocols to different members of the research team (University of Lethbridge protocol #1912, University of Ottawa protocol #3991, Government of British Columbia 236743-ACA, and University of British Columbia protocols A16-0210 and B16-0222). DNA extraction and sequencing Total genomic DNA was extracted from 25 mg of larval tail tissue or metamorphic leg muscle using Qiagen DNeasy Blood and Tissue Kits (Qiagen, Valencia, CA, USA). The manufacturer’s protocol was followed with the following modifications: We added a physical homogenization step, using a VWR bead mill with 2.8 mm ceramic beads to improve tissue digestion. Tissues were then digested for three hours at 56°C, after which, we added 6 μl of RNAase A solution (10 mg/ml). Samples were eluted in two steps, using 80 μl and 60 μl of buffer AE respectively. DNA concentration and quality were assessed using a Thermo Fisher Qubit 4 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and NanoDrop 1000 spectrophotometer (Nanodrop Technologies Inc. Wilmington, DE, USA). Samples were diluted and normalized to a concentration of 20 ng/μl in 20 μl (400 ng total) before library preparation. We obtained sequence data for each individual using double digest restriction site-associated DNA sequencing (ddRADseq; Petersen et al., 2012). Library preparation was done at the Plateforme d’Analyses Génomiques of the Institut de Biologie Intégrative et des Systèmes (IBIS, Université Laval, Quebec, Canada) and was based on the genotyping-by-sequencing procedures described by Poland et al., (2012) . Double-stranded DNA was digested using the restriction enzymes PstI (5-CTGCAG-3) and MspI (5’-CCGG-3’). A set of corresponding Illumina adapters were ligated to the ends of the fragmented DNA and the inline barcodes were ligated to the adapters. Prior to PCR amplification, a Blue Pippin (SAGE Sciences Inc. Beverly, MA, USA) was used to select fragment sizes between 100 and 300 bp using an elution time set between 50 and 65 minutes on a 2% gel. Individuals were sequenced in a single library that included both a technical replicate (replicated from DNA extraction) and two negative controls. A single library was sent to the Centre d’expertise et de service Génome Québec at McGill University in Montreal Quebec for paired-end sequencing on an Illumina NovaSeq 6000. There, multiple libraries were normalized, pooled, and then denatured in 0.02N NaOH and neutralized using HT1 buffer. Using the Xp protocol as per the manufacturer’s recommendations, pools were loaded at 200 pM on an Illumina NovaSeq S4 lane with 2 x 150 reads on a 300-cycle run in paired-end mode. Base calling was performed using RTA v3 and bcl2fastq2 v. 2.20 (Illumina Inc., 2019) was used to demultiplex the libraries. A total of 1,506,075,128 reads were generated across all individuals. Loci assembly and filtering Raw reads were demultiplexed using process_radtags in Stacks (version 2.66 Catchen et al., 2013 ). Removing reads with uncalled bases or adaptor contamination left 971,781,833 reads across all individuals. Fastp (version 0.23.4, Chen, 2023 ) was used to trim barcodes and adapters by removing five and three bases from the forward and reverse reads, respectively. We also used Fastp to remove reads with an average Phred quality score less than 20, resulting in 929,487,304 reads (1,853,538 to 56,194,548 reads per individual). We then used multiqc (version 1.21, Ewels et al., 2016 ) to verify sequencing quality and to ensure GC content along the reads was not excessive. Reads were aligned to the western toad reference genome (Jackson County, Colorado, USA; Trumbo et al., 2023 ) using the Burrows-Wheeler algorithm in bwa-mem (version 0.7.18, Li and Durbin, 2009 ). The resulting SAM files were sorted and converted to BAM format using Samtools (version 1.20, Li, 2011 ). We called single nucleotide polymorphisms (SNPs) from across mapped regions of the reference genome using the gstacks module in Stacks and filtered out stacks with a quality score less than 20. SNPs with high heterozygosity often result from paralogs or alignment errors. To remove SNPs with excess heterozygosity, we set the maximum observed heterozygosity (--max-obs-het) to 0.5 in the populations module in Stacks, treating all individuals as a single population. Further filtering was conducted in R (version 4.3.3, R core team, 2024 ), using the package SNPfiltR ( DeRaad, 2023 ). To ensure accurate genotype calls, we removed SNPs with a genotype depth of less than five or with more than 22 reads. We additionally removed genotype calls with quality of less than 20. The ratio of reads with reference versus alternative bases in heterozygous individuals is expected to be balanced, with large imbalances potentially signaling issues with coverage, multilocus contigs or other artifacts ( O’Leary et al., 2018 ). Thus, we filtered loci with an allele balance less than 0.25 or greater than 0.75 (i.e. the minimum number of reads to call a heterozygote was two). We then applied a minor allele count of three (i.e. with mac = 3, there must be at least one homozygote and one heterozygote for the alternative allele across all individuals to keep a SNP). We iteratively filtered for both missingness of individuals and SNPs (e.g. O’Leary et al. 2018 ). Given that most sites in this study are represented by a single individual, we prioritized the retention of individuals over SNPs during this process. Finally, we removed SNPs that are likely in linkage disequilibrium using the function --indep-pairwise in PLINK (version 1.90b, Chang et al., 2015 ). This function calculates the correlation coefficient between pairs of SNPs within a specified window and removes SNPs until all the correlations within the window are below a given threshold. We used a window size of 50 base pairs and tested step sizes of 5 or 10 base pairs and R 2 threshold values between 0.5 and 0.8 before settling on a final step size of 5 and R 2 of 0.8 (see Supplementary Methods). Additional details of loci assembly and filtering, including the number of SNPs retained at each step, can be found in the Supplementary Methods. Characterizing genetic structure across the range We used four approaches to assess genetic structure across the Canadian range of western toads. First, we visualized the distribution of individuals in multivariate genetic space using principal component analysis (PCA) and discriminant analysis of principal components (DAPC). PCA was run using PLINK ( Chang et al., 2015 ). DAPC was run using the R package Adegenet ( Jombart, 2008 ). Prior to running DAPC, we used the find.clusters function to determine the optimal number of genetic clusters based on the Bayesian Information Criterion (BIC) and to assign individuals to these clusters. We retained four principal components and two discriminant functions. In addition to ordination analyses, we assessed phylogenetic relationships between individuals by generating a maximum likelihood tree in IQTree (version 2.2.2.7, Minh et al., 2020 ), using the ascertainment bias correction (ASC) to account for the lack of invariant sites. Use of the ASC required removal of partially consistent sites (SNPs for which one variant is ambiguously called). Removing these sites left 9,648 SNPs for this analysis. IQTree was run with 1000 Ultrafast bootstrap replicates using the TMV+F+ASC+G4 model of sequence evolution, which was identified as the best model by ModelFinder ( Kalyaanamoorthy et al., 2017 ). Finally, we used the admixture model in STRUCTURE (version 2.3.4, Pritchard et al., 2000 ) to estimate the proportion of each individual’s ancestry derived from inferred genetic groups. STRUCTURE requires the user to specify the number of groups (K) a priori and we tested values of K from 1 to 10. Structure was run with a burn-in period of 300,000 steps, 1,000,000 Markov Chain Monte Carlo (MCMC) steps, and 10 replicate runs for each value of K. CLUMPAK ( Kopelman et al., 2015 ) was used to average the runs for each value of K, and ΔK ( Evanno et al., 2005 ) was used to determine the optimal value of K. All analyses revealed the same clustering of individuals (see Results), and we used the R package dartR ( Mijangos et al., 2022 ) to estimate F ST between genetic groups. To understand the role of IBD in structuring genetic diversity within the species’ range, we used Mantel tests to determine if there was a significant association between pairwise geographic and genetic distance. Because most sites in the dataset were represented by a single individual, we used the proportion of shared alleles between pairs of individuals as an estimate of genetic distance (rather than differences in allele frequencies). To avoid issues of non-independence associated with multiple individuals at some sites, for the few sites with multiple individuals, we randomly selected a single individual for this analysis. Genetic distances between individuals were calculated using the gl.dist.pop function from the dartR package and pairwise Euclidean geographic distances between individuals were calculated using the distm function in the R package geosphere (Hijmans, 2024a). Separate Mantel tests were run for: 1) all samples, 2) all samples excluding the three individuals from southeastern BC and southwestern Alberta that formed a separate genetic group (see Results), 3) all samples excluding two individuals from island sites in BC, 4) individuals from the Calling range and 5) individuals from the Non-Calling range. Mantel tests were run in R using the gl.ibd function from the dartR package. Although patterns of IBD are expected in widespread species with continuous distributions ( Turbek et al., 2023 ), if there is a sharp genetic break in the range, we expect estimates of pairwise genetic distances between groups to be higher than estimates within groups at all distance classes, including over relatively short geographic distances—that is we expect differences in the y-intercept in a plot of geographic distance versus genetic distance when pairs of samples are evaluated by comparison type. To test for a genetic break in IBD associated with the boundary between the Calling and Non-Calling populations, we focused on pairwise genetic distances between individuals that were separated by no more than 150 km (less than the average width of the CRMs), noting the comparison type in each case (Calling to Non-Calling: CN, Calling to Calling: CC, Non-calling to Non-Calling: NN). We then calculated the difference in mean pairwise genetic distances between the different comparison types (i.e. mean CC versus mean CN; mean NN versus mean CN). To test the statistical significance of these mean differences, we compared the observed mean difference in each case to a null distribution of differences in means generated by randomizing comparison type assignment across the pairwise genetic distance values 1000 times. Mean differences greater than the 95 th percentile of this distribution were considered to be statistically significant (i.e. one-sided test of significance). Concordance between genetic boundaries and call morphology Pauly (2008) previously examined 201 museum specimens of sexually mature, male, western toads from 84 sites across the Canadian portion of the range. Individuals were scored for the presence of vocal slits (from which vocal sacs originate). Specimens with one or both vocal slits were scored as Calling whereas those without vocal slits were scored as Non-Calling ( Pauly, 2008 ). We assigned the 39 individuals from that study with accurate GPS coordinate information and that were within 50 km of a genetic sample (a reasonable distance for life-time dispersal in this species: e.g. Bull, 2006 ; Thompson, 2019 ) to genetic clade based on Euclidean distance to the nearest genetic sample included in our phylogenetic tree. We then used a Fisher’s exact test to test the association between morphological classification and genetic clade assignment. Ecological niche models and niche overlap Ecological Niche Models (ENMs) have emerged as an important tool in molecular ecology for assessing ecological differences between genetic groups (e.g. Campbell et al., 2022 ; Gallego-Garcia et al., 2023; MacDonald et al., 2025b ) and the potential role of ecological barriers in maintaining groups (e.g. Glor and Warren, 2011 ; Lee-Yaw and Irwin, 2015 ). We generated separate ENMs for the Calling and Non-Calling populations to explore the climate variables that best predict presence and the distribution of suitable habitat for each population. Species occurrence records for the models were downloaded from the Global Biodiversity Information Facility (GBIF; December 26, 2024) and filtered to include occurrences from 1990 onward. To match the resolution of the environmental variables (below), we retained records with an accuracy of <1 km and retained one occurrence per grid cell of the environmental variables. Occurrence records were classified as belonging to the Calling (N = 215) or Non-Calling (N = 1523) population based on the recognized morphological boundaries ( Pauly 2008 ; Environment and Climate Change Canada, 2020 ). Thirty-one climate variables representing climate normals for the period between 1990 and 2020 were downloaded from ClimateNA ( Wang et al., 2016 ; Mahony et al., 2022 ) at a resolution of 30 arc seconds (1 km by 1 km) and were cropped to the boundaries of Alberta and BC. To avoid issues with collinearity in the models, we removed variables with a high variance inflation factor (VIF >10; calculated in R: usdm Naimi et al. 2014 ), retaining eight variables for the models: Julian date on which the frost-free period begins, degree-days below 18°C, mean temperature of the coldest month, precipitation as snow (mm), summer (June to August) precipitation (mm), spring (March to May) precipitation (mm), mean annual relative humidity (%), and summer heat moisture index. Niche models were generated in the R package predicts ( Hijmans, 2024b ) using Maxent (version 3.4.3, Phillips et al., 2006 ), which often outperforms other presence-only modeling algorithms ( Elith et al., 2011 ; Bradie and Leung, 2017 ). Separate models were generated for the Calling and Non-Calling population using the currently recognized range of each population in Canada ( Environment and Climate Change Canada, 2020 ). To reduce the impacts of sampling bias in our models, we sampled 5000 background points for each model from a bias grid corresponding to the density of observations of pond-breeding amphibians in each study extent (i.e. target-group background as per Barber et al., 2022 ; Baker et al. 2024 ). Records for the target-group (American bullfrog, boreal chorus frog, Canadian toad, Columbia spotted frog, long-toed salamander, northern leopard frog, northern pacific tree frog, northwestern salamander, western tiger salamander, wood frog) were retrieved from GBIF on October 21, 2021 and were filtered in a manner consistent with the occurrence records (above). The bias grid was generated by applying a Gaussian kernel density estimation to the target-group occurrences using the density.ppp function from the spatstat.core package in R (Baddeley et al., 2023). We used the ENMeval package ( Kass et al., 2021 ) in R to select the optimal regularization multiplier (RM) parameter and feature class setting for the Maxent models. Values of RM tested during tuning were: 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4. Feature classes tested were: L, LQ, H, LQH, LQP, LQT, LQHP, LQPT, and LQHPT where L= linear, Q = quadratic, H = hinge, P= product and T = threshold. The Akaike information criterion correction (AICc) was used to select the final combination of RM and feature classes with which to generate each model. Models were generated with the convergence threshold set to 0.001% ( Phillips et al., 2006 ) and the number of maximum iterations set to 5000 to allow models to converge ( Young et al., 2011 ). All other Maxent settings were kept as default values, including prevalence (𝜏 = 0.5), which has produced accurate models for this species in other contexts (e.g. Bergman et al., 2024 ). For each model, we determined the area under the receiver operating characteristic curve (AUC) based on average AUC from spatial block (4-fold) cross validation. Models were projected across BC and Alberta using Maxent’s cloglog output, resulting in a continuous prediction surface for each population. To avoid extrapolation, we used the mess function in the R package predicts to generate a Multivariate Environmental Similarity Surface (MESS; Elith et al., 2010 ) for each population and used this surface to mask areas in our prediction surfaces with conditions outside those used to calibrate the models. These areas can be treated as harboring environmental conditions that are “untested” by the population in question. To facilitate interpretation, we converted continuous model predictions into binary surfaces of suitable habitat using the 10% omission rate of the input locality records as a threshold for calling a cell suitable or not ( Vale et al., 2014 ; Rosner-Katz et al., 2020 ; Bergman et al., 2024 ). We used both the continuous and binary predictions of suitable climatic space to assess how climate may influence the boundaries of each population and the potential for gene flow between the two populations. We additionally visualized the amount of niche overlap between the Calling and Non-Calling populations based on all 31 ClimateNA variables using the PCAenv approach of Broennimann et al. (2012) as executed by the ecospat.plot.niche.dyn function in the R ecospat package ( Broennimann et al., 2025 ). We used the ecospat.niche.dyn.index function in ecospat to calculate niche stability ( Guisan et al., 2014 ). In our case, because both the accessible (background) and occupied environmental space for the Calling population were almost entirely nested within the much larger, occupied environmental space of the Non-Calling population (see Results), we used the test to ask what proportion of the Calling population’s occupied niche is also occupied by the Non-Calling population (i.e. test configured with the Calling population as the “native species” and D restricted to the shared environment space occupied by both populations). Finally, we used the ecospat.niche.similarity.test function in ecospat ( Warren et al., 2008 and Broennimann et al., 2012 ) to quantify niche overlap between the two populations using Schoener’s D ( Schoener, 1968 ). D ranges from 0 (no overlap) to 1 (complete overlap) and we tested whether the observed value of D was greater than a null distribution of values generated by fixing the observed niche of one population and randomly selecting a niche for second population within its range 1000 times. This test asks whether the niche of the population being randomized is more similar to that of the other population than expected by chance, given its accessible environment. We ran the test twice, alternating which niche was held constant. Results Sequencing and filtering After filtering and removing individuals with large amounts of missing data, our final genomic dataset consisted of 40 individuals from 34 sites distributed from coastal BC to the eastern edge of the species’ range in Alberta and inclusive of populations in and adjacent to the CRMs ( Fig. 1B ). Individuals were genotyped for 11,950 SNPs. Mean depth of coverage was 19.0X, with no more than 5% missing data per SNP and less than 33% missing data per individual. Pronounced genetic structure across the Canadian range of western toads All analyses assessing genetic structure clearly identified two groups of western toads corresponding to what are referred to as the Calling population to the east of the CRMs and the Non-Calling population to the west of the CRMs (hereafter we will use Calling and Non-Calling populations to refer to both the current Designatable Units and the genetic groups interchangeably). Specifically, Calling and Non-Calling toads separated along PC1, which explained 27.6% of the variation in the principal components analysis ( Fig. 2A ). Similar separation between the two groups was confirmed by the DAPC analysis (Fig. S1). Calling and Non-Calling individuals also formed distinct clades in the maximum likelihood tree, with bootstrap support values of 100 and 96 respectively ( Fig. 2B ). Finally, K of 2 in the STRUCTURE analysis received the highest support (ΔK value of 5146; Fig. S2), with Calling and Non-Calling individuals forming separate groups ( Fig. 2C ). Although there were some admixed individuals along the boundary between groups ( Fig. 2D ), average assignment of Calling and Non-Calling individuals to their respective groups was high (Q Non-Calling = 0.94; Q Calling = 0.98). Download figure Open in new tab Figure 2. Population structure within western toads ( Anaxyrus boreas) across the Canadian portion of the species’ range. Principal components analysis revealed three discrete genetic clusters corresponding to Calling and Non-Calling individuals, as well as an additional group (CotC) associated with the Crown of the Continent Ecosystem (A). These groups formed monophyletic clades in the unrooted maximum likelihood tree, with high bootstrap support (B). Two groups (K = 2) were optimal in the STRUCTURE analysis, representing Calling and Non-Calling individuals (C), which are separated by the Canadian Rocky Mountains (D). In all panels, orange symbols and shading represent the Calling population, blue symbols and shading represent the Non-Calling population, and black symbols and shading represent individuals from the CotC. Symbols in D represent locations from which individuals were sampled (one individual per location) and show the proportion of each individual’s ancestry to the two groups identified by the STRUCTURE analysis at K = 2. The colours of the thick borders around all points represent clade membership based on the maximum likelihood phylogeny in B. The extent of the Crown of the Continent Ecosystem in the south is indicated by the black polygon. In addition to the differences between the Calling and Non-Calling populations, our results revealed a third genetic group in Canada. Specifically, individuals from three sites in southwestern Alberta and southeastern BC separated from the Calling and Non-Calling populations along PC2, which explained 13.2% of the variation in the principal components analysis ( Fig. 2A ). These sites are located in what is referred to as the Crown of the Continent Ecosystem ( UNESCO, 2025 ; Fig. 2D ) and are hereafter referred to as the CotC group. The DAPC analysis with the lowest BIC (value of 231) supported K = 3, with the CotC group also forming a distinct group in this analysis (Fig. S1). Likewise, the CotC individuals formed a separate clade in the phylogenetic tree with bootstrap support of 100 ( Fig. 2B ). Although these individuals appear admixed in the STRUCTURE analysis when K = 2, they form a separate group at K = 4 (K=3 pulled out the Calling and Non-calling populations and a group for which the biological significance was unclear; Fig S3). Genetic differentiation between all three genetic groups was high, with F ST between the CotC group and the Calling (F ST = 0.52) and Non-Calling populations (F ST = 0.32) similar to or higher than between the Calling and Non-Calling populations (F ST = 0.43). Across all sites there was a significant pattern of IBD (Mantel R 2 = 0.60, p < 1.0e-4; Fig. 3A ). Significant IBD was also observed when focusing on individuals from either the Calling or Non-Calling populations (i.e. excluding CotC individuals: Mantel R 2 = 0.61, p < 1.0e-4) and when removing individuals from island sites (Mantel R 2 = 0.56, p < 1.0e-4). We also found a significant pattern of IBD within the ranges of both the Calling and Non-Calling populations (Calling: Mantel R 2 =0.36, p < 0.05; Non-Calling: Mantel R 0.69, p < 1.0e-4; Fig. 3B ). Download figure Open in new tab Figure 3. Pairwise comparisons of geographic distance (km) versus genetic distance (1 – proportion of shared alleles) for different subsets of the data (A and B). (A) The solid black line is the Mantel test regression line of the relationship between geographic and genetic distance (IBD) for all individuals. Green squares are pairwise comparisons of “CotC” to “CotC” individuals, and red points are pairwise comparisons between “CotC” individuals and either Calling or Non-Calling individuals. (B) The same plot from A but with comparisons involving the CotC group removed and points color-coded by comparison type for clarity (orange: pairs of individuals within the Calling population; blue: pairs of individuals within the Non-Calling population; purple: pairs of individuals from different populations). The orange and blue dashed lines are the Mantel test regression lines from the isolation-by-distance analyses for the Calling and Non-Calling populations respectively. The regression line through Calling vs Non-Calling comparisons (purple dashed line) is for illustration purposes only. (C) Mean pairwise genetic distances for different comparison types for those individuals within 150 km of each other. Asterisks above the bars connecting between-population comparisons indicate significance based on a randomization test. We found evidence for a sharp genetic break across the boundary between Calling and Non-Calling individuals. Specifically, across all geographic distances, pairwise genetic distance when comparing Calling to Non-Calling individuals was higher than when comparing Calling to Calling individuals or Non-Calling to Non-Calling individuals ( Fig. 3B ). Importantly, the mean genetic difference between Calling and Non-Calling individuals was significantly greater than the mean genetic differences between individuals from the same clade over relatively small distances (up to 150 km), including between sites spanning the boundary between genetic groups (mean difference of proportion of shared alleles CC vs CN = 0.137, p-value <0.05; mean difference of proportion of shared alleles NN vs CN = 0.074, p-value <0.05; Fig S4; Fig 3C ). There were too few individuals in the third CotC group to test for IBD within this group. However, we note that pairwise genetic distances between the three individuals in the CotC group (green squares: Fig. 3A ) were all smaller than pairwise genetic distances between each of these individuals and individuals in either the Calling or Non-Calling populations (red points: Fig. 3A ). Agreement between genetic boundaries and call morphology We found a significant association between call morphology and genetic assignment (Fisher’s exact test, p-value = 2.5e-07). All 14 individuals that were morphologically Calling (i.e. presence of one or more vocal slits) were geographically closest to sites that fell into the Calling genetic clade in the phylogenetic tree ( Table 1 ). Of the 25 individuals that were morphologically Non-Calling (i.e. absence of vocal slits), the majority (17) were geographically closest to individuals belonging to the Non-Calling genetic clade ( Table 1 ). Of the remaining eight morphologically Non-Calling individuals, four were from one site, which was geographically closest to an individual in the Calling genetic clade (i.e. mismatch between call morphology and genetic assignment). Notably, this site also included a morphologically Calling individual and was the only mixed phenotype site included in the analysis. This site falls along the genetic boundary between the Calling and Non-Calling groups and is an area of genetic admixture based on the STRUCTURE analysis ( Fig. 4 ). The remaining four morphologically Non-Calling individuals were geographically closest to the CotC genetic group ( Table 1 ). Download figure Open in new tab Figure 4. Analysis of niche overlap between Calling and Non-Calling populations of western toads in Canada. Comparison of continuous (A, B) and binary (C, D) predictions of the distribution of suitable habitat for the Calling (left panels) and Non-Calling (right panels) populations based on Maxent models. Black points in A and B represent the input occurrence records for each model. Warmer colors in these panels represent areas of higher predicted suitability; cooler colors represent areas of lower predicted suitability. Binary surfaces (insets in C and D) were generated from the continuous surfaces using the 10% omission threshold based on the corresponding occurrences, with green shading showing suitable habitat and light pink shading showing unsuitable habitat. A zoomed-in, arbitrary 400 x 650 km area around the boundary between genetic groups is shown in C and D to facilitate comparison of suitable habitat towards this boundary. Triangles in C and D show locations where individuals were assessed for call morphology, with black triangles showing sites where all individuals had vocal slits (Calling) and white triangles showing sites where all individuals lacked vocal slits (Non-Calling; data originally presented in Pauly 2008 ). Gray triangles show site along the boundary between groups where both phenotypes were found in the same location (see main text). Larger triangles represent locations used to test the concordance between genetic and morphological groups. Circles in the zoomed-in region show the genetic sites and the proportion of each individual’s ancestry assigned to the two populations based on the STRUCTURE analysis (see Fig. 2D ). In all surfaces, gray areas show places where conditions were outside of the range of conditions used to calibrate the models (i.e. untested climatic space) as identified through a Multivariate Environmental Similarity Surface (MESS) analysis. Panels E and F show the available (polygon outlines) and occupied (shading) climatic space for the Calling (orange) and Non-Calling (blue) populations respectively based on 31 climatic variables. Darker shading within the occupied climatic space represents the density of occurrence records for the Non-Calling (E) and Calling (F) population, with dark blue pixels in both panels showing the overlap between populations. View this table: View inline View popup Download powerpoint Table 1. The association between call morphology (presence or absence of vocal slits) and genetic group membership based on the geographically closest genetic site. Only morphologically-characterized specimens within 50 km of a genetic sample were used. Distribution of suitable habitat and niche overlap The Maxent models performed well at predicting presence based on spatial block cross-validation (Non-Calling: AUC = 0.73, Calling: AUC = 0.73; Table S2). Both populations responded similarly to most climatic variables (Fig. S5) and both models suggested that the number of degree-days below 18°C strongly influences the probability of presence (permutation importance 18.4% and 33.7% for Calling and Non-Calling models respectively), with the probability of presence decreasing with a greater number of cooler days in the year for both populations. However, most of the other variables deemed to be important for predicting presence by the models were different for the two populations (Table S3). For example, the variable with the highest percent contribution to the Non-Calling population model was the amount of precipitation as snow (44%), which showed a threshold-like effect in terms of probability of presence for this population, but did not have much of an influence on the probability of presence for the Calling population (Fig. S5). In contrast, the Julian date on which the frost-free period begins had the highest percent contribution to the model for the Calling population (38%) and showed a positive relationship with probability of presence with this population, but did little to explain the probability of presence for the Non-Calling population (Fig. S5). Areas of high climate suitability for the Non-Calling population were predicted throughout the interior plateau of BC, on Haida Gwaii, Vancouver Island, and in parts of Alberta ( Fig. 4A , 4C). Areas of high suitability for the Calling population were predicted in central Alberta, in parts of the CRMs, and in the center of BC ( Fig. 4B , 4D). There were parts of both BC and Alberta that were predicted to be suitable for both populations, including an area to the north of our transects along the BC-Alberta border ( Fig. 4A and B ). However, much of the southern CRMs and the general area associated with the genetic and phenotypic transition between the Calling and Non-Calling populations was predicted to be unsuitable for the Non-Calling population ( Fig. 4C ) and was outside the range of conditions experienced by the Calling population within its geographic range (i.e. “untested by this population”; Fig. 4D ). Both the available and the occupied climatic space for the Calling population within its range were almost entirely nested within the much broader, climatic space occupied by the Non-Calling population ( Fig. 4E , 4F). Within the set of conditions found in the Calling population’s range, almost all conditions occupied by the Calling population were also occupied by the Non-Calling population (niche stability index = 0.89). Schoener’s D was very low (D = 0.005) within the climatic space available to both populations and was not significantly greater than when drawing the Non-Calling population’s niche at random from the much larger set of conditions available to it within its range (p-value = 0.309; Fig. S6). In contrast, D was significantly greater than values obtained by randomizing the niche of the Calling population (p-value = 0.048; Fig. S7)—that is, the niche of the Calling population is significantly concentrated within that of the Non-Calling population within the limited environmental space available to both populations. Discussion We aimed to understand population structure across the Canadian range of western toads in Canada. Our results revealed two genetic groups corresponding to Calling and Non-Calling populations, as well as a third, currently undescribed group of western toads associated with the northern part of the Crown of the Continent Ecosystem. Our analysis of the distribution of variation in call morphology and suitable climate space for the Calling and Non-Calling populations respectively suggests that the genetic boundary between Calling and Non-Calling individuals aligns closely with a sharp transition in a major reproductive trait as well as gaps in suitable climate conditions that may limit the extent to which these groups interact, at least along southern parts of the CRMs. This work is timely as the status of western toads in Canada, including whether the two western toad populations should continue to be recognized as legal units (i.e. Designatable Units; DUs), is due for re-assessment by the Committee on the Status of Endangered Wildlife in Canada (COSEWIC; species are reassessed ∼10 years and the species was last assessed in 2012). We thus discuss our results in light of the specific criteria used by COSEWIC to guide these decisions. Support for the continued recognition of two populations of western toads in Canada In Canada, Designatable Units (DUs) are established by the Committee on the Status of Endangered Wildlife in Canada (COSEWIC) using two main criteria. The first criterion focusses on the “discreteness” of populations ( COSEWIC, 2023 ). Populations are considered to be discrete if there are heritable traits or markers that distinguish them and/or there are natural barriers to the transmission of heritable traits (i.e. barriers to gene flow: COSEWIC, 2023 ). Several of the results presented here speak to the discreteness of the Calling and Non-Calling populations. Our genomic results suggest that Calling and Non-Calling populations are distinct across much of their genome, with levels of genetic differentiation similar to, or even higher than, those observed among conservation units within other species (e.g. Moore et al., 2014 ; Jorge et al., 2022 ; Gallego-Garcia et al., 2023; MacDonald et al., 2025b ) and even among what are considered distinct species (e.g. Baumsteiger et al., 2017 ; Campbell et al., 2022 ). These genetic differences are associated with differences in terms of the presence of a pronounced advertisement call, which may influence conspecific attraction (e.g. Buxton et al., 2015 ) and thus opportunities for the transmission of genetic information. Finally, although our niche models and assessment of niche overlap suggest that the two populations have similar responses to variation in climatic conditions and that the occupied niche of the Calling population is entirely nested within that of the Non-Calling population, much of the area in and immediately adjacent to the CRMs is predicted to be unsuitable for the Non-Calling population and is outside the range of conditions experienced by the Calling population within its range. Along with the limited amount of genetic admixture observed in this region, these results suggest that the CRMs may act as at least a partial barrier to gene flow between the two populations. Together, these results support the continued recognition of the Calling and Non-Calling populations of western toads as discrete groups. Isolation by distance (IBD) is expected in widespread species that are continuously distributed across the landscape ( Wright, 1943 ; Turbek et al., 2023 ) and one challenge when assessing the discreteness of potential groups in widespread species is to ensure that such groups are not simply an artifact of IBD ( Chambers and Hillis, 2020 ; MacGuigan et al., 2022 ; Turbek et al., 2023 ). For this reason, we sampled individuals from the west coast, through to the eastern edge of the species’ range, ensuring that individuals on either side of the putative boundary between the groups were included. Although we observed strong IBD across the full extent of the study area and within the ranges of each of the populations, IBD alone is unlikely to account for the observed genetic structure. Specifically, genetic distances are larger on average between Calling and Non-Calling individuals than between individuals within either population at each distance class. Importantly, at relatively small geographic distances (<150 km), including between adjacent samples in and around the transition between groups, we found that mean pairwise genetic distance was significantly higher than expected when comparing samples across populations than when comparing samples within populations. Thus, while IBD plays a role in shaping patterns of genetic variation in western toads, Calling and Non-Calling populations show genetic differences above and beyond what is expected to arise from IBD alone. The second criteria used to delineate DUs in Canada focuses on the evolutionary significance of populations ( COSEWIC, 2023 ). This criterion often emphasizes the distinct evolutionary histories of populations ( COSEWIC, 2023 ). Our focus on the Canadian portion of the species’ range means our understanding of the biogeographic history of western toads remains incomplete. However, it is likely that western toads, like other species (e.g. Hoffman and Blouin, 2004 ; Klütsch et al., 2012 ; Cheng et al., 2014 ), have undergone periods of isolation in separate refugia during the Pleistocene glaciations. Goebel et al. (2009) previously suggested that western toads may have maintained refugial populations in the Klamath-Siskiyou Mountains, near the border of Oregon and California, and expanded north along coastal BC and Alaska. More recently, Gamlen-Greene (2022) hypothesized that Haida Gwaii may have potentially served as an additional western refugia for the species ( Gamlen-Greene, 2022 ). However, neither study included sites from further east in the species’ range in Canada. The eastern half of the western toad’s range overlaps with other species in eastern BC and Alberta that are thought to have persisted in ice-free areas of the Rocky Mountains (reviewed by Brunsfeld et al., 2001 ; Shafer et al., 2010 ) and interior of Idaho (e.g. Lait et al., 2012 ) during the Pleistocene glaciations, as well as in cryptic refugia within and between the Laurentide and Cordilleran ice sheets in Canada ( Shafer et al., 2010 ). Thus, it seems likely that western toads had access to multiple glacial refugia and that periods of isolation among different populations have shaped the observed genetic structure within the species’ range. Differences in adaptive traits are also considered by COSEWIC when evaluating the evolutionary significance of different groups (i.e. “DU-wide differences in adaptive traits not found elsewhere in Canada”: COSEWIC, 2023 ). The presence of a male vocal sac and pronounced advertisement call in a geographically-discrete part of the western toad’s range and absence of these traits elsewhere is highly unique amongst anurans ( Pauly, 2008 ; Elias-Costa and Faivovich, 2025 ). The adaptive significance of these differences in calling in western toads is unclear, but in other anurans, variation in the amount of calling that males do and the characteristics of their calls (e.g. frequency) is associated with differences in predation pressure (e.g. Dapper et al., 2010; Bonachea and Ryan, 2011 ) and the acoustic environment (e.g. Goutte et al., 2018 ; Zhao et al., 2021 ). Furthermore, differences in call behaviour and characteristics have promoted divergence in some anuran systems ( Pröhl, 2003 ; Funk et al., 2009 ) and are thought to have driven speciation in some cases ( Boul et al., 2007 ). Future work is needed to assess whether the differences in call behaviour of Calling and Non-Calling populations of western toads are adaptive and/or contribute to reproductive isolation between the populations. However, DU-wide differences in a key reproductive trait make it likely that Calling and Non-Calling populations are on separate evolutionary trajectories. Broader implications for the suturing of genetic diversity in western Canada The places where different species, or populations within species, meet and interact (i.e. “suture zones”) has long captured the interests of biogeographers and evolutionary ecologists ( Remington, 1968 ; Swenson and Howard, 2004 ; Swenson and Howard, 2005 ; Wait and Peñalba, 2025 ). The CRMs are a major geographic barrier that has partitioned biodiversity at multiple scales ( Malanson et al., 2018 ). For instance, the CRMs represent the transition between three ecozones: the Montane Cordillera in BC and the Boreal foothills and Prairies in Alberta. This transition defines the boundaries between different faunal provinces ( COSEWIC, 2023 ) and is where the western range limits of many species meet the eastern range limits of close congeners. The CRMs are also associated with genetic breaks in several taxa that maintain populations on both sides of the mountains (e.g. silky pocket mouse, North American red squirrel, little brown myotis: Jensen et al., 2024 ; Canadian Lynx: Watt et al., 2021 ; American goldfinch: Cloutier et al., 2024 ; Black-capped chickadees, Boreal chickadees: Lait et al., 2024 ). Our study adds the first amphibian species to this list, further suggesting that the CRMs are an important suture zone in the north. Our results also revealed a third genetic group located in southeastern BC and southwestern Alberta. A recent survey of genetic differences across the range of western toads in the United States similarly found a distinct genetic group associated with the Rocky Mountains in Glacier National Park, Montana (Colorado US Fish and Wildlife, 2017), which is also in the Crown of the Continent Ecosystem. Although that study did not include samples in Canada, the close geographic proximity of that group to the samples included in this study, makes it likely that these patterns are reflective of the same genetic substructure. Toads in this part of the Rocky Mountains may be unique with respect to their calling behaviour. Although the specimens examined by Pauly (2008) lacked vocal sacs (both those in Canada and in a broader survey which included sites in Montana), there are reports of toads in southeastern BC and southwestern Alberta producing consistent calls (per. communication L-A. Isaac, BC Ministry of Water, Land and Resource Stewardship; Waterton Lakes National Parks Ecological Monitoring program, Parks Canada Waterton Lakes Field Unit 2021), and Pauly (2008) reported males producing quieter, long, pulsed calls at a single location in northern Montana despite the absence of vocal sacs. Formal analysis of recorded calls suggested these vocalizations differ from the calls of individuals in Alberta and may be an “extra-long release call” ( Pauly, 2008 ) but additional acoustic and behavioural analyses are needed to clarify the nature and function of vocalizations in this part of the species’ range. The CotC group highlights the potential for the Crown of the Continent Ecosystem to harbour unique genotypes. The Crown of the Continent Ecosystem, spanning mountainous regions of the southwest corner of Alberta, southeast corner of BC, and northern Montana, is characterized by large elevational differences, includes a range of different ecosystems from alpine forests to prairie grasslands, and is one of the most biologically diverse places in North America ( Prato and Fagre, 2010 ). This region is thought to have harbored ice-free areas during the last glacial maximum and may have been an important northern refugia for western species ( Brunsfeld et al., 2001 ; Shafer et al., 2010 ). Although studies specifically testing the phylogeographic importance of this region are lacking, genetic differences between populations in the Crown of the Continent Ecosystem and elsewhere have been reported in other groups (e.g. Long-toed salamanders: Lee-Yaw and Irwin, 2012 ; Half-moon hairstreak butterflies: Macdonald et al., 2025b ). Additional studies focusing on this region would enhance our understanding of the extent to which the Crown of the Continent Ecosystem is a hotspot of genetic diversity. Conclusions and Future Directions This study is one of the first studies to evaluate genomic differences in a widespread amphibian species in western Canada. We found evidence to support distinct genetic groups associated with Calling and Non-Calling populations of western toads in Canada. Paired with phenotypic and ecological differences between these populations, the continued recognition of these populations as separate conservation units (i.e. DUs) seems warranted. At the same time, we note that we did not have samples from the most northern parts of the species range. Importantly, our niche models revealed areas of high climatic suitability for both populations to the east of the CRMs in areas north of our sampling. Although there are gaps in available morphological information from this part of the species’ range, there is at least one known location in the Grande Prairie region of Alberta where both phenotypes are found ( Fig. 4 C,D; Pauly 2008 ). At this latitude, the CRMs are less steep, valleys are wider, and the foothills extend further east than in southern parts of the CRMs. Thus, it’s possible that there are more opportunities for contact between the two populations in this region. Additional samples from this region and from across more northern parts of the species’ range are thus necessary to fully characterize the genetic boundaries between populations. An outstanding question at this point is the extent to which the different genetic groups observed are reproductively isolated. Fine-scale sampling across the boundaries between groups to quantify levels of introgression and hybridization, paired with tests of the role of differences in calling behaviour in mate recognition would shed light on this question. Comparison of levels of reproductive isolation on either side of the southern CRMs to levels in those northern areas where climatic suitability is predicted to be high for both the Calling and Non-Calling populations would be particularly interesting, providing the ultimate test of whether any barriers to gene flow established along putatively different colonization routes around a geographic barrier are maintained upon secondary contact (i.e. akin to Moritz et al., 1992 ; Helbig, 2005 ). Regardless of whether Calling and Non-Calling populations are reproductively isolated, the presence of a highly distinct genetic group that is geographically contained entirely within Canada (i.e. the Calling population) is an unusual situation (see MacDonald et al. 2025a for a similarly “curious” example). From the perspective of both understanding and managing genetic diversity, our results thus highlight the importance of filling gaps in our assessments of intraspecific genetic variation in previously glaciated areas. Data Accessibility All code is available in an online GitHub repository ( https://github.com/Jaynabergman/WETO_ddRAD ) Raw reads from all samples will be available upon publication on the Short Read Archive except for the sample from Haida Gwaii. Data from the Haida Gwaii sample are available on request from the BC Ministry of Forests, subject to review by the Council of the Haida Nation. Author Contributions Jayna Bergman : Data curation, Formal analysis, Validation, Software, Visualization, Writing – Original Draft. Juan Enciso-Romero : Software, Formal analysis. Gregory B. Pauly : Data curation, Writing – Review & Editing. Roseanna Gamlen-Greene : Data curation, Writing – Review & Editing. Melissa Todd : Data curation, Writing – Review & Editing. Julie A. Lee-Yaw : Conceptualization, Methodology, Project administration, Funding acquisition, Supervision, Writing – Review & Editing Data Availability All code is available in an online GitHub repository ( https://github.com/Jaynabergman/WETO_ddRAD ) Raw reads from all samples will be available upon publication on the Short Read Archive except for the sample from Haida Gwaii. Data from the Haida Gwaii sample are available on request from the BC Ministry of Forests, subject to review by the Council of the Haida Nation. Funding Statement This work was supported by a NSERC Canada Graduate Scholarship to JCB, Ontario Graduate Scholarship to JCB, Alberta Conservation Association grant in biodiversity to JAL (ACA-030-00-90-321), and a NSERC Discovery Grant to JAL (RGPIN-2020-04611). The British Columbia Ministry of Forests Research Program supported the collection of some of the samples in British Columbia. Conflict of interest statement The authors have no conflict of interest to report. Acknowledgements The Western Toads included in this study are found on the traditional Lands of many First Nations and we recognize the ongoing stewardship of the Indigenous Peoples who have cared for this species and their habitats. We thank the following for their support of the field efforts and sample collection: S. Marks, H. Meikle, P. Monterio, B. Cadsand, E. Wind, M. Thomspon, G. Skinner, G. Dare, C. Gill, K. Pearson, T. Einfeldt, and C. Bergeron. We wish to acknowledge C. Cullingham and J. Martin for comments on an earlier thesis version of this manuscript, and N. Kester for discussion on the niche overlap analysis. Funder Information Declared Natural Sciences and Engineering Research Council Alberta Conservation Association, https://ror.org/030rn8443 British Columbia Ministry of Forests Research Program Footnotes https://github.com/Jaynabergman/WETO_ddRAD References ↵ Addis , B. R. , Lowe , W. H. , Hossack , B. R. , & Allendorf , F. W . ( 2015 ). Population genetic structure and disease in montane boreal toads: More heterozygous individuals are more likely to be infected with amphibian chytrid . Conservation Genetics , 16 ( 4 ), 833 – 844 . doi: 10.1007/s10592-015-0704-6 OpenUrl CrossRef Baddeley A , Turner R ( 2005 ). “ spatstat: An R Package for Analyzing Spatial Point Patterns. _Journal of Statistical Software_ , *12*(6), 1-42. doi: 10.18637/jss.v012.i06 . OpenUrl CrossRef PubMed ↵ Baker , D. J. , Maclean , I. M. D. , & Gaston , K. J . ( 2024 ). Effective strategies for correcting spatial sampling bias in species distribution models without independent test data . Diversity and Distributions , 30 ( 3 ), e13802 . doi: 10.1111/ddi.13802 OpenUrl CrossRef ↵ Ballard , J. W. O. , & Whitlock , M. C . ( 2004 ). The incomplete natural history of mitochondria . Molecular Ecology , 13 ( 4 ), 729 – 744 . doi: 10.1046/j.1365-294X.2003.02063.x OpenUrl CrossRef PubMed Web of Science ↵ Barber , R. A. , Ball , S. G. , Morris , R. K. A. , & Gilbert , F . ( 2022 ). Target-group backgrounds prove effective at correcting sampling bias in Maxent models . Diversity and Distributions , 28 ( 1 ), 128 – 141 . doi: 10.1111/ddi.13442 OpenUrl CrossRef ↵ Barbosa , S. , Mestre , F. , White , T. A. , Paupério , J. , Alves , P. C. , & Searle , J. B . ( 2018 ). Integrative approaches to guide conservation decisions: Using genomics to define conservation units and functional corridors . Molecular Ecology , 27 ( 17 ), 3452 – 3465 . doi: 10.1111/mec.14806 OpenUrl CrossRef ↵ Baumsteiger , J. , Moyle , P. B. , Aguilar , A. , O’Rourke , S. M. , & Miller , M. R . ( 2017 ). Genomics clarifies taxonomic boundaries in a difficult species complex . PLoS ONE , 12 ( 12 ). doi: 10.1371/journal.pone.0189417 OpenUrl CrossRef ↵ Beebee , T. J. C . ( 2005 ). Conservation genetics of amphibians . Heredity , 95 ( 6 ), 423 – 427 . doi: 10.1038/sj.hdy.6800736 OpenUrl CrossRef PubMed Web of Science ↵ Benefer , C. M. , Van Herk , W. G. , Ellis , J. S. , Blackshaw , R. P. , Vernon , R. S. , & Knight , M. E. ( 2013 ). The molecular identification and genetic diversity of economically important wireworm species (Coleoptera: Elateridae) in Canada . Journal of Pest Science , 86 ( 1 ), 19 – 27 . doi: 10.1007/s10340-012-0454-x OpenUrl CrossRef ↵ Bergman , J. C. , Finn , K. J. , & Lee-Yaw , J. A . ( 2024 ). Study extent influences the predictions and performance of species distribution models: A case study of six amphibian species at the edge of their geographic distributions in western Canada . Biodiversity and Conservation , 33 ( 1 ), 4295 – 4318 . doi: 10.1007/s10531-024-02953-3 OpenUrl CrossRef ↵ Bickford , D. , Lohman , D. J. , Sodhi , N. S. , Ng , P. K. L. , Meier , R. , Winker , K. , Ingram , K. K. , & Das , I . ( 2007 ). Cryptic species as a window on diversity and conservation . Trends in Ecology & Evolution , 22 ( 3 ), 148 – 155 . doi: 10.1016/j.tree.2006.11.004 OpenUrl CrossRef PubMed Web of Science ↵ Bonachea , L. A. , & Ryan , M. J . ( 2011 ). Simulated Predation Risk Influences Female Choice in Túngara Frogs, Physalaemus pustulosus . Ethology , 117 ( 5 ), 400 – 407 . doi: 10.1111/j.1439-0310.2011.01889.x OpenUrl CrossRef ↵ Boul , K. E. , Chris Funk , W. , Darst , C. R. , Cannatella , D. C. , & Ryan , M. J . ( 2007 ). Sexual selection drives speciation in an Amazonian frog . Proceedings of the Royal Society B: Biological Sciences , 274 ( 1608 ), 399 – 406 . doi: 10.1098/rspb.2006.3736 OpenUrl CrossRef PubMed Web of Science ↵ Bradie , J. , & Leung , B . ( 2017 ). A quantitative synthesis of the importance of variables used in MaxEnt species distribution models . Journal of Biogeography , 44 ( 6 ), 1344 – 1361 . doi: 10.1111/jbi.12894 OpenUrl CrossRef ↵ Broennimann , O. , Fitzpatrick , M. C. , Pearman , P. B. , Petitpierre , B. , Pellissier , L. , Yoccoz , N. G. , Thuiller , W. , Fortin , M.-J. , Randin , C. , Zimmermann , N. E. , Graham , C. H. , & Guisan , A . ( 2012 ). Measuring ecological niche overlap from occurrence and spatial environmental data: Measuring niche overlap . Global Ecology and Biogeography , 21 ( 4 ), 481 – 497 . doi: 10.1111/j.1466-8238.2011.00698.x OpenUrl CrossRef ↵ J Silvertown , J Antonovics Brunsfeld , S. J. , Sullivan , J. , Soltis , D. E. , & Soltis , P. S . ( 2001 ). Comparative phylogeography of North-western North America: a synthesis . In: J Silvertown , J Antonovics , eds. Integrating Ecological and Evolutionary Process in a Spatial Context . Oxford UK : Blackwell Science , 319 – 339 . ↵ Broennimann O , Di Cola V , Guisan A. ( 2025 ). _ecospat: Spatial Ecology Miscellaneous Methods_. R package version 4.1.2 , ↵ Bull , E. L . ( 2006 ). Sexual Differences in the Ecology and Habitat Selection of Western Toads ( Bufo boreas ) in Northeastern Oregon . Herpetological Conservation and Biology 1 ( 1 ), 27 – 38 . OpenUrl ↵ Buxton , V. L. , Ward , M. P. , & Sperry , J. H . ( 2015 ). Use of chorus sounds for location of breeding habitat in 2 species of anuran amphibians . Behavioral Ecology , 26 ( 4 ), 1111 – 1118 . doi: 10.1093/beheco/arv059 OpenUrl CrossRef ↵ Cairns , N. A. , Cicchino , A. S. , Stewart , K. A. , Austin , J. D. , & Lougheed , S. C . ( 2021 ). Cytonuclear discordance, reticulation and cryptic diversity in one of North America’s most common frogs . Molecular Phylogenetics and Evolution , 156 , 107042 . doi: 10.1016/j.ympev.2020.107042 OpenUrl CrossRef PubMed ↵ Campbell , E. O. , MacDonald , Z. G. , Gage , E. V. , Gage , R. V. , & Sperling , F. A. H . ( 2022 ). Genomics and ecological modelling clarify species integrity in a confusing group of butterflies . Molecular Ecology , 31 ( 8 ), 2400 – 2417 . doi: 10.1111/mec.16407 OpenUrl CrossRef ↵ Catchen , J. , Hohenlohe , P. A. , Bassham , S. , Amores , A. , & Cresko , W. A . ( 2013 ). Stacks: An analysis tool set for population genomics . Molecular Ecology , 22 ( 11 ), 3124 – 3140 . doi: 10.1111/mec.12354 OpenUrl CrossRef PubMed Web of Science ↵ Chambers , E. A. , Lara-Tufiño , J. D. , Campillo-García , G. , Cisneros-Bernal , A. Y. , Dudek , D. J. , León-Règagnon , V. , Townsend , J. H. , Flores-Villela , O. , & Hillis , D. M . ( 2025 ). Distinguishing species boundaries from geographic variation . Proceedings of the National Academy of Sciences , 122 ( 19 ), e2423688122 . doi: 10.1073/pnas.2423688122 OpenUrl CrossRef PubMed ↵ Chambers , E. A. , Marshall , T. L. , & Hillis , D. M . ( 2022 ). The Importance of Contact Zones for Distinguishing Interspecific from Intraspecific Geographic Variation. Systematic Biology , syac 056 . doi: 10.1093/sysbio/syac056 OpenUrl CrossRef PubMed ↵ Chambers , E. , & Hillis , D . ( 2020 ). The Multispecies Coalescent Over-Splits Species in the Case of Geographically Widespread Taxa . Systematic Biology , 69 ( 1 ), 184 – 193 . doi: 10.1093/sysbio/syz042 OpenUrl CrossRef ↵ Chang , C. C. , Chow , C. C. , Tellier , L. C. , Vattikuti , S. , Purcell , S. M. , & Lee , J. J . ( 2015 ). Second-generation PLINK: Rising to the challenge of larger and richer datasets . Gigascience , 4 ( 1 ), s13742 – 015 -0047–0048. doi: 10.1186/s13742-015-0047-8 OpenUrl CrossRef ↵ Chen , S . ( 2023 ). Ultrafast one-pass FASTQ data preprocessing, quality control, and deduplication using fastp . iMeta , 2 . doi: 10.1002/imt2.107 OpenUrl CrossRef ↵ Cheng , E. , Hodges , K. E. , Melo-Ferreira , J. , Alves , P. C. , & Mills , L. S. ( 2014 ). Conservation implications of the evolutionary history and genetic diversity hotspots of the snowshoe hare . Molecular Ecology , 23 ( 12 ), 2929 – 2942 . doi: 10.1111/mec.12790 OpenUrl CrossRef PubMed ↵ Cloutier , A. , Chan , D. T. C. , Poon , E. S. K. , & Sin , S. Y. W . ( 2024 ). The genetic consequences of historic climate change on the contemporary population structure of a widespread temperate North American songbird . Molecular Phylogenetics and Evolution , 201 , 108216 . doi: 10.1016/j.ympev.2024.108216 OpenUrl CrossRef PubMed ↵ Coates , D. J. , Byrne , M. , & Moritz , C . ( 2018 ). Genetic Diversity and Conservation Units: Dealing With the Species-Population Continuum in the Age of Genomics . Frontiers in Ecology and Evolution , 6 , 165 . doi: 10.3389/fevo.2018.00165 OpenUrl CrossRef ↵ COSEWIC . ( 2012 ). Western toad (Anaxyrus boreas): COSEWIC assessment and status report. Committee on the Status of Endangered Wildlife in Canada . ↵ COSEWIC . ( 2023 ). COSEWIC guidelines for recognizing designatable units . https://cosewic.ca/index.php/en/reports/preparing-status-reports/guidelines-recognizing-designatable-units.html Dapper , A. L. , Baugh , A. T. , & Ryan , M. J . ( 2011 ). The Sounds of Silence as an Alarm Cue in Túngara Frogs, Physalaemus pustulosus . Biotropica , 43 ( 3 ), 380 – 385 . doi: 10.1111/j.1744-7429.2010.00707.x OpenUrl CrossRef ↵ DeRaad D. ( 2023 ). _ SNPfiltR: Interactively Filter SNP Datasets_. R package version 1.0.1 , . Dincă , V. , Lee , K. M. , Vila , R. , & Mutanen , M . ( 2019 ). The conundrum of species delimitation: A genomic perspective on a mitogenetically super-variable butterfly . Proceedings of the Royal Society B: Biological Sciences , 286 ( 1911 ), 20191311 . doi: 10.1098/rspb.2019.1311 OpenUrl CrossRef PubMed ↵ Elias-Costa , A. J. , & Faivovich , J . ( 2025 ). Evolution of Vocal Sacs in Anura . Bulletin of the American Museum of Natural History , 2025 ( 470 ). doi: 10.1206/0003-0090.470.1.1 OpenUrl CrossRef ↵ Elith , J. , Kearney , M. R. , & Phillips , S. J . ( 2010 ). The art of modelling range-shifting species . Methods in Ecology and Evolution , 1 , 330 – 342 . OpenUrl CrossRef ↵ Elith , J. , Phillips , S. J. , Hastie , T. , Dudık , M. , Chee , Y. E. , & Yates , C. J . ( 2011 ). A statistical explanation of MaxEnt for ecologists . Diversity and Distributions , 15 . ↵ Environment and Climate Change Canada . ( 2020 ). Management Plan for the Western Toad (Anaxyrus boreas), Calling and Non-calling Populations, in Canada ( Species at Risk Act Management Plan Series ; v + 39 pp.). Ottawa, ON : Environment and Climate Change Canada . ↵ Evanno , G. , Regnaut , S. , & Goudet , J . ( 2005 ). Detecting the number of clusters of individuals using the software structure: A simulation study . Molecular Ecology , 14 ( 8 ), 2611 – 2620 . doi: 10.1111/j.1365-294X.2005.02553.x OpenUrl CrossRef PubMed Web of Science ↵ Everson , K. M. , Gray , L. N. , Jones , A. G. , Lawrence , N. M. , Foley , M. E. , Sovacool , K. L. , Kratovil , J. D. , Hotaling , S. , Hime , P. M. , Storfer , A. , Parra-Olea , G. , Percino-Daniel , R. , Aguilar-Miguel , X. , O’Neill , E. M. , Zambrano , L. , Shaffer , H. B. , & Weisrock , D. W . ( 2021 ). Geography is more important than life history in the recent diversification of the tiger salamander complex . Proceedings of the National Academy of Sciences , 118 ( 17 ), e2014719118 . doi: 10.1073/pnas.2014719118 OpenUrl Abstract / FREE Full Text ↵ Ewels , P. , Magnusson , M. , Lundin , S. , & Kaller , M . ( 2016 ). MultiQC: summarize analysis results for multiple tools and samples in a single report . Bioinformatics , 32 ( 19 ), 3047 – 3048 . OpenUrl CrossRef PubMed Fišer , C. , Robinson , C. T. , & Malard , F . ( 2017 ). Cryptic species as a window into the paradigm shift of the species concept . Molecular Ecology , 27 ( 3 ), 613 – 635 . doi: 10.1111/mec.14486 OpenUrl CrossRef ↵ Forester , B. R. , Murphy , M. A. , Mellison , C. , Petersen , J. , Pilliod , D. S. , Van Horne , R. , Harvey , J. , & Funk , W. C. ( 2022 ). Genomics-informed delineation of conservation units in a desert amphibian . Molecular Ecology , 1 – 21 . ↵ Fu , J. , & Wen , L . ( 2023 ). Impacts of Quaternary glaciation, geological history and geography on animal species history in continental East Asia: A phylogeographic review . Molecular Ecology , 32 ( 16 ), 4497 – 4514 . doi: 10.1111/mec.17053 OpenUrl CrossRef ↵ Funk , W. C. , Cannatella , D. C. , & Ryan , M. J . ( 2009 ). Genetic divergence is more tightly related to call variation than landscape features in the Amazonian frogs Physalaemus petersi and P. freibergi . Journal of Evolutionary Biology , 22 ( 9 ), 1839 – 1853 . doi: 10.1111/j.1420-9101.2009.01795.x OpenUrl CrossRef PubMed ↵ Funk , W. C. , Forester , B. R. , Converse , S. J. , Darst , C. , & Morey , S . ( 2019 ). Improving conservation policy with genomics: A guide to integrating adaptive potential into U.S. Endangered Species Act decisions for conservation practitioners and geneticists . Conservation Genetics , 20 ( 1 ), 115 – 134 . doi: 10.1007/s10592-018-1096-1 OpenUrl CrossRef ↵ Funk , W. C. , McKay , J. K. , Hohenlohe , P. A. , & Allendorf , F. W . ( 2012 ). Harnessing genomics for delineating conservation units . Trends in Ecology & Evolution , 27 ( 9 ), 489 – 496 . doi: 10.1016/j.tree.2012.05.012 OpenUrl CrossRef PubMed Web of Science ↵ Gallego-García , N. , Caballero , S. , & Shaffer , H. B . ( 2021 ). Are Genomic Updates of Well-Studied Species Worth the Investment for Conservation? A Case Study of the Critically Endangered Magdalena River Turtle . Journal of Heredity , 112 ( 7 ), 575 – 589 . doi: 10.1093/jhered/esab063 OpenUrl CrossRef PubMed Gallego-García , N. , Vargas-Ramírez , M. , & Shaffer , H. B . ( 2023 ). The importance of cryptic diversity in the conservation of wide-ranging species: The red-footed tortoise Chelonoidis carbonarius in Colombia . Molecular Ecology , 32 , 4531 – 4545 . Scopus. doi: 10.1111/mec.17052 OpenUrl CrossRef ↵ Gamlen-Greene , R . ( 2022 ). The ecology, distribution and population genetics of amphibians on Haida Gwaii, British Columbia [PhD] . The University of British Columbia . García-Rodríguez , O. , Hardouin , E. A. , Pedreschi , D. , Richards , M. B. , Stafford , R. , Searle , J. B. , & Stewart , J. R . ( 2024 ). Contrasting Patterns of Genetic Diversity in European Mammals in the Context of Glacial Refugia . Diversity , 16 ( 10 ), 611 . doi: 10.3390/d16100611 OpenUrl CrossRef ↵ Glor , R. E. , & Warren , D . ( 2011 ). Testing ecological explanations for biogeographic boundaries . Evolution , 65 ( 3 ), 673 – 683 . doi: 10.1111/j.1558-5646.2010.01177.x OpenUrl CrossRef PubMed ↵ Goebel , A. M. , Ranker , T. A. , Corn , P. S. , & Olmstead , R. G . ( 2009 ). Mitochondrial DNA evolution in the Anaxyrus boreas species group . Molecular Phylogenetics and Evolution , 50 ( 2 ), 209 – 225 . doi: 10.1016/j.ympev.2008.06.019 OpenUrl CrossRef PubMed Web of Science ↵ Gordon , M. R. ( 2017 ). Three new bufonid (Bufo (Anaxyrus)) species discovered within the Great Basin and the consequences of taxonomic crypsis [MSc] . University of Nevada . ↵ Gordon , M. R. , Simandle , E. T. , Sandmeier , F. C. , & Tracy , C. R . ( 2020 ). Two New Cryptic Endemic Toads of Bufo Discovered in Central Nevada, Western United States (Amphibia: Bufonidae : Bufo [ Anaxyrus ]) . Copeia , 108 ( 1 ), 166 . doi: 10.1643/CH-18-086 OpenUrl CrossRef ↵ Goutte , S. , Dubois , A. , Howard , S. D. , Márquez , R. , Rowley , J. J. L. , Dehling , J. M. , Grandcolas , P. , Xiong , R. C. , & Legendre , F . ( 2018 ). How the environment shapes animal signals: A test of the acoustic adaptation hypothesis in frogs . Journal of Evolutionary Biology , 31 ( 1 ), 148 – 158 . doi: 10.1111/jeb.13210 OpenUrl CrossRef PubMed ↵ Graham , B. A. , & Burg , T. M . ( 2012 ). Molecular markers provide insights into contemporary and historic gene flow for a non-migratory species . Journal of Avian Biology , 43 ( 3 ), 198 – 214 . doi: 10.1111/j.1600-048X.2012.05604.x OpenUrl CrossRef ↵ Guisan , A. , Petitpierre , B. , Broennimann , O. , Daehler , C. , & Kueffer , C . ( 2014 ). Unifying niche shift studies: Insights from biological invasions . Trends in Ecology & Evolution , 29 ( 5 ), 260 – 269 . doi: 10.1016/j.tree.2014.02.009 OpenUrl CrossRef PubMed Web of Science ↵ Helbig , A. J . ( 2005 ). Evolutionary genetics: A ring of species . Heredity , 95 , 113 – 114 . OpenUrl CrossRef PubMed ↵ Hewitt , G . ( 2000 ). The genetic legacy of the Quaternary ice ages . Nature , 405 ( 6789 ), 907 – 913 . doi: 10.1038/35016000 OpenUrl CrossRef PubMed Web of Science ↵ Hewitt , G. M . ( 2004 ). Genetic consequences of climatic oscillations in the Quaternary . Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences , 359 ( 1442 ), 183 – 195 . doi: 10.1098/rstb.2003.1388 OpenUrl CrossRef PubMed Hijmans R , ( 2024 ). _geosphere: Spherical Trigonometry_. R package version 1.5-20 , . ↵ Hijmans R , ( 2024b ). _predicts: Spatial Prediction Tools_. R package version 0.1-17 , . ↵ Hoffman , E. A. , & Blouin , M. S . ( 2004 ). Evolutionary history of the northern leopard frog: Reconstruction of phylogeny, phylogeography, and historical changes in population demography from mitochondrial DNA . Evolution , 58 ( 1 ), 145 – 159 . doi: 10.1111/j.0014-3820.2004.tb01581.x OpenUrl CrossRef PubMed Web of Science ↵ Jensen , A. J. , Cove , M. V. , Goldstein , B. R. , Kays , R. , McShea , W. , Pacifici , K. , Rooney , B. , & Kierepka , E . ( 2024 ). Geographic barriers but not life history traits shape the phylogeography of North American mammals . Global Ecology and Biogeography , 33 ( 8 ), e13875 . doi: 10.1111/geb.13875 OpenUrl CrossRef ↵ Jombart , T . ( 2008 ). adegenet: A R package for the multivariate analysis of genetic markers . Bioinformatics , 24 ( 11 ), 1403 – 1405 . OpenUrl CrossRef PubMed Web of Science ↵ Jorge , R. F. , Lima , A. P. , & Stow , A. J . ( 2022 ). Selection and localised genetic structure in the threatened Manauense Harlequin Frog ( Bufonidae: Atelopus manauensis ) . Conservation Genetics , 23 ( 3 ), 559 – 574 . Scopus. doi: 10.1007/s10592-022-01436-8 OpenUrl CrossRef ↵ Kalyaanamoorthy , S. , Minh , B. Q. , Wong , T. K. F. , Von Haeseler , A. , & Jermiin , L. S. ( 2017 ). ModelFinder: Fast model selection for accurate phylogenetic estimates . Nature Methods , 14 ( 6 ), 587 – 589 . doi: 10.1038/nmeth.4285 OpenUrl CrossRef PubMed ↵ Kass , J. M. , Muscarella , R. , Galante , P. J. , Bohl , C. L. , Pinilla-Buitrago , G. E. , Boria , R. A. , Soley-Guardia , M. , & Anderson , R. P. ( 2021 ). ENMeval 2.0: Redesigned for customizable and reproducible modeling of species’ niches and distributions . Methods in Ecology and Evolution , 12 ( 9 ), 1602 – 1608 . doi: 10.1111/2041-210X.13628 OpenUrl CrossRef ↵ Klütsch , C. F. C. , Manseau , M. , & Wilson , P. J . ( 2012 ). Phylogeographical Analysis of mtDNA Data Indicates Postglacial Expansion from Multiple Glacial Refugia in Woodland Caribou ( Rangifer tarandus caribou ) . PLoS ONE , 7 ( 12 ), e52661 . doi: 10.1371/journal.pone.0052661 OpenUrl CrossRef PubMed ↵ Kopelman , N. M. , Mayzel , J. , Jakobsson , M. , Rosenberg , N. A. , & Mayrose , I . ( 2015 ). Clumpak: A program for identifying clustering modes and packaging population structure inferences across K . Molecular Ecology Resources , 15 ( 5 ), 1179 – 1191 . doi: 10.1111/1755-0998.12387 OpenUrl CrossRef PubMed ↵ Lait , L. A. , Enciso-Romero , J. , Lekamlage , T. T. M. , Veale , A. , Abeyrama , D. K. , & Burg , T. M . ( 2024 ). RADseq data reveal widespread historical introgression in four familiar North American songbirds . Evolution , 78 ( 5 ), 860 – 878 . doi: 10.1093/evolut/qpae014 OpenUrl CrossRef ↵ Lait , L. A. , Friesen , Vicki. L. , Gaston , A. J. , & Burg , T. M. ( 2012 ). The post-Pleistocene population genetic structure of a western North American passerine: The chestnut-backed chickadee Poecile rufescens . Journal of Avian Biology , 43 , 541 – 552 . OpenUrl CrossRef ↵ Lee-Yaw , J. A. , & Irwin , D. E . ( 2012 ). Large geographic range size reflects a patchwork of divergent lineages in the long-toed salamander ( Ambystoma macrodactylum ) . Journal of Evolutionary Biology , 25 ( 11 ), 2276 – 2287 . doi: 10.1111/j.1420-9101.2012.02604.x OpenUrl CrossRef PubMed ↵ Lee-Yaw , J. A. , & Irwin , D. E . ( 2015 ). The importance (or lack thereof) of niche divergence to the maintenance of a northern species complex: The case of the long-toed salamander ( Ambystoma macrodactylum ) . Journal of Evolutionary Biology , 28 ( 4 ), 917 – 930 . doi: 10.1111/jeb.12619 OpenUrl CrossRef PubMed ↵ Lee-Yaw , J. A. , Irwin , J. T. , & Green , D. M. ( 2008 ). Postglacial range expansion from northern refugia by the wood frog, Rana sylvatica . Molecular Ecology , 17 ( 3 ), 867 – 884 . doi: 10.1111/j.1365-294X.2007.03611.x OpenUrl CrossRef PubMed ↵ Lehnert , S. J. , Bradbury , I. R. , Wringe , B. F. , Van Wyngaarden , M. , & Bentzen , P. ( 2023 ). Multifaceted framework for defining conservation units: An example from Atlantic salmon ( Salmo salar ) in Canada . Evolutionary Applications , 16 ( 9 ), 1568 – 1585 . doi: 10.1111/eva.13587 OpenUrl CrossRef ↵ Lesbarrères , D. , Ashpole , S. L. , Bishop , C. A. , Blouin-Demers , G. , Brooks , R. J. , Echaubard , P. , Govindarajulu , P. , Green , D. M. , Hecnar , S. J. , Herman , T. , Houlahan , J. , Litzgus , J. D. , Mazerolle , M. J. , Paszkowski , C. A. , Rutherford , P. , Schock , D. M. , Storey , K. B. , & Lougheed , S. C . ( 2014 ). Conservation of herpetofauna in northern landscapes: Threats and challenges from a Canadian perspective . Biological Conservation , 170 , 48 – 55 . doi: 10.1016/j.biocon.2013.12.030 OpenUrl CrossRef ↵ Li , H . ( 2011 ). A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data . Bioinformatics , 27 ( 21 ), 2987 – 2993 . doi: 10.1093/bioinformatics/btr509 OpenUrl CrossRef PubMed Web of Science ↵ Li , H. , & Durbin , R . ( 2009 ). Fast and accurate short read alignment with Burrows–Wheeler transform . Bioinformatics , 25 ( 14 ), 1754 – 1760 . doi: 10.1093/bioinformatics/btp324 OpenUrl CrossRef PubMed Web of Science ↵ Lyman , R. A. , & Edwards , C. E . ( 2022 ). Revisiting the comparative phylogeography of unglaciated eastern North America: 15 years of patterns and progress . Ecology and Evolution , 12 ( 4 ). doi: 10.1002/ece3.8827 OpenUrl CrossRef ↵ MacDonald , Z. G. , Dupuis , J. R. , Glasier , J. R. N. , Sissons , R. , Moehrenschlager , A. , Shaffer , H. B. , & Sperling , F. A. H . ( 2025a ). Genomic and ecological divergence support recognition of a new species of endangered Satyrium butterfly (Lepidoptera, Lycaenidae) . ZooKeys , 1234 , 291 – 307 . doi: 10.3897/zookeys.1234.143893 OpenUrl CrossRef PubMed ↵ MacDonald , Z. G. , Dupuis , J. R. , Glasier , J. R. N. , Sissons , R. , Moehrenschlager , A. , Shaffer , H. B. , & Sperling , F. A. H . ( 2025b ). Whole-Genome Evaluation of Genetic Rescue: The Case of a Curiously Isolated and Endangered Butterfly . Molecular Ecology , 34 ( 4 ). OpenUrl ↵ MacGuigan , D. J. , Mount , G. G. , Watkins-Colwell , G. J. , Near , T. J. , & Lambert , M. R . ( 2022 ). Genomic Data Clarify Aquarana Systematics and Reveal Isolation-by-Distance Dominates Phylogeography of the Wide-Ranging Frog Rana clamitans . Ichthyology & Herpetology , 110 ( 3 ). doi: 10.1643/h2021129 OpenUrl CrossRef ↵ Mahony , C. R. , Wang , T. , Hamann , A. , & Cannon , A. J . ( 2022 ). A global climate model ensemble for downscaled monthly climate normals over North America . International Journal of Climatology , 42 ( 11 ), 5871 – 5891 . doi: 10.1002/joc.7566 OpenUrl CrossRef ↵ Malanson , G. P. , Fagre , D. B. , & Zimmerman , D. L . ( 2018 ). Scale dependence of diversity in alpine tundra, Rocky Mountains, USA . Plant Ecology , 219 ( 8 ), 999 – 1008 . doi: 10.1007/s11258-018-0852-0 OpenUrl CrossRef ↵ Mijangos , J. L. , Gruber , B. , Berry , O. , Pacioni , C. , & Georges , A . ( 2022 ). dartR v2: An accessible genetic analysis platform for conservation, ecology and agriculture . Methods in Ecology and Evolution , 13 ( 10 ), 2150 – 2158 . doi: 10.1111/2041-210X.13918 OpenUrl CrossRef ↵ Miller , C. V. , Bossu , C. M. , Sarraco , J. F. , Toews , D. P. L. , Rushing , C. S. , Roberto-Charron , A. , Tremblay , J. A. , Chandler , R. B. , DeSaix , M. G. , Fiss , C. J. , Larkin , J. L. , Haché , S. , Nebel , S. , & Ruegg , K. C . ( 2024 ). Genomics-informed conservation units reveal spatial variation in climate vulnerability in a migratory bird . Molecular Ecology , 33 ( 1 ). Scopus. doi: 10.1111/mec.17199 OpenUrl CrossRef ↵ Minh , B. Q. , Schmidt , H. A. , Chernomor , O. , Schrempf , D. , Woodhams , M. D. , Von Haeseler , A. , & Lanfear , R. ( 2020 ). IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era . Molecular Biology and Evolution , 37 ( 5 ), 1530 – 1534 . doi: 10.1093/molbev/msaa015 OpenUrl CrossRef PubMed ↵ Moore , J.-S. , Bourret , V. , Dionne , M. , Bradbury , I. , O’Reilly , P. , Kent , M. , Chaput , G. , & Bernatchez , L . ( 2014 ). Conservation genomics of anadromous Atlantic salmon across its North American range: Outlier loci identify the same patterns of population structure as neutral loci . Molecular Ecology , 23 ( 23 ), 5680 – 5697 . Scopus. doi: 10.1111/mec.12972 OpenUrl CrossRef ↵ Moritz , C. , Schneider , C. J. , & Wake , D. B . ( 1992 ). Evolutionary relationships within the Ensatina eschscholtzii complex confirm the ring species interpretation . Systematic Biology , 41 ( 3 ), 273 – 291 . OpenUrl CrossRef Web of Science ↵ Naimi , B. , Hamm , N. A. S. , Groen , T. A. , Skidmore , A. K. , & Toxopeus , A. G . ( 2014 ). Where is positional uncertainty a problem for species distribution modelling? Ecography , 37 ( 2 ), 191 – 203 . doi: 10.1111/j.1600-0587.2013.00205.x OpenUrl CrossRef PubMed ↵ Nobarinezhad , M. H. , Challagundla , L. , & Wallace , L. E . ( 2020 ). Small-Scale Population Connectivity and Genetic Structure in Canada Thistle ( Cirsium arvense ) . International Journal of Plant Sciences , 181 ( 4 ), 473 – 484 . doi: 10.1086/706882 OpenUrl CrossRef ↵ O’Leary , S. J. , Puritz , J. B. , Willis , S. C. , Hollenbeck , C. M. , & Portnoy , D. S . ( 2018 ). These aren’t the loci you’e looking for: Principles of effective SNP filtering for molecular ecologists . Molecular Ecology , 27 ( 16 ), 3193 – 3206 . doi: 10.1111/mec.14792 OpenUrl CrossRef ↵ Pauly , G. ( 2008 ). Phylogenetic systematics, historical biogeography, and the evolution of vocalization in Nearctic Toads (Bufo) [PhD] . The University of Texas at Austin . ↵ Pauly , G. B. , Piskurek , O. , & Shaffer , H. B . ( 2007 ). Phylogeographic concordance in the southeastern United States: The flatwoods salamander, Ambystoma cingulatum , as a test case . Molecular Ecology , 16 ( 2 ), 415 – 429 . doi: 10.1111/j.1365-294X.2006.03149.x OpenUrl CrossRef PubMed Peterson , B. K. , Weber , J. N. , Kay , E. H. , Fisher , H. S. , & Hoekstra , H. E . ( 2012 ). Double Digest RADseq: An Inexpensive Method for De Novo SNP Discovery and Genotyping in Model and Non-Model Species. PLoS ONE , 7 ( 5 ), e37135 . doi: 10.1371/journal.pone.0037135 OpenUrl CrossRef PubMed ↵ Phillips , S. J. , Anderson , R. P. , & Schapire , R. E . ( 2006 ). Maximum entropy modeling of species geographic distributions . Ecological Modelling , 190 ( 3 ), 231 – 259 . doi: 10.1016/j.ecolmodel.2005.03.026 OpenUrl CrossRef PubMed ↵ Poland , J. A. , Brown , P. J. , Sorrells , M. E. , & Jannink , J.-L . ( 2012 ). Development of High-Density Genetic Maps for Barley and Wheat Using a Novel Two-Enzyme Genotyping-by-Sequencing Approach . PLoS ONE , 7 ( 2 ), e32253 . doi: 10.1371/journal.pone.0032253 OpenUrl CrossRef PubMed ↵ Prato , T. , & Fagre , D. B . ( 2010 ). Sustainable Management of the Crown of the Continent Ecosystem . The George Wright Forum , 27 ( 1 ), 77 – 93 . OpenUrl ↵ Pritchard , J. K. , Stephens , M. , & Donnelly , P . ( 2000 ). Inference of Population Structure Using Multilocus Genotype Data . Genetics , 155 ( 2 ), 945 – 959 . doi: 10.1093/genetics/155.2.945 OpenUrl Abstract / FREE Full Text ↵ Pröhl , H . ( 2003 ). Variation in Male Calling Behaviour and Relation to Male Mating Success in the Strawberry Poison Frog ( Dendrobates pumilio ) . Ethology , 109 ( 4 ), 273 – 290 . doi: 10.1046/j.1439-0310.2003.00863.x OpenUrl CrossRef Web of Science ↵ R Core Team . ( 2024 ). R: A language and environment for statistical computing. R Foundation for Statistical Computing https://www.R-project.org/ ↵ T. Dobzhansky , M. K. Hecht , and W. C. Steere Remington , C. L . ( 1968 ). Suture-zones of hybrid interaction between recently joined biotas. Pp. 321–428 . in T. Dobzhansky , M. K. Hecht , and W. C. Steere , eds. Evolutionary biology . Plenum Press , New York . ↵ Rosner-Katz , H. , McCune , J. L. , & Bennett , J. R . ( 2020 ). Using stacked SDMs with accuracy and rarity weighting to optimize surveys for rare plant species . Biodiversity and Conservation , 29 ( 11–12 ), 3209 – 3225 . doi: 10.1007/s10531-020-02018-1 OpenUrl CrossRef ↵ Schoener , T. W . ( 1968 ). The Anolis Lizards of Bimini: Resource Partitioning in a Complex Fauna . Ecology , 49 ( 4 ), 704 – 726 . doi: 10.2307/1935534 OpenUrl CrossRef Web of Science ↵ Shafer , A. B. A. , Cullingham , C. I. , Côté , S. D. , & Coltman , D. W . ( 2010 ). Of glaciers and refugia: A decade of study sheds new light on the phylogeography of northwestern North America . Molecular Ecology , 19 ( 21 ), 4589 – 4621 . doi: 10.1111/j.1365-294X.2010.04828.x OpenUrl CrossRef PubMed Web of Science ↵ Soltis , D. E. , Gitzendanner , M. A. , Strenge , D. D. , & Soltis , P. S . ( 1997 ). Chloroplast DNA intraspecific phylogeography of plants from the Pacific Northwest of North America . Plant Systematics and Evolution , 206 ( 1–4 ), 353 – 373 . doi: 10.1007/BF00987957 OpenUrl CrossRef Web of Science ↵ Soltis , D. E. , Morris , A. B. , McLACHLAN , J. S. , Manos , P. S. , & Soltis , P. S . ( 2006 ). Comparative phylogeography of unglaciated eastern North America: Phylogeography of unglaciated eastern North America . Molecular Ecology , 15 ( 14 ), 4261 – 4293 . doi: 10.1111/j.1365-294X.2006.03061.x OpenUrl CrossRef PubMed Web of Science ↵ Sommer , R. S. , & Zachos , F. E . ( 2009 ). Fossil evidence and phylogeography of temperate species: ‘Glacial refugia’ and post-glacial recolonization . Journal of Biogeography , 36 ( 11 ), 2013 – 2020 . doi: 10.1111/j.1365-2699.2009.02187.x OpenUrl CrossRef ↵ Stacy , E. M. , Robards , M. D. , Jung , T. S. , Kukka , P. M. , Sullivan , J. , Hohenlohe , P. A. , & Waits , L. P . ( 2025 ). Comparing microsatellites and single nucleotide polymorphisms to evaluate genetic structure and diversity in wolverines ( Gulo gulo ) across Alaska and western Canada . Journal of Mammalogy , 106 ( 3 ), 561 – 575 . doi: 10.1093/jmammal/gyae151 OpenUrl CrossRef PubMed Suárez , D. , Arribas , P. , Jiménez-García , E. , & Emerson , B. C . ( 2022 ). Dispersal ability and its consequences for population genetic differentiation and diversification . Proceedings of the Royal Society B: Biological Sciences , 289 , 20220489 . doi: 10.1098/rspb.2022.0489 OpenUrl CrossRef PubMed ↵ Sukumaran , J. , & Knowles , L. L . ( 2017 ). Multispecies coalescent delimits structure, not species . Proceedings of the National Academy of Sciences , 114 ( 7 ), 1607 – 1612 . doi: 10.1073/pnas.1607921114 OpenUrl Abstract / FREE Full Text ↵ Swenson , N. G. , & Howard , D. J . ( 2004 ). Do suture zones exist? Evolution , 58 ( 11 ), 2391 – 2397 . doi: 10.1111/j.0014-3820.2004.tb00869.x OpenUrl CrossRef PubMed Web of Science ↵ Swenson , N. G. , & Howard , D. J . ( 2005 ). Clustering of Contact Zones, Hybrid Zones, and Phylogeographic Breaks in North America . The American Naturalist , 166 , 581 – 591 . OpenUrl CrossRef PubMed Web of Science ↵ Switzer , J. F. , Johnson , R. , Lubinski , B. A. , & King , T. L . ( 2009 ). Genetic structure in the Anazyrus boreas species group (Anura, Bufonidae): An evaluation of the Southern Rocky Mountain population (p. 70 ) [ U.S. Fish and Wildlife Service ]. ↵ Thompson , P. D . ( 2019 ). Translocation of Boreal Toad ( Anaxyrus boreas boreas ) into Two Springs in the Grouse Creek Mountains, Utah, Including Demographic Observations . Western North American Naturalist , 79 ( 1 ), 24 . doi: 10.3398/064.079.0103 OpenUrl CrossRef ↵ Trumbo , D. R. , Hardy , B. M. , Crockett , H. J. , Muths , E. , Forester , B. R. , Cheek , R. G. , Zimmerman , S. J. , Corey-Rivas , S. , Bailey , L. L. , & Funk , W. C . ( 2023 ). Conservation genomics of an endangered montane amphibian reveals low population structure, low genomic diversity and selection pressure from disease . Molecular Ecology, n/a(n/a ). doi: 10.1111/mec.17175 OpenUrl CrossRef ↵ Turbek , S. P. , Funk , W. C. , & Ruegg , K. C . ( 2023 ). Where to draw the line? Expanding the delineation of conservation units to highly mobile taxa . Journal of Heredity , 114 ( 4 ), 300 – 311 . doi: 10.1093/jhered/esad011 OpenUrl CrossRef PubMed ↵ UNESCO . ( 2025 ). Crown of the Continent Biosphere Reserve . Man and the Biosphere Programme, UNESCO. Retrieved May, 2025, from UNESCO website: https://www.unesco.org/en/mab/crown-continent . U.S. Fish and Wildlife Service . ( 2017 ). Species Status Assessment Report for the Eastern Population of the Boreal Toad (Anaxyrus boreas boreas) . ↵ Vale , C. G. , Tarroso , P. , & Brito , J. C . ( 2014 ). Predicting species distribution at range margins: Testing the effects of study area extent, resolution and threshold selection in the Sahara-Sahel transition zone . Diversity and Distributions , 20 ( 1 ), 20 – 33 . doi: 10.1111/ddi.12115 OpenUrl CrossRef ↵ Wait , D. R. , & Peñalba , J. V . ( 2025 ). Suture zones, speciation, and evolution . Evolution , 79 ( 3 ), 329 – 341 . doi: 10.1093/evolut/qpae184 OpenUrl CrossRef PubMed ↵ Wang , T. , Hamann , A. , Spittlehouse , D. , & Carroll , C . ( 2016 ). Locally Downscaled and Spatially Customizable Climate Data for Historical and Future Periods for North America . PLOS ONE , 11 ( 6 ), e0156720 . doi: 10.1371/journal.pone.0156720 OpenUrl CrossRef PubMed ↵ Waples , R. S. , & Lindley , S. T . ( 2018 ). Genomics and conservation units: The genetic basis of adult migration timing in Pacific salmonids . Evolutionary Applications , 11 ( 9 ), 1518 – 1526 . doi: 10.1111/eva.12687 OpenUrl CrossRef PubMed ↵ Warren , D. L. , Glor , R. E. , & Turelli , M . ( 2008 ). Environmental niche equivalency versus conservatism: Quantitative approaches to niche evolution . Evolution , 62 ( 11 ), 2868 – 2883 . doi: 10.1111/j.1558-5646.2008.00482.x OpenUrl CrossRef PubMed Web of Science ↵ Waters , J. M. , Emerson , B. C. , Arribas , P. , & McCulloch , G. A . ( 2020 ). Dispersal Reduction: Causes, Genomic Mechanisms, and Evolutionary Consequences . Trends in Ecology & Evolution , 35 ( 6 ), 512 – 522 . OpenUrl CrossRef PubMed ↵ Watt , C. M. , Kierepka , E. M. , Ferreira , C. C. , Koen , E. L. , Row , J. R. , Bowman , J. , Wilson , P. J. , & Murray , D. L . ( 2021 ). Canada lynx ( Lynx canadensis ) gene flow across a mountain transition zone in western North America . Canadian Journal of Zoology , 99 ( 2 ), 131 – 140 . doi: 10.1139/cjz-2019-0247 OpenUrl CrossRef ↵ Wilson , G. A. , Fulton , T. L. , Kendell , K. , Scrimgeour , G. , Paszkowski , C. A. , & Coltman , D. W . ( 2008 ). Genetic diversity and structure in Canadian northern leopard frog ( Rana pipiens ) populations: Implications for reintroduction programs . Canadian Journal of Zoology , 86 ( 8 ), 863 – 874 . doi: 10.1139/Z08-062 OpenUrl CrossRef ↵ Wright S , 1943 . Isolation by distance . Genetics , 28 ( 2 ), 114 – 138 . OpenUrl FREE Full Text ↵ Young , N. , Carter , L. , Evangelista , P. , & Jarnevich , C. S. ( 2011 ). A Maxent Model v3.3.3e Tutorial . ↵ Zeisset , I. , & Beebee , T. J. C . ( 2008 ). Amphibian phylogeography: A model for understanding historical aspects of species distributions . Heredity , 101 ( 2 ), 109 – 119 . doi: 10.1038/hdy.2008.30 OpenUrl CrossRef PubMed Web of Science ↵ Zhao , L. , Santos , J. C. , Wang , J. , Ran , J. , Tang , Y. , & Cui , J . ( 2021 ). Noise constrains the evolution of call frequency contours in flowing water frogs: A comparative analysis in two clades . Frontiers in Zoology , 18 ( 1 ), 37 . doi: 10.1186/s12983-021-00423-y OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted June 23, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Pronounced genetic structure associated with differences in a reproductive trait and climatic barriers in Canadian populations of the western toad (Anaxyrus boreas) Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Pronounced genetic structure associated with differences in a reproductive trait and climatic barriers in Canadian populations of the western toad ( Anaxyrus boreas ) Jayna C. Bergman , Juan Enciso-Romero , Gregory B. Pauly , Roseanna Gamlen-Greene , Melissa Todd , Julie A. Lee-Yaw bioRxiv 2025.06.17.660222; doi: https://doi.org/10.1101/2025.06.17.660222 Share This Article: Copy Citation Tools Pronounced genetic structure associated with differences in a reproductive trait and climatic barriers in Canadian populations of the western toad ( Anaxyrus boreas ) Jayna C. Bergman , Juan Enciso-Romero , Gregory B. Pauly , Roseanna Gamlen-Greene , Melissa Todd , Julie A. Lee-Yaw bioRxiv 2025.06.17.660222; doi: https://doi.org/10.1101/2025.06.17.660222 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7621) Biochemistry (17645) Bioengineering (13867) Bioinformatics (41873) Biophysics (21420) Cancer Biology (18550) Cell Biology (25443) Clinical Trials (138) Developmental Biology (13360) Ecology (19866) Epidemiology (2067) Evolutionary Biology (24289) Genetics (15587) Genomics (22472) Immunology (17707) Microbiology (40318) Molecular Biology (17142) Neuroscience (88457) Paleontology (666) Pathology (2826) Pharmacology and Toxicology (4815) Physiology (7634) Plant Biology (15111) Scientific Communication and Education (2042) Synthetic Biology (4285) Systems Biology (9813) Zoology (2268)