DOI: 10.5061/dryad.9cnp5hqsr
title: Future sea level rise in northwest Mexico is projected to decrease the distribution and habitat quality of the endangered Calidris canutus roselaari (Red Knot)

description: Sea level rise (SLR) is one of the most unequivocal consequences of
 climate change, yet the implications for shorebirds and their coastal
 habitats is not well understood, especially outside of the north temperate
 zone. Here, we show that by the year 2050, SLR has the potential to cause
 significant habitat loss and reduce the quality of the remaining coastal
 wetlands in Northwest Mexico—one of the most important regions for
 Nearctic breeding migratory shorebirds. Specifically, we used species
 distribution modelling and a moderate SLR static inundation scenario to
 assess the effects of future SLR on coastal wetlands in Northwest Mexico
 and the potential distribution of Calidris canutus roselaari (Red Knot), a
 threatened long-distance migratory shorebird. Our results suggest that
 under a moderate SLR scenario, 55% of the current coastal wetland extent
 in northwest Mexico will be at risk of permanent submergence by 2050, and
 the high-quality habitat areas that remain will be 20% less suitable for
 C. c. roselaari. What is more, 8 out of the 10 wetlands currently
 supporting the largest numbers of C. c. roselaari are predicted to lose —
 on average — 17.8% of their highly suitable habitat areas, with two sites
 completely losing all their highly suitable habitat. In combination with
 increasing levels of coastal development and anthropogenic disturbance in
 Northwest Mexico, these predicted changes suggest that the potential
 future distribution of C. c. roselaari (and other shorebirds) will likely
 contract, exacerbating their ongoing population declines. Our results also
 make clear that SLR will likely have profound effects on ecosystems
 outside the north temperate zones, providing a clarion call to natural
 resource managers. Urgent action is required to begin securing sufficient
 space to accommodate the natural capacity of wetlands to migrate inland
 and implement local-scale solutions that strengthen the resilience of
 wetlands and human populations to SLR.
