Modelling spatial biodiversity in the world’s largest mangrove ecosystem—The Bangladesh Sundarbans: A baseline for conservation

Mangrove forests are among the most threatened and rapidly vanishing, but poorly understood ecosystems. We aim to uncover the variables driving mangrove biodiversity and produce baseline biodiversity maps for the Sundarbans world heritage site—the Earth's largest contiguous mangrove ecosystem.


| INTRODUC TI ON
Tropical and subtropical mangrove forests (between 30° N and 30° S) provide numerous ecosystem services and support coastal livelihoods worldwide (Lee et al., 2014). However, they are among the most threatened and rapidly vanishing habitats on Earth (Polidoro et al., 2010;Richards & Friess, 2016). The mangrove biome has already lost about 50% of its coverage since the 1950s , and IUCN has listed 40% of mangrove tree species as threatened (Polidoro et al., 2010). Increasing anthropogenic pressures and anticipated sea level rise are likely to alter the structure and functions of the remaining endangered mangrove forests (Duke et al., 2007), in particular, the Sundarbans UNESCO world heritage site-the Earth's largest contiguous mangrove ecosystem.
Making spatial predictions of biodiversity is important for pinpointing the locations or communities requiring immediate or longterm protection and conservation actions, in evaluating threats to those communities and in monitoring spatial distributions and temporal dynamics in biodiversity (Socolar, Gilroy, Kunin, & Edwards, 2015). A variety of biodiversity modelling approaches (e.g., stacked species distribution models, macroecological models, ordination and stochastic models -Ferrier & Guisan, 2006;Mateo, Mokany, & Guisan, 2017) have been applied to understand the spatial patterns of species richness and composition in different forest ecosystems (e.g., neotropical, boreal and temperate forests). However, their application to mangrove forests is limited (but see Record, Charney, Zakaria, & Ellison, 2013) due to the scarcity of field data (Ellison, 2001), thus resulting in poor understanding of mangrove biogeography.
Each of the three established components of biodiversity (alpha, beta and gamma -Whittaker, 1960) characterizes different fundamental attributes of natural communities, and therefore has specific conservation implications. For example, spatial maps of alpha diversity can help in specifying the most species-rich habitats while beta diversity maps can determine the most heterogeneous communities, where protecting larger areas will encompass more biodiversity.
Similarly, gamma diversity measures can identify the overall areas with the highest biodiversity. Thus far, mangrove biodiversity studies have mostly relied on alpha diversity, and in particular species richness (Ellison, 2001;Osland et al., 2017;Record et al., 2013) which, by ignoring the variability in species relative abundances, has known weaknesses in identifying areas for prioritization (Veach, Minin, Pouzols, & Moilanen, 2017). At a regional scale, mangrove plant communities may look spatially homogeneous because mangrove forests are relatively species-poor compared to the upland tropical forests.
However, at finer scales, considerable heterogeneity in vegetation structure becomes apparent (Farnsworth, 1998). Therefore, looking at how the components of biodiversity respond to biotic and abiotic variables is important for constructing more informative and practically useful biodiversity maps.
Mapping biodiversity indices is important in order to investigate spatio-temporal variations in natural communities, to locate habitats or communities or species that require immediate protection and to support spatially explicit conservation planning (Devictor et al., 2010). Both habitat-based and covariate-free (direct interpolation methods such as Kriging) approaches have been used for mapping biodiversity indices. Although covariate-free approaches have been criticized for low predictive ability (Granger et al., 2015), the relative performance of the approaches has rarely been tested using field data.
Testing the "zonation" hypothesis (i.e., the distinct ordering of tree species along the shore-to-inland gradient, Ellison, Mukherjee, & Karim, 2000) and explaining the "biodiversity anomaly" (i.e., why mangrove plant species richness drops along the latitudinal gradient, Ricklefs, Schwarzbach, & Renner, 2006) have been the key agendas dominating the mangrove biodiversity literature in the last two decades. While such studies have substantially improved our insight into species sorting and richness, limited attention has been paid to understanding how abiotic, biotic and historical anthropogenic pressures have contributed to spatial variations in mangrove diversity and composition. Such knowledge gaps have obstructed the success of conservation initiatives in many tropical coastal regions (Lewis, 2005) such as the Sundarbans.
This study focused on the threatened mangrove plant communities of the Sundarbans which are under severe threat from historical forest exploitation, habitat degradation and future climate change impacts (Sarker, Reeve, Thompson, Paul, & Matthiopoulos, 2016). Using a newly introduced abundance-based framework for biodiversity partitioning  and a habitat-based biodiversity modelling approach, our overarching goal was to uncover the influences of fine-scale habitat conditions and historical events in shaping the current spatial distributions of alpha, beta and gamma diversity. Our more specific questions include the following: What are the key drivers of mangrove biodiversity? How do the predictive abilities of covariate-driven habitat models compare with those of covariate-free direct interpolation approaches? Where are the biodiversity hotspots in the Sundarbans currently located? Are these hotspots well protected? Finally, we demonstrate and discuss the demonstrates the utility of direct interpolation approaches when environmental data are unavailable.

