As part of the Maine Beach Mapping Program (MBMAP), MGS surveys annual alongshore shoreline positions (see Beach_Mapping_Shorelines). Using these shoreline positions and guidance from the USGS Digital Shoreline Analysis System (DSAS). DSAS is referenced as Thieler, E.R., Himmelstoss, E.A., Zichichi, J.L., and Ergul, Ayhan, 2009, Digital Shoreline Analysis System (DSAS) version 4.0— An ArcGIS extension for calculating shoreline change: U.S. Geological Survey Open-File Report 2008-1278. For more information on DSAS and the methodology DSAS employs, please see: https://woodshole.er.usgs.gov/project-pages/DSAS/. The supporting DSAS User Guide which describes how DSAS works and how statistics are calculated is available here: http://www.maine.gov/dacf/mgs/hazards/beach_mapping/DSAS_manual.pdf. MGS wrote a database procedure following protocols outlined in DSAS that allows for the calculation of different shoreline change rates and supporting statistics. This was done so that MGS no longer needed to depend on USGS updates to the DSAS software to keep current with ArcGIS software updates. The script casts shoreline-perpendicular transects at a set spacing (in this case, 10-m intervals along the shoreline), from a preset baseline (located landward of the monitored shorelines), and calculates a range of shoreline change statistics, including: Process Time: The time when the statistics were calculated. TransectID: The ID of the transect (including the group or line section ID; for example, 1-1, is line 1, transect 1) SCE: Shoreline Change Envelope. The distance, in meters, between the shoreline farthest from and closests to the baseline at each transect. NSM: Net Shoreline Movement. The distance, in meters, between the oldest and youngest shorelines for each tranect. EPR: End Point Rate. A shoreline change rate, in meters/year, calculated by dividing the NSM by the time elapsed between the oldest and youngest shorelines at each transect. LRR: Linear Regression Rate. A shoreline change rate, in meters/year, calculated by fitting a least-squares regression line to all of the shoreline points for a particular transect. The distance from the baseline, in meters, is plotted against the shoreline date, and slope of the line that provides the best fit is the LRR. LR2: The R-squared statistic, or coefficient of determination. The percentage of variance in the data that is explained by a regression, or in this case, the LRR value. It is a dimensionless index that ranges from 1.0 (a perfect fit, with the best fit line explaining all variation) to 0.0 (a bad fit, with the best fit line explaining little to no variation) and measures how successfully the best fit line (LRR) accounts for variation in the data. LCI95: Standard error of the slope at the 95% confidence interval. Calculated by muliplying the standard error, or standard deviation, of the slope by the two-tailed test statistic at the user-specified confidence percentage. For example if a reported LRR is 1.34 m/yr and a calculated LCI95 is 0.50, the band of confidence around the LRR is +/- 0.50. In other words, you can be 95% confidence that the true rate of change is between 0.84 and 1.84 m/yr. LRR_ft: The Linear Regression Rate, converted to feet/year. LCI95_ft: The LCI95, converted to feet. EPR_ft: The End Point Rate converted to feet.

Predicted standared error of GIS geostistical interpolation of a surface that models the depth to bedrock. Derived from known georeferenced locations where depths to bedrock have been observed. These primarily include bedrock outcrops and well or boring locations. Interpolation method was ordinary kriging, using a lag size of 448.6 ft. A maximum of 20 neighbors and a minimum of 8 neighbors were used in interpolation. The predicted standard error is the standard deviation of the predicted surface, and is a function of distance from the nearest data point.

1. Biantun_toponyms_points.shp: Point data of Biantun toponyms in Yunnan. Used for kernel density analysis, Geographically weighted regression, standard deviation ellipse, etc.2. Ethnic_minority_points.shp: Point data of minority toponyms in Yunnan. Used to generate the distribution map of minority toponyms.3. dem_Yunnan.tif: Extract elevation, slope.4. Population distribution: Including the early Ming, mid-Ming dynasty, late Ming dynasty, and Qing dynasty population distribution. Used to extract the number of people corresponding to Biantun toponym points and generate Geo-informatic Tupu.5. rivers.shp: Extract the nearest distance from the Biantun toponym point to the water system.

Aim

To understand the representativeness and accuracy of expert range maps, and explore alternate methods for accurately mapping species distributions.

Location

Global

Time period

Contemporary

Major taxa studied

Terrestrial vertebrates, and Odonata

Methods

We analyzed the biases in 50,768 animal IUCN, GARD and BirdLife species maps, assessed the links between these maps and existing political and various non-ecological boundaries to assess their accuracy for certain types of analysis. We cross-referenced each species map with data from GBIF to assess if maps captured the whole range of a species, and what percentage of occurrence points fall within the species’ assessed ranges. In addition, we use a number of alternate methods to map diversity patterns and compare these to high resolution models of distribution patterns.

Results

