Frontiers of Biogeography Quantifying amphibian range fragmentation in the southeastern

. An often overlooked component of research on factors that drive amphibian geographic distributions is description of species range shape. Broad-scale range disjunction has implications for phylogeography, ecology, and conservation, but descriptions of fragmentation are usually based on subjective visual assessment of range maps. Here, we describe a method for objectively quantifying range fragmentation and use this method to describe the patterns of amphibian species range shapes in the southeastern United States, home to the highest amphibian species richness in North America. Species ranges varied widely in degree of fragmentation, from completely contiguous to highly fragmented, and degree of isolation of range fragments. Incorporating ecological niche models added information about finer-scale fragmentation. We also demonstrate that this method can add objectivity to studies that use ecological niche modeling to assess change in range fragmentation through time, enhancing research in conservation and biogeography.


Introduction
Geographic ranges of species are fundamental to the study of biogeography. The size and shape of species ranges are related to the biotic and abiotic factors and historical processes that influence the distribution and abundance of organisms (Brown et al. 1996). Climate change and human-mediated habitat loss or alteration often leads to fragmentation, which, in addition to loss of habitat, also involves an increase in the number of disjunct patches in the species range, a decrease in the size of the patches, and an increase in isolation of patches (Fahrig 2003). Depending on the dispersal capability of a species, habitat fragmentation can reduce or eliminate dispersal between populations, as well as the ability to migrate in response to environmental change, increasing the risk of extinction (Lande 1988, Cushman 2006. Quantifying the extent of range disjunction for species present in a given area would aid in identifying broad-scale historical biogeographic patterns and selecting species for comparative phylogeographic study. Similarly, phylogeographic studies that use ecological niche modeling and paleodistribution modeling to explore species distribution shifts through time generally rely on qualitative visual assessment of model differences and would benefit from objective measurement of distribution shape (e.g., Newman and Austin 2015).
Species range maps are all some form of extrapolation from point occurrences, as the most basic description of a distribution is simply a collection of points in space at a particular time. However, broader descriptions of the area where a species occurs are often more useful. To this end, many attempts have been made to define and measure the extent of species ranges (Reaka 1980, Schoener 1987, Spitzer and Lepš 1988, Ford 1990, Gaston 1994. Gaston (1991) describes two ways of defining a species' geographic range: extent of occurrence and area of occupancy. Extent of occurrence in general broadly encompasses the entire area where a species is found -essentially the minimum convex polygon that includes all known occurrences. Areas of occupancy exclude regions within the wider extent of occurrence where the species is not found, such as regions with unsuitable habitat, and depends on the scale at which it is measured. Range maps in field guides usually depict extent of occurrence, sometimes with varying degrees of area of occupancy taken into account, depending on scale.
Missing from the biogeography literature is an explicit method for quantifying the degree of disjunction of a species range defined by any means. Here, we propose a method for assessing the extent of fragmentation

Quantifying species range disjunction
We propose a method for objectively assessing the degree of fragmentation of any species range. Stevens and Enquist (1996) define a range fragment as "a cohesive cluster of sightings that is represented as a continuous blob on a distribution map." We retain this definition while emphasizing their caveat that bias is unavoidable in map creation -bias in observation as well as determination of cohesiveness. Degree of range disjunction and population isolation is a direct matter of scale (Erickson 1945), and small ranges tend to be mapped with greater detail and magnification than larger ranges (Brown et al. 1996). The method we propose can thus either be used at any scale as an initial examination of differences in range shape across species, or at multiple scales for a single species to explore the effects of scale, habitat, etc., on descriptions of range fragmentation.
Species range shape is analyzed in Fragstats v.4.2.1.603 (McGarigal et al. 2012) using two metrics: the landscape DIVISION index (D) and Euclidean nearest neighbor distance (NN). These metrics are useful in particular because they can be applied to binary range maps (i.e., presence/absence or suitable/unsuitable), as well as to the land cover or habitat rasters more commonly used in ecological studies. DIVISION measures the probability that two points placed randomly on the landscape (= entire species range) will be on the same undissected patch. Fragstats considers a "patch" to be a contiguous group of raster cells connected horizontally, vertically, or diagonally. Mathematically, D is based on patch area and total landscape area, where the landscape is simply the combination of all patches and excluding all raster cells outside of the species range. D is calculated by Fragstats using the following formula: where a is the area (m 2 ) of patch ij, and A is the total landscape area (m 2 ). D ranges from 0 (a single contiguous patch, no fragmentation) to 1 (each patch is a single raster cell, highest fragmentation). Because D is uninformative about the degree of geographic isolation of patches in a species range, we also calculate the average distance between patches in each range. NN is the shortest straight-line distance between a focal patch and its nearest neighboring patch, measured from the centers of the closest two cells of the neighboring patches, averaged across all patches in a range.

County-based range maps
For county-based range maps formatted as ESRI polygon shapefiles: each shapefile is projected to the North American Albers Equal Area Conic coordinate system and then converted to a presence/absence binary raster grid with cell size equal to the shortest of the width or height of the shapefile, divided by 250. This is the default option when performing the same conversion in ArcGIS (ESRI). Rasterizing vector polygons of range maps can be problematic if an inappropriate grid cell resolution is used. If the resolution is too large, patches that were disjunct on the original (polygon) range map may inadvertently be joined on the raster. In our analyses (see Example below), we visually inspected all raster maps, and none showed problems from resolution. In addition, the use of a unique resolution for each species has no qualitative effect on results (see Example below). Because unnecessarily high raster resolution for species with large ranges produces extremely large file sizes that can be problematic in downstream analysis, raster resolution based on range shapefile dimensions is most efficient.
Another issue of rasterization is the presence of extraneous grid cells outside the range but near the edge; this is especially common along coastlines, where the range may be highly irregular in shape and the original map may include islands that are part of counties where the species is present. A grid cell that is not connected in any direction to other patches is considered by Fragstats to be a separate patch and can thus cause misleading results. Extraneous cells should therefore be removed from rasters. Connected cells include the four cells adjacent to the focal cell, as well as the four diagonal neighbor cells.

Ecological niche models
Ecological niche models (ENMs; also referred to as species distribution models) can provide a higher resolution map of a species range than the county scale. ENMs combine species occurrence data with environmental data (e.g., climate layers) to generate a probability surface of environmental suitability extrapolated from environmental conditions at locations where the species is known to occur. For a given species, an ENM is generated using climate (temperature and precipitation) and land cover layers. In our example detailed below, we downloaded natural history collection specimen occurrence records from the GBIF online database 1 . Occurrence records outside of the county-based range for that species were discarded to minimize potential error from misidentified specimens or incorrect georeferencing. We used 19 bioclimatic layers downloaded from Worldclim (Hijmans et al. 2005) at a resolution of 1 km 2 and clipped to an extent that encompassed all species ranges. A land cover layer was downloaded from the National Map 2 database Frontiers of Biogeography 2019, 11.1, e37772 © the authors, CC-BY 4.0 license 3 of the U.S. Geological Survey at a scale of 100 m and resampled to match the resolution and extent of the climate layers. An ENM for each species was generated in Maxent v.3.3.3k (Phillips et al. 2006), as this method has been shown to model realized species distributions better than county-based range maps (Phillips et al. 2006). ENMs were converted to binary (suitable/unsuitable) using the threshold of maximum sum of test sensitivity and specificity (Liu et al. 2013).
Especially when comparing fragmentation metrics using both a county-based range map and an ENM for a given species, it is important to minimize overprediction of the ENM into areas where the species probably does not actually occur due to various factors other than climate. To minimize overprediction, the resulting ENM raster grids are clipped using one of two masks: if the county-based range map is completely contiguous (D = 0), the respective ENM is clipped to the extent of the range polygon. If the county-based range map consists of multiple patches, the ENM is clipped to the extent of the minimum convex polygon encompassing all patches of the range. Clipping to the minimum convex polygon allows for the possibility of species with disjunct ranges but more contiguous ENMs. All ENMs should then be resampled (using nearest-neighbor resampling) to match the resolution of the county-based range rasters.
An important caveat is that NN should only be used with polygon-based ranges, such as the county range maps. Because Fragstats considers each isolated cell a separate patch, and NN is averaged over all patches, extraneous raster cells have a large impact on NN. Thus, highly patchy rasters such as ENMs would have biased NN values. Fragstats contains a wide variety of metrics at different scales to describe size, shape, and composition of a landscape. Future extensions of the DIVISION index applied to quantifying species range disjunction should attempt to incorporate degree of patch isolation that is appropriate for all types of species ranges.

Example: amphibians of the southeastern United States
To demonstrate the utility of our proposed method, we quantify fragmentation of species ranges of amphibians in the southeastern United States (hereafter, "the Southeast"). Amphibians are particularly vulnerable to even subtle changes in climate (Duellman and Trueb 1994) and face a growing threat from pathogens such as Batrachochytrium (Olson et al. 2013) and Ranavirus (Price et al. 2014). Of course, the extent of population vulnerability depends on a variety of intrinsic and extrinsic factors such as microhabitat and genetic variation, but it is nevertheless clear that there is an urgent need to understand factors driving population decline.
The Southeast is a hotspot for amphibian biodiversity (Rissler and Smith 2010), in large part because the region was not directly affected by Pleistocene glaciation. However, alternating high and low sea levels throughout the glacial cycle caused river drainage fragmentation and fusion (Kozak et al. 2006) and seawater flooding of coastal streams (Wright and Frey 1965). Rivers have been shown to be biogeographic barriers for amphibians (e.g., Kozak et al. 2006, Pauly et al. 2007, Shepard and Burbrink 2011, Herman and Bouzat 2016, and eight major southeastern rivers drain into the Gulf of Mexico: Mississippi, Pearl, Pascagoula, Mobile, Escambia-Conecuh, Choctawhatchee, Apalachicola, and Suwannee, as well as two major tributaries of the Mobile River: Upper Tombigbee and Sipsey (Ward et al. 2005). In addition, multiple physiographic provinces come into contact in the Southeast, providing potential opportunities for selection based on habitat; the Coastal Plain meets foothills in all but three states (Louisiana, Mississippi, Florida), along a boundary known as the Fall Line. In the Southeast, habitat loss is the primary factor driving amphibian population decline (Stuart et al. 2004).

Study species
We followed Mitchell and Gibbons (2010) in defining the southeastern United States as including the following states: Alabama, Florida, Georgia, Kentucky, Louisiana, Mississippi, North Carolina, South Carolina, Tennessee, and Virginia. The selection of this area was based on its unique topographical and environmental features, such as physiographic provinces, habitat type, and climate. Most terrestrial and semi-aquatic amphibian species with geographic ranges either partially or entirely within at least one of the southeastern states were initially included in this study, but a few species were excluded from analyses due to exceptional difficulty in defining range extents. Our data set included 74 salamander species and 36 frog species (Supplementary Information, Table S1).

Quantification of amphibian range disjunction
We compared our fragmentation metric across two different methods of range mapping: county-based ranges and ENMs. For each of the 110 species, we downloaded a georeferenced map of the county-level geographic range in ESRI shapefile format from the IUCN Red List (IUCN 2016) or other sources for a few species with outdated IUCN maps (Supplementary Information, Table S1). Manipulation and reformatting of range maps were completed in R v.3.3.2 (R Core Team 2016). We calculated D for all 110 species and NN for all species with D > 0, and we classified the species ranges into three groups -low, moderate, and high fragmentation or isolation -for each metric based on visual clustering of histogram bars.
For the 110 amphibian species included in our county-based range analysis, D ranged from 0 -0.6732 (Fig. 1). Four species were classified as having the highest degree of range disjunction: Plethodon serratus (D = 0.6732), Rana sevosa (D = 0.6667), Plethodon websteri (D = 0.5876), and Hyla andersonii (D = 0.522). Each of these species ranges consists of multiple patches of similar size (Fig. 2). Eight species ranges were classified as moderately disjunct, with D values of 0.4056 -0.2471. The remaining 98 species ranges were classified as having low or no disjunction  Of those 98 ranges, 65 were completely contiguous and thus had a D of 0. NN was calculated for 44 species and ranged from 5 -481 km (Fig. 1). Three species were classified as having a high NN; two of those species (Hyla andersonii, Plethodon serratus) also had a high D, while the third species (Hyla chrysoscelis) had a low D. Species ranges with moderate or high D and moderate or high NN are characterized by multiple fragments of similar size isolated by moderate to large distances (Fig. 2). Similarly, species with moderate or high D and low NN are characterized by multiple fragments of similar size but isolated by smaller distances. Non-contiguous ranges with low D and all values of NN are characterized by a large primary patch with one or more much smaller peripheral patches.
To examine the effects of finer-scale range mapping on D, we selected a subset of 65 species for further analysis using ENMs. Species were selected based on availability and geographic coverage of georeferenced locality data for specimen records for a species (see below). For example, a species was omitted if it had fewer than 20 georeferenced specimen records or a large portion of its range was missing specimen records (e.g., >50% of patches of a patchy range missing records or greater than approximately 50% of a contiguous range missing records). We calculated D for each ENM in Fragstats. ENMs were classified into high, moderate, and low fragmentation based on visual clustering of histogram bars. We tested the effect of range mapping method (county-based versus ENM) on D using a Wilcoxon signed rank test (Wilcoxon 1945).
Of those species, one (Plethodon websteri) had a highly disjunct county-based range, and three had a moderately disjunct range. All range maps and ENMs, with associated values for D and NN, are shown in Supplementary Information, Figs. S1 and S2. We tested for an effect of raster resolution on D and NN by creating new county-based range rasters for all species with D > 0 and new ENM rasters for all species in the original ENM subset; these new rasters all had a resolution of 400 m. We chose this high resolution because very small ranges cannot be scaled down to much lower resolution, and resolutions higher than approximately 400 m would have generated unmanageably large files for large ranges. We calculated D for all new rasters and NN for all new county-based range rasters. We then conducted Wilcoxon signed rank tests to compare these results with the results from rasters with differing resolutions. Using a single resolution across all rasters -versus the original data set of rasters with resolutions based on extent of range shapefile -did not qualitatively affect D or NN values. D did not significantly differ between the two methods for county-based ranges (p = 0.88) or ENMs (p = 0.067). Raster resolution did statistically affect NN for ranges (p < 0.001); however, we argue that this difference is not meaningful in this case because the data sets are highly similar qualitatively (Fig. S3), and we present NN as a useful metric in comparative analysis, rather than an exact description of fragment isolation.

Effects of mapping method
In comparing D for county-based ranges and ENMs, the primary pattern that emerges involves contiguous or nearly-contiguous ranges that show moderate or high ENM fragmentation (Fig. 4). This is not surprising, as range maps constructed at the county scale obviously generally overpredict the actual species distribution. However, specific cases of difference or similarity between ranges and ENMs can hint at processes potentially driving these patterns. In some cases, conflicting D between a species' range and its ENM is due to range maps not taking into account major geographic features dividing populations. For example, in our analyses, the species with high ENM fragmentation but low or moderate range fragmentation are restricted to high elevation mountaintops (Desmognathus santeetlah, Plethodon welleri) or either side of a river (Desmognathus wrighti). In other cases, though, the range map does incorporate such unsuitable habitat; for example, the range of Plethodon montanus is moderately disjunct, depicting several small isolated patches representing mountaintops. Another source of conflicting D values for a range and ENM is the extent to which putative ENM overprediction (predicted presence outside of range boundary) is excluded, as well as the threshold used to convert the ENM from a continuous scale to binary. For example, Plethodon shermani has a moderately disjunct range with two patches, but its ENM is nearly contiguous because we allowed for overprediction between disjunct range patches. In addition, there are no strict guidelines for selecting ENM thresholds, but the threshold used in our analyses -maximizing the sum of sensitivity and specificity -has been shown to be the most consistently high-performing across data sets (Liu et al. 2013). Applying a higher threshold causes ENM patches to shrink and, if the threshold is high enough, become more disjunct, increasing ENM fragmentation. Our study focuses on quantifying fragmentation and uses ENMs to illustrate the often dramatic effects of the definition of a species range on an objective metric. Discussion of methods to improve individual ENMs is beyond the scope of this study.
A third source of conflict between level of fragmentation of a species' range and its ENM is inherent to the organism and/or its environment, rather than an artifact of the methods used to depict ranges and ENMs. In one case, Plethodon websteri, the species range was highly disjunct, but the ENM was essentially contiguous, and it remained contiguous until the threshold was set unreasonably high (data not shown). This suggests that factors other than those included in the model may be driving the highly fragmented shape of the range of Plethodon websteri.

Change in fragmentation through time
An important aspect of many studies at all scales concerns change in fragmentation through time.
To demonstrate the utility of D and NN for assessing change in range fragmentation through time, we tested the assertion made in a previous study (Newman and Austin 2015) that the distribution of Plethodon serratus was more contiguous during the last glacial maximum (LGM; ~21,000 YBP) than today. That conclusion was drawn from visual comparison of the species' ENM and paleodistribution models generated by projecting the ENM onto climate layers from the LGM. Here, we tested the validity of that assessment by using D to quantify the fragmentation of the ENM and paleodistribution model. First, an ENM was generated using only the climate layers. To generate the paleodistribution model, LGM climate layers based on the Community Climate System Model (CCSM) simulation were downloaded from Worldclim at a spatial resolution of 5 km 2 . The paleodistribution model was converted to binary as above and resampled to match the resolution of the ENM. Because the paleodistribution model is a projection of the ENM, we did not clip its extent. D was calculated for the paleodistribution model and compared to D for the species' climate-based ENM. Unlike its climate-based ENM (D = 0.6322), the paleodistribution model for Plethodon serratus had a low degree of fragmentation (D = 0.0556; Fig. S4), confirming the previous assessment that the range of Plethodon serratus was much more contiguous during the LGM than today.

Discussion
Studies of historical biogeography and comparative phylogeography often examine patterns of species range shape in a particular region, yet no simple metric exists to objectively quantify these patterns. We have demonstrated that the DIVISION index is useful for quantifying broad-scale geographic range disjunction using different types of range maps. D coupled with a measure of NN provides additional information by incorporating degree of isolation. The biological effects of range disjunction and habitat fragmentation depend heavily on the specific ecology of the species in question. The method we present here does not address, for example, the permeability of the matrix between habitat patches, nor does it consider the dispersal ability of a species. Instead, this is a broadly applicable method for simply quantifying patchiness and isolation of patches, but the biological meaning of the D and NN values obtained must be interpreted in the context of the species of study.
Population allopatry is driven by a variety of factors (Raven 1972), and the scale of allopatry considered will depend on the questions of study. In many cases, broad-scale range disjunction, such as the range of Plethodon serratus and ranges of other species with high D and high or medium NN, is driven largely by historical factors over tens of thousands of years, such as glacial cycles, and is of primary interest to many researchers in phylogeography and historical biogeography (Avise 2000). Conversely, smaller-scale habitat fragmentation is often of high interest to researchers in conservation. While this method presented here can easily be applied to assessments of finer-scale habitat fragmentation by, for example, using higher resolution land cover and land use raster layers to generate ENMs, the focus of our study on broad-scale range shape addresses a gap in the biogeographic literature that is generally outside the scope of ecological studies that describe and quantify local habitat fragmentation. But cases such as Hyla andersonii, which has a highly disjunct range and also a highly specialized habitat (Warwick et al. 2015), point to the need for methods that are able to quantify allopatry at any scale.
The example presented here for Plethodon serratus shows how the DIVISION index can be used to objectively quantify change in fragmentation. As global temperatures are predicted to continue rising over the next 100 years (IPCC 2013), studies using modeling to predict the future distribution of a species will become more critical to conservation decisions and management efforts. DIVISION and NN can quantify the change in range fragmentation over time that otherwise would be assessed visually on a map. This method, and future extensions that further incorporate degree of isolation, provides a foundation for studies in fields encompassing historical biogeography, phylogeography, conservation, and ecology.