Projecting global mangrove species and community distributions under climate change

Given the multitude of ecosystem services provided by mangroves, it is important to understand their potential responses to global climate change. Extensive reviews of the literature and manipulative experiments suggest that mangroves will be impacted by climate change, but few studies have tested these predictions over large scales using statistical models. We provide the first example of applying species and community distribution models (SDMs and CDMs, respectively) to coastal mangroves worldwide. Species distributions were modeled as ensemble forecasts using BIOMOD. Distributions of mangrove communities with high species richness were modeled in three ways: as the sum of the separate SDM outputs, as binary hotspots (with >3 species) using a generalized linear model, and continuously using a general boosted model. Individual SDMs were projected for 12 species with sufficient data and CDMs were projected for 30 species into 2080 using global climate model outputs and a range of sea-level rise projections. Species projected to shift their ranges polewards by at least 2 degrees of latitude consistently experience a decrease in the amount of suitable coastal area available to them. Central America and the Caribbean are forecast to lose more mangrove species than other parts of the world. We found that the extent and grain size, at which continuous CDM outputs are examined, independent of the grain size at which the models operate, can dramatically influence the number of pseudo-absences needed for optimal parameterization. The SDMs and CDMs presented here provide a first approximation of how mangroves will respond to climate change given simple correlative relationships between occurrence records and environmental data. Additional, precise georeferenced data on mangrove localities and concerted efforts to collect data on ecological processes across large-scale climatic gradients will enable future research to improve upon these correlative models.


Introduction 37
Sea-level rise and altered weather patterns resulting from global climate change have 38 impacted and will continue to impact coastal systems, altering the ecological and economic 39 services that they offer (Nicholls et al. 2007). In coastal tropical and sub-tropical areas and Farnsworth 1997, Ye et al. 2003). These reviews and experiments suggest that individual 59 5 mangrove species' distributions may contract and local species richness and productivity may 60 decrease in regions where climate-change scenarios forecast that precipitation and run-off will 61 decrease while salinity soil sulfides increase (Snedaker 1995, Ellison 1994. In contrast, where 62 precipitation and run-off increase, upland nutrients will be deposited, salinity will be reduced, 63 and acid-sulfide soils will be moderated, leading to increased productivity, opportunities for 64 range expansion of individual mangrove species, and potential for increases in local species have 50% more (most of Southeast Asia). Although we recognize that there also is variation 151 among GCMs (IPCC 2007), it was beyond the scope of this study to run different GCMs on the 152 SRES A1b scenario. 153 As we did for the mangrove occurrence data, we assigned to each coastal grid cell the 154 nearest value (within a 40-km radius) of each of the current and future environmental variables. 155 To account for possible spatial error in the river discharge layer to coastal cells, this layer was 156 first resampled at a 14 km grid size, taking the maximum value within that larger region before 157 assigning values to the coastal cells. to determine the oceans in which the species occur, and then set projected probabilities to zero at 228 longitudes beyond these regions (Table 4). After selecting crop lines for each species, we 229 examined global projected distributions to ensure that the crop lines did not intersect areas 230 predicted to have continuous occurrences. Thus, summary statistics of model outputs should not 231 be very sensitive to the location of crop lines. 232

