Mapping evolutionary process: a multi-taxa approach to conservation prioritization

Human-induced land use changes are causing extensive habitat fragmentation. As a result, many species are not able to shift their ranges in response to climate change and will likely need to adapt in situ to changing climate conditions. Consequently, a prudent strategy to maintain the ability of populations to adapt is to focus conservation efforts on areas where levels of intraspecific variation are high. By doing so, the potential for an evolutionary response to environmental change is maximized. Here, we use modeling approaches in conjunction with environmental variables to model species distributions and patterns of genetic and morphological variation in seven Ecuadorian amphibian, bird, and mammal species. We then used reserve selection software to prioritize areas for conservation based on intraspecific variation or species-level diversity. Reserves selected using species richness and complementarity showed little overlap with those based on genetic and morphological variation. Priority areas for intraspecific variation were mainly located along the slopes of the Andes and were largely concordant among species, but were not well represented in existing reserves. Our results imply that in order to maximize representation of intraspecific variation in reserves, genetic and morphological variation should be included in conservation prioritization.


Introduction
Human-induced land use changes are leading to rapid habitat fragmentation. As a result, local populations are increasingly unable to track ongoing climate change (Travis 2003) and will have to adapt to the new environmental conditions or go extinct. The ability to adapt genetically depends on the amount of heritable variation currently present in a population. Thus, one way to mitigate the potential impacts of climate change on biodiversity is to ensure that the evolutionary processes underlying the generation and maintenance of such variation are preserved (Smith et al. 1993). Evolutionary processes have long been regarded of major importance in reserve design (e.g. Frankel 1974;Moritz 2002;Gregory et al. 2006;Smith and Grether 2008), but in practice have seen limited use in conservation planning. Nevertheless, given the current rate of anthropogenic climate change, there is an urgent need for approaches that enable conservation biologists to take into account evolutionary processes in prioritization schemes.
Recently, we proposed a novel, integrative approach to incorporate evolutionary processes in conservation prioritization based on spatially explicit models of intraspecific genetic and phenotypic variation, in combination with species-level data (Thomassen et al. 2010). The approach is based on the assumption that retaining the highest possible levels of adaptive variation will be a useful strategy to maximize a species' adaptive potential (e.g. Vandergast et al. 2008;Tymchuk et al. 2010;Sgrò et al. 2011). In contrast to species richness, adaptive intraspecific variation is likely to be the result of relatively recent microevolutionary processes and represents options for an evolutionary response to current and future environmental conditions. To accomplish its inclusion in conservation prioritization, ecological modeling approaches used to map intraspecific variation should be integrated with existing conservation strategies. To this end, we suggested a framework that includes the modeling of species distributions and corresponding environmentally associated intraspecific variation using population-level genetic and phenotypic variation and the subsequent integration of the resulting models with species-level information and socio-economic data to finalize prioritization (Thomassen et al. 2010). As a result of recent advances in spatially explicit modeling methods and the availability of climate and high-resolution remotely sensed environmental information, such an integrative approach is now feasible. Climate and remote sensing data capture many habitat properties that are important to characterize a species' ecological niche and to identify selective processes that result in adaptive variation. Because environmental variables often cover most of the global land surface, it is possible to predict environmentally associated intraspecific genetic and phenotypic variation to areas that have not been previously sampled.
Here, we extend our model to include multiple amphibian, bird, and mammal species. Our aim is to assess the utility of environmentally associated intraspecific variation in conservation prioritization. We focus on the association between intraspecific variation and environmental heterogeneity, because it is reasonable to suspect that this association is the combined result of direct selective forces on linked genes and of decreased gene flow because of reduced dispersal ability caused by habitat differences. Thus, intraspecific variation insofar explained by spatial heterogeneity in environmental variables may represent adaptive variation, even if the markers examined are evolving neutrally. Moreover, the morphological characters studied are fitness-related in many species. We investigate two important issues. First, if historical evolutionary processes that have resulted in speciation remain important today, traditional methods of reserve design based on species-level data may also capture relevant intraspecific variation. However, if this is not the case, valuable intraspecific variation may not be conserved with traditional methods. Second, because it is not feasible to gather data for all species present in a region under consideration for conservation, reserve design must rely on a limited subset of species that may or may not be representative of the entire community. Our objectives in this paper are to (i) model the distributions of a set of target species; (ii) model genetic and morphological variation in those species; (iii) design reserves based on patterns of intraspecific variation; (iv) design reserves based on species-level data; (v) compare the overlap between priority areas based on intraspecific variation and species-level data; and (vi) assess whether priority areas based on intraspecific variation are consistent among species. Our goal is not to produce a reserve plan, which would require the inclusion of socioeconomic, cost, and threat criteria, as well as the input of stakeholders. Rather, we aim to examine the utility of our proposed strategy from a purely biological perspective, focusing on maximizing different aspects of biodiversity. Our framework is also largely independent of the choice of modeling and reserve selection software, and the methods that we use here might readily be replaced by others developed in the future.
We conducted our research in continental Ecuador, a region that has received considerable conservation focus, because of its high species richness, endemism, and level of threat (e.g. Orme et al. 2005;Finer et al. 2010). Ecuador is geographically highly diverse, harboring steep environmental gradients along the slopes of the Andes, with the Pacific coastal region on the west, and the Amazon basin on the east. The topographic and climatologic heterogeneity of the country has likely resulted in diversifying selection, and divergence through isolation, resulting in high biodiversity (e.g. Graham et al. 2004;Chaves et al. 2007;Milá et al. 2009). Two previously defined biodiversity hotspots occur in Ecuador: the Chocó /Darien/West Ecuador and the Tropical Andes ecoregions (Myers et al. 2000). Lowland, eastern Ecuador is part of the Amazon basin ecoregion. Ecuador is highly suitable for our study, because our results can be compared to previous prioritization plans and existing protected areas. Further, its variety of environmental gradients and opportunities for isolation make it an ideal test of our approach for complex regions harboring much biodiversity.

