Conservation genomics of an endangered montane amphibian reveals low population structure, low genomic diversity and selection pressure from disease

Wildlife diseases are a major global threat to biodiversity. Boreal toads (Anaxyrus [Bufo] boreas) are a state‐endangered species in the southern Rocky Mountains of Colorado and New Mexico, and a species of concern in Wyoming, largely due to lethal skin infections caused by the amphibian chytrid fungus Batrachochytrium dendrobatidis (Bd). We performed conservation and landscape genomic analyses using single nucleotide polymorphisms from double‐digest, restriction site‐associated DNA sequencing in combination with the development of the first boreal toad (and first North American toad) reference genome to investigate population structure, genomic diversity, landscape connectivity and adaptive divergence. Genomic diversity (π = 0.00034–0.00040) and effective population sizes (Ne = 8.9–38.4) were low, likely due to post‐Pleistocene founder effects and Bd‐related population crashes over the last three decades. Population structure was also low, likely due to formerly high connectivity among a higher density of geographically proximate populations. Boreal toad gene flow was facilitated by low precipitation, cold minimum temperatures, less tree canopy, low heat load and less urbanization. We found >8X more putatively adaptive loci related to Bd intensity than to all other environmental factors combined, and evidence for genes under selection related to immune response, heart development and regulation and skin function. These data suggest boreal toads in habitats with Bd have experienced stronger selection pressure from disease than from other, broad‐scale environmental variations. These findings can be used by managers to conserve and recover the species through actions including reintroduction and supplementation of populations that have declined due to Bd.