description: <strong>Study area and model
 species</strong> <em>C. c.
 roselaari</em> are long-distance migratory shorebirds that breed in
 Alaska (USA) and on Wrangel Island (Russia) during the boreal summer and
 migrate along the Pacific coast to spend their nonbreeding season in NW
 Mexico (Carmona et al., 2013). Their nonbreeding habitats are restricted
 to beaches, coastal lagoons, and deltas spanning 32 - 21° N, although
 incidental records have been recorded further south. The subspecies has an
 estimated population size of just 21,700 individuals (Lyons et al., 2016)
 and is thought to be declining, leading them to be listed as threatened
 and endangered in Canada and Mexico, respectively (COSEWIC, 2007;
 SEMARNAT, 2010). Due to their small population size, restricted range,
 high site fidelity, and specialist habitat requirements, <em>C. c.
 roselaari</em> are expected to be highly vulnerable to SLR and can
 serve as an umbrella species for other shorebirds that also spend the
 nonbreeding season in NW Mexico (Muñoz-Salas et al., 2023).
 Our study area extended from 32 - 21° N and included the states
 of Sonora, Sinaloa, and Nayarit, which we have termed the “mainland
 region,” as well as the Baja California Peninsula, which we have called
 the “peninsular region.” <em>C. c. roselaari</em> nonbreeding
 distribution is divided into three separate climatic zones: Mediterranean,
 in the northwestern tip of Baja California; arid, along the rest of Baja
 California and in western Sonora; and humid/dry tropical throughout the
 rest of the study area (Vidal-Zepeda, 2005). Except for Nayarit and parts
 of Sinaloa, coastal wetlands have a discontinuous distribution along the
 coastline, with discrete, well-defined wetlands separated by large
 expanses of arid-xeric ecosystems.
 <strong>Species Data Collection</strong>
 We assembled a dataset to model the current potential
 distribution of <em>C. c. roselaari</em> using ‘presence’
 locations collected from 2000 – 2020 in NW Mexico from three different
 sources. (1) We gathered <em>C. c. roselaari</em> sightings
 from within this time period using the <em>eBird Basic
 Dataset</em> (eBird, 2020). (2) We tracked the annual movements of
 58 <em>C. c. roselaari</em> using 3-g GPS satellite
 transmitters (PinPoint 75 Argos; Lotek Inc.)  attached to the backs of
 adults with cyanoacrolyte glue at Grays Harbor, WA, during April–May 2017
 and 2018. Tags recorded locations accurate to ±10 m every two days during
 the post-breeding period until early to mid-October when they fell off of
 molting birds or their batteries were depleted. (3) We carried out
 on-the-ground surveys from 7-15 December 2019 and 2-12 January 2020 along
 the coastlines of Sonora, Sinaloa, and Nayarit. We chose 110 survey
 locations following a stratified random sampling design that included
 sites within 4 km of the shoreline and the nearest road access. Sites were
 chosen using the land use and vegetation dataset for Mexico (INEGI, 2019),
 and consisted of habitats types that are known to be used by shorebirds.
 We conducted 5-minute point counts within 3 hours of high tide during
 which we counted every shorebird within a 400-m radius of a random site.
 Field surveys ensured that the <em>C. c. roselaari</em>
 nonbreeding distribution was well represented, helped reduce biases
 associated with eBird observations concentrated in popular birding areas,
 and added areas that could have been visited by tagged <em>C. c.
 roselaari</em> but not recorded after October — the point at which
 <em>C. c. roselaari</em> shed their satellite
 transmitters. <strong>Presence data
 preparation</strong> Records from eBird were
 manually inspected for positional uncertainty and only those that
 presented enough detail to ensure an accurate location were kept
 (<em>n</em> = 1,100). We filtered the transmitter locations by
 location class and kept only those in classes 3D and A3 (10- and
 &lt;250-m accuracy, respectively; <em>n</em> = 980). Field
 surveys revealed 10 sites in NW Mexico where <em>C. c.
 roselaari</em> presence could be detected; each site was treated as
 an individual presence record. To account for the autocorrelation typical
 of transmitter datasets, and to reduce records close to each other, we
 pooled together the transmitter, eBird, and survey locations and, then,
 rarified the dataset by eliminating points closer together than 2 km. This
 procedure reduced our final dataset to a sample size of 112
 presences.   <strong>Environmental data
 preparation </strong> We selected 10
 environmental predictor variables related to climatic conditions and the
 static physical environment that are thought to define the ecological
 niche of <em>C. c. roselaari</em> (Table 1). Climatic
 variables vary seasonally; we therefore selected mean values (e.g., mean
 temperature) from Dec and Jan, the months during which <em>C. c.
 roselaari</em> populations are thought to be most stable in their
 distribution across their non-breeding range (Carmona et al., 2013). In
 contrast, we included variables related to the static physical environment
 (e.g., elevation and distance to the nearest wetland) collected at any
 time of year within the past 10 years. All environmental predictors were
 transformed, cropped, and resampled to match their extents, projections,
 and resolutions (~0.8 x 0.8 km), respectively.  To reduce
 multicollinearity, layers with variance inflation factors (USDM package;
 Naimi <em>et</em> al., 2013) &gt;12 were excluded (Table
 1).  For instance, we included three metrics of the distance to the
 nearest wetland (Euclidean distance to intertidal habitat, Euclidean
 distance to the coastline, and Euclidean distance to the polygon of each
 coastal wetland) but only distance to the polygon of each coastal wetland
 was retained as all metrics were highly correlated and the latter provided
 the best results (Table 1). <strong>Table
 1.</strong> Environmental predictor variables used to develop
 species distribution models for <em>C. c. roselaari</em> in NW
 Mexico 
 <strong>Environmental predictor</strong>
 
 <strong>Description</strong> 
 <strong>Data Source</strong>
 Elevation
 Global elevation in meters
 WorldClim (Fick and Hijmans,
 2017) 
 Temperature Mean
 December and January temperature 
 WorldClim (Fick and Hijmans, 2017) 
 Distance to wetland
 Euclidean distance to the nearest polygon
 of coastal wetlands Own
 data 
 Distance to urban areas 
 Euclidean distance to the nearest polygon of  urban
 areas Own data based on Inegi
 (2019) 
 Land cover and vegetation 
 MCD12Q1 MODIS/Terra+Aqua Land Cover Global 500m
 Friedl & Sulla-Menashe
 (2019) https://lpdaac.usgs.gov/products/mcd12q1v006/ Nightlights Proxy for disturbance and development intensity NOAA (2020); https://ngdc.noaa.gov/eog/download.html Silt Global soil silt content 0-5 cm depth https://soilgrids.org/; Poggio et al. 2021) Clay Global soil clay content 0-5 cm depth https://soilgrids.org/; Poggio et al. 2021) Soil Ph Global soil pH x 10 in H20 0-5 cm depth https://soilgrids.org/; Poggio et al. 2021) Soil Coarseness Global coarse fragments 0-5 cm depth https://soilgrids.org/; Poggio et al. 2021)   <strong>Sea level rise scenario for 2050 </strong> To assess the inundation risk in NW Mexico in 2050, we used a static ‘bathtub’ model obtained from Climate Central (2020). The model is based on global-scale datasets for elevation, tide, and coastal flooding likelihoods for the year 2050. The model parameters included the <em>CoastalDEM</em> elevation dataset (v. 1.1; Kulp & Strauss, 2018), an RCP scenario of 4.5, and medium ‘luck’ based on mid-range results from the sea-level projection range (50th percentile). The output of this product is a spatial layer with areas identified as vulnerable to permanent SLR alone and to minor floods that may rise and fall slowly. While static models can overestimate the extent of inundation and ignore bio-geomorphological feedbacks known to buffer the effects of SLR (e.g., marsh migration; Klingbeil et al., 2021), they are better suited for large-scale inundation assessments where a lack of high-resolution data — such as in the Global South — and high computational costs limit the use of dynamic models. As such, the static models used here should be considered a “worst-case” scenario to identify areas deserving further investigation with dynamic models. <strong>Coastal wetland inventory </strong> To map all available coastal wetlands in NW Mexico, we used the <em>Land Use and Vegetation Series V</em> habitat classification (INEGI, 2019)—the most comprehensive landcover dataset available for Mexico. We retained seven classes associated with coastal wetlands and <em>C. c. roselaari</em> habitat which included: halophilic xerophyte vegetation, hydrophilic halophilic vegetation, aquaculture, bare, water body, dune vegetation and mangrove. We clipped and retained all wetlands within 50 km of the coast. After visual inspection using satellite imagery, we further refined the perimeter of each wetland to include a buffer of 500 m above the high tide line. When wetlands had well-defined limits, the perimeter was manually digitized. When limits were ambiguous, such as for the wetlands spanning the coasts of Sinaloa and Nayarit, the extent of each wetland was defined by considering the watershed to which it belonged. These fine-scale perimeter delimitations included habitats that appeared suitable for <em>C. c. roselaari</em> (e.g., shallow water, intertidal mudflats, etc.) and excluded — as best as satellite imagery allowed — non-suitable habitats, such as deep water (e.g., estuarine tidal channels). Our inventory resulted in 296 polygons representing individual wetlands or fragments of continuous wetlands. <strong>2.7 Ecological Niche Modelling</strong> We followed the modelling framework of Naimi & Araújo (2016) to estimate the current and future ecological niche of <em>C. c. roselaari</em> under SLR. This involved three steps: processing (model fitting), post processing, and model transfer. Processing and post-processing the data involved using our <em>C. c. roselaari</em> presence information to extract a data frame of values from our set of environmental predictor variables. With this data frame, we ran four replicates of three modelling algorithms within the <em>SDM</em> package (Naimi & Araújo, 2016) in the R Programming Environment (v3.6.2; R Core Development Team, 2021) that included generalized linear, Maxent, and MARS models and resulted in a set of 12 models. After fitting these models, we used a K-fold cross-validation (30% of observations set aside for testing) along with threshold-dependent (TSS) and independent statistics such as area under the curve (AUC) and correlation coefficients (COR) for model evaluation. In presence-only studies, AUC measures the probability that a model will rank a randomly selected site where the species is present higher than a random background site, also known as pseudo-absence (Phillips et al. 2006; Phillips and Dudik, 2008). In contrast, COR values range from 0 – 1 and measure whether large differences in prediction values correspond to a higher presence probability for a species (Phillips and Dudik, 2008). We used the variable of importance percentage as a measure to evaluate the role of each predictor variable in explaining <em>C. c. roselaari</em> distribution. We then used the <em>predict</em> function for each model, with the output of this process being a raster of the current and future habitat suitability probability for <em>C. c. roselaari</em> (hereafter, ‘habitat suitability’). Habitat suitability values in this framework range from 0 (not suitable) to 1 (highly suitable) and can be used as a proxy for a species’ potential distribution. Finally, we created an ensemble of all 12 model predictions by stacking all raster models and calculating a weighted average across all model values. This ensemble represented the current habitat suitability for <em>C. c. roselaari</em>. To assess the habitat suitability for <em>C. c. roselaari</em> in the future, we used the SLR scenario for 2050 in combination with the current habitat suitability for <em>C. c. roselaari</em> to mask out those areas that are predicted to be inundated by 2050 and, thus, unavailable to <em>C. c. roselaari</em>. We then used this layer to train, process, and transfer a new model that would project future <em>C. c. roselaari</em> habitat suitability in 2050. <strong>2.8 Effects of sea level rise on wetland habitat and <em>C. c. roselaari</em> distribution. </strong> We conducted three analyses to determine whether our results supported our distributional contraction hypothesis and to identify areas with losses or gains in suitability. We defined a distributional contraction as a reduction in wetland extent by 2050 and first measured habitat contractions by calculating how much wetland area is at risk of permanent inundation. To do this, we subtracted areas from each wetland polygon predicted to be unavailable under a 2050 SLR scenario. Second, we mapped areas where <em>C. c. roselaari</em> habitat suitability was gained or lost by using a raster calculator that subtracted the future and current habitat suitability of <em>C. c. roselaari</em>. Third, we calculated a ‘highly suitable habitat’ (HSH) index for the current and future distribution of <em>C. c. roselaari</em> within each wetland. The HSH index consists of the median value of all pixels considered to be prime habitat (i.e., pixel values &gt;0.6). Fourth, we quantified the percentage loss and gain of habitat suitability by 2050 per wetland polygon by assessing how much suitability — on average — a given wetland would gain or lose under a 2050 SLR scenario. We then took the HSH index for 10 of the most important wetlands for <em>C. c. roselaari</em> (based on historic counts at the sites) and assessed whether habitat quality at these wetlands is likely to increase or decrease by 2050 compared to the current baseline. <strong> REFERENCES</strong> Carmona, R., N. Arce, V. Ayala-Perez, A. Hernández-Alvarez, J.B. Buchanan, L. J. Salzer, P. S. Tomkovich, J. A. Johnson, R. E. Gill, B. J. McCaffery, J. E. Lyons, L. J. Niles, and D. Newstead (2013). Red Knot Calidris canutus roselaari migration connectivity, abundance and non-breeding distribution along the Pacific coast of the Americas. Wader Study, 120(3), 168–180. Climate Central (2020). Land projected to be below tideline in 2050. https://www.climatecentral.org/sea-level-rise COSEWIC. (2007). COSEWIC Annual Report (p. 120). The Committee on the Status of Endangered Wildlife in Canada. https://wildlife-species.az.ec.gc.ca/species-risk-registry/virtual_sara/files//cosewic/car_COSEWIC_annual_report_2007_e.pdf Klingbeil, B. T., J.B. Cohen, M.D. Correll, C.R. Field, T.P. Hodgman, A. I. Kovach, E. E. Lentz, B.J. Olsen, W.G. Shriver, W. A. Wiest and C.S. Elphick (2021). High uncertainty over the future of tidal marsh birds under current sea-level rise projections. Biodiversity and Conservation, 30(2), 431–443. https://doi.org/10.1007/s10531-020-02098-z Kulp, S. A. and B. H. Strauss (2018). CoastalDEM: A global coastal digital elevation model improved from SRTM using a neural network. Remote Sensing of Environment, 206, 231–239. https://doi.org/10.1016/j.rse.2017.12.026 Lyons, J. E., W.L. Kendall, J. A. Royle, S. J. Converse, B. A. Andres and J.B. Buchanan (2016). Population size and stopover duration estimation using mark–resight data and Bayesian analysis of a superpopulation model. Biometrics, 72(1), 262–271. https://doi.org/10.1111/biom.12393 Muñoz-Salas, E. I., E. Palacios, L. Alfaro and M. E. Reiter (2023). Winter population trends and environmental drivers for three species of temperate shorebirds. Global Ecology and Conservation, 46, e02557. https://doi.org/10.1016/j.gecco.2023.e02557 Naimi, B. and M.B. Araújo (2016). Sdm: A reproducible and extensible R platform for species distribution modelling. Ecography, 39(4), 368–375. https://doi.org/10.1111/ecog.01881 R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Phillips, S. J., R. P. Anderson, and R. E. Schapire (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190(3-4), 231–259. https://doi.org/10.1016/j.ecolmodel.2005.03.026 Phillips, S. J. and M. Dudík (2008). Modeling of Species Distributions with Maxent: New Extensions and a Comprehensive Evaluation. Ecography, 31(2), 161–175. http://www.jstor.org/stable/30244563 SEMARNAT. (2010). NORMA Oficial Mexicana NOM-059-SEMARNAT-2010. https://www.profepa.gob.mx/innovaportal/file/435/1/NOM_059_SEMARNAT_2010.pdf Vidal-Zepeda, R. (2005). Las regiones climáticas de México. 213. https://doi.org/10.14350/rig.30041 INEGI. (2019). Uso de suelo y vegetación 2016. https://www.inegi.org.mx/temas/usosuelo/#descargas    