Community distribution models and species richness 233
Mangroves tend to occur in association with multiple mangrove species, each of which 234 may occur at specific tidal elevations (Macnae 1968). At the coarse scale of this study, we are 235 interested primarily in identifying areas where multi-species mangrove assemblages are likely to 236 occur, rather than distinguishing between different types of mangrove communities. We modeled 237 local species richness ("alpha diversity") because we had inadequate data to model species 238 turnover ("beta diversity"). 239 We modeled mangrove species richness using three different approaches: a composite 240 model, a continuous-response model, and a binary-response model. For the composite model, we 241 13 combined the independent projections of the 12 individual SDMs by summing the predicted 242 occurrences within each coastal cell. For the continuous and binary models, we calculated the 243 current species richness within each coastal cell based upon all 30 major mangrove species in our 244 GBIF data set (Table 1). 245 In the binary model, we sought to identify those cells where multi-species mangrove 246 communities are most likely to exist. To do this, we assigned each cell with three or more species 247 out of the 30 total species in our GBIF data a value of one and each cell with less than three 248 species was assigned a value of zero. In this analysis, we modeled presence of cells with high 249 species richness relative to the other cells in our data set. This process yielded 355 presences of 250 high richness cells. We used three species as the threshold because this was the highest value that 251 would yield enough presences of these high richness cells for sufficient predictor-to-response 252 ratios in the models. The presence of three species may not indicate a true hotspot of mangrove 253 diversity in the field. However, this threshold is appropriate within the context of the GBIF data 254 set, and allows us to confidently weed out cells where only one or two mangrove species exist. 255 We considered using different thresholds for defining high richness in the eastern and western 256 hemispheres, because one might expect greater overall richness in the eastern hemisphere. 257 However, we only see more high-richness cells in the east when the threshold is set at four or 258 five species per cell (Fig. 1), at which levels there are insufficient sample sizes. We further felt 259 that it was more appropriate to treat all of the data uniformly in the model, rather than imposing 260 further rules that may introduce more potential for bias. 261 We ran the binary richness data through the same BIOMOD modeling process that we did 262 for each of the individual species. For weights in the binary model, we used the actual number of 263 species observed in each cell (Fig. 1). 264 For the continuous CDMs, we did not have access to a comprehensive software package 265 for ensemble distribution model selection and prediction based on non-binary data (e.g., 266 BIOMOD does not model abundance). Instead, we fit the full GBM and GLM models using all 267 21 environmental variables as predictors and the number of mangrove species reported within 268 each grid cell as the response. We compared models with the full suite of predictor variables to 269 those fit using subsets of variables: the five variables with the greatest influence; or by iteratively 270 discarding the least influential variable between pairs of variables with greater than 0. univariate GLMs comparing observed species richness in the holdout data sets to predicted 282 species richness. Because our ultimate aim was to examine large-scale patterns in mangrove 283 species diversity, we also tested predictive performance of the full GBM and full GLM models at 284 a coarser resolution. In the coarse-resolution tests, we aggregated the predicted and observed 285 data in the holdout regions to a 500-km grid cell size before comparing predicted and observed 286 species densities. For the final selected model, we fit the subset of predictor variables to the 287 entire world, and then projected forward using the environmental variables in the 2080 3m sea 288 level rise scenario because the results of the SDMs we ran previously were not sensitive to the 289 different sea-level rise scenarios. 290 For the composite, binary, and continuous CDMs, we generated 500-km grid cell maps of 291 forecasted change in species richness between current conditions and future scenarios. We also 292 calculated means of latitude in each cell weighted by the fitted species richness in current and 293 future scenarios for the three models. The GBM model with the full suite of variables had the 294 best predictive performance in most scenarios (Table 5), and so we used this model for our future 295 projections. As with the SDMs, model evaluation with holdout data suggested that models 296 trained with the least pseudo-absences had the best predictive performance when tested against 297 the data with the original ~4 km (i.e., 2.5 minute) grid size. However, coarse scale maps 298 produced by these models exhibited many nonsensical predictions for current mangrove 299 occurrences, including high species richness in high latitude regions. When examining 300 predictions that had first been re-scaled to a 500 km grid size, inclusion of pseudo-absences 301 improved model likelihoods, and produced maps of current fitted distributions that better 302 matched our expectations. Because our study is focused on global changes in mangrove 303 distributions, we opted for including 2000 pseudo-absences in the final model. This yielded an 304 approximate presence to absence ratio of 1:1, similar to that used in the individual SDMs with 305 500 pseudo-absences. Code for all SDMs and CDMS performed using R statistical software 306 version 14.0 are included in Supplementary Material. 307

Evaluation of SDM and CDM outputs 308
We evaluated model outputs by generating summary maps at a coarser resolution in order 309 to generalize patterns across regions. We generated these maps with 500 km grid cells and 1000 310 km grid cells. Within each of the larger cells, we summed the predicted species richness in all of 311 the 4 km grid cells. The result is a mangrove species density value for each of the measured cells. 312 This density is different from the mean species richness, because it incorporates both species 313 richness and the number of occupied cells. A 500-km cell centered on Panama has much more 314 coastline than a 500-km cell centered on the coast of Peru. Thus, even if every 4 km coastal cell 315 had the same number of species, the species density measured in the 500 km grid cells would be 316 higher in Panama than Peru. 317

Coastline versus Latitude 318
Our study analyzes latitudinal shifts in coastal species. To frame our results, we also 319 needed to understand how the world's coastlines are distributed with respect to latitude. To this 320 end, we summed the total number of ~4 km grid cells within each 2-degree latitudinal bin. We 321 also performed a separate analysis using ArcMap in which we compared the total length of our 322 coastline vector data within 15° of the equator, and between 15° and 30° from the equator. The 323 vector data was generated at a 1.7 km resolution. 324

