DOI: 10.5061/dryad.08kprr55b
title: Data from: Hindcast-validated species distribution models reveal future vulnerabilities of mangroves and salt marsh species

description: Rapid climate change threatens biodiversity via habitat loss, range
 shifts, increases in invasive species, novel species interactions, and
 other unforeseen changes. Coastal and estuarine species are especially
 vulnerable to the impacts of climate change due to sea level rise and may
 be severely impacted in the next several decades. Species distribution
 modeling can project the potential future distributions of species under
 scenarios of climate change using bioclimatic data and georeferenced
 occurrence data. However, models projecting suitable habitat into the
 future are impossible to ground truth. One solution is to develop species
 distribution models for the present and project them to periods in the
 recent past where distributions are known to test model performance before
 making projections into the future. Here, we develop models using abiotic
 environmental variables to quantify the current suitable habitat available
 to eight Neotropical coastal species: four mangrove species and four salt
 marsh species. Using a novel model validation approach that leverages
 newly available monthly climatic data from 1960-2018, we project these
 niche models into two time periods in the recent past (i.e., within the
 past half-century) when either mangrove or salt marsh dominance was
 documented via other data sources. Models were hindcast-validated and then
 used to project the suitable habitat of all species at four time periods
 in the future under a model of climate change. For all future time
 periods, the projected suitable habitat of mangrove species decreased, and
 suitable habitat declined more severely in salt marsh species.