K E Y W O R D S
biodiversity conservation, endangered species, generalized additive models, habitat rehabilitation, protected area, sea level rise potential applications of these novel insights and biodiversity maps for future mangrove research, biodiversity protection, monitoring and spatial conservation planning.

| Study area
The Sundarbans (10,017 km 2 ), a part of Earth's largest delta, the Ganges-Brahmaputra, is distributed in Bangladesh and India. Due to its outstanding universal ecological and economic value, the Bangladesh part of the Sundarbans (21°30′-22°30′N, 89° 00′-89°55′E, 6,017 km 2 ) was declared a UNESCO world heritage site in 1997 (Gopal & Chauhan, 2006). It was also declared a Ramsar wetland ecosystem under the Ramsar Convention in 1992 (Chowdhury, Kitin, Ridder, Delvaux, & Beeckman, 2016). The Sundarbans is washed by the tide twice a day, and freshwater flowing from the Ganges and the opposing saltwater influx from the Bay of Bengal together control its hydrology (Wahid, Babel, & Bhuiyan, 2007).
The climate is humid tropical with four main seasons as follows: pre-monsoon (March-May), monsoon (June-September), post-monsoon (October-November) and the dry winter season (December-February). The average annual precipitation is 1700 mm, and the mean temperatures in pre-monsoon, monsoon, post-monsoon, and dry winter are 29, 30, 26 and 20°C, respectively .
The Bangladesh Forest Department (BFD) established these PSPs ( Figure 1) in 1986. As part of the 2008-2014 surveys, our team together with the BFD tagged every tree with stem diameter ≥4.6 cm (because mangroves grow very slowly, this threshold value has been used in all previous forest inventories in the Sundarbans since the 19th century, Iftekhar & Saenger, 2008), at 1.3 m from the ground with a unique tree number and recorded tree counts for the PSPs.
In total, we recorded 49,409 trees from 20 mangrove species (see Appendix S1 in Supporting Information).
In 2014 (January-June), we collected nine soil samples from each PSP (soil depth = 15 cm) adopting a soil sampling design (see Appendix S1 in Supporting Information) to account for the withinplot variations in soil variables. We then determined soil sand, silt F I G U R E 1 Sampling sites (triangles) in the Sundarbans, Bangladesh. Blue areas represent water bodies and clay percentages, salinity, pH, oxidation reduction potential, NH 4 , P, K, Mg, Fe, Zn, Cu and sulphide concentrations. For determining sand, slit and clay percentages, we used the hydrometer method (Gee & Bauder, 1986). We measured soil salinity (as electrical conductivity) in a 1:5 distilled water:soil dilution (Hardie & Doyle, 2012) using a conductivity metre. Soil pH and oxidation reduction potential were measured in the field using soil pH and oxidation reduction potential metres. We followed the Kjeldahl method (Bremner & Breitenbeck, 1983) to determine soil NH 4 and the molybdovanadate method (Ueda & Wada, 1970) to determine total P concentrations.
Soil K, Mg, Fe and Zn concentrations were measured using an atomic absorption spectrophotometer. For each soil variable, we recorded the average reading from nine soil samples.
We retrieved five elevation readings (above-average sea level) from each PSP using the available digital elevation model (accuracy at pixel level = ±1 m) (IWM, 2003) and then averaged them to account for sampling error. We also calculated the "upriver position" (URP), the straight-line distance of each PSP from the river-sea interface (Duke, Ball, & Ellison, 1998) and classified their position as (i) "downstream," representing the lower third (0%-33% upriver from the sea-Bay of Bengal), (ii) "intermediate," representing the middle third (34%-66% upriver from the sea) and (iii) "upstream," representing the upper third (67%-100% upriver from the sea) of the estuarine system. This classification system is useful for understanding variability in diversity and species compositions along the downstream (saltwater-dominated river system)-upstream (freshwater-dominated river system) gradient.

