An integrated framework to identify wildlife populations under threat from climate change

Abstract Climate change is a major threat to global biodiversity that will produce a range of new selection pressures. Understanding species responses to climate change requires an interdisciplinary perspective, combining ecological, molecular and environmental approaches. We propose an applied integrated framework to identify populations under threat from climate change based on their extent of exposure, inherent sensitivity due to adaptive and neutral genetic variation and range shift potential. We consider intraspecific vulnerability and population‐level responses, an important but often neglected conservation research priority. We demonstrate how this framework can be applied to vertebrates with limited dispersal abilities using empirical data for the bat Plecotus austriacus. We use ecological niche modelling and environmental dissimilarity analysis to locate areas at high risk of exposure to future changes. Combining outlier tests with genotype–environment association analysis, we identify potential climate‐adaptive SNPs in our genomic data set and differences in the frequency of adaptive and neutral variation between populations. We assess landscape connectivity and show that changing environmental suitability may limit the future movement of individuals, thus affecting both the ability of populations to shift their distribution to climatically suitable areas and the probability of evolutionary rescue through the spread of adaptive genetic variation among populations. Therefore, a better understanding of movement ecology and landscape connectivity is needed for predicting population persistence under climate change. Our study highlights the importance of incorporating genomic data to determine sensitivity, adaptive potential and range shift potential, instead of relying solely on exposure to guide species vulnerability assessments and conservation planning.


Appendix 4 -Supplemental Figures
Pages 21-24 The ddRAD library preparation protocol was based on the methodology originally reported by Peterson et al. (2012), with modifications / refinements as described in Manousaki et al. (2016).
Briefly, each of the 95 DNA samples was simultaneously digested by two high fidelity restriction enzymes (RE): SbfI (CCTGCA|GG recognition site), and SphI (GCATG|C recognition site), both sourced from New England Biolabs (NEB, UK). Digestions were incubated at 37°C for 50 min, using 10 U of each enzyme per microgram DNA in 1× CutSmart Buffer (NEB), in a 6 µL total reaction volume. After cooling the reactions to room temperature, 3 µL of a premade barcode / adapter mix was added to the digested DNA, and incubated at room temperature for 10 min.
This adapter mix comprised individual-specific barcoded combinations of P1 (SbfI-compatible) and P2 (SphI-compatible) adapters at 6 nM and 72 nM concentrations respectively, in 1× reaction buffer 2 (NEB). Adapters were compatible with Illumina sequencing chemistry (see Peterson et al., 2012). The barcoded adapters were designed such that adapter-genomic DNA ligations did not reconstitute RE sites, while residual RE activity limited concatemerization of genomic fragments. The adapters included an inline five-or seven-base barcode for sample identification. Ligation was performed over 3 hr at 22°C by addition of a further 3 µL of a

Bioinformatics for genomic data analysis
The MiSeq generated reads were processed using a software pipeline designed specifically for RAD analysis, Stacks (v.1. 17;Catchen et al., 2013). First, the 'process_radtags' function was used to demultiplex the individual samples. During this step sequence reads with quality scores below 10, missing either restriction site or with ambiguous barcodes were discarded. Barcodes were removed and all sequences trimmed to be no greater than 148 bases long. For the purposes of this analysis paired-end reads were treated as separate loci, read 2 sequences being appended to read 1 sequence files. These sequences were assigned to RAD loci and genotypes using the 'denovo_map.pl' component of Stacks.
The key parameter values employed in identifying RAD loci were; a minimum stack depth of 10 (m=10), a maximum of 2 mismatches allowed in a locus (M=2) in an individual and up to 1 mismatch between loci when building the catalog (n=1). Finally the 'populations' component of Stacks was used to export filtered data (polymorphic loci containing 1-3 SNPs and present in at least 70% of samples for each population) in PLINK file format (PED and MAP files).

Population structure analysis methods
We ran assignment tests in STRUCTURE (Pritchard et al. 2000) for population cluster values of K=1-9. We performed five independent runs for each K with 100,000 burn-in iterations followed by 500,000 MCMC steps. We assumed the admixture model and did not include any prior information on populations. The number of distinct clusters was determined using STRUCTURE HARVESTER (Earl & von Holdt 2012) based on the conservative Evanno's method (Evanno et al. 2005). We re-ran the analysis for each identified population cluster to look for sub-structuring.

Outlier scan methods
Bayescan (Foll & Gaggiotti 2008) was run with 1,000,000 iteration, 50,000 burn-in and 20 pilot runs. Results were visualised in R using the script provided with the Bayescan download package with no modifications, setting false discovery rate (FDR) to 0.05.
LOSITAN (Antao et al. 2008) was run with 1,000,000 simulations under the Infinite Alleles mutation model, using the 'Neutral' mean Fst and forcing mean Fst options, and setting FDR to 0.05 and confidence intervals to 0.99.