description: <em>Data acquisition</em>
             We obtained specimen-based occurrence data for each
 species from iDigBio (Integrated Digitized Biocollections; idigbio.org)
 and GBIF (Global Biodiversity Information Facility; gbif.org) and
 supplemented these data with locality data from personal collections for
 three mangrove species (<em>Avicennia germinans</em>,
 <em>Laguncularia racemosa</em>, <em>Rhizophora
 mangle</em>). Four of the species included in the analysis are
 mangroves (<em>Avicennia germinans</em>, black mangrove;
 <em>Laguncularia racemosa</em>, white mangrove; and
 <em>Rhizophora mangle</em>, red mangrove) or
 mangrove-associated species (<em>Conocarpus erectus</em>,
 buttonwood). For simplicity, these four species will hereafter be
 collectively referred to as ‘mangroves,’ even though <em>Conocarpus
 erectus</em> is not considered a true mangrove (Tomlinson, 2016). We
 also selected four salt marsh species (<em>Batis
 maritima</em>, turtleweed; <em>Sesuvium
 portulacastrum</em>, sea purslane; <em>Spartina
 alterniflora</em>, smooth cordgrass; and <em>Sporobolus
 virginicus</em>, seashore dropseed) for analyses. These four species
 were selected because they occur in close proximity to one
 another—indicating the presence of salt marsh habitat—and because of their
 broad and exclusively coastal distributions in the Neotropics. We used SDM
 to investigate changes in suitable habitat for all eight species. The raw
 data were cleaned using standard approaches and R scripts (e.g., Marchant
 et al., 2017); duplicates and incorrect data (e.g., latitude and longitude
 of 0) were removed from the data set (all scripts used in this paper were
 deposited in GitHub (github.com/richiehodel/coastal_ENM), and all cleaned
 occurrence data, layers, and models were deposited in Dryad. We included
 species that had exclusively coastal or estuarine distributions, and only
 species with at least 50 occurrence points (after cleaning) were used in
 the analyses. Given the complexities of the modeling approach, we focused
 on the Neotropics as opposed to a global analysis; only mangrove and salt
 marsh species with native ranges in the Americas were used (i.e.,
 cosmopolitan species were excluded). Certain species that inhabit salt
 marshes, but that have extensive inland distributions, including
 freshwater wetlands, were excluded (e.g., <em>Distichlis
 spicata</em>).             We acquired
 bioclimatic environmental layers from Worldclim 2.1 (worldclim.org; Fick
 & Hijmans, 2017) for multiple time periods. The bioclimatic
 layers, which contain temperature and precipitation data for every
 continent except Antarctica, have been used extensively and successfully
 in SDM studies (Booth, 2018). In Worldclim 2.1, annual precipitation,
 maximum temperature, and minimum temperature data are available for every
 month from 1960-2018 at 2.5 arc minute resolution; these three variables
 can be used to calculate values of all 19 bioclimatic variables (Harris et
 al., 2014; Fick & Hijmans, 2017; Hijmans et al., 2017). We
 considered the present to be 2013-2018, the 1980s salt marsh dominance
 period to be 1984-1989, and the early 2000s mangrove dominance to be
 2001-2006. These time periods were selected to capture the optimal amount
 of either mangrove or salt marsh dominance during each documented
 oscillation (Cavanaugh et al., 2019), and we selected these windows of
 time so that the present and past time periods were all six years.
 Although many of the study species may be longer-lived than each of the
 time periods (i.e., six years), we prioritized using time periods that
 captured either mangrove or salt marsh dominance. Due to the oscillations
 of mangrove versus salt marsh dominance, many individual plants were
 likely exterminated on short time scales. We used all occurrence data to
 construct an SDM for each species for our defined present time (2013-2018)
 regardless of when the specimens were collected. It would be ideal to use
 separate occurrence specimens from each time period to assess SDM
 performance, but this was not possible with the temporal distribution of
 georeferenced data points. For each six-year time period, we averaged the
 annual precipitation, maximum temperature, and minimum temperature for
 each month (e.g., average values of these three variables were calculated
 across the six January months, six February months, etc. in each time
 period) and used the resulting 12 monthly averages to calculate the
 standard 19 bioclimatic variable values using the ‘biovars’ function in
 the ‘dismo’ R package for each six-year time period (Hijmans et al.,
 2017). The standard 19 bioclimatic variables are not available on a
 monthly basis because some of them incorporate seasonality and require
 data for at least one year. By using monthly data for annual
 precipitation, maximum temperature, and minimum temperature variables, all
 of the 19 bioclimatic variables can be calculated (Hijmans et al.,
 2017). All layers were then trimmed so that the extent
 of the study area was between -120 and -32 degrees longitude, and -36 and
 36 degrees latitude using custom scripts and the R package ‘raster’
 (Hijmans et al., 2015) and exported in ASCII format (Fig. 1). This study
 area was selected because it included subtropical and tropical regions of
 both the Northern and Southern Hemispheres, captured the ecotone between
 mangrove and salt marsh species in both Hemispheres, and allowed for an
 expansion zone as some species may expand their ranges in the future as
 the climate changes. Regions such as Hawaii, where some Neotropical
 mangrove species have been introduced, were not included in the study. We
 used an R script and the R package ‘raster’ (Hijmans et al., 2015) to
 measure the pairwise correlation of the 19 bioclimatic variables. When
 variables were correlated with one another (<em>r</em>
 &gt; 0.7), only one of the layers was retained for subsequent analyses
 (Dormann et al., 2013). After removing correlated layers, we had a data
 set of six bioclimatic variables (BIO2, mean diurnal temperature range;
 BIO5, maximum temperature of the warmest month; BIO6, minimum temperature
 of the coldest month; BIO12, annual precipitation; BIO15, precipitation
 seasonality; BIO18, precipitation of warmest quarter). BIO6 and BIO1 were
 highly correlated (<em>r</em> = 0.956), and BIO1 (mean annual
 temperature) was excluded even though it is frequently included in SDM
 analyses because BIO6 has been identified as an important variable shaping
 range limits of coastal species (Tomlinson, 2016). All layers were clipped
 using the ‘mask’ function in the ‘raster’ R package (Hijmans et al., 2015)
 such that all cells with an elevation greater than 10m were considered ‘no
 data’ cells. This was done to ensure that the SDM analyses were not
 trained on inland regions representing areas where these coastal species
 do not occur.   <em>Species
 Distribution Modeling</em> The occurrence data
 obtained from digitized herbaria records and the six environmental layers
 were used as input for the SDM analyses. SDM uses the occurrence data for
 each species in the present to identify pixels that have suitable habitat
 for the species of interest based on environmental data. We used the
 maximum entropy algorithm implemented in MAXENT v3.4.1 (Phillips et al.,
 2006; Phillips et al., 2017) to conduct SDM analyses. The maximum entropy
 algorithm uses presence data and random background sampling to develop the
 model, and it has been shown to perform well with presence-only data
 (Elith et al., 2006; Wisz et al., 2008). Optimal settings for MAXENT model
 fit were determined using the ‘ENMevaluate’ function in the ENMeval R
 package (Muscarella et al., 2014). We investigated regularization
 multipliers from 0.5 to 4 at intervals of 0.5 and the following
 features/combinations of features: linear, linear/quadratic,
 linear/quadratic/hinge, linear/quadratic/hinge/product,
 linear/quadratic/product/threshold, and
 linear/quadratic/hinge/product/threshold. The ‘ENMevaluate’ function was
 run for each species, using the same 10,000 background points, occurrence
 data for the species of interest, and the ‘maxnet’ algorithm with the
 ‘checkerboard2’ method. The DAICc scores for all models tested for each
 species were compared to determine the optimal model to be inputted into
 MAXENT. Other non-default settings used include five-fold
 cross-validation, a minimum training presence threshold, and fading by
 clamping. Cloglog output was used because it produces an estimate for each
 pixel between 0 and 1 that represents the probability of presence
 (Phillips et al., 2017). We assessed each model’s
 prediction ability by using partial receiver operating characteristic
 (pROC), which measures the ratio of the area under the receiver operating
 characteristic curve (AUC). AUC ranges from zero to one and measures the
 model’s ability to predict suitable habitat, with one indicating perfect
 discrimination between suitable and unsuitable habitat. The pROC is the
 ratio of the partial AUC divided by random expectation, and it can range
 from 0 to 2, with 1 representing random model performance (Escobar et al.,
 2018). For independent occurrence points, this metric measures the
 relationship of omission error and proportion of suitable area under
 conditions of low omission errors (Peterson et al., 2008). Jackknife tests
 of regularized training gain were used to measure the relative
 contribution of each bioclimatic variable to the model. Average habitat
 suitability values for each pixel were modeled for the present, past, and
 future of each species, and these values were used in downstream analyses.
 For a given region, the sum of habitat suitability scores was considered
 the total suitable habitat.  
 <em>SDM validation and projection into
 future</em> We measured each species’ suitable
 habitat and how it was projected to change from the present to the past.
 First, we defined an area representing the northeastern Florida ecotone
 used in previous studies to use for hindcast validations: between -82 and
 -80 degrees longitude, and 28 and 31 degrees latitude (Cavanaugh et al.,
 2019). For convenience, we hereafter refer to this region as ‘NE Florida.’
 We considered the SDM to be properly fit when it accurately inferred the
 anticipated relative change in suitable habitat between the average
 mangrove species and the average salt marsh species for all past time
 periods in the NE Florida validation region. We also used a larger
 geographic region (between -87 and -79 degrees longitude, and 24 and 31
 degrees latitude) to test if the hindcast validations were consistent when
 a larger region was used; we hereafter refer to this region as ‘Florida’
 for simplicity. Once the SDMs were hindcast-validated,
 the same approach was used to infer the projected change in suitable
 habitat for the three time periods in the future. To project future values
 of environmental variables, we used a widely used and well-validated
 climate model—the CNRM-CM6-1 model, which is a fully-coupled
 atmosphere-ocean general circulation model developed by Centre National de
 Recherches Météorologiques (CNRM) for the sixth generation of the IPCC
 Coupled Model Intercomparison Project 6 (CMIP6), and with the shared
 socioeconomic pathway 245 (Eyring et al., 2016). This climate model was
 selected because it is one of 49 used in the most recent IPCC CMIP6, is
 compatible with the WorldClim 2.1 data used for hindcast analyses
 (https://worldclim.org/data/cmip6/cmip6climate.html), and the shared
 socioeconomic pathway was selected because it represents a central part of
 the range of plausible future pathways. SDMs for the full geographic study
 region were projected into both the past and future time periods. Future
 time periods were determined by the availability of WorldClim 2.1
 data.
description: The programx MAXENT v3.4.1, which was used for SDM analyses, is
 freely available. ASC layers can be viewed using the open-source program
 QGIS. R and RStudio (both freely available) and associated open-source
 packages were used to process and analyze data.