| Covariate selection
We followed Twilley and Rivera-Monroy's (2005) mangrove-centric conceptual framework to construct a biologically informative variable set for our mangrove biodiversity models. This framework integrates abiotic and biotic constraints to explain vegetation structure and productivity at local and regional scales. The abiotic constraints comprise resources, regulators and hydroperiod. Resources (i.e., nutrients) are assimilated by trees. Here, we selected three essential plant macro-nutrients-soil NH 4 , P and K-for their critical roles in mangrove growth and development (Reef, Feller, & Lovelock, 2010).
Regulators are non-resource variables that control tree eco-physiology (Guisan & Thuiller, 2005). Here, we selected soil salinity, pH and silt. Hydroperiod (i.e., inundation frequency, duration and depth) controls the regional and local hydrology that in turn influence species distributions in coastal environments (Crase, Liedloff, Vesk, Burgman, & Wintle, 2013). PSP-level hydroperiod data were unavailable, so we used elevation as a proxy of the likely variation in hydroperiod across the region.
Biotic interactions (e.g., competition or facilitation) between plants can influence species composition at a local scale (Howard et al., 2015). Competitive exclusion of weak competitors in stressed mangrove habitats may lead to species-poor mangrove communities dominated by a single or few opportunistic species (Saenger, 2002). To account for such influences, initially, we considered two candidate biotic variables: (i) "community size" (CS)-total number of individuals in each PSP, and (ii) total basal area in each PSP. Diversity models using basal area as a covariate had lower explanatory powers, compared to models with "CS." Therefore, we selected CS as a proxy of biotic interactions.
We incorporated URP of each PSP in our covariate set to account for the influence of the river systems on species composition along the downstream-upstream gradient. In riverine estuaries, tidal inundation levels, soil physical and chemical properties can significantly vary along the riverbank-inner forest gradient, which influences colonization success and survival of mangrove plants (Berger, 2008).
To account for such variations, we included the straight-line distance of each PSP from the nearest riverbank (henceforth DR).
Natural disturbances (such as tree disease and mortality) and anthropogenic disturbances (such as tree harvesting) offer opportunities for tree recruitment through gap creation, thus influencing vegetation composition (Duke, 2001). To account for the influences of natural and human disturbances on current diversity and species composition, we incorporated historical harvesting (HH) and disease prevalence (DP) as covariates in our models. Here, HH and DP represent the total number of illegally harvested and diseased (e.g., "topdying" disease of Heritiera fomes) trees in each PSP from historical records . Finally, using Variance Inflation Factors (VIF, Robinson & Schumacker, 2009), we checked for multicollinearity in our covariates (see Appendix S2) and removed covariates leading to VIF greater than 2.5. This led to the removal of oxidation reduction potential from our covariate set (see Appendix S2).
Implemented in Reeve et al.'s (2016) framework, this allows us to partition the alpha, beta and gamma diversity of an ecosystem (called a metacommunity [MC]) into its subcommunity (SC) components, thus allowing comprehensive and consistent quantification and modelling of all biodiversity components in a spatial context.
In this study, each PSP represents a SC, and the combined PSPs form the MC. This approach allows us to understand and easily compare the species diversity and composition in every single SC in relation to the MC (the whole Sundarbans ecosystem). We measured SC alpha, beta and gamma diversity. Here, the normalized alpha diversity index (denoted ̄) represents the diversity of a single SC (PSP) in isolation. The normalized beta diversity index (denoted ̄) measures representativeness and assesses how well a SC represents the species composition of its MC. It is maximized (1) when the MC is homogenous, and a SC's species composition is identical to that of the MC and therefore represents it perfectly. Low ̄ therefore suggests high spatial heterogeneity in species composition within the MC, and high ̄ suggests spatial homogeneity.
The gamma diversity (denoted γ) is the conventional gamma diversity (Hill, 1973;Jost, 2006;Leinster & Cobbold, 2012) at the MC level that reflects the total species diversity in an unpartitioned ecosystem. Reeve et al.'s (2016) framework partitions the MC gamma diversity into SC gamma diversity that measures each PSP's average contribution to (or influence on) the MC diversity per tree. This diversity measure combines the alpha diversity of a SC with its beta diversity to form an assessment of the overall contribution of the PSP to the MC.
Following Hill (1973), Jost (2006Jost ( , 2007 and Leinster and Cobbold's (2012), the values of all the biodiversity measures are moderated by a viewpoint parameter, q, taking a value between 0 and ∞ representing how conservative the measure is in accounting for species abundance. For ̄ and , the diversity at q = 0 measures species richness; at q = 1 measures the exponential of Shannon entropy (Shannon, 1948); and at q = 2 measures the inverse of Simpson's concentration index (Simpson, 1948). For all analyses, we present the results using the above three q values (0, 1, and 2), writing them as 0̄, 1̄, 2 , etc.

| Biodiversity modelling
We constructed generalized additive models (GAMs, Wood, 2011) to quantify how the different biodiversity components responded to different variables. Guided by data and using nonparametric smoothing functions, GAMs can capture response-predictors relationships without a priori knowledge of the functional form of these relationships (Guisan & Thuiller, 2005). These advantageous features of GAMs are well suited for uncovering unknown biodiversity-environment linkages in dynamic ecosystems such as the Sundarbans where multiple environmental gradients have interactive effects on species distributions (Sarker et al., 2016). All analyses were done in R version 3.2.3 (R Core Team, 2016). Biodiversity GAMs were built using cubic basis splines with the Gamma error distribution using the "mgcv" package version 1.8 -7 (Wood, 2011). Model selection and model averaging were carried out using the "MuMIn" package version 1.15.1 (Barton, 2015). Biodiversity measures were calculated using the "rdiversity" package version 1.0 (Mitchell & Reeve, 2017).
We exhaustively fitted GAMs for each diversity index with all possible combinations of covariates. Then, we ranked the fitted GAMs using the second-order AIC (AICc) because the ratio between sample size and the number of covariates was <40 (Burnham & Anderson, 2002). Models whose AIC c had values less than 2 units from the best model (∆AIC c < 2) were retained as competing models (Burnham & Anderson, 2002). The relative support for each of the competing models was then determined using their Akaike weights (AIC cw vary between 0 and 1, and the sum of all AIC cw across the competing models is 1). To reduce model selection uncertainty and bias, we then conducted model averaging to predict the diversity indices. To determine the strength of the covariates, we ranked them based on their relative importance (RI) values. RI of each covariate was calculated by totalling the AIC cw of the models in which the covariate was included. RI values vary between 0 and 1, where 0 specifies that the target covariate is not included in any of the competing models while 1 means that the covariate is included in all competing models. We measured goodness of fit of the biodiversity models using the R 2 (coefficient of determination) statistic between the observed and estimated values of the diversity indices.

| Biodiversity mapping
We applied two different approaches to make spatial biodiversity predictions. First, we used our habitat-based models (GAMs) and interpolated covariate surfaces to produce model-averaged predictions. Second, we used a direct interpolation method-ordinary kriging-to make purely spatial predictions. We compared these two approaches because environmental data collection is challenging, whereas tree surveys are conducted annually at the PSPs. Hence, it is useful to know how close the predictions of the habitat-based biodiversity models were compared to direct interpolation methods.
The size of each grid cell of the interpolated surfaces was 625 m 2 .
We compared the predictive abilities of GAMs with ordinary kriging, using the normalized root mean square error (NRMSE) statistic derived from a leave-one-out cross-validation procedure. For normalization, the root mean square error statistic was divided by the range of the actual diversity values. Ordinary kriging was performed using the "gstat" package version 1.0 -26 (Pebesma, 2004) in R.
A protected area network comprising three Wildlife Sanctuaries has been operational in the Sundarbans since the 1970s. To evaluate its capacity to support the remaining biodiversity hotspots in the Sundarbans, we superimposed this onto our biodiversity maps. All the biodiversity maps were constructed using the "raster" package version 2.4 -18 (Hijmans, 2017) in R.

| Habitat-based biodiversity models
The explanatory power and the goodness of fit of the alpha, beta and gamma diversity GAMs varied when we increased weight on species relative abundances (q = 0, 1 and 2) in the subcommunities (SCs). 1̄ (Shannon entropy) GAM explained more deviance (DE = 71%) and showed a better fit to the data (Adj. R 2 = 0.71) compared to those for 0̄ (species richness) and 2̄ (Simpson's concentration) ( Table 1), suggesting that, for alpha diversity, the model with a moderate focus on species relative abundances in the SCs (i.e., q = 1) could capture more signal compared to the models that only considered species presence-absence (q = 0) or offered more importance to the more dominant species (q = 2) in the SCs. Like 1̄, the 1̄ GAM could capture more signal than 0̄ and 2̄ GAMs. In contrast, for beta diversity, with DE = 65% and Adj. R 2 = 0.70, the 2̄ GAM captured more signal than the 0̄ and 1̄ GAMs, implying that our covariates could more successfully explain the variability in species composition across the SCs when the variability was mostly contributed by more dominant species.

| Drivers and responses of biodiversity components
The RI of the covariates in influencing biodiversity indices also varied when we changed weight on species relative abundances in the SCs. For example, while HH had no influence on 0̄ (possibly due to high number of shared species between SCs or HH did not lead to species extirpation), it had stronger effects on 1̄ and 2̄, indicating that the influence of past tree harvesting in shaping current community composition becomes clearer when we account for the variability in species relative abundances across the SCs.
The partial response plots of the best alpha, beta and gamma diversity GAMs (for q = 0, 1 and 2) showed similar relationships across the models (Figure 2, Appendix S3). While alpha diversity (for 1̄) increased with increasing DR (>1,500 m) and URP (>80%), it decreased with increasing HH (>175 tree cuts/0.2 ha), silt (>20%), CS (>450 trees/0.2 ha) and pH (>7.25). The response of alpha diversity varied for different nutrients. The K concentration that maximized 1̄ was 5.5 gm/Kg while increasing soil P (>35 mg/Kg) was related to decreasing 1̄. Mangrove communities showed increasing representativeness (for 2̄) , that is, homogeneity in species composition with increasing HH (>150 tree cuts/0.2 ha), silt TA B L E 1 Results of generalized additive models (GAMs) for nine diversity measures. Summaries of model fit in rightmost three columns are only shown for the best model (DE = deviance explained). Numbers in the main part of the table (enclosed in box) represent the Relative Importance (RI) of each covariate. Dark-shaded cells highlight covariates that were retained in the best model for each biodiversity index. Light-shaded cells represent covariates retained in other models within the candidate set. Dashed boxes indicate no participation of that covariate in any of the candidate models. The covariate shorthands are community size (CS), upriver position (URP), salinity, distance to riverbank (DR), historical harvesting (HH), acidity (pH), silt concentration, disease prevalence (DP), soil total phosphorus (P), soil potassium (K), elevation above-average sea level (ELE) and soil NH 4 . F I G U R E 2 Effects of covariates inferred from our best generalized additive models fitted to the biodiversity indices for q = 1. The solid line in each plot is the estimated spline function (on the scale of the linear predictor), and shaded areas represent the 95% confidence intervals. Estimated degrees of freedom are provided for each smooth following the covariate names. Zero on the y-axis indicates no effect of the covariate on diversity index values. Covariate units: CS = total number of individuals in each plot, URP = % upriver, soil salinity = dS/m, DR = distance (m) of each PSP from the riverbank, historical harvesting (HH) = total number of harvested trees in each plot since 1986, silt (%), disease prevalence (DP) = total number of diseased trees in each plot since 1986, P = mg/Kg and K = gm/Kg (>20%), DP (>25 diseased trees/0.2 ha) and CS (>450 trees/0.2 ha).
F I G U R E 3 Spatial distributions of SC alpha, beta and gamma diversities (for q = 0-2) over the entire Sundarbans generated through generalized additive models. Higher values of ̄ and γ indicate greater species diversity and community contribution to the overall diversity of the ecosystem. Lower values of ̄ indicate greater heterogeneity in species composition (i.e., community distinctness from the metacommunity), and higher values of ̄ represent greater representativeness (i.e., homogeneity) in species composition. The black contours represent the three protected areas

| Biodiversity maps
Spatial diversity maps are presented in Figure 3. Alpha diversity maps (first row) uncovered that hotspots in species richness (q = 0), Shannon entropy (q = 1) and Simpson's concentration (q = 2) were restricted to the northern (specifically, the Kalabogi region) and eastern (specifically the Sharankhola region) Sundarbans. Beta (second row) and gamma (third row) diversity maps revealed that the entire Sundarbans looks homogeneous when we only looked at species presence or absence (q = 0), that is, not accounting for the between-species variability in relative abundances. Allowing increasing weight on species abundance (q = 1 and 2) revealed that the most heterogeneous mangrove communities and the communities that contributed most to the overall biodiversity of the ecosystem were restricted to the northern upstream habitat. Additionally, our maps indicated that the most diverse (i.e., biodiversity hotspots) and heterogeneous mangrove communities are situated outside the established protected area network. Prediction error was always reduced by the use of environmental covariates, but particularly for predictions of alpha and gamma diversity. In case of beta diversity, while the predictive ability of the GAM was better than that of kriging for 0̄ and 1̄, both approaches had almost similar prediction error for 2̄ (Table 2).

| D ISCUSS I ON
This study provides a baseline quantification and habitat-based modelling of alpha, beta and gamma diversity of threatened mangrove communities. Contrary to the common assumption that one or two straightforward environmental gradients (salinity and inundation) control mangrove biodiversity (Ellison, 2001), our results revealed that several environmental drivers, biotic interactions and historical events contribute to the emergence of observed spatial patterns of mangrove biodiversity. The high explanatory power and predictive power of our biodiversity models confirm their usefulness in making spatially explicit predictions of species diversity and composition in dynamic mangrove ecosystem. The ability of the models to reveal previously unknown linkages between the biodiversity components and abiotic, biotic and disturbance variables have yielded novel biological insights and thus now prompt many ecological questions for future studies. Our analyses uncovered a strong impact of HH and DP in shaping current distributions of the biodiversity components in the Sundarbans, implying the importance of integrating past disturbance events in habitat-based models for more accurate predictions. We detect a significant negative effect of HH on alpha and gamma diversities, although DP has no visible effect. This discrepancy may be related to local extinction of many rare endemics during past formal and informal logging activities and high DP (top-dying and heart rot diseases) in the specialists (i.e., H. fomes and X. mekongensis) (Banerjee, Gatti, & Mitra, 2017) that might not lead to their extirpation but reduced their relative abundances in a higher amount TA B L E 2 Comparison of predictive accuracy (through leaveone-out cross-validation) of the habitat-based (GAMs) and Kriged diversity models using normalized root mean square error (NRMSE) of the predicted versus the actual diversity values. NRMSE is expressed here as a percentage, where lower values indicate less residual variance. compared to the generalists. However, for beta diversity, both HH and DP contributed to increasing homogeneity in species composition across the SCs (Figure 2). This again indicates that the diseases have not infected all trees equally rather, they have only infected and removed a few specialists such as H. fomes and X. mekongensis which have resulted in increasing homogeneity in the mangrove communities. Therefore, by using the approach of Reeve et al. (2016) to look at how DP simultaneously affects alpha, beta and gamma diversity, we are now able to get indications of the pathogenicity of the disease (i.e., whether it is a generalist and infects and removes all species equally or it is specialized on specific host species).
Harvesting-and disease-induced tree mortalities have created many large as well as small forest gaps in the Sundarbans. Intriguingly, the large diameter tree species (i.e., H. fomes and X. mekongensis) that still dominate the less saline habitats, recruit poorly in the forest gaps (Iftekhar & Islam, 2004). Instead, these forest gaps are increasingly colonized by the disturbance specialists (e.g., C. decandra) (Mukhopadhyay et al., 2015). Therefore, increasing colonization and dominance of disturbance specialists in the historically disturbed billion tons of sediments (Mitra & Zaman, 2016). Therefore, future research is required to understand species-specific sensitivities and adaptations (e.g., modified rooting architecture) to siltation because this will help to forecast which species may colonize the newly formed islands and which are compatible for replanting under future siltation scenarios.
In their pioneering work, Ellison et al. (2000) found no evidence for "zonation" in the Sundarbans. In contrast, we detect a clear pattern of increasing alpha and gamma diversities along the riverbank/ forest interior gradient. Communities that are at least 1,500 m away from the riverbank have higher alpha diversity and 800 m away have higher gamma diversity compared to the near-bank communities ( Figure 2), implying late-successional forest interior communities are more diverse than the early successional riverbank communities.
Salinity has been considered a key constraint limiting species richness in coastal ecosystems ). It appears from our analyses that salinity has no effect on species richness although the importance of salinity slightly increased for Shannon entropy and Simpson concentration, implying the role of salinity becomes clearer when we account for between-species variability in relative abundance. Considering beta diversity, increasing salinity contributes to increasing compositional heterogeneity among the SCs (Figure 2).
This pattern suggests high plot-to-plot variation in composition in the degraded saline soils for population declines and range contraction of many salt-intolerant specialists (e.g., H. fomes) and increasing colonization success of the salt-tolerant generalists such as E. agallocha and C. decandra (Aziz & Paul, 2015;Iftekhar & Saenger, 2008;Mukhopadhyay et al., 2015).
Nitrogen, phosphorus and potassium were found to be the important soil nutrients limiting mangrove forest structure in coastal areas in Brazil, Florida and South Africa (Da Cruz et al., 2013;Lovelock, Ball, Feller, Engelbrecht, & Ling, 2006;Naidoo, 2009).
Interestingly, these resource variables received less support in our biodiversity models, reconfirming the high importance of regulators and historical disturbances in structuring mangrove communities (Twilley & Rivera-Monroy, 2005).

| Mangrove biodiversity maps
Our biodiversity maps for the Sundarbans (Figure 3) reveal that currently the most species-rich ( 0̄) mangrove communities are confined to the northern (specifically, Kalabogi) and eastern (specifically, Sarankhola) regions. Due to the proximity of Baleshwar and Posur rivers, these areas receive greater amount of freshwater than the rest of the ecosystem, thus securing suitable conditions for many salt-intolerant and rare plant species. The remaining ecosystem is relatively species-poor. 1̄ and 2̄ maps not only show similar patterns but also pinpoint the areas-the north-western and south-western Sundarbans-where the super-dominance of generalists has resulted in lower alpha diversity. These areas are prone to regular saltwater flooding and high salinity fluctuation which together were found to inhibit regeneration and growth of many species (Ghosh, Kumar, & Roy, 2016). Spatial variability in beta diversity becomes clearer when more weight was put on the dominant species ( 1̄, 2̄) , compared to the rare species ( 0̄) . In general, the most heterogenous communities and the communities that contribute most to the overall biodiversity of the whole ecosystem ( 0 , 1 , 2 ) are currently restricted to the northern upstream habitats supporting tree species facing the risk of local (X. mekongensis) and global (H. fomes) extinction (Sarker et al., 2016).
Restricted distributions of diverse and distinct mangrove communities in a few specific areas clearly indicate for historical pressures on Sundarbans's floral composition, as reported by many (Aziz & Paul, 2015;Ghosh et al., 2016;Gopal & Chauhan, 2006).
The freshwater supply from the transboundary rivers into the Sundarbans has substantially declined (3,700 m 3 /s to 364 m 3 /s) since the construction of the Farakka dam (1974) in India (Mirza, 1998). The average soil salinity has already increased by 60% since 1980 (Aziz & Paul, 2015). Opportunistic harvesting of trees and heavy siltation in the internal channels are ongoing (Rahaman et al., 2015). Therefore, our findings lead us to conclude that additional harvesting, siltation, cuts in freshwater supply and range expansions of the generalists under projected sea level rise (Karim & Mimura, 2008) may convert the whole Sundarbans into a species-poor homogeneous ecosystem.
The existing approaches for biodiversity mapping without including environmental data are shown to produce inaccurate spatial predictions of diversity indices (Granger et al., 2015). In this study, in general, the environmental data-driven GAMs showed better predictive ability than the covariate-free direct interpolation method (Table 2), thus, supporting the inclusion of fine-scale environmental, biotic and historical disturbance data for more accurate mapping of biodiversity indices when these data are available. However, similar performances of these approaches in predicting 2̄, and small differences in prediction error for 0̄ and 0 , indicates the utility of direct interpolation methods when environmental data are not available.

| Conservation applications
Sea level rise is likely to have drastic impacts on mangrove forests worldwide, particularly, the Sundarbans. Under the projected range of sea level rise by 2,100 (30-100 cm) (Karim & Mimura, 2008), the Sundarbans is likely to lose 10%-23% of its present area (Payo et al., 2016) with alteration to soil biogeochemistry (Banerjee et al., 2017) and estuarine hydrology (Wahid et al., 2007). Given the severity of these future environmental impacts on Sundarbans, identifying the existing and future environmental stressors of mangrove biodiversity is important. We detect siltation, salinity and pH as the dominant environmental stressors responsible for decreasing mangrove diversity (  The country has also ratified the "Bangladesh Biodiversity Act 2017" to stop illegal trade of forest flora and fauna. It has also adopted a SMART (Spatial Monitoring and Reporting Tool) patrol management system since 2015 to expand the scope of its current mangrove protection efforts. Our baseline biodiversity maps can guide these valuable biodiversity protection and conservation initiatives. In addition, these maps can contribute to successful implementation of the REDD+ (Gardner et al., 2012) initiatives for enhancing carbon stock (through biodiversity conservation) as well as financial returns.

| CON CLUS IONS
This study provides the first comprehensive and coherent quantification and habitat-based modelling of alpha, beta and gamma diversity in threatened mangrove communities of the Sundarbans.
We find that several environmental drivers, biotic interactions and historical events have combined effects on the biodiversity components. Specifically, salinity intrusion, HH, increasing CS, siltation, disease and soil alkalinity are the dominant stressors responsible for reducing mangrove diversity. Although habitat-based models showed better predictive ability than the covariate-free approach, the small margin of differences between the approaches demonstrates the utility of direct interpolation approaches when environmental data are unavailable. Our biodiversity maps uncover that the most diverse and distinct mangrove communities (biodiversity hotspots) have restricted distributions in the freshwater-dominated northern and eastern regions. Although these biodiversity hotspots are susceptible to human exploitation, they are not included in the existing protected area network, thus suggesting for an immediate reconfiguration of the protected area network. We believe details on the environmental stressors, and our biodiversity maps, collectively, will contribute to designing and implementing climate-smart mangrove enhancement, restoration and reforestation initiatives.
In addition, our maps can guide the existing and future mangrove biodiversity protection, monitoring and REDD+ initiatives.

ACK N OWLED G EM ENTS
We thank the editor and two anonymous referees for their constructive suggestions. We gratefully acknowledge the numer-