Drivers of global mangrove loss and gain in social-ecological systems
Mangrove cover change variables. We used the Global Mangrove Watch (GMW) v2.0 dataset from 1996 to 201656 to calculate four response variables across landscape mangrove geomorphic units24 over two time periods, 1996–2007 and 2007–2016: (1) percent net loss (units that had a net change in mangrove cover of 0), (3) percent gross loss (units that had a decrease in mangrove cover, not accounting for any increase), and (4) percent gross gain (units that had an increase in mangrove cover, not accounting for any decrease). Percent variables were calculated relative to the area at the start of the time period and were log transformed to meet the assumptions of the statistical models. We initially also considered 5 primary response variables (Supplementary Table 3), including net change in mangrove area ranging from negative (loss) to zero (no change) to positive (gain), however, the data did not meet model assumptions of equal variance (Supplementary Table 9). It was therefore necessary to separate areas of net loss and net gain and areas of gross loss and gross gain to remove zeros and log-transform to achieve normal distribution. Area of mangrove change was correlated with size of the mangrove geomorphic unit (higher area of mangrove loss or gain in bigger units), therefore we included geomorphic unit size as an explanatory variable in the models with primary response variables. We selected the transformations of these primary variables – percent net loss, percent net gain, percent gross loss, and percent gross gain to include in the analysis, because the percent changes control for differences in relative sizes of geomorphic units and because net change alone can underestimate the extent of change57.Examining mangrove change across geomorphic settings is likely to be relevant to socioeconomic and environmental conditions. Mangroves occur in the intertidal zone in diverse coastal geomorphic settings (e.g., deltas, estuaries, lagoons) shaped by rivers, tides, and waves58,59. The distribution, structure, and productivity of mangroves varies spatially with regional climate and local geomorphological processes (e.g., river discharge, tidal range, hydroperiod, and wave activity) that control soil biogeochemistry60,61,62,63. These geomorphic settings are defined by natural landscape boundaries (e.g., catchments/bays) which also often delineate boundaries of human settlements. A global mangrove biophysical typology v2.2 dataset64 was used for the delineation of landscape mangrove geomorphic units, which used a composite of the GMW dataset from the 1996, 2007, 2010, and 2016 timesteps to classify the maximal extent of mangrove cover into 4394 units (classified as delta, estuarine, lagoon or open coast). The mangrove geomorphic units do not include non-mangrove patches, unless they have been lost from the unit over time. The mean size of geomorphic units was 33.63 ha. Some splits of geomorphic units were undertaken to reduce size and divide by country boundaries. The four largest deltas (northern Brazil Delta ID 70000, Sundarbans Delta ID 70004, Niger Delta ID 70009, and Papua coast Delta ID 70013) were split into 4, 5, 4, and 2 units, respectively to aid with data processing. Mangrove geomorphic units that overlapped two countries (Peru/Ecuador, Singapore/Malaysia, and Papua New Guinea/Australia) were split by the national boundary.The country governing each geomorphic unit was assigned to match national-level variables to geomorphic units. To capture mangroves that are mapped outside of country coastline boundaries, we did a union of the GADM country shapefile v3.665 and the Exclusive Economic Zones (EEZs) v1166. The following manual country designations were made to resolve overlapping claims in the EEZs: (1) Hong Kong was merged with China as Hong Kong does not have a mapped EEZ; (2) The overlapping claim of Sudan/Egypt was maintained as a joint Sudan/Egypt designation, as this is an area of disputed land called the Halayib Triangle. However, for this study, mangrove units within this area were assigned to Egypt because Egypt currently has military control over the area; (3) Mayotte (claimed by France and Comoros) was assigned to Mayotte as it is a separate overseas territory of France recognised in GADM that has different socioeconomic variables; (4) The protected zone established under the Torres Strait Treaty was assigned to Australia as these islands are Australian territory.Areas of mangrove cover in 1996, 2007, and 2016, and gross losses and gains in each geomorphic units over the two time periods were assessed in ArcMap 10.867. Percent losses and gains were calculated in R 4.0.268. In using the GMW mapping, a minimum mapping unit of 1 ha is recommended for reliable results5, therefore we removed all geomorphic units less than 1 ha from the analysis, which reduced the available sample size from 4394 across 108 countries to 4235 units across 108 countries. In calculating percent net gains, 11 and 12 of the units returned infinity values for 1996–2007 and 2007–2016, respectively, because there was no initial mangrove cover. In these instances, 100% gain was assigned to these units.Socioeconomic variables (Supplementary Table 4)Economic growthPrevious global analyses of mangroves have been limited by data availability on economic activity to national metrics, such as a country’s Gross Domestic Product (GDP)12,18. Night-time lights satellite data provide local measures of economic activity that are comparable through time and available globally9,69. The data improve estimates of GDP in low to middle income countries69 and are strongly correlated with local indicators of human development70 and electricity consumption and GDP at the national-level71. We used the Night-time Lights Time Series v472 stable lights data, where transient lights that are deemed ephemeral, e.g., fires, have been filtered out and non-lit areas set to zero73, choosing the newer satellites where applicable70. As a proxy for local economic growth, we calculated the change in annual average stable lights within a 100 km buffer of the centroid of each geomorphic unit from 1996 to 2007 and 2007 to 2013 (no data available past 2013) using the ‘raster’ package in R74. The 100 km buffer was chosen to account for pressures from human activity within and surrounding the mangrove area, and to avoid bias with larger spatial units70.Market accessibilityTravel time to the nearest major market (national or provincial capital, landmark city, or major population centre) has been shown to be a stronger predictor of fish biomass on coral reefs than population density or linear distance to markets27. We used the global map of travel time to cities for 201575 to estimate the average travel time from each geomorphic unit to the nearest city via surface transport using the ‘raster’ package in R74, as an indicator of access to markets to trade commodities (e.g., rice, shrimp, palm oil).Economic complexityPrevious studies have examined the effect of GDP on mangrove change18, however, this is a blunt measure of country capability. Measuring a country’s economic complexity, that is the diversified capability of a nation’s economy, is preferable. For example, a country with high GDP but low economic complexity can be prone to regulatory capture by high-value natural resource industries and resource corruption26. Therefore, we used the Economic Complexity Index (ECI)76 for countries as an indicator of regulatory independence. The ECI had better coverage of countries in later years (Supplementary Table 4), therefore the ECI for the end of the time periods was used (2007 and 2016), although we recognise this may reduce the detection of trends because of potential time lags in impacts.DemocracyWe used the Varieties of Democracy (VDEM) index v10 which measures a country’s degree of freedom of association, clean elections, freedom of expression, elected executives, and suffrage77, and has been indicated to influence NDC ambition in countries to address climate change78. We adopted the VDEM index for the start of the time periods (1996 and 2007) to account for potential time lags in impacts.Community forestry supportWe determined the extent that community forestry (CF) is implemented across countries through a systematic review of articles returned in the Web of Science database (Core collection; Thomson Reuters, New York, U.S.A.). We used the search terms: TS = (“community forestry” OR “community-based forestry” OR “social forestry”) AND (TI = ”country” OR AB = ”country”) to identify how many CF case studies were reported in each country, and whether any were in mangroves. As scientific literature is biased towards particular regions, we also reviewed relevant FAO global studies79,80,81 and online databases (ICCA registry82 and REDD projects database83) to identify additional case studies (Supplementary Fig. 5). We then generated scores of 0–3 for each country based on summing values assessed using these criteria: +1 (1–50 CF case studies); +2 ( >50 CF case studies); +1 (CF case study in mangroves). There may have been some double counting as we counted the number of case studies in each article, and we will have missed CF projects not published or communicated in English. However, this is likely to have had a limited impact on the scoring method.Indigenous landThe proportion of Indigenous peoples’ land versus other land per country was calculated from national-level data84. Whilst this study involved Indigenous peoples’ land mapping at a global scale, the spatial data was not published, and thus we could only evaluate the influence of Indigenous land at the national level rather than local level.Restoration effortThe number of mangrove restoration sites per country was calculated from combining two datasets collated by C. Lovelock (2020) and Y.M. Gatt and T.A. Worthington (2020) identifying mangrove restoration project locations from web searches in English and for scientific and grey literature using Google Scholar. Duplications were removed and the number of sites was used as an indicator of effort. This will underrepresent effort in countries with few, large sites, and where restoration projects are not published or communicated in English.Climate commitmentsThe Paris Agreement is a global programme for countries to commit to climate action by submitting Nationally Determined Contributions (NDCs) to the United Nations Framework for the Convention of Climate Change (UNFCCC). First, we reviewed NDCs for mangrove-holding nations from the NDC Registry85 submitted as of 07/01/2021 to determine the extent that mangroves or coastal ecosystems were included in national climate policy (scoring method in Supplementary Table 4). We hypothesised that countries with mangrove or coastal ecosystem NDCs may be more likely to promote mangrove conservation or restoration. While the first NDCs were submitted around 2015, at the end of our time series, we suspected higher commitments would point towards a stronger baseline in environmental governance. Most countries submitted updated or second NDCs during 2021 however these were not considered relevant to the time periods assessed. Google Translate was used to interpret NDCs in languages other than English.Ramsar wetlandsThe ecological character of Ramsar wetlands have been found to be significantly better than those of wetlands generally86. The area of Ramsar coastal and marine wetlands from the Ramsar Sites Information Service87 was calculated per country. Thirty-eight mangrove-holding countries are not signatories to the Ramsar Convention, and these countries were assigned a value of 0. The area of Ramsar wetlands per country was scaled by dividing by the country’s area of mangroves in 1996.Environmental governanceWe assessed the Environmental Performance Index (EPI)88 as an indicator of a country’s effectiveness in environmental governance. The biodiversity and habitat (BDH) issue category assesses countries’ actions toward retaining natural ecosystems and protecting the full range of biodiversity within their borders. We took the BDH score for 2020 for the 2007–2016 time period and the BDH score for 2010 for the 1996–2007 time period (calculated by subtracting the ten-year change from BDH 2020). However, due to collinearity with other variables this index was excluded from the analysis (see statistical analysis).Protected area managementWe also assessed Marine Protected Area (MPA) staff capacity as an indicator of the effectiveness of management of protected areas for countries. We used published global marine protected area (MPA) management data14 which is based on the Management Effectiveness Tracking Tool (METT), the World Bank MPA Score Card, and the NOAA Coral Reef Conservation Programme’s MPA Management Assessment Checklist. Adequate staff capacity was the most important factor in explaining fish responses to MPA management globally, followed by budget capacity, but they were significantly correlated14. We, therefore, calculated the mean staff capacity across MPAs per country as our indicator. Mangroves can be included in terrestrial protected areas, which are not represented in this dataset, however, this measure provides an indicator of national governance of protected areas. However, due to collinearity with other variables this indicator was excluded from the analysis (see statistical analysis). The extent of protected areas was not included in the analysis because it has already been found to influence mangrove loss18.Biophysical variables (Supplementary Table 5)Coastal geomorphic typeMangrove extent change likely varies among different coastal geomorphic settings because human activities or environmental changes occur more commonly in some geomorphic settings than others. For example, losses of lagoonal mangroves were nearly twice as large as those in other geomorphic types24. Landscape geomorphic units from the global mangrove typology dataset v2.264 were classified as delta, estuary, lagoon or open coast.Sediment availabilityMangrove expansion and retreat are driven by sediment deposition and erosion, which are influenced by sediment availability from rivers and wave action, and alterations in hydrodynamic regimes47,89. We used the sediment trapping index from the global free-flowing rivers (FFR) dataset90 to indicate sediment availability from rivers within different geomorphic units. A mangrove catchment dataset was created based on the HydroSHEDS database91. River networks that intersected with mangrove geomorphic units were linked to that unit’s ID. Where rivers intersected multiple units, they were manually assigned by visual inspection. River basins that intersected either with the geomorphic units directly or the river networks were also linked to that unit’s ID. The FFR dataset90 was then spatially joined to the mangrove catchment dataset to identify the most downstream (i.e., the coastal outlet) segment of each FFR and its associated sediment trapping index. Not all geomorphic units (n = 3475) were linked to an FFR, however, an individual unit could be linked with several FFRs. Therefore, the unit sediment trapping index was the weighted mean of the river values, with weighting based on each FFR’s average long-term (1971–2000) naturalised discharge (m3s−1), with discharge set to the minimum value for segments with zero flow. Geomorphic units without connecting FFRs were given an index of zero (no sediment trapping). The sediment trapping index represents the percentage of the potential sediment load trapped by anthropogenic barriers along the river section. The focus on river barriers may obscure larger scale oceanic patterns that influence mangrove losses and gains (e.g., movement of mud banks from the Amazon River over 1000’s of kilometres92) or increases in sediment that could be coming from soils with catchment deforestation and erosion.Habitat fragmentationMany countries with high mangrove loss have been associated with elevated fragmentation of mangrove forests, although the relationship is not consistent at the global scale93. We calculated the clumpiness index of mangrove patches within geomorphic units within each time period, as this habitat fragmentation metric is independent of areal extent93. Whilst habitat fragmentation can be human-driven, clumpiness measures the patchy distribution of mangroves, which can also be due to natural factors inducing edge effects. We used a similar approach to Bryan-Brown, et al.86 to quantify the clumpiness index. The ‘landscape’ was defined as the combined extent of the mangrove geomorphic units across four timesteps (1996, 2007, 2010, and 2016) from the GMW dataset56. For the three focal years in this study (1996, 2007, and 2016) each geomorphic unit (n = 4394) was converted into a two-class polygon, where class one represented mangroves present during that time step and class two mangroves present in the other time steps (i.e., areas of mangrove loss). The polygons were transformed to a projected coordinate system (World Cylindrical Equal Area) and converted to rasters with a resolution of 25 m. Each raster was imported into R version 3.6.394, with clumpiness calculated using the package ‘landscapemetrics’ v1.5.095.Clumpiness describes how patches are dispersed across the landscape and ranges between minus one, where patches are maximally disaggregated, to one, where patches are maximally aggregated, a value of zero represents a case whereby patches are randomly distributed across the landscape. The clumpiness index requires that both classes are present in the landscape, therefore a no data value (NA) was returned for units where no loss of mangroves had occurred, or where there was 100% gain of mangroves in a later time period. The number of directions in which patches were connected was set to eight. The following manual fixes were conducted for NA values returned: 1) Where NA was returned for units where no loss of mangroves had occurred in another time period, i.e., class 1 (mangrove present) = 1 and class 2 (mangrove loss) = 0, assume +1 (maximally clumped); and 2) Where NA was returned for units where there was 100% gain of mangroves in a later time period, i.e., class 1 (mangrove present) = 0, class 2 (mangrove present) = 1 (100% gain), assume −1 (maximally disaggregated).Tidal amplitudeIn settings of low tidal range, mangrove vertical accretion is less likely to keep pace with rapid sea level rise3. However, in settings of high tidal range, mangroves may be more extensive and vulnerable to conversion to aquaculture or agriculture because of larger tidal flat extents. The Finite Element Solution global tide model (FES2014)96 is considered one of the most accurate tide models for shallow coastal areas97 and was selected to estimate the mean tidal amplitude within each geomorphic unit using the principal lunar semi-diurnal or M2 tidal amplitude as this is this most dominant tidal constituent98. To account for potential variation in the tidal amplitude across large geomorphic units, the raster pixel value for M2 tidal amplitude96 closest to the centroid of each mangrove patch within each unit was calculated, with the smallest value set at 0.01 m. For each geomorphic unit, the tidal amplitude was calculated as the weighted mean of the patch values, with weighting based on the patch area relative to the total unit area.Antecedent sea-level riseThe distribution of mangroves on shorelines changes over time with sediment accretion, erosion, subsidence, and sea-level rise (SLR)99, and periods of low sea level can cause mangrove dieback100. We used regional mean sea-level trends between January 1993 and December 2015 from the global sea level Essential Climate Variable (ECV) product v.2101,102 to estimate the mean antecedent SLR for each geomorphic unit. Spatial variation in regional sea-level trends generally range between −5 and +5 mm yr−1 (global mean of 3 mm yr−1)13. Extreme values ( >5 mm yr−1) observed in the dataset are subject to high levels of uncertainty (Sea Level CCI team, pers. comm.), and were therefore truncated to 5 mm yr−1. The raster pixel value for SLR102 closest to the centroid of each mangrove patch within each geomorphic unit was calculated. The geomorphic unit antecedent SLR values was calculated as the weighted mean of the patch values within the unit.DroughtWhilst long-term precipitation and temperature influence mangrove distribution globally62, periods of low rainfall have been reported to cause extensive mangrove dieback at regional scales, particularly when combined with high temperatures and low sea levels103. We used the Standardized Precipitation-Evapotranspiration Index (SPEI) from the global SPEI database v.2.6104 as an index of drought severity. SPEI is derived from precipitation and temperature and is considered an improved drought index that allows spatial and temporal comparability105,106. The mean SPEI raster pixel value was calculated for each time period and then averaged across the geomorphic units using the ‘ncdf4’107 and ‘raster’ packages74 in R.Tropical storm frequencyLarge-scale destruction of mangroves across regions have been reported from strong winds, high energy waves, and storm surges associated with tropical storms108. We used the International Best Track Archive for Climate Stewardship (IBTrACS) dataset since 1980 v4109 to calculate the number of tropical cyclone occurrences (points along their paths) within a 200 km buffer of the centroid of geomorphic units within each time period using the sf package110 in R. Maximum wind velocity and surface pressures are likely experienced within 100 km of a cyclone’s eye111, therefore the 200 km buffer zone was selected to cover the average size of geomorphic units (33.63 ha), and all tropical storms potentially influencing mangrove growth. Whilst tropical storms affect only 42% of the world’s mangroves60, they are likely to be important stressors within cyclone-impacted countries.Minimum temperatureExtreme low temperature events were a driver of mangrove loss in subtropical regions, such as Florida and Louisianan of the US, and China28,112. We used the WorldClim bioclimatic variable 6 (minimum temperature of the coldest month averaged for the years 1970–2000)113 to calculate the mean minimum temperature across the geomorphic units using the ‘sf’110 and ‘raster’ packages74 in R. Where NAs were returned due to no overlapping raster layer, the value of the closest raster pixel to the centroid of the geomorphic unit was assigned.Statistical analysisWe used multi-level linear modelling to investigate relationships between mangrove cover change variables and socioeconomic and biophysical variables to consider landscape (level 1) and country (level 2) predictors in a hierarchical approach114. For each response variable, we modelled the response for 1996–2007 and 2007–2016, using explanatory variables specific to the time-period where available. Data inspection revealed that high percent loss or gain was concentrated in small geomorphic units, therefore to avoid bias in our results, we removed geomorphic units less than 100 ha from the analysis, which further reduced the available sample size to 3134 units across 95 countries. Statistical analysis was undertaken in R 4.0.268.The response variables were log-transformed to fit normal distribution. We tested for collinearity between our explanatory variables using Pearson’s correlation coefficient (r > 0.5) (Supplementary Tables 6 and 7). MPA staff capacity and EPI were excluded from our models because MPA staff capacity was correlated with ECI 2007 and ECI 2016 (both r = 0.54), and EPI 2020 was correlated with VDEM 2016 (r = 0.63). To improve model fit, travel time to the nearest city, mangrove restoration effort and Ramsar wetland area (relative) were log+1-transformed, and tidal amplitude was log-transformed.Two linear multi-level (mixed-effects) models were fitted for each response variable using the lme function in the ‘lme4’ package115 (Supplementary Table 8). First, a random intercept model with intercepts of landscape-level predictors varying by country was fitted. Then a random intercept and slope (coefficients) model with intercepts of landscape-level predictors varying by country, as well as slopes for socioeconomic predictors considered to have between-country variation (travel time to nearest city and night-time lights growth) was fitted, as we expect that mangrove cover change may respond to economic growth and market accessibility depending on national governance. A likelihood ratio test between the null linear model and the null random intercept model for each response variable showed that effects varied across countries and therefore we included country as a random effect (Supplementary Table 9). We also conducted likelihood ratio tests between the random intercept model and the random coefficient model to test whether the effect of travel time and night-time lights on mangrove change varies across countries. If significant, the model including random slopes for travel time and night-time lights was used (Supplementary Table 9). Mixed-effects models were fitted by maximum likelihood and model fit was validated by inspection of residual plots for the four response variables included in the analysis; percent net loss, percent net gain, percent gross loss, and percent gross gain (Supplementary Table 9).To test for spatial autocorrelation we performed spatial autoregressive (SAR) models using the errorsarlm function in the ‘spatialreg’ package116. SAR models were first fitted using a range of neighbourhood distances (50, 500, and 1000 km in 100 km intervals) for the net change variable117. Distance of 500 km showed the smallest AIC and was therefore adopted for all response variables. Neighbourhood lists of the centroid coordinates of the geomorphic units were defined with the row-standardised (‘W’) coding using the ‘spdep’ package118. We then produced Moran’s I correlograms using the correlog function in the ‘ncf’ package119 and the centroid coordinates of the geomorphic units. Correlograms for the multi-level model and SAR model were compared for each response variable (Supplementary Fig. 4). The SAR models did not improve spatial autocorrelation for any of the mangrove cover change variables and therefore the multi-level models were adopted.Hotspot estimatesWe defined hotspots as geomorphic units where raw values of percent net and gross loss and gain between 2007 and 2016 ((gamma)) differed by more than two standard deviations (sd) from the country average ((mu)).$${{{{{{rm{More}}}}}}},{{{{{{rm{loss}}}}}}}/{{{{{{rm{more}}}}}}},{{{{{{rm{gain}}}}}}}=left(gamma -mu right) , > , (2,times {{{{{{rm{sd}}}}}}})$$
(1)
$${{{{{{rm{Less}}}}}}},{{{{{{rm{loss}}}}}}},/,{{{{{{rm{less}}}}}}},{{{{{{rm{gain}}}}}}}=left(gamma -mu right) , < , -(2,times {{{{{{rm{sd}}}}}}})$$
(2)
We excluded countries with only one geomorphic unit. Large deviations of the raw value from the country average were found for small units at a threshold below 50 km2, therefore we removed all units smaller than 50 km2 to overcome bias of hotspots towards smaller sites. This likely removed the identification of several hotspots. For example, Myanmar has had some large gains due to river sediments in the Gulf of Martaban (net gain of 100 % in Estuary 5834 and 39 % in Open Coast 62244), however, these areas were small (8 and 2 km2, respectively) and were therefore removed from the hotspot estimates.We analysed the factors contributing to hotspots by spatial investigation of satellite imagery in Google Earth with mangrove specialists from those countries. The hotspots were also assessed against protected area datasets for those countries120,121,122,123.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