description: # Data to model potential and future distribution of red knots under a
 2050 SLR scenario
 [https://doi.org/10.5061/dryad.9cnp5hqsr](https://doi.org/10.5061/dryad.9cnp5hqsr) ## README The analysis are divided in three broad sections: **Species distribution modelling.** The R script to create the current and future habitat suitability of red knots is called “SpeciesDistributionModelling.R”. The location data of C. canutus roselaari from which the habitat suitability analysis were derived is called 1._Red_knot_locations.csv. To protect endangered species, coordinate decimal points were rounded up to two decimals. The environmental predictor rasters used to model the current and future habitat suitability under a sea level scenario are described in the table below: | **Environmental predictor rasters** | **Description** | **File name** | | :---------------------------------- | :------------------------------------------------------------ | :---------------------------------- | | Elevation | Global elevation in meters | Elev.tif and Elev2050.tif | | Temperature | Mean December and January temperature | Temp.tif and Temp2050.tif | | Distance to wetland | Euclidean distance to the nearest polygon of coastal wetlands | DistWet.tif and Distwet2050.tif | | Distance to urban areas | Euclidean distance to the nearest polygon of  urban areas | DistUrban.tif and DistUrban2050.tif | | Land cover and vegetation | MCD12Q1 MODIS/Terra+Aqua Land Cover Global 500m | Veg.tif and Veg2050.tif | | Nightlights | Proxy for disturbance and development intensity | NightL.tif and NightL2050.tif | | Silt | Global soil silt content 0-5 cm depth | Soilsilt.tif and Soilsilt2050.tif | | Clay | Global soil clay content 0-5 cm depth | Soilclay.tif and Soilclay2050.tif | | Soil Ph | Global soil pH x 10 in H20 0-5 cm depth | Soilph.tif and Soilph2050.tif | | Soil Coarseness | Global coarse fragments 0-5 cm depth | SoilCoar.tif and SoilCoar2050.tif | The final outputs of this analysis are described below: * A model ensemble of the current Red Knot Distribution. The values of this raster file represent habitat suitability ranging from 0 (not suitable) to 1 (very suitable). * A model ensemble Ensemble of the Future 2050 Red Knot Distribution. The values of this raster file represent habitat suitability ranging from 0 (not suitable) to 1 (very suitable). * The variables of importance for the Current and Future distribution to create Figure 3C. **Data handling in ArcGis Pro** The model ensembles rasters created using the species distribution modelling were loaded into ArcGis pro. We used the "raster to polygon tool" to convert these rasters to polygon. Then, we performed an “Intersection” between ensembles of red knot distribution and the wetland inventory (file WetlandInventory2.shp). This allowed us to isolated only the pixels found within the coastal wetland polygons for the current and future distribution of Red knots. After this process, each pixel polygon and group of pixel polygons within a wetland would have a unique identifier, an area calculation in hectares and a associated habitat suitability. Then, we extracted the centroid of each coastal wetland and obtained its latitude and longitude. We appended the wetland centroid coordinates to each pixel polygon within each wetland. Finally, we categorized each pixel according to the region (Pacific, Gulf or Mainland) where they were found. The resulting file is called: WetlandSuitabilityTodayFuture.csv The structure of WetlandSuitabilityTodayFuture.csv is explained below: | **Column** | **Description** | | :----------- | :------------------------------------------------------------- | | OID\_ | Unique identifier made during the intersection | | ID | Unique ID of each wetland | | ID\_region | Unique ID of each wetland within a region | | Region | Region where the wetland is located | | TodayArea | The current extent of the wetland | | FutureArea | The future extent of the wetland under a SLR scenario | | HabLoss | The amount of habitat potentially lost to SLR by 2050 | | HabRemaining | The amount of habitat potentially remaining due to SLR by 2050 | | Name | The name of the wetland (incomplete) | | x | Longitude centroid coordinate | | y | Latitude centroid coordinate | **Data analysis and figures in R** * The following figures use the R script DataAnalysisAndFigures.R and the datasets “WetlandSuitabilityTodayFuture.csv” and “TodayandFutureSuitabilityV2.csv” * Figure 2 uses “WetlandSuitabilityTodayFuture.csv” * Figure 4 requires to merge, the WetlandSuitabilityTodayFuture.csv with TodayandFutureSuitabilityV2.csv—a dataset that includes the raw habitat suitability values of each single pixel on a given wetland from the current and 2050 scenario. The merging process is included in the script. *  Figure 3 is a multi-panel. Panel A and B are the ensembles from the species distribution modelling. Panel C is made using the variables of importance. Panel D is a raster subtraction between the ensembles. The code is found in the script. * Figure 5 is made by filtering dataset by the 10 priority coastal wetland in NW Mexico for red knots and by habitat suitability &gt; 0.6. The median habitat suitability is then calculated for each of these wetlands now and in the future. The file structure for WetlandSuitabilityTodayFuture.csv is explained below: | **Column** | **Description** | | :---------- | :--------------------------------------------------------------------- | | OBJECT ID | Unique identifier created in ArcGis for each row | | ID | Unique identifier for each coastal wetland | | ID2 | Identifier for each pixel polygon | | Name | Name of each coastal wetland (incomplete) | | Scenario | Today (current habitat suitability); Future (2050 habitat suitability) | | Suitability | Values ranging from 0 (not suitable) to 1 (highly suitable) | | PixelHa | Pixel size in hectares | **Usage notes** R is required to run the following scripts: * SpeciesDistributionModelling.R * DataAnalysisAndFigures.R ArcGis Pro or other GIS software is recommended but not needed to visualize and/or perform intersections. Any GIS software and R can be used to open, visualize and analyze .tif files and shapefiles. WetlandInventory2 file requires three files with the extension .shp, shx and .dbf located in the same workspace or folder to be opened and visualized using GIS. Excel or R can be used to visualize .csv files