Results 325
Species distribution models 326 The current distribution of each of the most common 12 mangrove species was best 327 predicted by a different set of five environmental variables (Tables 2 and 3); precipitation in the 328 warmest and coldest quarters appeared in the list of top five predictors for more than half of the 329 mangrove species. In the variable selection process, river discharge and horizontal tide were 330 identified as important environmental predictors only for Rhizophora apiculata, R. racemosa, 331 and R. stylosa (Tables 2 and 3). The predictive performance of the models was high: TSS values 332 for the twelve species averaged 0.97 (range 0.950 -0.988), but in a few instances the SDMs 333 predicted current mangrove distributions outside of their current known latitudinal range (Fig. 2). 334 Rather than focusing only on minimum and maximum latitudes, we therefore also examined the 335 mean and standard deviations of the absolute values of latitude. 336 All 12 common mangrove species were forecast to change their absolute mean latitude 337 and total suitable coastal area relative to current climatic conditions (Fig. 2). Half of the modeled 338 species were projected to have a poleward shift of two degrees of latitude or more in the absolute 339 mean latitudes of their distributions under the future climate scenario (Fig. 2). These six species 340 also were forecast to suffer losses in the total area of suitable coastal habitat available within 341 their expanded ranges (Fig. 2). This loss of the amount of suitable coastal habitat available for 342 species with poleward range shifts could be due to the lower amount of total coastline in higher 343 tropical latitudes compared to equatorial areas (Fig. 3). All of the species that did not experience 344 a poleward shift in the absolute mean values of their distributions gained total suitable coastal 345 habitat under the future scenario regardless of the amount of sea-level rise. 346 The four species with current ranges limited to the Americas, western and central Africa, 347 and the western Pacific islands -Avicennia germinans, Laguncularia racemosa, Rhizophora 348 mangle, Rhizophora racemosa -were all forecast to experience overall losses in total suitable 349 coastal habitat and poleward shifts under the future climate scenario compared to current 350 climatic conditions (Figs. 2, and 4-27). The NCAR-CCSM3 GCM forecasts that the annual 351 precipitation in these regions will decrease by at least 50% and that annual temperature will 352 increase by at least 2°C. Our forecasts of mangrove loss in these areas supports previous 353 hypotheses that individual mangrove species' distributions will contract and richness will decline 354 as rainfall and runoff decrease while salinity and extent of acid-sulfide soils increase (Snedaker 355 1995, Ellison 1994). 356

18
The remaining eight species, with current ranges limited to eastern Africa, Asia, and 357 Australia, had more variable forecasts. Lumnitzera littorea and Rhizophora mucronata were 358 projected to shift polewards and lose suitable coastal habitat, while Avicennia marina, Ceriops 359 tagal, Lumnitzera racemosa, and Rhizophora apiculata were forecasted to gain potential coastal 360 area with absolute mean latitudinal gains of less than two degrees. Sonneratia alba and 361 Rhizophora stylosa were projected to gain coastal habitat and experience decreases in absolute 362 mean latitude (i.e., equatorial range contractions). With forecasted gains in suitable coastal area 363 of 260% to 290% of its current projected distribution, R. stylosa was forecast to gain a 364 remarkable 110 to 185% additional habitat relative to its current distribution. 365

Community distribution models and species richness 366
The means of the absolute value of latitude weighted by fitted current species density 367 were 14.5º, 14.3º and 17.0º for the composite model, the binary model, and the continuous 368 model, respectively. The projected mean latitudes for the 3m sea-level rise were 14.6º, 14.2º, and 369 15.7º for the same three models. The projected maps of change in species density differed 370 between the three model types, although there were a few areas of overlap (Fig. 28). All three 371 models predicted gains in mangrove species density across much of southeastern Asia, southern 372 Brazil, northern Chile, eastern Australia, southeastern Africa, parts of northern Africa, and parts 373 of northwestern Mexico. All three models also predicted losses of mangrove species density in 374 the Caribbean Islands, parts of Central America and parts of northern Australia (Fig. 28). 375

Coastline versus latitude 376
In summing the ~4 km coastal cells vs. latitude, we found that the total length of coastline 377 between the equator and ±15° was 42% greater than the length of coastline between 15° and 378 30°N or S (i.e., 182,000 km versus 129,000 km, respectively; Fig. 3). The vector analysis 379 similarly showed 43% more coastline within 15° of the equator than between 15° and 30° from 380 the equator. dataset, suggesting that our forecasts for these species are more robust than those for species with 411 sparser occurrence records, such as many species in the Indo-West Pacific (Table 1) For researchers hoping to advance techniques for distribution models based on 425 continuous data, our model selection process offers a further lesson in considering spatial scale. 426 We found that the extent and grain size at which continuous model outputs are examined, 427 independent of the grain size at which the models operate, can dramatically influence the number 428 of pseudo-absences needed for optimal parameterization. That small scales are best modeled 429 without pseudo-absences, but large-scale models are benefited by pseudo-absence is somewhat 430 intuitive. Without pseudo-absences, the models evaluate finer scale differences within sites 431 occupied by mangroves, whereas with many pseudo-absences, the models can better evaluate the 432 coarser scale differences between areas with and without mangroves. This issue should only 433 apply to continuous data where all presences are not identical, unlike in binary data. 434 The SDMs and CDMs presented here provide a first approximation of how mangroves 435 will respond to climate change given simple correlative relationships between occurrence records