On average 20-30% of species’ non-coastal range boundaries overlapped with administrative national boundaries. In total, 60% of areas with the highest spatial turnover in species (high densities of species range boundaries marking high levels of shift in the community of species present) occurred at political boundaries, especially commonly in Southeast Asia. Different biases existed for different taxa, with gridded analysis in reptiles, river-basins in Odonata (except the Americas) and county-boundaries for Amphibians in the US. On average, up to half (25-46%) species recorded range points fall outside their mapped distributions. Filtered Minimum-convex polygons performed better than expert range maps in reproducing modeled diversity patterns.

Main conclusions

Expert range maps showed high bias at administrative borders in all taxa, but this was highest at the transition from tropical to subtropical regions. Methods used were inconsistent across space, time and taxa, and ranges mapped did not match species distribution data. Alternate approaches can better reconstruct patterns of distribution than expert maps, and data driven approaches are needed to provide reliable alternatives to better understand species distributions.

Methods Materials and methods

We use a combination of approaches to explore the relationship between species range maps and geopolitical boundaries and a subset of geographic features. In some cases we used the density of species range boundaries to explore the relationship between these and various features (i.e. administrative boundaries, river basin boundaries etc.). Additionally, species richness and spatial turnover are used to explore changes in richness over short geographic distances. Analyses were conducted in R statistical software unless noted otherwise. All code scripts are available at https://github.com/qiaohj/iucn_fix. Workflows are shown in Figure S1a-c with associated scripts listed.

Species ranges and boundary density maps

ERMs (Expert range maps) were downloaded from the IUCN RedList website for mammals (5,709 species), odonates (2,239 species) and amphibians (6,684 species; https://www.iucnredlist.org/resources/grid/spatial-data). Shapefile maps for birds were downloaded from BirdLife (10,423 species, http://datazone.birdlife.org/species/requestdis), and for reptiles from the Global Assessment of Reptile Distributions (GARD) (10,064 species; Roll et al., 2017). Each species’ polygon boundaries were converted to a polylines to show the boundary of each species range (Figure S1a-II; codes are lines 7 – 18 in line2raster_xxxx.r ; xxxx varies based on the taxa). The associated shapefile was then split to produce independent polyline files for each species within each taxon (see Figure S1a-I, codes are lines 29 to 83 in the same file above.).

To generate species boundary density maps, species range boundaries were rasterized at 1km spatial resolution with an equal area projection (Eckert-IV), and stacked to form a single raster for each taxon (at the level of amphibians, odonates, etc.). This represented the number of species in each group and their overlapping range boundaries (Figure S1b-II, codes are in line2raster_all.r). Each cell value indicated the number of species whose distribution boundaries overlapped with each cell, enabling us to overlay this rasterized information with other features (i.e. administrative boundaries) so that the overlaps between them can be calculated in R. These species boundary density maps underlie most subsequent analyses. R code and caveats are given in the supplements, links are provided in text and Figure S1.

Geographic boundaries

Spatial exploration of species range boundaries in ArcGIS suggested that numerous geographic datasets (i.e. political and in few cases geographic features such as river basins) were used to delineate the species ranges for different regions and taxa (this is sometimes part of the methodology in developing ERMs as detailed by Ficetola et al., 2014). Thus in addition to analyzing the administrative bias and the percentage of occurrence records within each species’ ERM for all taxa, additional analyses were conducted when other biases were evident in any given taxa or region (detailed later in methods on a case-by-case basis).

For all taxa, we assessed the percentage of overlap between species range boundaries and national and provincial boundaries by digitizing each to 1km (equivalent to buffering thie polyline by 500m), both with and without coastal boundaries. An international map was used because international (Western) assessors use them, and does not necessarily denote agreed country boundaries (https://gadm.org/). The different buffers (500m, 1000m, 2500m, 5000m) were added to these administrative boundaries in ArcMap to account for potential, insignificant deviations from political boundaries (Figure S1b). An R script for the same function is provided in “country_line_buffer.r”.

To establish where multiple species shared range boundaries we reclassified the species range boundary density rasters for each taxa into richness classes using the ArcMap quartile function (Figure S1). From these ten classes the percentage of the top-two, and top-three quartiles of range densities within different buffers (500m, 1000m, 2500m, 5000m) was calculated per country to determine what percentage of highest range boundary density approximately followed administrative borders. This was done because people drawing ERMs may use detailed administrative maps or generalize near political borders, or may use political shapefiles that deviate slightly. It is consequently useful to include varying distances from administrative features to assess how range boundary densities vary in relation to administrative boundaries. Analyses of relationships between individual species range boundaries and administrative boundaries (coastal, non-coastal) were made in R and scripts provided (quantile_country_buffer_overlap.r).

Spatial turnover and administrative boundaries

Heatmaps of species richness were generated by summing entire sets of compiled species ranges for each taxon in polygonal form (Figure 1; Figure S1b-I). To assess abrupt diversity changes, standard deviations for 10km blocks were calculated using the block statistics function in ArcMap. Abrupt changes in diversity were signified by high standard deviations based on the cell statistics function in ArcGIS, which represented rapid changes in the number of species present. Maps were then classified into ten categories using the quartile function. Given the high variation in maximum diversity and taxonomic representation, only the top two –three richness categories were retained per taxon. This was then extracted using 1km buffers of national administrative boundaries to assess percentages of administrative boundaries overlapping turnover hotspots by assessing what proportion of political boundaries were covered by these turnover hotspots.

Taxon-specific analyses

Data exploration and mapping exposed taxon and regional-specific biases requiring additional analysis. Where other biases and irregularities were clear from visual inspection of the range boundary density maps for each taxa, the possible causes of biases were assessed by comparing range boundary density maps to high-resolution imagery and administrative maps via the ArcGIS server (AGOL). Standardized overlay of the taxon boundary sets with administrative or geophysical features from the image-server revealed three types of bias which were either spatially or taxonomically limited between: 1) amphibians with county borders in the United States, 2) dragonflies and river basins globally and 3) gridding of distributions of reptiles. In these cases, species boundary density maps were used as a basis to identify potential biases which were then explored empirically using appropriate methods.