Materials and methods
An overview of the workflow used in this study is provided in Fig. 1. First, geographic distribution models for each of seven target species were created using Maxent (Phillips et al. 2006) in conjunction with species-presence point localities and environmental variables. Second, using the modeled species distributions to define the spatial extent, we predicted genetic and morphological intraspecific variation in a generalized dissimilarity modeling framework (GDM; Ferrier et al. 2007). GDM modeled beta-diversity by evaluating how genetic and morphological variation vary as a function of geographic distance and environmental heterogeneity. For each trait, these analyses resulted in a GIS layer with values representing 50 classes of similarity in the respective traits, each class representing an equal amount of trait variation. These classes were separated as individual layers in ArcGIS (version 9.3) and subsequently used in reserve selection software (ResNet; Sarkar et al. 2009) that aimed at representing each class at a set target to prioritize areas for conservation based on genetic and morphological variation. In addition, we identified areas harboring the highest levels of intraspecific variation, thus taking into account both alpha-and beta-diversity. We also prioritized areas for conservation based on species richness using species distributions of birds, mammals, and amphibians, available from public databases.

Environmental data
We used a set of moderately high-resolution (1 km) climate and satellite remote sensing variables to characterize the often sharp habitat transitions in Ecuador (Data S1; Table S1). A detailed description of these variables can be found elsewhere (Buermann et al. 2008; see the Supporting Information in Thomassen et al. 2010 for maps of the variables used). Briefly, we used bioclimatic metrics (WorldClim; Hijmans et al. 2005), which are derived from meteorological station data , and capture variations in the annual mean, extreme, and seasonality of temperature and precipitation. To characterize ecosystem function and structure, we used optical passive and microwave active satellite data and derived products. The optical data were derived from MODIS (Justice et al. 1998) and include Leaf Area Index (LAI; Myneni et al. 2002) and percent tree cover (Hansen et al. 2002). Microwave data include QuikScat (QSCAT; Long et al. 2001) and Shuttle Radar Topography Mission (SRTM) elevation data. All raw satellite data were collected during 2001, and information on corresponding ecological inference can be found in Table S1.

Species-level data
Species distributions of birds (949 species; Ridgely et al. 2007) and mammals (336 species;Patterson et al. 2007) (Data S1) were obtained from the NatureServe database. Amphibian distribution maps were obtained from the Global Amphibian Assessment, accessed through the IUCN Red List database (467 species; IUCN, Conservation International, and NatureServe, 2008). These data were compiled based on field surveys and expert opinion and have seen widespread use in conservation planning (e.g. Hess et al. 2006;Darst et al. 2009 Maps of pattern and process Figure 1 Schematic of our workflow for this paper. (i) Species distribution models for the seven target taxa were created using Maxent (Phillips et al. 2006), species-presence point localities, and a set of environmental variables. The modeled species distributions served to delimit our study area for modeling intraspecific variation. (ii) We predicted the patterns of genetic and morphological diversity using in situ data and the environmental variables in a generalized dissimilarity modeling framework (Ferrier et al. 2007). For each genetic or morphological trait, these analyses resulted in a GIS layer with values representing 50 similarity classes. These classes were separated as individual layers in ArcGIS (version 9.3) and used in (iii) subsequent ResNet reserve selection software (Sarkar et al. 2009) to prioritize areas for conservation based on genetic and morphological variation. (iv) We also used species distributions of birds, mammals, and amphibians, available from public databases, to prioritize areas for conservation based on species richness and complementarity. Finally, we combined the results from steps (iii) and (iv) and compared all results with each other and with currently protected areas.
pelagic species were not included in our analyses. Downloaded data consisted of GIS feature layers for each species, which were clipped to the extent of continental Ecuador, resulting in a total of 1752 species for further analyses, subdivided as follows: 467 amphibians, 336 mammals, and 949 birds. Feature layers were subsequently converted to raster grid files in ASCII format at a 2-km spatial resolution, which was found to result in files that were small enough for computational handling in the ResNet reserve selection software. Data layers were processed using ArcGIS 9.3 (ESRI, Redlands, CA, USA).

Genetic and morphological data
To map intraspecific variation, we focused on seven species: three bird species: wedge-billed woodcreeper (Glyphorynchus spirurus), masked flowerpiercer (Diglossa cyanea), and streak-necked flycatcher (Mionectes striaticollis); three mammal species: silky, chestnut, and Seba's short-tailed bats (Carollia brevicauda, C. castanea, and C. perspicillata, respectively), and one amphibian species [the frog Pristimantis (formerly Eleutherodactylus) w-nigrum]. These species were chosen, because they are abundant, easily sampled, and represent a range of different vagilities and niches. They have different habitat preferences with respect to elevation and vegetation cover, occupy a range of vegetation strata, and show a range of food habits, including insectivores, frugivores, and nectarinivores. Genetic data were available for all species, except for two of the three bat species (C. brevicauda and C. castanea). Amplified fragment length polymorphism (AFLP) loci were used for the wedge-billed woodcreeper (136 loci, 178 individuals, 15 sites) and for the bat C. perspicillata (311 loci, 83 individuals, 9 sites). For the flowerpiercer and streak-necked flycatcher, we selected 10 microsatellite loci each and genotyped them for 102 individuals among 12 sampling locations and 106 individuals from nine sampling locations, respectively. For the frog P. w-nigrum, we sequenced two anonymous loci and an 840-bp region of the recombination activating gene 1 (RAG1) (n = 76-127 per locus, nine sampling sites). We calculated appropriate genetic population pairwise distances (Nei's D, Fst, and Fst) for subsequent use in GDM. Whereas the genetic markers used are different across taxa and may evolve at different rates, it seems justified to assume that the variation at this intraspecific level is the result of relatively recent evolutionary processes, which are the ones most likely to be relevant given future environmental change.
Morphological variation among populations was assessed through the measurement of fitness-related traits. For birds, these include wing, tail, tarsus, and bill lengths (the latter is defined as exposed culmen or the length of the bill from tip to the onset where the feathers start); bill width; and bill depth (wedge-billed woodcreeper n = 195, 15 sites; flowerpiercer n = 90, 10 sites; flycatcher n = 122, eight sites). For the bats, the following measurements were taken or calculated from landmarks of skulls using geometric morphometrics or inter-landmark measurements from the body on a total of 413 individuals (C. brevicauda: n = 167, 43 sites; C. castanea: n = 86, 25 sites; C. perspicillata: n = 160, 44 sites): centroid size for the size of the skull (hereafter referred to as skull size), angle of curvature of the zygomatic arch (hereafter referred to as zygomatic arch), and forearm length (Jarrín-V. et al. 2010). For P. w-nigrum, we measured the following morphological variables in 224 adult males among 16 sites: snout-vent length (SVL), gape width and the lengths of metacarpal phalanges (from the distal margin of the carpals to the tip of the third finger), radioulna, metatarsal phalanges (from the distal margin of the tarsal to the tip of the fourth toe), tarsal, tibio-fibula, femur (from the distal margin of the femur to the anterior margin of the urostyle), and lower jaw (from the posterior margin of the lower jaw to the tip of the snout). We calculated pairwise Euclidean distances from site averages as a measure of dissimilarity between sampling sites and divided these by the sum of the standard deviations within each site to take into account within-population variation. Further details on the genetic and morphological data used can be found in the Data S1.