| INTRODUC TI ON
Wildlife diseases are a major threat to species of conservation concern around the world (Brearley et al., 2013;Daszak et al., 2000;Smith et al., 2006Smith et al., , 2009)).These diseases can cause population declines and extirpations, particularly in concert with other environmental stressors like global climate change, by invading new geographic regions and establishing in historically naïve wildlife populations (Buttke et al., 2021;Price et al., 2019;Rohr et al., 2011).
While the consequences of disease on wildlife population dynamics have been of interest to ecologists and conservation biologists for some time (Anderson & May, 1979;Hudson et al., 1998), host genetic information can supplement traditional approaches in disease ecology, enhance our understanding of host-pathogen interactions and inform species conservation (Fitak et al., 2019;McKnight et al., 2017;Storfer et al., 2021).This is especially true with the recent advances in genomic sequencing methods that can reveal patterns of genome-wide variation and allow managers to incorporate both neutral and adaptive diversity in conservation planning and decision-making (Allendorf et al., 2010;Flanagan et al., 2018;Funk et al., 2019).
Disease-driven population declines can have a variety of genetic consequences that affect a population's ability to persist with disease and their ability to recover.Declines are associated with local extirpations, creating increased isolation between remaining populations which leads to loss of gene flow and connectivity.In turn, isolated populations are often subjected to losses of genetic diversity due to genetic drift, increased probabilities of inbreeding and decreased effective population sizes (Gilpin & Soule, 1986).The loss of genetic diversity can also exacerbate the deterministic effects driving population extinction (e.g.disease; Fagan & Holmes, 2006, Spielman et al., 2004), resulting in a negative feedback loop that has been termed the 'extinction-vortex' (Forester et al., 2022;Frankham et al., 2002;Gilpin & Soule, 1986;Soulé & Mills, 1998).Therefore, understanding the neutral genetic parameters of residual populations post-decline provides crucial information about extinction risk (Schwartz et al., 2007).
Conservation managers are often interested in quantifying levels of genetic diversity and inbreeding in remnant populations, in addition to describing the neutral genetic structure and landscape connectivity of populations to aid conservation planning.For example, integrating population genetic estimates of diversity, gene flow and inbreeding depression within and among remaining populations can improve the success of management actions such as translocation and reintroductions (Armstrong & Seddon, 2008;Rhodes Jr & Latch, 2010).Genetic diversity metrics can inform translocation plans by aiding the selection and identification of source and recipient populations to promote genetic rescue or establish new self-sustaining populations during reintroductions (Whiteley et al., 2015).In addition, understanding how the environment structures, facilitates or hinders connectivity and gene flow among populations can aid in the selection of source populations based on the surrounding habitat matrix and enhance the probability of future gene flow between reintroduced and remnant populations (Manel et al., 2003(Manel et al., , 2010;;Segelbacher et al., 2008).
With the development of genomic approaches such as reduced-representation sequencing (e.g.double-digest restriction site-associated DNA sequencing; ddRAD), managers can assess large numbers of genetic loci across the genome and incorporate information about the adaptive diversity and the adaptive potential of remnant populations into conservation planning (Flanagan et al., 2018;Forester et al., 2022;Funk et al., 2019).Identifying populations that are locally adapted to specific environmental parameters or a given stressor could improve the success of reintroductions or translocations by matching specific genotypes to current or future environmental conditions.Information on the adaptive diversity of populations can also be used by managers to inform captive breeding strategies and select source populations with beneficial genotypes to augment poorly adapted recipient populations.In addition, new reference genomes developed for non-model wildlife species now allow researchers to identify portions of the genome that may be under selection due to environmental factors, including disease (DeCandia et al., 2018;Formenti et al., 2022;Zamudio et al., 2020).
For example, Margres et al. (2018) were able to identify genes involved in the immune response of Tasmanian devils (Sarcophilus harrisii) to lethal devil facial tumour disease.By combining new wildlife reference genomes with cost-effective, reduced-representation sequencing, wildlife researchers can identify genes that may help an organism adapt to novel diseases and other contemporary or future stressors (Balkenhol et al., 2015;Joost et al., 2013;Manel et al., 2010).
Among the various wildlife taxa affected by disease, amphibians have suffered dramatic worldwide declines and extinctions (Scheele et al., 2019;Skerratt et al., 2007;Wake & Vredenburg, 2008).Chytridiomycosis, an amphibian disease caused by two fungi, Batrachochytrium dendrobatidis (Bd) and Batrachochytrium salamandrivorans (Bsal), is 'the worst infectious disease ever recorded among vertebrates in terms of the number of species impacted, and its propensity to drive them to extinction' (Gascon, 2007).To date, this disease has caused the greatest known loss of biodiversity due to a pathogen ever documented (Scheele et al., 2019).Managers and conservation organizations tasked with restoring and recovering affected amphibian species can benefit from the integration of large numbers of genetic loci from across the genome in the planning of reintroductions, translocations and captive breeding efforts.For example, this information can help identify source populations for translocation that are not only too divergent but also not too similar to recipient populations (Fitzpatrick & Funk, 2020).More specifically, outbreeding depression is most likely to occur between populations when they are distinct species, have fixed genomic differences, have had no gene flow within the last 500 years or inhabit very different environments (Frankham et al., 2011).

The boreal toad (Anaxyrus [Bufo] boreas) in the southern Rocky
Mountains (SRM) is an example of a species that has experienced alarming declines over the last four decades, and most of the declines are linked to Bd (Hardy et al., 2022;Mosher et al., 2018;Muths et al., 2003Muths et al., , 2008;;Scherer et al., 2005).While the geographic range of the boreal toad species complex is vast, spanning the western USA, Canada and Mexico (Stebbins, 2003), populations in the SRM have experienced the most severe declines, prompting wildlife agencies (Colorado Parks and Wildlife [CPW], New Mexico Department of Game and Fish and Wyoming Game and Fish Department) to designate this species as endangered or of special concern.In 2018, representatives from state and federal agencies comprising the Boreal Toad Conservation Team for the SRM, along with other regional species experts, employed a structured decision-making process to evaluate all possible conservation strategies (Converse et al., 2017;Gerber et al., 2018).The best management strategies include translocations to promote genetic rescue in remnant populations, or reintroductions in formerly occupied mountain ranges (Gerber et al., 2018).Maintaining captive assurance and breeding colonies in the hopes of preserving historical genetic lineages and promoting beneficial genotypes through crossing is also a priority for the con-

| Study system
The southern Rocky Mountain ecoregion (SRM) extends from southcentral Wyoming through Colorado and into north-central New Mexico.It is characterized by cool summers, severe winters, a relatively short frost-free period (25-150 days, depending on latitude and elevation) and deep snowpack at higher elevations.The boreal toad's historical distribution within the SRM is not known with certainty, but it was considered common to abundant in most higherelevation regions of Colorado and south-central Wyoming and was documented in north-central New Mexico until the mid-1980s (Carey et al., 2005;Hammerson, 1999;Muths & Nanjappa, 2005).

| Field sampling
In 2018 and 2019, state and federal biologists sampled all known active breeding sites across the SRM and collected genomic and Bd samples from each individual toad captured (Figure 1).Approximately half (N = 16) of the sites sampled had recorded Bd presence between 2003 and 2019, while the remaining 20 sites had never recorded Bd present (Table S1).Genomic (host) samples consisted of buccal swabs and tissues (i.e.toe and tail clips from adults and tadpoles, respectively), and Bd (disease) samples consisted of skin swabs from adults, juveniles and metamorphs (<2 cm snout-vent length [SVL]).Buccal swabbing was performed for 30-60 s for each terrestrial-stage toad, and two swabs were collected per individual.Skin swabbing was performed at least 20 times on the ventral surface for adults and juveniles (25 times for metamorphs), and 5 times on the rear feet and webbing of all individuals, and two skin swabs were also collected per individual.All genomic and Bd swabs and tissue samples were immediately placed in 95% ethanol and stored on ice in the field.
Ethanol was dried from the swabs in a fume hood, and DNA was extracted from tissues and buccal swabs at Colorado State University using Qiagen DNEasy kits, with a modified protocol for the lowerconcentration buccal swab sampling (Goldberg et al., 2003).

| Reference genome
We developed a boreal toad reference genome -the first North American toad [family Bufonidae] genome -in collaboration with Dovetail Genomics, Inc. (dovet ailge nomics.com).High-quality, genomic DNA was isolated from flash-frozen tissues of a male boreal toad, originally captured from a breeding site in Jackson County, Colorado, and subsequently part of the captive assurance colony at CPW's Native Aquatic Species Restoration Facility.A draft de novo assembly was developed from long-read sequencing of genomic DNA on a PacBio HiFi platform, as well as Omni-C sequencing using chromatin proximity to help span long-distance repetitive regions in the genome (Wenger et al., 2019).The aligned de novo assembly and Omni-C library were used as input data for HiRise, a software pipeline designed specifically for using proximity ligation data to scaffold large genome assemblies (Putnam et al., 2016).Omni-C library sequences were aligned to the draft input assembly using bwa (Li & Durbin, 2010).The separation of Omni-C read pairs mapped within draft scaffolds was analysed by HiRise to produce a likelihood model for genomic distances between read pairs, and the model was used to identify and break putative mis-joins, to score prospective joins and to make high-confidence joins across long-distance scaffolds.
For each Omni-C library, chromatin was fixed in place with formaldehyde in the nucleus and then extracted.Fixed chromatin was digested with DNAse I, and chromatin ends were repaired and ligated to a biotinylated bridge adapter followed by proximity ligation of adapter-containing ends.After proximity ligation, crosslinks were reversed and the DNA was purified.Purified DNA was treated to remove biotin that was not internal to ligated fragments.Sequencing libraries were then generated using NEBNext Ultra enzymes and Illumina-compatible adapters.Biotin-containing fragments were isolated using streptavidin beads before polymerase chain reaction (PCR) enrichment of each library.Libraries were sequenced on an Illumina HiSeq 4000 to produce >30× mean sequencing coverage.
HiRise was then implemented on the 141.9 million reads with mapping quality (MQ) > 50 for scaffolding.Wtdbg2 (Ruan & Li, 2020) was used to assemble and combine scaffolds, with the following parameters: -genome_size 5.0 g -read_type sq -min_read_len 5000.Blobtools v1.1.1 (Laetsch & Blaxter, 2017) was used to identify potential contamination in the assembly based on the National Center for Biotechnology Information's (NCBI) Basic Local Alignment Search Tool (BLAST) v.2.9 results of the assembly against the nucleotide database.Scaffolds that were identified as contaminants were removed from the assembly.The filtered assembly was then used as an input to run purge_dups v1.1.2(Guan et al., 2020), which allowed potential haplotypic duplications to be removed, resulting in the final reference genome assembly.
Annotation of the final reference genome was performed using reference transcriptome data developed from ribonucleic acid sequencing (RNA-seq) in an ongoing challenge experiment study by Corey-Rivas et al. at New Mexico Highlands University.In addition, BLAST was used to find homologous genes across other invertebrate and vertebrate organisms with published reference genomes.

| ddRAD sequencing and genotyping
We performed double-digest, restriction site-associated DNA sequencing (ddRAD; Peterson et al., 2012) to develop homologous markers spread randomly across the genome for subsequent landscape and population genomic analyses.DNA was digested using NlaIII (4-cutter) and EcoR1 (6-cutter) restriction enzymes, followed by dual indexing with 48 uniquely barcoded P1 adapters, and P2 primer indices for demultiplexing (Peterson et al., 2012).P2 adapters contained 10 base pair (bp) degenerate indices for removing PCR duplicates (Schweyen et al., 2014).We size selected using a Sage Science Blue Pippin for 376-476 bp fragments (i.e.300-400 bp fragments, plus 76 bp of P1 and P2 adapters), performed final PCRs using 12 cycles and assessed final libraries' fragment size distributions and concentrations on an Agilent Tapestation 2200 and Qubit 2.0 fluorometer.University of Oregon Genomics Core (gc3f.uoreg on.edu) performed 150 bp paired-end sequencing on a HiSeq 4000.
Our bioinformatics pipeline consisted of FastQC (www.bioin forma tics.babra ham.ac.uk/ proje cts/ fastqc/ ) to assess the read quality and verify the presence of barcodes and cut sites, followed by removal of PCR duplicates using Stacks V.2.53 clone_filter (Catchen et al., 2013).CutAdapt was used to trim off the excess P2 degenerate barcodes and adapters up to the cut sites (Martin, 2011), followed by Stacks process_radtags to demultiplex individuals using the unique P1 barcodes (Catchen et al., 2013).We aligned all ddRAD reads to the new boreal toad reference genome using the Burrows-Wheeler algorithm (i.e.bwa-mem V.0.7.17; Li & Durbin, 2010), which is designed to efficiently align short-sequence reads to large-reference genomes.Finally, we called our mapped single nucleotide polymorphism (SNP) loci across the boreal toad reference genome using the Stacks ref_map pipeline (Catchen et al., 2013).
After developing the SNP catalogue, we filtered SNP loci using Stacks populations module to remove low-coverage SNPs by sites (SNPs present < 75% of sites across the study area) and by individuals (SNPs present < 75% of individuals), likely paralogues (SNPs with heterozygosity > 50%), singletons (rare SNPs with a minor allele count < 2) and SNPs located on the same sequencing read that were expected to have excessively high-linkage disequilibrium (LD; selected one random SNP per locus in loci with multiple SNPs; Catchen et al., 2013).We then used Plink V.1.07to filter out low-coverage individuals (i.e.individuals with <75% of the total SNPs in the Stacks catalogue; Purcell et al., 2007) and calculate the total genotyping rate across all SNPs to assess SNP missingness, followed by the R package RADiator to calculate the mean depth of coverage across all SNP loci (Gosselin et al., 2017).In total, we sequenced 5 libraries of 62 individual toads per library, for a total of 310 individuals.

| Population structure and genomic diversity
We extracted at least 50 ng of DNA from buccal swabs and tissues (i.e.toe and tail clips from adults and tadpoles, respectively) from 310 individual boreal toads across 36 sites within 15 regions across the study area in the SRM (Figure 1).We summarized boreal toad population genomic measures by region in the SRM by combining the 36 breeding sites into 15 regions based on geographic proximity (<20 km apart), lack of obvious landscape barriers (e.g.no major highways, high-elevation ridgelines or other obvious barriers to dispersal between them) and a population of K = 1 supported using Admixture (Alexander et al., 2009; Figure 1).By combining breeding sites in close proximity, we were able to reduce the number of sites with small sample sizes (e.g.<5 individuals) due to recent Bd-related population crashes (Figure 1), which can result in imprecise estimates of population-level genetic structure (e.g.F ST ) and genomic diversity (Willing et al., 2012).Population structure across the study area was assessed using discriminant analysis of principal components (DAPC) with Adegenet (Jombart, 2008, Jombart et al., 2010).In addition, the programs Admixture (Alexander et al., 2009) and SplitsTree (Huson & Bryant, 2006) were used to further explore population structure.
F ST values among all 15 regions were calculated using Stacks populations (Catchen et al., 2013).Patterns of isolation by distance (IBD; Wright, 1943) were explored using Mantel tests with the ecodist package in R (Goslee, 2009) to compare pairwise genetic distances to topographically corrected geographic distances.
To evaluate levels of genomic diversity and potential inbreeding within regions, we used Stacks populations to calculate observed and expected heterozygosity (H OBS and H EXP ), nucleotide diversity (π), inbreeding coefficients (F IS ) and the total number of polymorphic loci (Catchen et al., 2013).We estimated effective population sizes (N e ) with the LD method as implemented in NeEstimator V.2.1, using SNPs with a minor allele frequency (MAF) > 0.1 since rare alleles may upwardly bias N e estimates (Do et al., 2014), as well as correcting for bias due to the number of chromosomes in boreal toads (2n = 22; Waples et al., 2016).Finally, we inferred the demographic history of boreal toads across the SRM using the site frequency spectrum (SFS) and Stairway Plot 2 (Liu & Fu, 2015, 2020).We used an estimated generation time of 5 years (Muths & Nanjappa, 2005), and an estimated mutation rate of 2.52e-9 (Crawford, 2003).

| Landscape connectivity
We investigated landscape connectivity using spatial autoregressive models (SAM; Peterson et al., 2019) and generalized Wishart modelling (GWM; Hanks & Hooten, 2013) using the rwc package in R (Hanks, 2018) as specified in Zimmerman et al. (2022).Landscape surfaces were developed in ArcGIS V.10.7 and consisted of 11 environmental factors hypothesized to be important for boreal toad movement and gene flow across the SRM (Table 1 and Figure 1).We first fit a set of single covariate models and an intercept-only (i.e.null IBD) model to our data over a range of raster grain sizes (0.7-25.3 km 2 ) to identify the raster grain size that best-reduced computation time, improved covariance matrix estimation and resulted in stable model inferences based on comparison of model rank and coefficient estimates.We resampled all rasters to the same geographic extent with bilinear interpolation and used the aggregate function in the raster R package to decrease raster covariate resolution.We found that model fitting failed for some variables below a grain size of 1.7 km 2 and model rank and coefficient estimates were comparable up to a grain size of 8.4 km 2 .Therefore, we fit all subsequent multi-covariate GWM models using a grain size of 8.4 km 2 .
We used a minimally correlated (Pearson's r < 0.5) set of variables for GWM multivariate modelling by retaining variables with TA B L E 1 Landscape variables used for connectivity analyses of boreal toads in the southern Rocky Mountains (SRM), including data source, spatial resolution and ecological justification.

| Adaptive divergence
We performed a multivariate genotype-environment association (GEA) test using redundancy analysis (RDA; Capblancq & Forester, 2021;Capblancq et al., 2018;Forester et al., 2016) to investigate whether toad populations showed evidence of local adaptation to environmental variation and Bd exposure.The environmental RDA included the 11 environmental variables used in connectivity analyses (Figure 1 and Table 1) after pruning highly correlated predictors (Pearson's r > 0.9).The Bd exposure RDA included two predictors: Bd prevalence measured in the field and time since first Bd detection obtained from previous monitoring efforts (Boreal Toad Conservation Team data; correlation between exposure measures was Pearson's r < 0.9).We analysed 18 study sites with ≥5 toads sampled per site (186 toads total) and extracted environmental values at each site for temperature and precipitation variables, and within a 10-km-radius buffer (representing a typical maximum dispersal distance of an adult boreal toad; Lucid et al., 2021;Murphy et al., 2010;Muths, 2003, Schmetterling & Young, 2008) for per cent tree canopy cover, heat load index, compound topographic index of wetness, streams and rivers, roads and per cent impervious surfaces.
Finally, we imputed missing SNP data using the most common genotype across all individuals, as using the mean value for imputation may upwardly bias heterozygosity (Capblancq & Forester, 2021).For both RDAs, we identified outliers along the first two RDA axes (i.e. axes explaining the most genomic variation) using a >2.5 standard deviation (SD) cutoff (Forester et al., 2016).

| SNP annotation
We used our newly developed boreal toad reference genome to determine whether SNPs identified as putatively adaptive from RDA analyses were located in or near (i.e. in LD with) genes that may confer resistance to Bd or local adaptation to environmental factors (Table 1 and Figure 1).First, we visually inspected LD plots within a 1 million bp window using VCFtools 0.1.17with a min and max alleles setting of 2 (Danecek et al., 2011) to determine the distance along the boreal toad genome where LD was high (e.g.r 2 > 0.7).We then input our outlier SNPs into BEDOPS v. 2.4.39 (Neph et al., 2012) and output genes that showed high LD with outlier SNP site coordinate in the annotated boreal toad genome to search for genes potentially under selection.Gene names from the boreal toad genome annotation were converted to the western clawed frog (Xenopus tropicalis) ensemble namespace using gprofiler2 v.0.2.1 (Kolberg et al., 2020), and Entrez gene names were identified for all output genes using mygene v.1.24.0 (Mark et al., 2020).We used the R Bioconductor package biomart v.2.45.8 (Durinck et al., 2005(Durinck et al., , 2009) ) to output gene ontology (GO) information for each gene found within our query sequences using the X. tropicalis Ensembl database.We repeated the above procedure using the human (Homo sapiens) Ensemble database due to higher quality of GO annotations for H. sapiens compared to X. tropicalis and expected homology across vertebrates.
We then searched the annotation reports for GO terms related to traits hypothesized to confer Bd resistance, as well as physiological responses to environmental variation (Figure 1).

| Field sampling, reference genome and sequencing
The final boreal toad reference genome was 4.57 billion base pairs

| Population structure and genomic diversity
We detected very low levels of population differentiation among regions (mean F ST = 0.03, range = 0.01-0.13and standard error = 0.002; Table S2).Little to no population structure was detected using Splitstree or Admixture (Figures S1 and S2a), but the first 2 DAPC axes showed some population divergence of five regions that were geographically peripheral to the rest of the range, consisting of southern Wyoming (WY), Elkhead Mountains (Elkhead), Northern Front Range (FR_N), Northern Sawatch Range (Saw_N) and Mosquito Range (Mosq) (Figure 2).A subtle pattern of IBD across the SRM was also detected using a Mantel test (Mantel's r = 0.08), although it was not significant (p = .48;Figure S2b).
Low overall genomic diversity (π, H OBS , H EXP ) was found within all regions of the SRM (Table 2).We also found small effective population sizes (N e ) and evidence of departures from HW proportions (e.g.F IS , H OBS < H EXP ) at many sites (Table 2), suggesting inbreeding

| Landscape connectivity
We performed landscape genomic analyses using GWM across 18 study sites with ≥5 toads per site (total of 186 individuals) and retained 5 minimally correlated (Pearson's r < 0.5) landscape variables for modelling consisting of annual precipitation, minimum temperature of the coldest month, per cent tree canopy cover, heat load index and per cent impervious surface (Table 3, Figures 1 and 3).All landscape models ranked above the intercept-only, null IBD model, with 23 multivariate models ranking higher than the top univariate model (minimum temperature of the coldest month; M2, Table 3).
The top multivariate model (M37, Table 3) predicted the highest gene flow among regions that were separated by areas of lower annual precipitation, colder minimum temperatures, less tree canopy cover, lower heat loads and lower per cent impervious surfaces.We also developed a conductance surface using the top-modelled GWM relationships, which showed relatively high boreal toad gene flow potential range wide, with populations in the northern portion of the study area surrounded by more landscape features likely to impede movement and gene flow than southern populations (Figure 3).

| Adaptive divergence
For the environmental GEA analyses, we first removed elevation and topographic roughness due to strong multicollinearity (Pearson's r > 0.9) with maximum temperature of the warmest month and compound topographic index of wetness (Figure S4).The environmental RDA identified 36 outliers related to minimum and maximum temperatures, precipitation, tree cover, streams and rivers, wetness, heat load, roads and impervious surfaces (Figure 4).By contrast, the Bd exposure RDA identified 310 outliers related to Bd prevalence and time since first detection (Figure 4).

| SNP annotation
We tested whether SNPs identified as outliers by GEA analyses were in or linked to (i.e. in high LD with) annotated genes in our new boreal toad reference genome.Visual inspection of LD plots calculated within a 1 million bp window suggested high linkage (e.g.r 2 < 0.7) between biallelic loci separated by <25 Kbp (Figure S5).We identified 88 genes (Table S3) within 25 Kbp of our outlier SNPs associated with Bd (Figure 4) out of 41,174 genes total (i.e.proportion of genes under selection = 0.002), of which 45 had assigned GO terms according to the X. tropicalis or H. sapiens ensemble databases (Table S4).Our GO term search of the annotation reports identified several genes linked to heart development (HHEX, DUSP6, PHRF1 and NIPBL) and regulation (ATR and DSG2), immune response (NOX1, CNIH1, FADD and REL), skin development (GJB3, ADAM9 and ESRP1) and cell membrane transport and regulation (ST3GAL5, CYBRD1, CHST1, RAB7A and TESK2), which could potentially be related to boreal toad resistance or tolerance to Bd in the SRM (Table S4).Seventeen genes were linked to SNPs related to the nine other landscape factor RDAs, one of which (PSTK) may be associated with temperature tolerance of boreal toads to varying climates across the SRM.

| DISCUSS ION
Genomic techniques are powerful tools for conservation planning and management, particularly for species affected by disease (Hohenlohe et al., 2021;Joost et al., 2013;Manel et al., 2010).Here, we addressed five objectives to inform managers about the neutral and adaptive genomics of an amphibian threatened by disease.
Using a newly developed annotated reference genome, and reduced representation sequencing, we identified patterns of high historic connectivity, low genetic variation and evidence of selection.Our results will help guide the conservation and recovery of this species by informing captive breeding, reintroduction and genetic rescue efforts to mitigate declines and extirpations due to Bd, thereby improving the resiliency of the species into the future.Below, we discuss how our findings contribute to understanding the genomic impacts of disease on wildlife populations, and the potential management implications for the conservation of boreal toads in the southern Rocky Mountains.

| Population structure and genomic diversity
We found low genetic structure and differentiation among toad populations in the SRM (i.e.low F ST and little population genomic structure; Table S2, Figures 2, S1, and S2).We originally expected to find higher differentiation and population structure within the SRM due to the precipitous declines that have increased isolation among remaining populations (Mosher et al., 2018) and the mountainous terrain that can hinder geneflow in many amphibians (e.g.Funk et al., 2005;Giordano et al., 2007).However, it is likely that the low population structure and differentiation we detect is a reflection of historical connectivity among densely distributed populations, many of which are now extirpated and an insufficient number of generations for recent Bd-related declines in the SRM to increase genetic structure (Bd arrival ~1980s; Carey, 1993;Corn et al., 1989).Anecdotal accounts suggest A. boreas was widespread and abundant in highelevation wetlands prior to the arrival of Bd (Hammerson, 1999;Muths & Nanjappa, 2005) and is known to commonly disperse overland up to 6 km (Bartelt et al., 2004;Bull, 2006;Muths, 2003 , 2003, 2018;Thompson, 2019).Boreal toads are also capable of making less-frequent long-distance movements of up to 13 km (Carey et al., 2005;Schmetterling & Young, 2008;Thompson, 2019).
Their historical abundance and dispersal abilities could have facilitated considerable gene flow historically and likely contributed to the low population structure we observed.Indeed, prior genetic studies on boreal toads in the northern Rocky Mountains portion of their range have also shown high dispersal and low population genetic structure (Lucid et al., 2021;Murphy et al., 2010) and gene flow was potentially high historically in the SRM before recent declines (Goebel, Oyler-McCance in USFWS, 2017).The lack of structure also supports previous work that suggests boreal toads throughout the SRM fall within a single genetic clade, without further divergence within the SRM (Goebel et al., 2009).
We also found low levels of genomic diversity (π), small effective population sizes (N e ) and evidence of departures from HW proportions within populations (i.e.positive F IS values, H OBS < H EXP ; Table 2) that may be linked to inbreeding or population substructure across the SRM.We believe three processes may contribute to this low genomic diversity.First, low diversity in the SRM may be a product of recent founder events.Boreal toads colonized their current high-elevation (2300-3700 metres) montane habitats in the SRM relatively recently, after glacial retreat at the end of the last glacial maximum approximately 12.3 ka ago (Brugger et al., 2019;Guido et al., 2007).Founder effects typically generate low genetic diversity in the founding population, whether occurring through natural dispersal (Arauco-Shapiro et al., 2020), intentional translocations (Hedrick et al., 2001) or unintentional invasion (Bai et al., 2012).In support of this scenario, we found genomic evidence from SFS data of boreal toad declines associated with ice ages from 10 to 100 K years ago (Figure S3).Second, SRM boreal toads likely display low genomic diversity due to their position at the southern range edge of their distribution and being isolated from populations to the West and North by hundreds of miles of inhospitable terrain (Goebel et al., 2009).
This isolation may effectively render the SRM a terrestrial 'island', and island populations are more prone to genetic drift and loss of diversity than populations on a 'mainland' (Cardoso et al., 2009;Eldridge et al., 1999;Frankham, 1997;Mills et al., 2004).
Third, it is possible that the low diversity we observed is also related to Bd-induced population declines.There is mixed evidence in the literature of the effects of Bd on genetic diversity, with some studies showing lower genetic diversity in frog populations impacted by Bd (e.g.Byrne et al., 2021;Torres-Sánchez & Longo, 2022), while others showed stable or even increased genetic diversity in frog populations infected with Bd (e.g.Banks et al., 2020;Horner et al., 2017) including in a microsatellite study of boreal toads from Glacier National Park in Montana (Addis et al., 2015).In addition, we do not have genomic diversity estimates from these populations prior to Bd arrival and subsequent declines, or estimates from other parts of the toad's range where declines have not been observed.Therefore, we cannot be certain that the low diversity we observe is attributable to Bd-related declines.In addition, any erosion of genetic diversity from recent declines has likely been very small, as population genetic theory demonstrates that heterozygosity is lost at a slow rate per generation of ΔH = −1/2N e (Nei, 1975).With boreal toads representing a relatively long-lived species (>10 years; Thompson, 2019), only three to five generations have elapsed since declines began, resulting in little erosion of genetic diversity.Therefore, we hypothesize that historical processes of colonization and subsequent isolation from the rest of the species' range are the predominant processes contributing to the low genomic diversity we observed.
Regardless of the relative importance of these processes, the limited genomic diversity observed here may impact the ability of extant toad SRM populations to adapt and survive ongoing disease threats, as standing genetic variation is important for resilience in wildlife species (Kardos et al., 2021;Savage et al., 2015).
Taken together, our results point to an SRM boreal toad population that was founded with, and has maintained, low levels of genomic diversity by sharing the same small pool of genetic variants across a relatively unstructured population.These results achieve our first two objectives by describing the genomic population structure and diversity within the SRM.Low population structure within the SRM indicates that any translocations of individuals among regions could be considered with negligible risk of outbreeding depression, as regions show little genetic differentiation among them.This information could facilitate reintroductions and translocations, as well as more efficient management of captive breeding populations.Additionally, the low genomic diversity highlights the need to reduce the potential for inbreeding depression and increase adaptive potential (DeCandia et al., 2018;Funk et al., 2019Funk et al., , 2021;;Hohenlohe et al., 2021).Because diversity is low across the SRM, it may become necessary for managers to consider source populations from outside the SRM to provide an influx of standing genetic variation into extant boreal toad populations.This is particularly true if non-SRM boreal toads have higher fitness when challenged with Bd, as field studies suggest (Hardy et al., 2022;Hossack et al., 2020;Pilliod et al., 2010).Future laboratory infection trials using wild sources of SRM and non-SRM boreal toads will add valuable information for conservation decision-making for SRM boreal toads.

| Landscape connectivity
Landscape connectivity analyses revealed several habitat factors that may limit or enhance boreal toad movement and gene flow across the SRM.We found support for all five of the minimally correlated covariates we proposed, but three of those five effects were in a direction contrary to our predictions.High annual precipitation, high tree canopy cover, high heat load, warmer minimum winter temperatures and higher impervious surfaces all appeared to restrict connectivity.While high heat loads and higher impervious surfaces make intuitive sense for restricting gene flow of an amphibian, the other results are not as obvious.It is possible that higher annual precipitation restricts gene flow because the majority of precipitation falls as snow in the range of SRM boreal toads.Higher snow fall may persist longer in the spring and begin earlier in the fall, preventing or restricting toad movement in these areas.Alternatively, annual precipitation and higher canopy cover could be masking what we term a 'goldilocks effect' whereby boreal toads prefer areas of high precipitation and canopy cover, making it unnecessary to disperse from these areas.Similarly, the presence of abundant, favourable habitats did not promote gene flow of mountain goats (Oreamnos americanus) in the Cascades Mountains of Washington (Shirk et al., 2010).We note that our landscape genomic analyses may lack sufficient resolution to detect meaningful environmental impacts on gene flow considering the low genetic structure across our study area.In other words, there is so little genetic variation among regions that meaningfully attributing that small amount of variation to environmental factors is likely difficult in our case.We are confident, however, that these findings contribute to our third objective.Managers may now use information about environmental factors that enhance or restrict gene flow across the landscape to guide decisions about selecting sites for reintroductions.For example, conservation of habitats promoting natural gene flow, such as habitats at intermediate elevations with high tree cover and low impervious surface and heat load, would be expected to increase the probability of individuals colonizing other nearby suitable habitats and maintaining gene flow between reintroduced and extant populations.In fact, a prior decision-making process suggested that selecting sites where boreal toads had a greater chance of colonizing other areas yielded the highest expected number of boreal toad populations among various management scenarios (Gerber et al., 2018).

| Adaptive divergence and disease
We developed the first North American toad reference genome and analysed adaptive divergence of boreal toads across the SRM to test for evidence of selection from recent Bd-related population crashes, as well as other broad-scale environmental factors that may lead to longer-term local adaptation.Adaptive landscape genomics using GEA revealed almost an order of magnitude more putatively adaptive SNP outliers related to Bd intensity (i.e.310 outliers related to Bd prevalence and time since first detection) than to all other landscape factors combined (36 outliers related to temperature, precipitation, forest cover, riparian habitats, heat load and urbanization; Figure 4, Tables S2 and S3).These findings suggest Bd is a much stronger selective force on boreal toad populations in the southern Rocky Mountains over less than 10 generations than longer-term evolutionary responses to local habitat conditions.Similarly, recent conservation and landscape genomics work on Tasmanian devils found a strong signature of selection caused by disease (devil facial tumour disease), which largely erased the signature of selection related to environmental factors (Fraik et al., 2020).These findings demonstrate that wildlife diseases may act to shape overall population and landscape genomic patterns, particularly in terms of adaptive genomic variation, and may outweigh the effects of local adaptation to different habitat types over longer timescales (Brugger et al., 2019;Guido et al., 2007).
Our new boreal toad reference genome annotation also allowed us to find additional genomic evidence of the strong selective force Bd infection imposes on boreal toad populations in the SRM by investigating the genes and modifier regions involved in adaptation to disease.Bd kills its amphibian hosts by infecting the keratinized tissue of the adult animal's skin, which disrupts osmoregulation and eventually leads to cardiac arrest (Voyles et al., 2009).We identified outlier SNPs putatively under selection that are in or linked to genes and modifier regions involved in fighting skin infections and eventual cardiac failure caused by lethal Bd infections in amphibians.For example, SNPs associated with Bd prevalence and/or years since first detection were linked to genes and modifiers related to skin (GJB3) and keratin (ADAM9), heart development (HHEX, DUSP6, PHRF1 and NIPBL) and heartbeat regulation (ATR, DSG2 and BIN1).We also detected genes related to immune response (NOX1, CNIH1, FADD and REL), including genes involved in B-cell production, which have been shown to be associated with genomic immunity to Bd infections in other amphibians (Ellison et al., 2015(Ellison et al., , 2020)), and that was associated with outlier SNPs flagged by RDA (Figure 4, Tables S2 and S3).Further investigation of genes with potential links to immune response (e.g.major histocompatibility complex, MHC; Trujillo et al., 2021) using experimental and targeted genomic approaches such as genome-wide association studies (GWAS; Wray et al., 2007) may be used to assess whether local adaptation to Bd is occurring for populations of boreal toads both inside and outside the SRM.
Our findings related to our fourth objective which investigated adaptive divergence and the possibility of outbreeding depression indicate that selection due to the environment is swamped by the strong pressure of Bd.Combined with our findings of low differentiation among regions, local adaptation to environmental factors is likely minimal across the SRM, limiting the potential for significant unintended outbreeding depression (Frankham et al., 2011) in offspring from matings between individuals from different regions when using translocations or captive breeding as conservation management strategies.
Results stemming from our fifth objective to investigate signs of selection to Bd may prove to be important in future management decisions.We show that there is a strong signature of selection related to Bd and that these outliers are near regions of the genome that may be important for defence against disease.Therefore, it may be beneficial to choose translocation sources from regions that have experienced Bd without signs of population crashes, versus those that have had Bd and are at low numbers, or those that may have larger numbers but are Bd naïve.In other words, the only risk of outbreeding depression from translocations or captive breeding may be in the form of mixing toads that are poorly adapted to Bd with those that are well adapted (i.e.extrinsic outbreeding depression caused by introducing maladapted genotypes).Further research is needed to confirm the importance of specific genes related to Bd tolerance or resistance and to link specific variants to populations to gauge their susceptibility using GWAS.We view these results as a critical start for understanding the evolution of disease resistance or tolerance in SRM toads and hope that they will inform management to ultimately improve the resiliency of boreal toads to Bd with future research and associated management.

| CON CLUS ION
Wildlife disease has become a critical threat to global biodiversity, especially within some taxa (e.g.amphibians; Brearley et al., 2013, DeCandia et al., 2018, Fraik et al., 2020).We developed the first North American toad reference genome and used it in combination with reduced-representation, cost-effective, ddRAD sequencing of hundreds of boreal toads from declining populations across the SRM to identify portions of the genome likely to be under selection for disease resistance, as well as to assess levels of gene flow, genomic diversity and effective population sizes (Hohenlohe et al., 2021;Luikart et al., 2003).The lack of structure among extant populations across the SRM could provide managers with more options for source populations for translocation that extend the nearest-neighbour approach employed to date.These source individuals can be used to bolster existing populations (i.e.genetic and demographic rescue) or initiate populations at unoccupied sites that have landscape features likely to promote natural gene flow and colonization.Additionally, our findings will help improve the efficiency of captive breeding efforts and provide opportunities for experimental cross-breeding or servation team.Genomic information could significantly increase the success of these conservation efforts by describing population genomic parameters to guide source and recipient population selection for assisted gene flow and rescue.Additionally, identifying potential signals of adaptation to disease or other environmental factors can aid captive breeding programmes in maintaining more resilient lineages.Using a conservation and landscape genomics approach, we aimed to address five objectives to inform the management and recovery of imperilled boreal toad populations in the SRM.(1) We estimated the extent of population structure across the SRM to identify potential source populations for captive breeding, reintroduction and genetic rescue efforts.(2) We measured the extent of genomic variation and effective population sizes of remaining populations.(3) We investigated key landscape factors that may influence boreal toad connectivity in the SRM to assess the potential for natural gene flow and recolonization.(4) We tested for evidence of divergent selection caused by environmental heterogeneity to assess the potential for outbreeding depression during genetic rescue efforts due to local adaptation.(5) Finally, we tested for loci under divergent selection related to variation among sites in Bd exposure that may be related to susceptibility to infection.We developed an annotated boreal toad reference genome to investigate the function of loci identified as under selection in objectives 4 and 5.

F
Study area for boreal toad conservation and landscape genomics in the southern Rocky Mountains (SRM) of Colorado and southern Wyoming, USA.Landscape factors for connectivity and adaptive divergence analyses consisted of (a) elevation, (b) minimum temperature of the coldest month, (c) maximum temperature of the warmest month, (d) annual precipitation, (e) per cent tree canopy cover, (f) compound topographic index of wetness, (g) heat load index, (h) topographic roughness, (i) streams and rivers, (j) per cent impervious surface and (k) roads (warmer red colours = lower values; cooler green colours = higher values).A total of 231 toads were sampled and genotyped across 36 sites in 2018 and 2019, with yellow sites = 1-5 toads genotyped, orange = 6-10 toads genotyped, red => 10 toads genotyped and white = historical sites where no toads were found in 2018 or 2019.Sampling sites were grouped into 15 regions for population genomic measures, consisting of Wyoming (WY), Elkhead Mountains (Elkhead), Front Range north (FR_N), Zimmerman Lake (Zim), Rocky Mountain National Park (RMNP), Gore Range north and south (Gore_N and Gore_S), I-70 north and south (I70_N and I70_S), Grand Mesa (Grand), Elk Mountains (Elk) and Mosquito Range (Mosq); and Sawatch Range north, central and south (Saw_N, Saw_C and Saw_S).
) may be barriers to adult toad movement and dispersal TA B L E 1(Continued)   the lowest single-covariate deviance information criterion (DIC;Spiegelhalter et al., 2002) among correlated sets of variables(Dormann et al., 2013;Row et al., 2017;Zimmerman et al., 2022).Calculation of the covariance parameter of GWM probability distribution is complicated by non-obvious correlation among sampled locations of raster covariates and can prevent models from converging(Zimmerman et al., 2022).All model fitting was performed in R with two independent chains of 150,000 Markov chain Monte Carlo (MCMC) iterations, with the first 75,000 iterations discarded as a burn-in period, using a random-walk Metropolis-Hastings sampler.Convergence of chains was evaluated through visual inspection of trace and density plots and calculation of the Gelman-Rubin diagnostic(Gelman & Rubin, 1992).Model comparison was based on DIC scores, where lower DIC values indicated better model fit.

E 2
Discriminant analysis of principal components (DAPC) showed low overall population structure of boreal toads across the southern Rocky Mountains (SRM), but there were five peripheral regions for which all individuals could be assigned with high probability to their own groups: the Wyoming (WY), Elkhead (Elkhead) and Front Range north (FR_N) regions in the north; the Mosquito Range (Mosq) region in the south; and the Sawatch Range north (Saw_N) region in the centre.Colours represent the location of different regions (blue to red = regions from north to south).DAPC axes 1 and 2 explain 7.8% and 1.1% of the genomic variation respectively (i.e.8.9% total).

TA B L E 2
Genomic diversity of boreal toads in the southern Rocky Mountains (SRM), including sample sizes (N), number of polymorphic loci, observed and expected heterozygosity (H OBS and H EXP ), inbreeding coefficients (F IS ), nucleotide diversity (π) and effective population sizes (N e ) with 95% confidence intervals (CI) calculated using the linkage disequilibrium (LD) method with NeEstimator corrected for number of chromosomes.may be occurring within many boreal toad populations in the SRM of Colorado and southern Wyoming, consistent with recent population crashes due to Bd. Demographic analyses using the SFS also suggest boreal toad declines may have occurred over the last 10-100 K years (FigureS3), indicating that declines over the last ice ages have also contributed to low overall levels of genomic diversity.
; Muths TA B L E 3 Landscape connectivity analyses using spatial autoregressive models (SAM) and generalized Wishart modelling (GWM) of boreal toads in the southern Rocky Mountains (SRM).

F
I G U R E 4 Genotype-environment associations (GEA) tested using redundancy analysis (RDA) shows putatively adaptive outlier single nucleotide polymorphisms (SNPs) (red circles) consisting of (a) 310 outliers related to Batrachochytrium dendrobatidis (Bd) prevalence (PREV) and years since first detection (YSFD); and (b) 36 outliers related to 9 other landscape factors consisting of annual precipitation (AP), minimum temperature of the coldest month (MTCM), maximum temperature of the warmest month (MTWM), per cent tree canopy cover (CC), compound topographic index of wetness (CTI), heat load index (HLI), streams and rivers (STR), roads (RDS) and per cent impervious surface (IMP).Putatively adaptive outlier SNPs are >2.5 standard deviations from mean RDA loadings.
Landscape connectivity surface represented by a gene flow (log-conductance) surface inferred from the top-ranked spatial autoregressive model (SAM) and generalized Wishart modelling (GWM).Positive values correspond to landscape facilitation of movement and gene flow, while negative values correspond to landscape resistance.Black points represent 18 study sites used with ≥5 toads sampled per site (186 toads total).X and Y axes are universal transverse Mercator units (UTMs) in metres (m), Zone 13.

Region N # Poly- morphic loci H OBS H EXP F IS π N e (95% CI) a
a Estimates of N e were unreliable for N < 5.