Global change is impacting biodiversity across all habitats on earth. New selection pressures from changing climatic conditions and other anthropogenic activities are creating heterogeneous ecological and evolutionary responses across many species' geographic ranges. Yet we currently lack standardised and reproducible tools to effectively predict the resulting patterns in species vulnerability to declines or range changes. We developed an informatic toolbox that integrates ecological, environmental and genomic data and analyses (environmental dissimilarity, species distribution models, landscape connectivity, neutral and adaptive genetic diversity, genotype‐environment associations and genomic offset) to estimate population vulnerability. In our toolbox, functions and data structures are coded in a standardised way so that it is applicable to any species or geographic region where appropriate data are available, for example individual or population sampling and genomic datasets (e.g. RAD‐seq, ddRAD‐seq, whole genome sequencing data) representing environmental variation across the species geographic range. To demonstrate multi‐species applicability, we apply our toolbox to three georeferenced genomic datasets for co‐occurring East African spiny reed frogs ( Our framework sets the stage for large scale, multi‐species genomic datasets to be leveraged in a novel climate change vulnerability framework to quantify intraspecific differences in genetic diversity, local adaptation, range shifts and population vulnerability based on exposure, sensitivity and landscape barriers.
Global climate change is affecting biodiversity in unprecedented ways, compounded by other anthropogenic impacts such as habitat degradation, fragmentation and loss (IPBES, 2019). For example, increased temperatures and frequencies of extreme climatic events are predicted to create new selection pressures by rapidly altering resource availability, exposure to pathogens and the structure and functioning of trophic networks for many species (Bellard et al., 2012;Hoffmann & Sgrò, 2011;Pinsky et al., 2019). How species respond to these new selection pressures depends on their 'vulnerability' (IPCC, 2007), which is defined as the combination of the stress to which a system is exposed, its sensitivity and aspects of its adaptive capacity (Foden et al., 2019). Climate change vulnerability assessments first emerged in the 1990s, as a tool that accounts for aspects of natural hazard and disaster planning, climate change effects and endangered species research. In their early iterations, vulnerability was mainly focused on people and communities (IPCC, 2014), though this was later applied to species and ecosystems.
Until recently, accounting for differences in vulnerability among populations within species was largely ignored in climate change vulnerability assessment approaches. However, highlighting intraspecific populations that are most at risk of local extinction, or identifying those with pre-adapted genotypes that can be sources for assisted gene flow and evolutionary rescue (Bell & Gonzalez, 2009), could greatly improve biodiversity conservation management by safeguarding populations and genetic diversity beneficial for resilience to future environmental change (Hoban et al., 2021(Hoban et al., , 2022;;Laikre et al., 2010). Neutral genetic diversity is important in this respect as it provides the basis for future evolution and could become selected upon when environmental conditions or geographic distributions change.
Populations with higher neutral genetic diversity may therefore have a higher chance of supporting individuals with advantageous mutations or traits (Ørsted et al., 2019). To address the risk of local extinction we refer throughout this manuscript to 'Exposure' (i.e. the nature, magnitude and rate of environmental change), 'Sensitivity' (i.e. the underlying neutral and adaptive genetic diversity which may buffer against environmental change) and 'Landscape barriers' (i.e. limitation to track favourable environmental conditions (Parmesan, 2006;Pecl et al., 2017) and to potentially spread beneficial neutral and adaptive genomic variation (Razgour et al., 2018)). We emphasise here that 'Landscape barriers' is a proxy for potential spread of genomic variation and does not account for dispersal capacity or number of generations across the analysed time periods. During the past decades, the dominant approach to climate change vulnerability assessments were based on forecasts of how species ranges are predicted to change using species distribution models (SDMs; Barbet-Massin et al., 2012;Elith & Leathwick, 2009;Guisan & Thuiller, 2005;Pacifici et al., 2015;Urban, 2015), in some cases refined using genetic data to build SDMs independently for intraspecific populations that have divergent ecological niches (e.g. Bittencourt-Silva et al., 2017;Collart et al., 2021;Ikeda et al., 2017). However, even when accounting for neutral population structure, a major limitation of these approaches has been that intraspecific local adaptation and differential responses to climate change have been largely ignored, potentially leading to inaccurate predictions of future distributions and misplaced conservation efforts (Foden et al., 2019;Hällfors et al., 2016). Adaptation to local environmental conditions is widespread across the tree of life (Hereford, 2009), and the geographic distribution of adaptive variation likely plays a fundamental role in the ability of populations within species to respond to global change (Capblancq et al., 2020;Exposito-Alonso et al., 2018, 2022;Forester et al., 2022). Assessing SDMs together with local adaptation, neutral genetic diversity and landscape connectivity and potential gene flow (e.g. Brennan et al., 2022;McGuire et al., 2016;Parks et al., 2022) is therefore essential to understand geographic differences in vulnerability under future global change scenarios.
Recent calls were made in the emergent field of climate change genomics (Lancaster et al., 2022) for the integration of genomic data to improve the accuracy of climate change vulnerability assessments (Capblancq et al., 2020;Fitzpatrick & Keller, 2015;Nadeau & Urban, 2019;Pauls et al., 2013;Waldvogel, Feldmeyer, et al., 2020).
Conceptual and analytical developments enabling the incorporation of intraspecific adaptations across species ranges (e.g. Aguirre-Liguori et al., 2021;Bay et al., 2018;Forester et al., 2023;Razgour et al., 2018Razgour et al., , 2019;;Ruegg et al., 2018), and phenotypic plasticity (Benito Garzón et al., 2019) have led to major advances in our ability to assess how population vulnerability varies across species ranges. Despite these recent advances, we lack practical and integrative tools to implement analyses across multiple taxonomic groups and geographic regions (see Pinsky et al., 2022). Due to the high multidisciplinarity and diversity of analyses required for most integrated climate change vulnerability assessments, researchers often tailor their approach to their own study system, without creating standardised data structures and code that can be applied more widely to any system (see Johnston et al., 2023;Waldvogel, Schreiber, et al., 2020). et al. (2018, 2019) to leverage information obtained from the raw data to estimate 'Exposure' (estimated from the magnitude of predicted climate change), 'Sensitivity' (estimated from both neutral and adaptive genetic diversity) and 'Landscape barriers' (the estimated limitations for future distributional shifts and evolutionary rescue given predicted future climate change). Landscape barriers is not included in the IPCC 4th assessment but may be particularly relevant for populations without sufficient standing genetic variation to adapt in-situ quickly, or for long-lived species with long generation times where shifting their range is a more likely response to environmental change than rapid adaptation, unless a sufficiently large fraction of individuals already possess pre-adapted genotypes (Razgour et al., 2018). Together, exposure, sensitivity and landscape barriers are mapped separately to identify populations with lower or higher scores for each metric and also combined as an average or custom combination to predict population vulnerability to global change across a species' range. We define population vulnerability as the likelihood that a population will become locally extinct due to global change impacts rather than other anthropogenic pressures such as overharvesting. Our population vulnerability metric is an approximation given the available genomic, spatial and environmental data, but it should be understood that it does not include any kind of population viability analysis.
Our overarching goal was to build an informatic toolbox to predict population vulnerability to global change, integrating and generalising code to make analyses applicable to any suitable population genomic dataset. To build our toolbox we expanded upon two recently published conceptual and analytical frameworks (Razgour et al., 2018(Razgour et al., , 2019) ) for climate change vulnerability assessments. Each of the two frameworks integrates genomic and environmental data to assess climate change vulnerability by incorporating a combination of SDMs, landscape connectivity analyses (using electrical circuit theory) and genetic diversity (neutral and adaptive). Candidate genomic regions under selection are identified using genotypeenvironment association (GEA) methods, which may be validated in our simulation scripts (see Figure 1, Figure S1) following a similar approach to Salmón et al. ( 2021) using randomisations and permutation tests (see Section 2.3 for full details). This information is then used to quantify 'genomic offset' per population based on the predicted mismatch of locally adapted genotypes to the predicted future climates in their current location, indicating the potential levels of their future maladaptation (Fitzpatrick & Keller, 2015). Therefore, FIGURE 1 Conceptual and analytical framework for the Life on the edge toolbox, incorporating 'Exposure' (current and projected future species distribution models (SDM) and Species their dissimilarity), 'Sensitivity (adaptive and neutral sensitivity), 'Landscape barriers' (predicted population connectivity) to predict a final 'Population vulnerability' metric for each population (which is a weighted combination of the other metrics). Software packages used are denoted in blue text (LFMM, latent factor mixed models; RDA, redundancy analysis). Table 1, Table S1). We integrate SDMs to assess dissimilarity in future environmental conditions ('Exposure'), landscape connectivity ('Landscape barriers'), standing genetic diversity ('Neutral sensitivity') and adaptive genetic diversity ('Adaptive sensitivity'), to estimate population vulnerability for all unique geographic locations with samples across species ranges. We highlight areas where all vulnerability metrics are in the upper and lower quantiles of the species results and highlight the highest and lowest metric values per population to guide conservation priorities. Our highly flexible toolbox makes it possible to 'plug in' any species with suitable data so that standardised analyses and comparisons across different taxa and regions can be readily made. The LotE toolbox therefore establishes the backbone of a generalised framework that aims to stimulate a new wave of data synthesis, increase reproducibility and standardise reporting for population level and species-level climate change vulnerability assessments (Waldvogel, Feldmeyer, et al., 2020). For ease of interpretation, Figure 2 summarises the main inputs, analyses and outputs for the toolbox.
To assign how metrics are quantified for exposure, neutral sensitivity, adaptive sensitivity, landscape barriers and population vulnerability, user-defined thresholds should be specified to the params file in a comma separated list for each individual metric (see Table S1).
For example, to assign the neutral sensitivity metric in increments of 0.1 so that a genetically diverse population with a high nucleotide diversity value of 0.225 gives a low neutral sensitivity value (i.e. =1) and a low nucleotide diversity value of 0.025 gives a high neutral sensitivity value (i.e. =10), define the variable neutral_sensitivity_nu-cleotide_diversity_thresholds as '0, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225'.
Life on the edge integrates diverse genomic and spatial analyses to create metrics of exposure, neutral and adaptive sensitivity and population vulnerability. The run_life_on_the_edge.sh wrapper script can be used to run each part of the toolbox as required by calling the desired functions. If users wish to provide carefully curated inputs without parameterising all steps of the toolbox (e.g.
Summary of modelling outputs and methodologies incorporated for the LotE toolbox. Modelling output General methodology employed Specific steps and [functions] Exposure • Species distribution models (SDMs) • Downloading and preparing spatial [prepare_spatial_data()] and environmental data [prepare_environmental_data()] • Spatially rarefying presence data and generating background (pseudoabsence) data for SDMs [spatially_rarefy_presences_and_create_background_data()] • Building and evaluating ensemble SDMs for current and future conditions [sdms_biomod2] • Calculating "Exposure" for each population [exposure()] Sensitivity • Identify climate adaptive loci • Investigate the sensitivity and robustness of identified climate adaptive loci using simulations • Categorise sampled individuals in adaptive ordination space • Quantify adaptive genetic sensitivity • Quantify neutral genetic sensitivity • Rerun SDMs based on categorised individual • Impute missing data [impute_missing_data()], LFMM [gea_ lfmm()], RDA [gea_rda()] to identify adaptive loci • RDA to categorise climate adapted individuals [gea_rda_individual_categorisation()] • Count numbers of adapted individuals per population [quantify_local_adaptation()] • Validate adaptive signal in the genomic data using randomisation simulations [-00_parameter_exploration-.R, -01_empirical_ data-.R, -02_randomise_data-.R, -03_perform_simulations-.R, -04_evaluate_significance-.R] • Calculate neutral genetic sensitivity based on nucleotide diversity (default) or heterozygosity [neutral_sensitivity()] • Calculate adaptive sensitivity based on results from GEAs and local adaptation/genomic offset results [adaptive_sensitivity()] • Build refined SDMs using only individuals that are categorised in specific adaptive ordination space [adaptive_sdms()] Landscape barriers • Run circuitscape to assess change in movement potential between current and future environmental conditions • Create circuitscape input and parameter files [create_ circuitscape_inputs()], categorise and plot "Landscape barriers" per population results [landscape_barriers()] Population vulnerability • Create final population vulnerability scores and outputs • Assimilate Exposure, Sensitivity, and Landscape barriers results to create population vulnerability metric and map all metrics [population_vulnerability()] • Create final summary PDF of all results [summary_pdfs()]
Note: Specific steps for each method are detailed, along with toolbox functions used in each case.
SDMs, imputed genotype data, GEA analyses, recommended for non-expert users), then parameters can be modified so that the toolbox recognises this, and the parameterisation steps for these analyses will be skipped (see adjustable parameter list in Table S1 for details). In Supporting Information Text S1, we provide a detailed overview of the input data and structure required and the main methods adopted by the LotE toolbox at each step, how these are implemented and which aspects are modifiable by the user.
Table 1 provides further details on the modelling output, steps taken to achieve that output and which functions in the toolbox are used at each step. The functions and scripts used (written in R, bash and Julia) are available in the LotE website (h t t p s : / / c d -b a r r a t t . g i t h u b . i o / L i f e _ o n _ t h e _ e d g e . g i t h u b . i o / ) and an accompanying vignette on example usage is available at h t t p s : / / c d -b a r r a t t . g i t h u b .
i o / L i f e _ o n _ t h e _ e d g e . g i t h u b . i o / V i g n e t t e . The toolbox enables HPC parallelisation, which is particularly useful for computationally intensive steps. We provide example benchmarking times for analyses to complete (Table 2) to provide users with an overview of processing times.
Users may store the toolbox in any location but need to adhere to the directory structure shown in Figure 2. The -data-folder contains four main directories, one of which contains mapping data (world shapefile files available from h t t p : / / n a t u r a l e a r t h d a t a . c o m , already provided), another the environmental data, the other two containing folders for each species' genomic and spatial data. To run the toolbox, appropriate environmental predictor data must first be downloaded (e.g. global data in .tif raster format) and stored in the environmental_data folder, the params (Params.tsv) file must be populated with the relevant parameters for each step of the analyses (see Table S1 for a description of each parameter and the vignette).
Environmental data downloads can be automated if requested using the geodata R package (Hijmans et al., 2023) by setting the env_data_download parameter to 'yes' in the params file as well as providing the resolution requested (climate_res_download), the shared socioeconomic pathway scenario (ssp_scenario_download), general circulation model (gcm_download) and the time of the future projection (time_proj_download). The three input files required to run a species dataset are the spatial coordinates of the genomic samples (Sample name, decimal Longitude, Latitude, in .csv format), and two standard PLINK (Purcell et al., 2007) formatted files (.ped and .map) for the genomic data (Figure 2). Conversion from other formats for genomic data such as the widely used VCF (variant call format) may be converted to PLINK format readily using available tools (e.g. Danecek et al., 2011). Details of data naming conventions can be found in the vignette, along with example input data available in the DRYAD repository (h t t p s : / / d a t a d r y a d . o r g / s t a s h / d a t a s FIGURE 2 Main inputs and data (yellow boxes), analyses (blue box) and outputs (green boxes) of the LotE toolbox. 'Species_binomial' is the name of the analysis for any given species, using genus name followed by species name separated by an underscore. Directory names are highlighted in bold, 'Exposure', 'Sensitivity' (including neutral and adaptive sensitivity) and 'Landscape barriers' become populated with the relevant output files for each analysis upon running LotE, which are then used to calculate output metrics per population. Information on specific R functions within the blue box and how they interact with the output directories can be found in Figure S1. The -scripts-and R_functions folders contain all the toolbox scripts and functions, and the -outputs-folder stores all output files in relevant subdirectories when running the toolbox. Blue lines represent locations for input files, dotted blue lines represent locations for input files if the user wants to circumvent the full toolbox workflow with their own input data (e.g. pre-prepared SDMs, a list of adaptive SNPs so that GEA analysis is unnecessary, imputed missing genotype data, or an already parameterised circuitscape input layer). e t / d o i : 1 0 . 5 0 6 1 / d r y a d . 2 r b n z s 7 t 4 ). To provide reliable outputs using LotE, we advise that the genomic samples cover an adequate range of environmental conditions that the species as a whole experiences (i.e. a representative proportion of the species range, see Vajana et al., 2023), and the number of SNPs should be sufficient to detect accurate signals of local adaptation (minimising false positives), covering a large proportion of the genome, though this will depend on the genomic characteristics and the degree of local adaptation in each species. To explore and validate signals of local adaptation detected in datasets, the sensitivity and simulation scripts we provide follow the approach outlined in Salmón et al. (2021). In brief, 100 simulations of the empirical data are created with randomised genotype-environment relationships, GEA analyses are performed on each of these simulations, tracking the p-values of all SNPs and then the adaptive signal in the empirical data is determined using a significance threshold (using z-scores; e.g. a significance threshold of >0.95) for the empirical data against the simulations data to identify statistically significant SNPs that are above this. These statistically significant SNPs may then be used for all local adaptation analyses going forward (using the 'use_only_statistically_significant_snps' set to 'yes' and 'which_loci' option set to '0' in params). For more details see Supporting Information Text S2, which also gives details on how to perform sensitivity analyses using LotE given different parameter combinations of thresholds for determining candidate SNPs using our implemented GEA methods.
For environmental data, we provide global terrestrial Worldclim2 data (30 arc seconds, clipped to East Africa) for testing purposes, but as mentioned previously, this can be automatically downloaded (at lower resolutions using the geodata R package; Hijmans et al., 2023) if required. Many users will require different spatial extents or higher resolution data if available for building accurate SDMs and detecting fine-scale environmental variation and local adaptation across populations on land. Marine or freshwater data could also be used here if focal taxa are non-terrestrial, however it is important that the environmental data are standardised to be the same spatial resolution and extent. These data can be, for example, georeferenced .tif files (e.g. bioclim or elevation) downloaded from public databases (e.g. Worldclim2 or CHELSA (Fick & Hijmans, 2017;Karger et al., 2017) for current and selected future conditions and SSP scenarios), or any other predictor data relevant to the study species that is available in raster format (e.g. land cover data). The user should decide on which spatial resolution, future time period and scenario is needed. The georeferenced genomic data will serve as known presence points for SDMs, which can be integrated with GBIF data that are cleaned and finalised using the toolbox, as well as the generation of appropriate pseudoabsence data. All input environmental data needs to be at the same spatial resolution and extent, an R script (00_process_environmental_data.R) is provided to assist the user in setting up environmental predictor data in the correct format for Worldclim2 and CHELSA data (i.e. separating a multi-band file representing many predictors together into individual predictor .tif files for current and future conditions). We recommend that before inputting genomic data (.ped and .map files), best practices (Paris et al., 2017) are followed to maximise polymorphism in input data while reducing potential 'false' loci caused by over-or under merging SNP loci, as well removing poorly sequenced individuals causing high levels of allelic dropout and lower numbers of loci and SNPs (Cerca et al., 2021).
Species #pops #SNPs #cells SDM details Runtime Afrixalus delicatus 14 8961 329,222 9 predictors, 510 models (CTA, ANN, RF, GAM, Maxent) 8 h 49 min [+16 h 32 simulation validation] Afrixalus sylvaticus 20 12,842 62,972 9 predictors, 510 models (CTA, ANN, RF, GAM, Maxent) 3 h 51 min [+11 h 43 simulation validation] Afrixalus fornasini 32 7309 304,616 9 predictors, 510 models (CTA, ANN, RF, GAM, Maxent) 8 h 58 min [+16 h 14 simulation validation] Myotis escalerai 67 18,356 345,042 6 predictors, 204 models (Maxent, CTA) 94 h 27 min Myotis crypticus 41 20,750 122,958 6 predictors, 204 models (Maxent, CTA) 25 h 48 min Note: #pops refers to the number of unique geographic locations (i.e. populations), #SNPs refers to the number of bi-allelic SNPs in the dataset, #cells refers to the number of grid cells (i.e. pixels) in the rasters used for spatial analyses. Genomic data pre-processing times are not included, and the extra simulation validation runtimes are listed. All analyses were each run as separate jobs on a single HPC cluster core (8 threads, non-volatile memory express storage) with 80GB RAM allocated.
TA B L E 2 Example benchmarking time for completion of the toolbox on the published data utilised for the three focal species.
Several important decisions relevant to each species dataset are required throughout the LotE toolbox. If users are not comfortable with the parameterisation of specific steps within LotE-particularly building SDMs, performing GEA analyses to identify candidate SNPs under selection, imputation of missing genotype data or the parameterisation of input connectivity layers for circuitscape analysis, we recommend that these are performed with assistance from relevant expertise as part of multidisciplinary teams, possibly outside of the toolbox and supplied as inputs to LotE. Running analyses without support or full understanding of the conceptual backgrounds and potential pitfalls of each method used in LotE will almost certainly lead to unreliable results. If any of the skip_sdm, skip_impute_genotypes, skip_gea, skip_circuitscape_layer_parameterisation parameters are set to 'yes', the toolbox will not perform these analyses and instead expect the relevant input files in the correct directories (details in Supporting Information Text S1). First, regarding environmental data, the selection of environmental predictors for SDMs and GEAs should be ecologically relevant to the study species. If required, multiple environmental predictor variables may be condensed into principal components (PCs) and those can be used as predictors in the GEAs. Predictor variables should be at a suitable spatial grain to be able to detect signals of local adaptation if they exist, sufficiently variable between sampled populations and expected to influence the distribution and/or genetic diversity of the species in question. The threshold q-values to determine which SNPs are putatively 'adaptive' are recommended to be set at a very conservative (e.g. False Discovery Rate >0.01 for LFMM analyses) and standard deviation from the mean loading (SD >2.5 for RDA analyses) by default to minimise false positives, but this can be modified in the params file if required ('lfmm_FDR_threshold' and 'rda_SD_threshold').
Furthermore, by using the simulation scripts of LotE (built-in by default), it is possible to estimate an acceptable false discovery rate of SNPs for a given dataset and then validate the adaptive signal using our simulation approach. For GEAs, if categorising local adaptations in individuals and populations, we have restricted analyses to two predictors simultaneously to easily parse adaptation to different conditions (e.g. hot-dry and cold-wet). However, more complex scenarios of adaptation to multiple predictors can be assessed by analysing pairs of additional predictors separately (see Barratt, Preißler, et al., 2024). Second, spatial occurrence data (presences) should be checked thoroughly to ensure that incorrect or unrealistic presence data are not included for SDMs (i.e. outside of the native range, or inverted coordinates for example) and that correct taxonomy is followed (e.g. only confirmed species records are included). We have taken measures using the CoordinateCleaner R package (Zizka et al., 2019) to deal with these potential problems, but data should be carefully inspected before analysis and interpretation of results. Additionally, SDMs require consideration of the geographic modelling extent, selection of background (pseudoabsence) data, data partitioning (training vs. testing) strategy and model evaluation in order to follow best practices in the field (see Araújo et al., 2019;Merow et al., 2013;Zurell et al., 2020 for guidelines), and the SDM output itself should be inspected to confirm that it is a reasonable prediction for the species and thus suitable for further use. Similarly, genomic offset predictions can be clipped to a buffer around the known presences to avoid predicting maladaptation in geographic space that is most likely to be unreachable by the species (see Table S1). Third, if using LotE on RAD-seq/ ddRAD-seq type data, an understanding of the types of errors that are associated with these kinds of data and how to minimise them is fundamental-we strongly advise that datasets have been appropriately analysed and curated before performing LotE analyses. Additionally, when assessing neutral and adaptive sensitivity, including the imputation of missing data and accounting for neutral population structure for GEA analysis, decisions are required to test a reasonable number of genetic clusters (k) represented by the data. In the GEA analyses themselves, the thresholds for defining putatively adaptive SNPs are also flexible to enable decisions on how tolerant the user is of false positives (see Forester et al., 2018;François et al., 2016). If there are adaptations to opposing conditions in the same population, this could be a genuine biological signal as a result of local gene flow, or that adaptive equilibrium may not have been reached across the landscape and thus the GEA approaches may not be suitable. In a case such as this, exploring parameter variation for GEA analyses using the simulation scripts we provide may help to more thoroughly evaluate false positives.
Candidate SNPs with low statistical significance compared to simulations are automatically removed from the list of adaptive loci using the 'remove_low_significance_adaptive_SNPs' and 'SNP_sim-ulation_significance_threshold' parameters in the params file (see Table S1), and this will be reported in the summary PDFs to assist evaluation of the adaptive signal for a given dataset.
Users of the LotE toolbox should be proficient in R and bash. As the toolbox utilises several programming languages which in turn require dependencies, correct initial setup is essential. A working installation of PLINK (Purcell et al., 2007) and circuitscape (Anantharaman et al., 2019) is required as well as a recent version of R (4.1.3 or later), a bash shell and a version of Singularity (Kurtzer et al., 2017).
Installation of package dependencies from within R needs to be
The LotE toolbox is fully transparent and parameterisable, with standardised workflows following best practices for running species distribution models (Araújo et al., 2019), and genotype-environment association analyses, (Capblancq & Forester, 2021;Forester et al., 2018). Transparency allows all results to be traced back to the data and helps avoid the toolbox being a 'black box'. We recommend user feedback and sanity checks at several decision-making points of the workflow in order to follow best practices in the relevant subfields of ecology and evolution. Our toolbox can be used as a single pipeline (e.g. from processed sequence and spatial data through to predicting population vulnerability) or in a modular fashion using specific functions. Furthermore, the toolbox offers flexibility, so if users wish to supply their own data (e.g. environmental data, SDMs, imputed genotypes, list of adaptive SNPs or input files for circuitscape analysis to assess landscape connectivity) it is possible to circumvent steps within the toolbox by simply adding the relevant files to the appropriate directories (see vignette).
Empirical datasets to demonstrate the utility of the LotE toolbox 2.7.1 | Afrixalus fornasini, A. delicatus and A. sylvaticus-'Novel' LotE analysis (including genomic data processing) To demonstrate the utility of the LotE toolbox we ran the toolbox in its entirety for three co-occurring East African spiny reed frog species, Afrixalus fornasini, A. delicatus and A. sylvaticus. Field work licences and permits are described in Barratt et al. (2018). We processed georeferenced genome-wide RAD-seq data from Barratt et al. (2018) (SRA accession number: PRJNA472166, Table S2) in Stacks 2 (Rochette et al., 2019), optimising the parameters to maximise information (Paris et al., 2017) and remove poorly sequenced samples (Cerca et al., 2021) and collated spatial data including published data in Barratt et al. (2018) and cleaned records from the Global Biodiversity Information Facility (G B I F . o r g , 2023). Environmental data from Worldclim2 was used-bioclim layers 1-19 and slope (Fick & Hijmans, 2017) and land cover (Schipper et al., 2020) at 30s spatial resolution (recategorised into 9 classes following Razgour et al., 2019 to reduce complexity in the model ('landcov1') see Table S3). We defined our modelling extent to capture the known range of each of the three species across East Africa, encompassing sampled populations across Tanzania, Kenya, Mozambique and Malawi, used variance inflation factors (VIFs) to reduce spatial autocorrelation in input SDM predictor variables, and quantified local adaptation to maximum temperatures of the warmest month (bioclim_5) and precipitation of the warmest quarter (bioclim_18) as these are known to be important predictors of the species distributions (Barratt et al., 2017(Barratt et al., , 2018)). We also predicted genomic offset per population, with the resulting prediction clipped to a 2 degree buffer around sampled populations using the genomic_offset_buff_dist_degrees option in the params file to avoid predicting maladaptation in areas where the species is unable to disperse to. We retained only statistically significant SNPs using our simulation scripts to reduce the likelihood of false negative genotype-environment associations being included in our analyses. We parameterised a cumulative resistance surface using five ecologically relevant variables (current SDM, slope, land cover, bioclim_5 and bioclim_18), with respective weights of 0.25, 0.1, 0.25, 0.2, 0.2 which were defined based on our ecological knowledge of the species and predictor effects on their dispersal.
Exposure, neutral sensitivity, adaptive sensitivity, landscape barriers and population vulnerability were all quantified using the 'defined' option, reading defined thresholds for each variable from the params file to determine scores. Full parameter settings for the Afrixalus fornasini, A. delicatus and A. sylvaticus analyses can be found in Table S4.
Second, we conducted a 'partial' LotE toolbox run using data from
Razgour et al. (2019) (European Nucleotide Archive accession no. PRJEB29086, Table S2) for two European bat species, Myotis escalerai and M. crypticus. We collated spatial (including cleaned records from the Global Biodiversity Information Facility, G B I F . o r g , 2023) and Worldclim2 environmental data at 30s spatial resolution (bio-clim_1, bioclim_4, bioclim_7, bioclim_5, bioclim_6, slope, Fick & Hijmans, 2017), as well as land cover (Schipper et al., 2020) (recategorised into 9 classes following Razgour et al., 2019, Table S3).
We generated background points and built SDMs using biomod2 (Thuiller et al., 2009), setting SDM parameters to match those used in the original manuscript, namely the spatial modelling extent, predictor variables, future time projections and general circulation models, and SDM modelling algorithms as well as the evaluation criteria for SDMs (ROC > 0.8). Processed genomic data (.ped and .map files) from Razgour et al. (2019) were input for the GEA analyses (LFMM and RDA) to investigate local adaptation to hot-dry and coldwet conditions based on maximum temperatures of the warmest month (bioclim_5) and precipitation of the warmest quarter (bio-clim_18). Adaptive SDMs, which are not standard within a normal full LotE analysis but can be generated easily with the adaptive_sdms() 2041210x, 2024, 11, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.14429 by Univ Of California Santa Cruz -UCSC, Wiley Online Library on [26/03/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License function of the LotE toolbox, were built using individuals parsed into 'hot-dry', 'cold-wet' adaptive categories for present and future conditions, again following Razgour et al. (2019). Full parameter settings for the Myotis escalerai and M. crypticus analyses can be found in Table S4. Field work licences and permits are described in Razgour et al. (2019).
Below we provide details on results for each of the main analyses across the datasets we analysed. Results obtained using the LotE toolbox on empirical data here closely matched those of the published data for Afrixalus fornasini, A. delicatus, A. sylvaticus, Myotis escalerai and M. crypticus, demonstrating that our toolbox is robust as well as being able to integrate diverse analyses to assist predictions of population vulnerability to global change.
After well as a final summary PDF (Appendices S2-S4) can be found in the Supporting Information.
After inspecting p-value distributions and adjusting the genomic in the Iberian Peninsula, particularly in more arid regions (southern Spain and Portugal), and expansions predicted in northern Portugal and parts of the Pyrenees (Figure S3A). M. crypticus habitat suitability also decreased under future climate projections, with the Iberian region predicted to be largely unsuitable and future suitability limited to high elevation regions including parts of the Alps and the Pyrenees. Separating the different categories of individuals (hot-dry and cold-wet adapted), both species results broadly matched Razgour et al. (2019, Figure 5). In M. escalerai, suitable environmental conditions for both categories were predicted to shift northwards, substantially affecting the predicted ranges, in particular the habitat suitability for cold-wet genotypes found in northern Iberia, which substantially contracted under predicted future conditions (Figure 5a,b). M. crypticus showed high habitat suitability in northern parts of the Iberian Peninsula, the Pyrenees, which were predicted to contract in future conditions, and parts of southern Europe, which were predicted to shift northwards in the future (Figure S3B). As with M. escalerai, hot-dry and coldwet adapted individuals' habitat suitability was predicted to decrease slightly, particularly for the latter, whose potential suitable range throughout the Pyrenees and central and south-western France almost completely disappeared under future predictions (Figure 5a,b). Full log file outputs from the 'partial' LotE analysis for M. escalerai and M. crypticus can be found in the Supporting Information (Appendix S5).
Current climate change genomics approaches to assess population vulnerability lack practical and integrative tools to implement analyses across multiple taxonomic groups and geographic regions (Pinsky et al., 2022). Adaptive responses to global change are likely to be insufficient for many species (Quintero & Wiens, 2013; Aguirre -Liguori et al., 2021;Bay et al., 2018;Razgour et al., 2018Razgour et al., , 2019;;Ruegg et al., 2018), but these have not been widely applicable.
With Life on the edge, we present a novel, generalised and customisable toolbox which can perform sophisticated analyses in a straightforward and reproducible way. The toolbox enables large scale data synthesis across multiple species and geographic areas and the outputs address ecological and evolutionary responses to future global change and can be used to guide conservation action, including accounting for the differences in specific metrics following the principles of complementarity (e.g. Beger et al., 2015). We envisage LotE as a key tool in the emergent field of climate change vulnerability assessments using genomics, with scope for future development and expansion of the concepts presented in this manuscript. Below, we discuss applications of the LotE toolbox in its current form, its limitations, and future additions, analytical and conceptual developments.
Life on the edge, with its geographical and taxonomic flexibility, contributes to answering several fundamental questions in ecology, evolution and conservation. Outputs from LotE may therefore be utilised for informing the management of populations to maintain genetic diversity (e.g. nucleotide diversity and/or adaptive genetic diversity) within and between populations (e.g. Segelbacher et al., 2022). For example, suitable donor and recipient populations for translocations/evolutionary rescue may be selected based on a combination of nucleotide diversity estimates, local adaptations and spatial connectivity. With these strategies, conservation managers may avoid introducing maladapted individuals to unsuitable climatic conditions (Chen, Grossfurthner, et al., 2022;Chen, Jiang, et al., 2022), thus strengthening the overall individual fitness, genetic diversity and adaptive potential of populations (Frankham, 2015;Frankham et al., 2019). At a regional scale, spatial conservation planning can benefit from the climate change genomics perspective of the outputs generated by the LotE toolbox. For example, any of the metrics for neutral and adaptive sensitivity or population vulnerability could be interpolated spatially to be useful for prioritising areas that support particularly high levels of genetic diversity or vulnerable populations, using common software such as Zonation, Marxan and relatives (Ball et al., 2009;Lehtomäki & Moilanen, 2013;Moilanen et al., 2005).
However, care must be taken when interpolating the data if sampling is uneven, for example using kriging to account for uncertainty.
When data from multiple species in an ecological community are available, evaluating congruence in metrics across species would provide a representative measure of community-level genetic diversity and vulnerability based on genome-wide data (e.g. Schielzeth & Wolf, 2021;Stange et al., 2021). Alternatively, the complementarity between these metrics may be used to inform conservation decisions in a given scenario, for example a population with low landscape barriers and high exposure might benefit from assisted migration, whereas a population with high landscape barriers and low adaptive sensitivity/high neutral sensitivity might be suitable for protected area implementation and habitat restoration along connectivity corridors. This could be scaled up (i.e. across countries and continents) with sufficient data across taxonomic groups and can be used to identify hotspots of vulnerability across taxonomic groups.
Similarly, investigating the ecological, environmental and anthropogenic correlates of LotE outputs can identify causal drivers of observed vulnerability patterns for species and ecological communities (similar to approaches in Howard et al., 2020;Maxwell et al., 2016;Tilman et al., 2017), which can inform broad conservation actions as well as understanding species-specific drivers of declines more thoroughly with population level data. However, care must be used to avoid circularity (e.g. the environmental input layers are not meaningful as predictors of the vulnerability patterns). We further acknowledge that phenotypic plasticity is an important component of adaptive capacity (Foden et al., 2019;Fox et al., 2019;Merila & Hendry, 2014), which should be included for systems where this is known, though in our framework here we do not include this as this information is not available for our example species.
Finally, the adaptive SDMs within the LotE toolbox (which are an optional add-on) can contribute to more realistic estimates of shifts in range suitability under global change scenarios, providing improved predictions of future biodiversity losses that may be offset with appropriate conservation measures (Hoffmann & Sgrò, 2011). Overlooking intraspecific population variability, in particular local adaptation, can result in an overestimation future biodiversity losses (Razgour et al., 2019), and it is increasingly clear that predictive models informed by empirical genomic data provide a more realistic alternative to simplistic modelling approaches that do not account for local adaptation (Bay et 2017;Forester et al., 2023).
In addition to applications of the current toolbox, there are three future conceptual and analytical developments that could expand the framework and its long-term impact and benefit to the research and conservation community-monitoring biodiversity change over time, performing simulations and sensitivity analyses to validate findings, and integrating phenotypic plasticity and functional genomics data.
The toolbox has been purposely designed to be dynamic, so that additional 'modules' may be created and integrated in future versions, thus enabling it to evolve in tandem with the research community and adopt best practices and state-of-the-art tools. New methods or tools (for example a new and improved SDM or GEA package) that supersedes existing methods may be integrated relatively simply into updated versions of the toolbox by creating additional functions to supplement or update the existing modules. This also applies to the type of genomic data that can be analysed, which can be scaled up from short read (RAD-seq) type data to whole genome sequencing data with ease by using alternative tools (e.g.
Korneliussen et al., 2014; McKenna et al., 2010). A primary focus for expanding LotE in the future is for monitoring biodiversity change and population vulnerability over time, for example by linking outputs with global conservation efforts such as the Sustainable Development Goals as well as 'Essential Biodiversity Variables' (Hoban et al., 2022). By integrating time series data and multiple biological replicates, it would be possible to track changes in exposure, sensitivity (adaptive/neutral), landscape barriers and population vulnerability as new data become available. Although suitable population-level genomic and spatial sampling time-series replicates are currently rare (but see Pfenninger et al., 2023), ever improving advances in sequencing technologies and the exponential accumulation of these dat a in online repositories (e.g. NCBI Sequence Read Archive, ENI European Nucleotide Archive, DNA Databank of Japan) make the tantalising prospect of genomics-informed biodiversity monitoring an achievable target in the near future (Bálint et al., 2018; Pfenninger & Bálint, 2022; Taus et al., 2017). Ultimately, automated 'scraping' of public repositories for species datasets as new data become available could provide a near real-time assessment of species' biodiversity status based on the latest available genomic and spatial information, similar to the explosion of data generation, availability and automation to monitor recent epidemiological outbreaks (e.g. h t t p : / / n e x t s t r a i n . o r g ). Second, validating the outputs of LotE using analyses and simulations will substantially improve predictions and confidence intervals in results (Hoban, 2014). LotE is amenable to parallelisation, enabling it to run across a range of parameter settings in the params file, with which results can be harvested and parameter variation effects on results can be explored in detail. Our simulation scripts demonstrate how this may be used to explore how parameter variation affects results and to validate adaptive signals using simulated data, but population vulnerability to climate can be empirically validated with appropriate population trend (e.g. Bay et al., 2018) and common garden experiments data (e.g. Fitzpatrick et al., 2021) where available. The integrative analyses of LotE may complement and enhance similar frameworks investigating genomic offset (e.g. Smith et al., 2021) that on their own do not consider migration and gene flow (e.g. Capblancq et al., 2020; Fitzpatrick & Keller, 2015; Rellstab et al., 2021) and could be evaluated against outputs from comparable frameworks using the same underlying datasets. With respect to input genomic data, we acknowledge there has been significant debate on the power of reduced representation library (RAD/ddRAD-seq) datasets to detect sufficient signals of local adaptation (Lowry et al., 2017, but see Catchen et al., 2017). Ideally, high coverage whole genome sequencing data would be the gold standard for detecting local adaptation in populations, which will improve the ability to detect both weak and strong signals of local adaptations. Given that these kinds of datasets for range-wide population sampling are presently rare, reduced representation library datasets currently offer the most feasible approach to synthesise data across taxa and regions, but Benjelloun et al., 2019).
Third, integrating phenotypic plasticity and functional genomics data are a rich potential avenue for expansion for the LotE toolbox, for example in model systems where there is adequate knowledge on physiological limits for species, or the underlying genetic basis of their functional traits and reaction norms (e.g. Oomen & Hutchings, 2022). Expansion of the toolbox to incorporate this information, particularly to strengthen the 'sensitiv- Wollenberg Valero et al., 2022), especially when available from multiple individuals and populations across a species' range. From a spatial perspective, mechanistic SDMs that explicitly incorporate process that limit species distributions (Kearney et al., 2010;Mathewson et al., 2017) and joint distribution models that potentially incorporate species interactions (Ovaskainen et al., 2016, Poggiato et al., 2021) could be used instead of correlative SDM outputs in the toolbox. widely applicable datasets for our approach. As sequencing technologies continue to improve, we will undoubtedly move towards larger numbers of whole genome sequencing datasets which will provide higher resolution data for assessing genetic diversity and local adaptation in particular.
We introduce the 'LotE' conceptual and analytical framework, a toolbox that facilitates the integration of environmental, molecu- et al., 2018).
Our study brings together authors from a number of different countries but does not include scientists based in the countries where the samples were collected. The toolbox described in this manuscript is mainly methodological and the underlying data are already published. However, as well as building on and citing previous work by local scientists, the toolbox in this manuscript is currently being used for region and species-specific questions in those respective countries, which do integrate local scientists prominently.
2041210x, 2024, 11, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.14429 by Univ Of California Santa Cruz -UCSC, Wiley Online Library on [26/03/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License