For amphibians, counties in the United States (US) were digitized using a county map from the US (https://gadm.org/), then buffered by with 2.5km either side. Amphibian species range boundary density maps were reclassified showing where species range boundaries existed (with other non-range boundary areas reclassified as “no data,”) and all species boundaries numerically indicated (i.e. values of 1 indicates one species range boundary, values of 10 indicates ten species range boundaries). Percentages of species boundary areas falling on county and in the buffers, in addition to species range boundaries which did not overlap with county boundaries were calculated to give measures of what percentage of the species boundaries fell within 2.5km of county boundaries.

For Odonata, many species were mapped to river basin borders. We used river basins of levels 6-8 (sub-basin to basin) in the river hierarchy (https://hydrosheds.org) to assess the relationship between Odonata boundaries and river boundaries. Two IUCN datasets exist for Odonata; the IUCN Odonata specialist group spatial dataset

As part of the Maine Beach Mapping Program (MBMAP), MGS surveys annual alongshore shoreline positions (see Beach_Mapping_Shorelines). Using these shoreline positions and guidance from the USGS Digital Shoreline Analysis System (DSAS). DSAS is referenced as Thieler, E.R., Himmelstoss, E.A., Zichichi, J.L., and Ergul, Ayhan, 2009, Digital Shoreline Analysis System (DSAS) version 4.0— An ArcGIS extension for calculating shoreline change: U.S. Geological Survey Open-File Report 2008-1278. For more information on DSAS and the methodology DSAS employs, please see: https://woodshole.er.usgs.gov/project-pages/DSAS/. The supporting DSAS User Guide which describes how DSAS works and how statistics are calculated is available here: http://www.maine.gov/dacf/mgs/hazards/beach_mapping/DSAS_manual.pdf. MGS wrote a database procedure following protocols outlined in DSAS that allows for the calculation of different shoreline change rates and supporting statistics. This was done so that MGS no longer needed to depend on USGS updates to the DSAS software to keep current with ArcGIS software updates. The script casts shoreline-perpendicular transects at a set spacing (in this case, 10-m intervals along the shoreline), from a preset baseline (located landward of the monitored shorelines), and calculates a range of shoreline change statistics, including: Process Time: The time when the statistics were calculated. TransectID: The ID of the transect (including the group or line section ID; for example, 1-1, is line 1, transect 1) SCE: Shoreline Change Envelope. The distance, in meters, between the shoreline farthest from and closests to the baseline at each transect. NSM: Net Shoreline Movement. The distance, in meters, between the oldest and youngest shorelines for each tranect. EPR: End Point Rate. A shoreline change rate, in meters/year, calculated by dividing the NSM by the time elapsed between the oldest and youngest shorelines at each transect. LRR: Linear Regression Rate. A shoreline change rate, in meters/year, calculated by fitting a least-squares regression line to all of the shoreline points for a particular transect. The distance from the baseline, in meters, is plotted against the shoreline date, and slope of the line that provides the best fit is the LRR. LR2: The R-squared statistic, or coefficient of determination. The percentage of variance in the data that is explained by a regression, or in this case, the LRR value. It is a dimensionless index that ranges from 1.0 (a perfect fit, with the best fit line explaining all variation) to 0.0 (a bad fit, with the best fit line explaining little to no variation) and measures how successfully the best fit line (LRR) accounts for variation in the data. LCI95: Standard error of the slope at the 95% confidence interval. Calculated by muliplying the standard error, or standard deviation, of the slope by the two-tailed test statistic at the user-specified confidence percentage. For example if a reported LRR is 1.34 m/yr and a calculated LCI95 is 0.50, the band of confidence around the LRR is +/- 0.50. In other words, you can be 95% confidence that the true rate of change is between 0.84 and 1.84 m/yr. LRR_ft: The Linear Regression Rate, converted to feet/year. LCI95_ft: The LCI95, converted to feet. EPR_ft: The End Point Rate converted to feet.