Based upon code from the National Fishes Vulnerability Assessment Project
===============================
Sam Silknetter (SCS) Updated 04/21/2020 Traci P. DuBose Updated 09/30/2022
This repository contains R code for analyzing species rarity and climate sensitivity for 90 anuran species native to the conterminous United States. We use downloaded observations from Global Biodiversity Information Facility (GBIF) and HerpMapper to describe where these species occur. Once filtered, occurrence points are then used to determine area of occupancy (AOO) for two grain sizes: small occupied watersheds (USGS HUC-12 sub-basins [1] and HydroBASINS level 12 [2]) and buffered occurrence points (1km radius). Bioclimatic variables calculated from WorldClim [3] are then used to determine species climate sensitivity (CS) for each species' grain-specific AOO. Climate variables examined include mean annual temperature (bio1), maximum temperature of the warmest quarter (bio5), minimum temperature of the coldest quarter (bio6), annual precipitation (bio12), and precipitation seasonality (bio15). At each grain size, CS and AOO metrics are combined into a relative rarity and climate sensitivity index (RCS)[4]. We then investigate the relationships between intrinsic sensitivity, taxonomy, and three conservation lists.
The analyses herein are intended for use by wildlife managers and conservation practitioners to identify species that have high relative, intrinsic sensitivity to changes in climate. By using relative metrics that can be applied across different spatial scales, these assessments of geographic sensitivity allow for direct comparisons between species with variable data availability, including species that are both well- and poorly-studied. One note: as a relative metric, the ranking and index value are sensitive to the spatial extent considered and the taxa included in the analysis.
[1] U.S. Geological Survey [USGS]. 2020. National Watershed Boundary Dataset (ver. USGS National Watershed Boundary Dataset in FileGDB 10.1 format (published 20200701). Available from https://www.usgs.gov/core-science-systems/ngp/national-hydrography/access-national-hydrography-products (accessed July 14, 2020).
[2] HydroBASINS level 12 dataset - Lehner B, Grill G. 2013. Global river hydrography and network routing: baseline data and new approaches to study the world’s large river systems. Hydrological Processes 27:2171–2186. https://hydrosheds.org/page/hydrobasins
[3] WorldClim Datasets - Fick SE, Hijmans RJ. 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37:4302–4315. https://www.worldclim.org/data/worldclim21.html
[4] Mims, M. C., D. H. Olson, D. S. Pilliod, and J. B. Dunham. 2018. Functional and geographic components of risk for climate sensitive vertebrates in the Pacific Northwest, USA. Biological Conservation 228:183-194.
===============================
Create Watershed Boundaries .rds and climate raster files -- "0_Creating_HUC12_RDS.Rmd"
Code that creates the spatial rds files used to quantify area of occupancy. You can download the USGS National Watershed Database geodatabase here. Within this script, we also create the subset of HydroBasin level 12s used to descrbe the area outside of the Contiguous US. The HydroBASINS geodatabase can be found here. In order to still incorporate USGS HUC-12 watershed subbasins, we stiched together HydroBasins's level 12 watershed polygons from outside the contiguous US with USGS's HUC-12 polygons. Each polygon has a unique ID that relates it to the larger watershed area. This script 1) isolates USGS watersheds within the contiguous US, 2) removes overlapping polygons between the HydroBasins and HUC-12 datasets, and 3) creates HydroBasin clipped polygons that bridge the two datasets described above. It also outputs three RDS files: a projected map of the contiguous US, a projected map of all HUC12s within the contiguous US, and a projected map of HydroBasins within the North American continent but outside the contiguous US and our created clipped HydroBasin polygons. We also calculate the average bioclimatic variables from 1961-2018 WorldClim data by calculating the bioclim variables from 1961-1969 and 2000-2018 and combing it with the 1970-2000 average so that each year is weighted equally.
Download Occurrence data -- "0_GBIF_occurrence_download.Rmd"
We downloaded GBIF occurrences using the rgbif R package to query the GBIF API. Chloe Moore and Tess Alexander downloaded the HerpMapper occurrences in September 2019. These occurrence datasets are combined using lines 104 to 159. To access our GBIF downloaded occurrences, use the download codes found in GBIF_codes20210330.txt and the occ_download_get function from the rgbif package.
Record conservation status -- "1_Anuran_Conserv_Status.Rmd"
Using the rredlist package, information from the National Species of Greatest Conservation Need, and the federal Endangered Species Act list, we record the conservation status of the anurans in our species list for those three geopolitical groups. This script produces anuran_conservation_info.csv, which is used as a table key in later scripts.
- Table S2.1: Taxonomic and Conservation status for anurans analyzed in this study. Includes number of observations used in analysis
- Table 1: N of species within each conservation status for the three geopolitical groups
Create shapefiles of anuran range and calculate AOO -- "2_Area_of_Occurrence_Code_Anurans.Rmd"
Code to calculate area of occupancy (AOO) at two grains (small watersheds and 1km Point Buffer) for all species in the final list. Aloso has a for loop to count the observations lost because of the filter by the native range. Output files include shapefiles at both grains/species, and a CSV with AOO values and ranks for both grain sizes (AOO HUC12 Output_DATE.csv), and an intermediate file that describes the AOO within and outside the contiguous US for both grain sizes (AOOs_raw_DATE.csv).
- Unused Figure: Example of occupied watershed across US border and the IUCN native range
- Unused Figure: Plot to compare area of watershed polygons within the USGS HUC-12 and HydroBasin datasets
- part of Table S2.1: Observation before 1960 and outside the conterminous United States
Identify Climate Outliers -- "3b_Identify_Climate_Outliers.Rmd"
After first generating shapefiles using 2_Area of Occupancy Code_Anurans.R, we can then pair those shapefiles with climate data to identify points found in climates that are far outside the norm for that species. These points might be poorly georeferenced, and so we want to remove them from our analysis. See the anuran QA/QA document (Appendix 3) for further explanation. This code pairs each point with climate data to estimate the point's z score for three climate variables. Then, points with z scores that have an absolute value greater than 4 are removed (larger for more wide ranging species). After this script is run, 2_Area of Occupancy Code_Anurans.R should be run again to remove those outliers from AOO calculation.
Calculate the Climate Specificity of each species -- "3c_Climate_Sensitivity_Anurans.Rmd"
This code must be run on the super computer (need higher working memory) and exists to calculate standard deviations of climate variables for each species. This requires climate variables downloaded and calculated in Script 0_Creating_HUC12_RDS.Rmd. Two large outputs are the average climate values within each occupied polygon (taxa_ws_climate_values_DATE.csv) and six csvs for each species for the climate values within the buffered points (SPP_climate_raster.csv). It also creates the intermediate files that calculate climate sensitivity across the North American range for these species: buff_climate_sums_DATE.csv and range.wide.weighted.CS_DATE.csv.
Bring it all together to calculate the RCS index -- "4_RCS_Index_Code.Rmd"
This code calculates Climate Sensitivity (CS) from standard deviation values generated in Script 4. Area of Occupancy (AOO) data (generated in Script 2) and CS data are merged into a single RCS output for all species and grains. Finally, the Relative Climate Sensitivity index (RCS) is calculated for each species and range metric. Three files are created: RCS_Table_DATE.csv, RCS_table_L48_DATE.csv, and the combination of the two former files Anuran RCS values 20210405.csv.
- Table S2.2 Spearman correlations among spatial extent & grain size
- Unused Figure: How AOO varies across spatial extent
- Figure 1: RCS Dot Plot in ggplot
- Unused Figure: Sample Size effect on RCS calculation
- Figure S2.6: Distribution of Anuran AOO
- Figure S2.7: AOO dot plot to depict rank vulnerability
- Figure S2.8: Climate Sensitivity dot plot to depict rank vulnerability
- Unused Figure: Distribution of Climate SDs
- Figure S2.1: How spatial extent alters RCS calculation - not valid anymore as not all occurrences downloaded (20% cut off)
Assess alignment between conservation status and RCS values -- "5_RCS_Predict_Conser_Status.Rmd"
Combining the output of 1_Anuran_Conserv_Status and 4_RCS_Index_Code, we calculate differences in mean rank RCS among the different conservation statuses in different conservation geopolitical groups and spatial extents. We also build a logistic regression (using gams to account for non-linearity to assess if RCS can predict protection status. Overlooked species are identified based on different thresholds for each conservation list. We also assess how RCS values are distributed among anuran genera and families.
- Table S2.4: Differences in rank RCS value among conservation statuses of three lists
- Figure 3abc: Mean rank RCS values and conservation status boxplot
- Table 2: Alignment between conservation status and the RCS index
- Figure 3def: RCS poorly predicts protection status GAMs
- Figure 2: RCS values and protection status among anuran genera
- Table S2.3: Chi square tests to assess overrepresentation in listed genera
- Figure S2.9: Differences in RCS values and protenction status among anuran families
- Figure S2.4: Relationship between AOO and times listed on a SWAP
Quantify Results -- "6a_Results_paragraph.Rmd"
Code used to summarize and evaluate the results of the above analyses. These summarizations are covered in the results section of the manuscript.
Generate Maps to explain results -- "6b_Anuran_map_figures.Rmd"
This code aggregates species occurrence at the HUC-12 level to the HUC-8 level. We also calculate the mean RCS value at the HUC-8 level.
- Figure S2.5a : Anuran species richness at the HUC-8 grain size
- Figure S2.5b: Average RCS value at the HUC-8 grain size
- Figure S2.10: Example of climate specificity among scales with Ascaphus montanus
- Unused Figure: Example RCS calculation using A. crepitans
Quality Control and Assurrance -- "7_anuran_rcs_qaqc.Rmd"
An RMarkdown document recording quality assurance and quality control measures used to evaluate the above analyses. All figures and tables in Appendix 3 are found in this RMarkdown.
- Check that raster::extract correctly computing climate values
- Check that the units match between our two datasets PRISM & WorldClim
- Check for NA values extracted - points that land on the border
- Build a histogram of each species’ climate values (is the SD capturing climate breadth well)
- Investigate appropriate cut offs for extreme climate values (which might indicate poorly georeferenced points)
Deciding Point at Climate Extreme Cutoffs -- "climate_cutoff_desc.Rmd"
An RMarkdown document recording variation in the z scores of climate variables associated with points potentially used in this analysis. This Rmd made z score cut off decisions needed in 3b_Identify_Climate_Outliers.Rmd. Figures and tables in Appendix 3 are found in this RMarkdown.
===============================