Modeling of species distributions and intraspecific variation
Modeling morphological and genetic variation of the target species across Ecuador first required the delineation of each species' geographic range. Maps of continuous habitat suitability for each species were generated using Maxent 3.0 (Phillips et al. 2006), using known occurrences from throughout their ranges (Table 1), with particular focus on Ecuador. Rather than relying on both presence and absence data points, Maxent uses presence-only data. The input data consist of a set of environmental layers for the study region and the observed species-presence localities within that region. The program then uses these data to estimate the environmental niche space that most accurately describes the observed occurrences. To test for the importance of environmental predictor variables in determining the species' range, we used the variable jackknifing approach implemented in Maxent. Using a subset of the original variable set, new models were computed and their performance as measured by the regularized training gain (the average log probability of the presence samples, corrected for a uniform distribution with gain = 0) compared to that of the Conservation of pattern and process Thomassen et al. full model (Phillips et al. 2006). We tested whether models were significantly better than random predictions using a one-tailed binomial test on the proportion of test sites falling outside the prediction resulting from a model that used 60% of the data for training and 40% for testing (Data S1) (Anderson et al. 2002). In addition, we used the area under the receiver operator curve (AUC) to evaluate model performance (Phillips et al. 2006). To convert the continuous Maxent predictions into presence-absence maps defining the study area for subsequent GDM analyses (see below), we used thresholds for habitat suitability, which were centered within the range of a number of optimized Maxent thresholds, including 'Equal training sensitivity and specificity' and 'balance threshold'. Criteria for choosing this threshold level included verification whether the Ecuadorian field sampling sites were included in the suitable area and comparisons to published range estimates in field guides.
To predict the distribution of environmentally associated genetic and phenotypic variation across the landscape, we used GDM (Ferrier et al. 2007), a matrix regression technique that predicts biotic dissimilarity (i.e. beta-diversity) between sites based upon environmental dissimilarity and geographic distance. GDM explicitly considers the influence of geographic distance on explaining biological variation. This is carried out in a two-step method: first, dissimilarities of a set of environmental predictor variables are fitted to the pairwise genetic (e.g. Fst, Nei's D, or Ust) or phenotypic dissimilarities (Euclidean distance), the response variables. The contributions of predictor variables to explaining the observed response variation are tested by permutation, and only those variables that are significant are retained in the final model.
The relative importance of predictor variables in a GDM can be assessed by means of response curves. These procedures result in a function that describes the relationship between environmental and response variables. Second, using the function resulting from the first step, a spatial prediction is made of the response variable patterns. The output is a map of classes, in our case 50, that summarize the variation in different phenotypes or genotypes. These classes were used in subsequent reserve design algorithms as a way to identify areas of high conservation priority.
Geographic distance may not represent the distance an individual might travel among localities, because of variation in habitat suitability, or the presence of barriers, such as wide rivers or high mountain ranges. We therefore also included least-cost-paths (LCP; PathMatrix 1.1; Ray 2005) and resistance distances (RD; Circuitscape 2.2, McRae 2006) in our models (Data S1), which take into account the permeability of the habitat matrix. LCP and RD were computed using friction surfaces in which unsuitable habitat (defined as those areas predicted to be unsuitable by a species distribution model, after applying a threshold value for suitability) was assumed to be 10 times as difficult to penetrate as suitable habitat.
To further evaluate the extent to which distance is potentially correlated with environmental differences, for each region and for each dependent variable we ran independent tests with the following sets of predictor variables: (i) environmental variables and distance (geographic, LCP, or RD); (ii) only distance (geographic, LCP, or RD); (iii) only environmental variables. Comparison of the results from these three runs provided an indication of the correlation between geographic distance and environmental differences. In order to assess the significance of the level of variation that was explained by our models, we ran additional models in which the environmental layers were substituted by layers with random values for each grid cell. The resulting percentage of variation explained was compared to that of the full model. We considered the performance of the full model not significant if it explained an equal amount or less of the total variation than a model with random environmental variables.
To identify areas harboring the highest amounts of genetic and phenotypic alpha-diversity, where reserves could most efficiently protect a large percentage of the intraspecific variation, we calculated the standard deviation from GDM predictions in an area of 3 · 3 grid cells at a spatial resolution of 1 km. Thus, we first identified variation that is related to environmental variables, and which is likely to be representative of adaptive variation, and subsequently mapped areas showing high levels of such genetic and morphological variation within a 3 · 3 km area.

Reserve selection
We designed conservation areas to cover two types of measures for overall biodiversity: (i) species occurrences and (ii) intraspecific genetic and morphological variation. The same area prioritization procedure was used for both data types: a rarity-complementarity algorithm implemented in the ResNet software package (Sarkar et al. 2009). The design of a set of areas to protect features of biodiversity is typically formulated as a constrained optimization problem in which the objective is to establish protected areas that meet the representation targets for the features while taking up as little land as possible . A representation target of 10% of the total habitat suitable for a species has seen widespread use in conservation planning exercises and is included in the Convention on Biological Diversity (Justus et al. 2008). We examined a range of conservation targets for the species and genetic and phenotypic classes, specified as a single occurrence (one site of 2 · 2 km), and 5%, 10%, and 20% of the total suitable habitat.
The input data for ResNet consist of a set of sites and a list of the elements of biodiversity in each site. Using this information, ResNet employs a heuristic algorithm to prioritize areas for conservation. The iterative algorithm first selects the site that contains the rarest species or class of genetic or phenotypic similarity. ResNet breaks ties in rarity by selecting the site with the greatest complementarity, which is the one containing the largest number of elements of biodiversity that are not currently represented adequately to meet their targets. ResNet stops selecting sites when all elements of biodiversity have met their targets in the selected areas. Because of the large number of grid cells in combination with the amount of species or classes of genetic and morphological similarity, we used a heuristic algorithm implemented in ResNet, which has been shown to find near-optimal solutions to large conservation planning problems (Sarkar et al. 2009). We utilized ResNet because the large number of grid cells in combination with the amount of species or classes of genetic and morphological similarity was tractable only with a heuristic algorithm. Further, the rarity-complementarity algorithm is more effective for the rapid analysis of large-scale data sets than other techniques that have been used to select reserves such as simulated annealing or integer programming (Kelley et al. 2002;Sarkar et al. 2004).
To avoid putting specific weight on either species-level data or intraspecific variation, we chose not to combine these two types of data into a single ResNet run, but instead made a composite map of areas harboring high levels of intraspecific alpha-diversity, reserves based on species-level data and those based on classes of similar morphology and genotype identified using GDM. To create this composite map, we started by mapping areas comprising the highest levels of genetic and morphological alpha-diversity for each species (see above). This was carried out in ArcGIS using the following three steps: (i) for each species we calculated the highest 10% of variation in each trait (alpha-diversity) and compiled all traits in a single map; (ii) the resulting GIS layer was reclassified such that areas comprising the highest 10% of variation were coded 1, and all other areas 0; (iii) the layers from the previous step were added across species, resulting in a layer with values 0-7, with increasing values indicating more cross-species overlap in high levels of variation. On top of this map, we then superimposed reserves selected by ResNet using intraspecific variation or species-level data, because the algorithm capitalizes on unique variation that may not be captured in areas harboring high levels of intraspecific alpha-diversity. The composite map is thus likely to be the most comprehensive representation of variation at all levels.

Species distribution models
Threshold-dependent one-tailed binomial tests on the Maxent extrinsic omission rate and proportional predicted area of the species distribution models (Fig. S1) were highly significant (Table 1), suggesting that the models performed significantly better than random. Evaluation of the AUC values also suggested high model performance (AUC values ‡0.895 for all models; test AUC values ‡0.744; Table 1). The variable importance for each Conservation of pattern and process Thomassen et al. of the models is shown in Table 1. The most important variable in defining distribution models for the streaknecked flycatcher and flowerpiercer was elevation, whereas precipitation variables were important for defining the ranges of the wedge-billed woodcreeper (see also Buermann et al. 2008) and the three bat species (Table 1). Temperature variables were particularly informative in defining the range of the frog species. The wedge-billed woodcreeper and the three bat species were predicted to be present across much of the Ecuadorian lowland areas west and east of the Andes (Fig. S1). The masked flowerpiercer was restricted to mid to high elevations. Likewise, the streak-necked flycatcher was predicted to occur from mid to high elevations, but with the highest probability in the mid-elevations, whereas the frog species occurred along a fairly narrow band at mid-elevations. The models corresponded well with published range maps (Lynch 1979;Eisenberg and Redford 2000;Ridgely and Greenfield 2001) and assessments based upon personal observations in the field (TBS, BM, PJ, CMK, CJS).

Models of intraspecific variation
Generalized dissimilarity models varied in their performance across species and across genetic and morphological traits. A total of 36 out of 60 models were found to perform better than random (Table 2). For the majority of these 36 models, environmental variables far outperformed measures of geographic distance as predictors of intraspecific variation ( Table 2), suggesting that isolationby-distance was not an important driver of diversification in these traits. However, for some of the models, geographic distance explained a considerable proportion of the total intraspecific variation, which suggests that this part of the observed variation is potentially attributable to isolation-by-distance. In these cases, environmental heterogeneity and geographic distance are partially correlated, and it is not possible to distinguish between the effects of the two on patterns of variation. In the Data S1, we describe the models for each species, and predictive maps of the models are shown in Fig. S2 (for the wedge-billed woodcreeper, see Thomassen et al. 2010). In summary, differentiation in the flowerpiercer is driven by moisture levels and elevation differences and located along elevational gradients on both sides of the Andes. In the flycatcher, vegetation and temperature variables were most important in explaining the observed genetic and morphological variation. Differences were most apparent between the eastern and western sides of the Andes, as well as along the western elevational gradient for microsatellites and tarsus length, and the eastern elevational gradient for tail length. In the bats, variation in the zygomatic arch was mostly related to seasonal measures of precipitation and vegetation characteristics. Forearm length was best explained by vegetation variables (including QSCAT, which distinguishes different forest types), which is potentially interpreted as adaptive given the importance of this trait in maneuverability in relation to vegetation density. Skull size was related to climate and vegetation variables, but also to geographic distance. Differentiation on the east side of the Andes was most pronounced between the lowlands and mid-elevations, but for some variables also along a latitudinal gradient (Fig. S2). Similarly, the pattern of variation in the bats on the west side of the Andes was most apparent along a latitudinal axis, with additional differentiation between the lowlands and mid-elevations. Finally, temperature and vegetation variables were most important in explaining the observed genetic and morphological variation in frogs. The resulting patterns of variation are located along the elevational gradients on both sides of the Andes, with additional variation located along a latitudinal gradient (Fig. S2).

Levels of intraspecific variation
To identify areas harboring the highest levels of intraspecific variation per unit area (alpha-diversity), we mapped the top 10% of standard deviations for GDM classes in a 3 · 3 km area (Fig. 2). The highest levels of alpha-diversity in morphological and genetic traits were found mostly along the slopes of the Andes (Fig. 2). These results are in contrast with the location of areas comprising the highest levels of species richness, which are in the Amazon basin (Fig. S3). Thus, the spatial patterns of species-level data and intraspecific variation are different, suggesting that the former may not be a suitable representation of the latter, and that, if the goal is to maximize intraspecific variation, it should be considered in conservation prioritization. Even though the wedge-billed woodcreeper and the bats C. brevicauda, C. castanea, and C. perspicillata also showed higher levels of genetic and morphological alpha-diversity within the Amazon basin, as well as along the Andes, the lowland areas capture only a fraction of the intraspecific variation. Interestingly, areas of high intraspecific alpha-diversity are at least partly concordant across the species examined here, suggesting that a subset of species might accurately represent a broad array of taxa.

Reserve selection
To generate maps of conservation priority areas based on species richness and complementarity of amphibians, birds, mammals, all species, and based on genetic and morphological variation, we performed ResNet runs at representation targets of one occurrence and 5%, 10%, and 20% of the total area occupied by each species or each GDM class of similar genotypes or phenotypes. The percentage of total land area selected for conservation increased linearly with higher targets to a maximum of 22.7% at a 20% target for birds (Table 3). The area selected for conservation based on genetic and morphological variation was up to 5.2% smaller than that based on biodiversity pattern, represented by species-level data. The overlap between priority areas for each taxonomic group was relatively low, ranging from 10.4% to 12.8% at a 10% representation target (Table 4). Interestingly, priority areas at the 10% representation target based on all species captured only 17.8-32.8% of those based on each individual taxonomic group. These results suggest that none of the three taxonomic groups (amphibians, birds, and mammals) is a good surrogate for any of the others and that only their combined analysis can ensure good representation of all species in priority areas. In addition, the overlap between genetic/morphological variation and Conservation of pattern and process Thomassen et al. species-level data at a 10% representation target was also low (7.6-12.4%; Table 4), suggesting that reserves based on species-level data do not effectively represent genetic and morphological variation, or vice versa. The existing reserve network covers about 16% of continental Ecuador and is overlapping with a higher percentage of selected areas based on genetic and morphological variation than based on species-level information (Table 3). Overall, the percentage of priority areas that were located in existing reserves ranged between 11.1% and 23.8%. The composite map of reserves using algorithm-based selection of species richness and intraspecific variation, and regions harboring the highest levels of intraspecific alpha-diversity, suggests that high priority should be given to regions located along the western and eastern slopes of the Andes (Fig. 3). Although the areas of highest intraspecific alpha-diversity are in the Andean foothills (Figs 2 and 3), when we selected areas to represent intraspecific variation with ResNet (Fig. 3), some sites in the Amazon lowlands were also selected. We hypothesized that the Amazon lowlands were selected because there was high intraspecific alpha-diversity in bat species in these lowland areas (Fig. 2E-G), as well as unique variation not present elsewhere. We tested these hypotheses by removing bats from the analysis and repeating the ResNet prioritization. This resulted in a more pronounced concentration of selected sites in the Andean foothills (Fig. S4). Nevertheless, areas in the Amazon basin continued to be selected, suggesting that unique genetic and morphological variation in wedge-billed woodcreepers from the Amazon basin contributed to its importance in the ResNet solutions. To determine which sites were regarded most important by the ResNet algorithm, we plotted the order of site selection of the map of Ecuador (Fig. S4), which suggested that sites in the Andes were selected first. This gave us confidence that the Andean foothills were indeed the most important in reserve selection.

Discussion
Climate change poses an urgent and significant threat to biodiversity (Fischlin et al. 2007;Sinervo et al. 2010). Species may not be able to respond to climate change unless they have the capacity to shift their ranges. However, range shifts may often not be possible for many reasons, including disparities in vagility between mutualistic species, isolation of populations because of habitat fragmentation, new biotic interactions resulting from a different species composition, or elevational constraints. Hence, many populations will have no choice but to adapt to the changing climate conditions in situ. A prudent conservation strategy therefore is to design reserves that maximize intraspecific variation. Our results provide the first, country-wide assessments of the utility of intraspecific genetic and morphological variation for conservation prioritization. Future work is planned to predict future morphological changes under climate change based on current spatial relationships between morphology and environment, but this is beyond the scope of this paper. In considering evolutionary processes in reserve design, conservation planners have relied on a number of different measures of biodiversity as surrogates for evolutionary process. First, the phylogenetic species definition has been used to identify evolutionarily significant units (ESU) that represent different (genetic) lineages for conservation prioritization (Moritz 1994;Peterson and Navarro-Sigüenza 1998). An advantage of this approach is that some level of sub-specific variation is likely to be preserved. However, it is probably biased toward broad-scale vicariant events and unable to capture the finer-scale adaptive variation we identified with our method. Waples et al. (2001) used another approach to define ESUs based on a combination of habitat characteristics, genetic data, and life history traits. Although this approach may capture adaptive variation in the study species, its use has been limited to a single species complex, where detailed studies of life history were available, which is unlikely to be the case for the vast majority of species. A second approach to preserve evolutionary processes is the protection of habitat variability, primarily along environmental gradients (e.g. Cowling and Pressey 2001;Smith et al. 2001Smith et al. , 2005Rouget et al. 2006;Klein et al. 2009). However, it is unclear to what extent total adaptive diversity is protected by focusing exclusively on gradients. Finally, phylogenetic diversity (PD) -as measured by the branch lengths in a phylogenetic tree -has served as a measure of genetic variation (e.g. Vane-Wright et al. 1991;Faith 1992;Forest et al. 2007). Such species-level information is useful in determining the phylogenetic diversity of communities, but may fail to capture extant population-level variation and adaptive diversity. It will likely miss very recent episodes of local adaptation and genetic isolation where sufficient time has not elapsed for species-level phylogenetic trees to show distinct genetic partitions. To help alleviate this issue, several phylogeographic studies have focused on intraspecific PD (Moritz and Faith 1998;Rissler et al. 2006). However, the genetic markers used in these studies typically evolve slowly, and the timescale of differentiation between lineages may be unsuitable for the purposes of preserving extant evolutionary processes. Moreover, patterns of genetic lineages may be the result of neutral rather than selective processes and consequently of less relevance to understanding adaptive variation that might be germane to climate change.
In contrast to the approaches discussed earlier, our spatially explicit models focus on the relationship between environmental heterogeneity and genetic and phenotypic variation and emphasize local adaptation and genetic turnover on recent timescales. The combined use of two types of data -genetic and morphological -is likely a benefit to our approach (Crandall et al. 2000). The genetic markers used may or may not be linked to adaptive variation, but the fact that a large proportion of the variation was explained by environmental variables, and not by geographic distance, suggests that the markers may be a suitable proxy for adaptive variation. In addition, although there is no direct evidence for the Figures are provided for reserves based on each of the three different taxonomic groups (amphibians, birds, mammals, and the three groups combined) for which species-level data were available, as well as for intraspecific variation (gen/morph). *The set of selected sites was required to represent at least a certain target percentage of the total area occupied by a species, or by a generalized dissimilarity modeling class of similar genotypes or phenotypes. We specified targets of one occurrence (one site of 2 · 2 km), 5%, 10%, and 20%.

Conservation of pattern and process
Thomassen et al. Figures are provided for overlap between reserves based on each of the three different taxonomic groups (amphibians, birds, mammals, and the three groups combined) for which species-level data were available and based on intraspecific variation (gen/morph). Results are shown for representation targets of 10% and 20% of the total area occupied by each species or by each class of similar genotypes or phenotypes. The table reads as: percentage of 'row' overlapping with 'column'.

Existing reserves ResNet solution based on pattern ResNet solution based on process
Overlap in number of species for highest levels of intra-specific variation 1 -3 4 -7 0 40 80 120 160 km Figure 3 Composite map of selected reserves at a 10% representation target using species-level data (green) and genetic and morphological data (blue) in ResNet and areas harboring high levels of intraspecific alpha-diversity (orange and red). Areas delimited with solid black lines indicate existing reserves. Areas indicated with dashed lines are those where one or more of the generalized dissimilarity models show high uncertainty, because the associated environmental conditions are outside the range of those encountered at sampled locations.
adaptive context of the morphological traits in the species studied, there is evidence in other species that these traits may be fitness-related (e.g. Grant et al. 1976;Wassersug and Sperry 1977;Arnold and Wade 1984;Smith 1990;Norberg 1994;Sedinger et al. 1995;Heinen and Hammond 1997;Bardwell et al. 2001;Benkman 2003;Jensen et al. 2004;Milá et al. 2008;Sol 2008;Nogueira et al. 2009). Thus, these traits are likely to represent adaptive variation. The combined use of morphological and genetic markers may provide more confidence in the resulting maps of priority areas than either alone.

Between-species comparisons
Because reserve selection is often scale-dependent and our focal species comprise distributions of varying sizes and coverage, we did not select reserves for each focal species to compare specific priority areas, but instead focused on the areas with the highest levels of intraspecific alpha-diversity to assess whether these priority areas are concordant among different species. Such areas are expected to be the ones where a given reserve could most likely harbor as much intraspecific variation as possible. All seven target species showed high levels of intraspecific alpha-diversity along the slopes of the Andes, and some in the Amazon basin, most notably the three bat species (Fig. 2). This result is perhaps not surprising given that the steep elevational gradient from lowlands to highlands is correlated with steep gradients in many other environmental variables. Interestingly, however, the highest levels of genetic and morphological alpha-diversity are mostly located in the mid-elevations from approximately 1500-2500 m, extending down to about 700 m in the southeast. In addition, areas harboring high levels of genetic variation do not in all species overlap with those comprising high morphological alpha-diversity. This may be explained by the fact that it is unlikely that the genetic markers used are linked to the morphological traits measured, and the two data types may thus represent different levels of diversity. In addition, variation in the phenotypic traits may in part be the result of a plastic response, without underlying genetic differences. Plasticity in itself is likely important to the versatility of populations to respond to climate change. Thus, the reduction of levels of plasticity because of habitat loss may have negative impacts on the persistence of populations and species under changing environmental conditions. As a consequence, considering phenotypic trait variation in conservation prioritization is important irrespective of whether the variation is because of underlying genetic differences or because of plasticity.

Comparison between reserves based on species richness versus intraspecific variation
Priority areas based on species-level data were only partly concordant with those based on intraspecific variation (Table 4) and may thus not be sufficiently effective in preserving either species richness or intraspecific variation given ongoing and future climate change. This result suggests that the reverse is also true: reserves based on intraspecific variation do not adequately capture species-level variation. A comparison of reserves based on intraspecific variation using either the ResNet rarity-complementarity algorithm or the criterion that reserves should harbor high levels of intraspecific alpha-diversity also reveals discordant patterns (Fig. 3). This can partly be explained by the distribution of unique variation not found elsewhere and by the way ResNet selects sites important for conservation. Sites comprising high levels of variation are selected first, because these represent the observed intraspecific variation very efficiently (Fig. S4). Thus, the target for the full range of genetic and morphological variation found in areas with high intraspecific alpha-diversity across species is quickly met in relatively small areas (the small, scattered blue areas along the Andes in Fig. 3). However, areas harboring unique variation, but little overall intraspecific diversity, have at that point not been selected yet, and the targets for the unique variation have not been met. In subsequent iterations of the algorithm, these are the areas selected next, but as overlap in unique intraspecific variation among species decreases, larger areas will need to be selected in order to meet the respective targets of representation. Thus, the reserves in areas harboring high levels of variation are smaller than those in areas of low variation, because of the higher efficiency of protecting intraspecific variation in the former. Taken together, these results suggest that the spatial patterns of species richness are not concordant with those of intraspecific variation and emphasize the need for both to be taken into account in reserve design. Our results also suggest that the evolutionary processes underlying observed patterns of species richness may occur on different spatial scales or have been differently distributed geographically than those causing current levels of intraspecific variation. For example, intraspecific variation may evolve in response to habitat heterogeneity on a smaller spatial scale than that necessary for speciation (Losos and Schluter 2000).

Comparison with Ecuador's existing reserves
The existing reserve network in Ecuador does not effectively capture species richness or intraspecific variation, despite the fact that a relatively large proportion of the Conservation of pattern and process Thomassen et al. total landmass is protected (16%). In addition, large areas of high intraspecific or species diversity remain unprotected by the existing reserves (Fig. 3). This discrepancy may not be surprising given the fact that reserves are often established using socio-economic and political criteria rather than biological ones (Pressey and Tully 1994;Bass et al. 2010). In fact, a similar ineffectiveness of the existing reserves was found in previous studies using ecosystem-level data across Ecuador (Sierra et al. 2002), species-level data for the equatorial Pacific region (Peralvo et al. 2007), and a combination of the representation of ecosystems and Red List species in the Tropical Andes and Chocó (Sarkar et al. 2009). Our study thus supports the conclusions of previous studies regarding the inefficacy of the existing reserve network in Ecuador, but is novel in showing that this is also true when using intraspecific variation as a measure of biodiversity.

Comparison to other proposed reserves
One of the most comprehensive conservation prioritization studies for the entire country of Ecuador was carried out using differences in habitat type as a measure of diversity (Sierra et al. 2002). Our results agree with Sierra et al. (2002) in placing high priority on the Andean slopes in regions that are largely unprotected in existing reserves. In addition, our study placed low to medium priority on some areas in the Amazon basin (blue areas in Fig. 3), which corresponds to the prioritization by Sierra et al. (2002). However, only our study identifies these areas as having unique intraspecific variation found nowhere else in Ecuador. Consequently, even though the Amazon basin may not harbor high levels of intraspecific alpha-diversity, a prioritization scheme that aims to protect the entire range of intraspecific variation in the species studied should also include areas in the Amazon basin. Moreover, the Amazon basin comprises the highest species richness per unit area found in Ecuador (Fig. S3).
In contrast to our study, the prioritization scheme by Sierra et al. (2002) put more emphasis on the far northwestern and central coasts, likely as a result of the presence of unique ecosystems, comprised of coastal lowland evergreen forests and dry shrub, found nowhere else in the country. Two additional studies that prioritized areas for conservation in Ecuador did so only for parts of Ecuador (Peralvo et al. 2007;Sarkar et al. 2009). Because of this difference in scale, and because of likely scale dependence in reserve selection (e.g. Sierra et al. 2002;Warman et al. 2004), we refrain from making comparisons with our study.

Advantages and limitations
A potential limitation of our study is that the range maps used as input for the species-level analyses are at a broad scale and do not capture finer-scale nuances of the species distributions, which may influence the results of prioritization software. In addition, it would be desirable to add intraspecific variation of additional species, most notable plant and insect species, which may improve the representation of intraspecific variation of the entire community. Nevertheless, we are confident in concluding that reserves based on species-level data are not concordant with those based on intraspecific variation and that there is a need for both in conservation prioritization.
Over-or underprediction of the species' ranges may have consequences for the models of intraspecific variation. Over-and underprediction is likely to occur at the extremes of the suitable range of environmental conditions. Thus, prediction of a species' range into areas where it is not present may result in the false prediction of unique intraspecific variation. In contrast, failure to predict the full species range may result in neglecting unique intraspecific variation that might be important to conserve. Nevertheless, based on expert knowledge and our own experience in the field with the seven target species, we are confident about the predicted species distribution. While our target species are unequally distributed across Ecuador, the combined distributions of the target species span almost the entire country. In addition, it is important to note that the patterns of high intraspecific alpha-diversity are concordant among target species (see above and Figs 2 and 3).
Several questions remain with regard to the implementation of the approach proposed here. For instance, the intraspecific genetic and morphological variation considered may not be equally important to fitness. Thus, when such data are available, weighting of the traits may improve the relevance of the prioritization scheme to maximizing the evolutionary potential of populations in the face of climate change. In addition, the consequence of adding or removing traits is unknown. In general, adding species or traits will likely result in down weighting the relative influence of individual traits, which may be beneficial to the utility of the prioritization scheme for as many species as possible. We demonstrated that removing the three bat species affected the resulting prioritization scheme, but not in a way that it significantly altered the locations of the top priority reserves.
Finally, the methods presented here appear to be generally applicable to other areas insofar as there is reasonable overlap between the ranges of species targeted for modeling intraspecific variation. Conservation prioritization is scale-dependent, and our approach is no exception.
Careful consideration of the size of the study area is therefore warranted. The framework presented here could be tailored to any region subject to data availability. The data required to use this framework include the longitude and latitude of species' occurrences along with data on intraspecific phenotypic variation and population structure inferred from genetic markers. The software and GIS layers utilized in this study are freely available on the internet and could be deployed to develop conservation plans for anywhere in the world. As for Ecuador, it should be emphasized that we did not include socioeconomic criteria, levels of threat, and opportunities, all critical elements in reserve design ( Fig. 1) (reviewed in Moffett et al. 2006). Future work will attempt to incorporate this information in order to develop a final prioritization strategy for Ecuador.

Conclusions
We have shown that the locations of the highest levels of intraspecific genetic and morphological variation among a diverse group of vertebrate taxa are broadly concordant. This suggests that a subset of taxa could potentially serve as suitable representatives for community-wide levels of intraspecific variation. We also demonstrate that priority areas for conservation-based species-level data are not concordant with those based on intraspecific genetic and morphological variation. Because maximizing the level of environmentally associated genetic and morphological variation is likely to benefit populations in the face of climate change, our results emphasize the need for incorporating intraspecific variation in conservation prioritization. Finally, we have developed a new framework for considering multiple types of diversity for conservation planning. While we used Ecuador as a test case for these methods, they are broadly applicable to other regions and species and should serve to improve conservation planning targeted at maintaining all levels of diversity.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Data S1. Supporting materials and methods. Figure S1. Predicted species distributions based on Maxent (Phillips et al. 2006) models for the following species: (A) wedge-billed woodcreeper G. spirurus; (B) streak-necked flycatcher M. striaticollis; (C) masked flowerpiercer D. cyanea; (D) the frog species P. w-nigrum; (E) silky short-tailed bat C. brevicauda; (F) chestnut short-tailed bat C. castanea; and (G) Seba's short-tailed bat C. perspicillata. Colors indicate the predicted logistic probability of occurrence (see color bar). Figure S2. Maps of the generalized dissimilarity models that performed better than random.

Conservation of pattern and process
Thomassen et al. Figure S3. Total species richness maps for birds, mammals, and amphibians in Ecuador. Figure S4. ResNet solutions for intraspecific variation at 10% representation targets. Table S1. List of environmental variables used in our analyses. Table S2. Elevation and measures of genetic diversity per site for M. striaticollis microsatellite data. Table S3. Elevation and measures of genetic diversity per site for D. cyanea microsatellite data. Table S4. Elevation and measures of genetic diversity per site for C. perspicillata AFLP data. Table S5. Elevation and measures of genetic diversity per site for P. w-nigrum nuclear DNA data.
Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.