Solving the fourth-corner problem: forecasting ecosystem primary production from spatial multispecies trait-based models

. Forecasting productivity and stress across an ecosystem is complicated by the multiple interactions between competing species, the unknown levels of intra- and interspecific trait plasticity, and the dependencies between those traits within individuals. Integrating these features into a trait-based quantitative framework requires a conceptual and quantitative syn-thesis of how multiple species and their functional traits interact and respond to changing envi- ronments, a challenge known in community ecology as the “ fourth-corner problem. ” We propose such a novel synthesis, implemented as an integrated Bayesian hierarchical model. This allows us to (1) simultaneously model trait – trait and trait – environment relationships by explicitly accounting for both intra- and interspecific trait variabilities in a single analysis using all available data types, (2) quantify the strength of the trait – environment relationships, (3) identify trade-offs between multiple traits in multiple species, and (4) faithfully propagate our modeling uncertainties when making species-specific and community-wide trait predictions, reducing false confidence in our spatial prediction results. We apply this integrated approach to the world ’ s largest mangrove forest, the Sundarbans, a sentinel ecosystem impacted simultaneously by both climate change and multiple types of human exploitation. The Sundarbans presents extensive variability in environmental variables, such as salinity and siltation, driven by changing seawater levels from the south and freshwater damming from the north. We find that tree species growing under stress have a typical functional response to the environmental drivers with inter-specific variability around this average, and the amount of variability is fur- ther contingent upon the nature and magnitude of the environmental drivers. Our model cap-tures the retreat in traits related to resource acquisition and a plastic enhancement of traits related to resource conservation, both clear indications of stress. We predict that, if historical increases in salinity and siltation are maintained, one-third of whole-ecosystem productivity will be lost by 2050. Our integrated modeling approach bridges community and ecosystem ecology through simultaneously modeling trait – environment correlations and trait – trait trade-offs at organismal, community, and ecosystem levels; provides a generalizable foundation for powerful modeling of trait-environment linkages under changing environments to predict their consequences on ecosystem functioning and services; and is readily applicable across the Earth ’ s ecosystems.


INTRODUCTION
Understanding and predicting the effects of global environmental changes on the Earth's ecosystems is a core challenge confronting community ecology, ecosystem ecology, and the nascent field of functional biogeography. To address this challenge, the use of plant functional traits, rather than taxonomic identities, has gained momentum (McGill et al. 2006, Diaz et al. 2007, Kraft et al. 2008, 2015a,b, Lavorel et al. 2011, Soudzilovskaia et al. 2013, Kunstler et al. 2015, Cernansky 2017, Faucon et al. 2017, Bjorkman et al. 2018a, Schweiger et al. 2018, Anderegg et al. 2019. By capturing essential features of plants' ecophysiology, morphology and life history, functional traits deliver a mechanistic connection between fundamental biological processes and species performance, community dynamics, ecosystem functions, and services (McGill et al. 2006, Adler et al. 2014. Therefore, ecologists have started using trait-based approaches to predict changes in community assembly (Kraft et al. 2008, Spasojevic and Suding 2012, Adler et al. 2013, Swenson 2013, Lasky et al. 2014, He et al. 2019, ecosystem function and services (Lamarque et al. 2014, Wood et al. 2015, Faucon et al. 2017, Hanisch et al. 2020, and global vegetation dynamics (Van Bodegom et al. 2014, Sakschewski et al. 2015. Functional traits include any measurable morphological or physiological feature that may influence the overall fitness or performance (e.g., growth, survival) of organisms (Violle et al. 2007). In nature, the productivity of primary producers is affected by a complex causal network encompassing multiple functional traits, multiple species and multiple environmental drivers (Fig. 1). Changes in traits via natural selection and species sorting (species elimination under persistent stress), may alter plant performance, community structure and primary production (Webb et al. 2010) but natural selection and species sorting do not act independently on single traits, single species, or as a result of single environmental stressors (Verberk et al. 2013). Therefore, attempting to model productivity without integrating these complex dependencies between traits, species and their environment constricts our mechanistic insights and limits reliable predictions.
Determining trait-environment relationships (TER) has been a long-standing problem, known as "the fourth-corner problem" in community ecology (Legendre et al. 1997). For its resolution, ecologists have proposed trait-based approaches such as the fourth-corner correlation (Legendre et al. 1997, Dray andLegendre 2008), the multivariate RLQ (Dolédec et al. 1996, Dray et al. 2014, the community weighted trait means (CWM) (Lavorel et al. 2007), and regression methods (Jamil et al. 2012, Pollock et al. 2012, Brown et al. 2014, Warton et al. 2015. Historically, the fourth-corner correlation and RLQ approaches have received much attention . Using a combination of environmental (R), species (L), and trait (Q) data tables and permutation tests, these approaches look at pairwise associations between a single trait and single environmental variable using significance testing. They are routinely applied by ecologists for identifying the response traits, i.e., "traits that mediate the response of plant species to the environment" (Pakeman 2011). While such trait-based approaches are insightful (Funk et al. 2017), substantial methodological difficulties remain (Webb et al. 2010, Verberk et al. 2013, Violle et al. 2014, Braga et al. 2018, ter Braak et al. 2018) that limit their utility for community ecology, ecosystem ecology and functional biogeography. For example, they cannot uncover the effect of multiple environmental drivers acting on the same functional trait, and significance testing says little about the strength of TER.
Recently proposed regression approaches relax these limitations by supporting model selection, validation, and quantitative predictions. They model species abundance (or presence-absence) as a function of traits and environmental variables and use their interaction terms (a product of a single trait with a single environmental variable) to describe TER . However, such approaches have focused mostly on one species at a time (Jamil et al. 2012) and ignored intraspecific covariation between traits (Funk et al. 2017), hence not managing to account for the complex interactions between species, traits, and the environment (Levine et al. 2017, Staniczenko et al. 2017). This may be the cause of their limited predictive power (Funk et al. 2017. The CWM approach has been widely applied to directly correlate community-level trait values with environmental variables , Bruelheide et al. 2018, Wüest et al. 2018, Boonman et al. 2020. However, describing TER via communitylevel instead of individual species-level data, likely understates individual heterogeneity (aggregation bias; Webb et al. 2010), thus resulting in inaccurate TER quantification and inconsistent parameter estimation (Verberk et al. 2013, Miller et al. 2018. The shortcomings of existing trait-based approaches have raised concerns (Webb et al. 2010, Verberk et al. 2013, Violle et al. 2014, Funk et al. 2017, Schliep et al. 2018, ter Braak 2019) about how accurately they will forecast species, community, and ecosystem responses under future environmental scenarios.
Trait distributions in individual species, nested within communities and ecosystems, respond dynamically to multiple environmental gradients across space or time. Again, traits are known to vary both between and within species (Jung et al. 2010, Bolnick et al. 2011, Pérez-Harguindeguy 2013. The hierarchically structured and dynamic nature of the TER and the intra-and interspecific variabilities in trait values necessarily constrains FIG. 1. Productivity in plant species is governed by community composition, functional traits, and spatiotemporal environmental heterogeneity, but these three determinants also interact with each other in non-trivial ways. the statistical or mathematical structure of quantitative approaches that are appropriately suited to support rigorous hypothesis testing and prediction. Webb et al. (2010) and Verberk et al. (2013) suggested that "traitsbased ecology is at a critical juncture" and that, to advance the field, we need an integrative quantitative framework that can appropriately analyze the implicit structure of trait data to translate the TER (i.e., the performance filter; Webb et al. 2010) into forecasts of community and ecosystem dynamics, thus bringing stronger theoretical linkages between community ecology, ecosystem ecology, and functional biogeography. However, by looking at species, traits and the environment in a disconnected way, and by failing to simultaneously account for both intra-and interspecific trait variations in models, existing trait-based approaches offer limited scope for accurate and precise spatial predictions. We therefore need a synthetic and data-driven approach for coexisting species, their traits, responses to environmental variables and the prediction uncertainties arising from (or being constrained by) these interactions.
Here, we develop an integrated Bayesian approach that simultaneously models trait-trait relationships (TTR) and trait-environment relationships (TER) for multiple traits, multiple species, and multiple environmental drivers at organismal, community, and ecosystem levels. This approach offers the formulation of a range of ecological hypotheses (about trait-environment linkages) and allows us to test them in a single analysis using all available data. By recognizing covariation between traits originating from life-history trade-offs and by fully propagating parameter uncertainty (from fitting through to prediction), this approach provides robust TER quantification and delivers a coherent and integrative quantitative framework to translate the TER into forecasts of community and ecosystem dynamics under future environmental conditions. Physiologically informed priors based on sound biological first principles can enter our model (Kindsvater et al. 2018). Therefore, the underlying mechanistic relationships between the traits and the environmental drivers (i.e., the performance filter) can be ported into other purely mechanic modeling approaches, such as First Principles Dynamical Systems Models and First Principles Bayesian Multilevel Models (Webb et al. 2010), that have not yet been extended to make trait-based predictions in ecology.
We apply our approach to the world's largest mangrove forest, the Sundarbans World Heritage Ecosystem (6,017 km 2 ), a sentinel ecosystem that is being impacted simultaneously by both climate change and multiple environmental and human pressures. Floristically, this ecosystem belongs to the Indo-Andaman mangrove province within the most species-rich Indo-West Pacific group and accounts for one-third of the global mangrove tree species count (Ghosh et al. 2015, Mahmood 2015. This sea-dominated, dynamic ecosystem is ideal for testing mechanistic, trait-based approaches because many of the highly productive tree species (including several globally endangered tree species), well-recognized for their contributions to mangrove ecosystem functioning and services (Rahman et al. 2015, Khan et al. 2020, need to maintain their fitness under constant environmental stress. We focus on nine prominent tree species (constituting 99% of the total tree populations in the Sundarbans; Iftekhar and Saenger 2008), eight prominent environmental drivers, and four key traits (canopy height [height], specific leaf area [SLA], wood density [WD], and leaf succulence [LS]) covering the "leaf economic spectrum" (Wright et al. 2004) and the "wood economic spectrum" (Chave et al. 2009) and reflecting the "acquisitive-conservative continuum" (Grime et al. 1997, Reich et al. 1997, Reich 2014, Díaz et al. 2015. Height and SLA represent a plant's acquisitive resource-use strategy, hence growth, while WD and LS describe a plant's conservative resource-use strategy, hence stress-tolerance ability (Pérez-Harguindeguy et al. 2013). Using our integrated approach, we asked (1) Which set of assumptions about trait responses to the environment are best supported by the data? (2) How do the different traits of each species respond to a range of environmental drivers? (3) Are there trait-trait tradeoffs that are detectable across the species? The Sundarbans mangrove forest has already become a highly stressed ecosystem due to a historical and ongoing reduction in freshwater flows into the river system, siltation, and salinity intrusion due to sea-level rise (Wahid et al. 2007, Aziz andPaul 2015). Further deterioration of the Sundarbans is likely to significantly affect plant functions, community composition and overall productivity of the forest. Therefore, we modeled the fate of the Sundarbans by developing productivity maps under present and future environmental scenarios, taking care to quantify the relative contribution of changes in community composition and species traits.

Study system
The Bangladesh Sundarbans (21°30 0 -22°30 0 N, 89°0 0 0 -89°55 0 E, 6,017 km 2 ) originated following the fragmentation of the Gondwanaland in the early Cretaceous. It is located on the Earth's largest delta, the Ganges-Brahmaputra, in South Asia. It has a humid, maritime, and tropical climate with mean annual rainfall of 1,700 mm and mean maximum annual temperature between 29.4°and 31.3°C (Chowdhury et al. 2016b). The soil type is silty-clay-loam. Sea water from the Bay of Bengal inundates most of the forest twice a day. This globally important Ramsar wetland ecosystem (Chowdhury et al. 2016b) not only supports the livelihoods of 3.5 million coastal inhabitants and protects them from cyclones/tsunamis (Aziz and Paul 2015) but also acts as a repository of many globally endangered plant and animal species (Iftekhar and Islam 2004). The Sundarbans is under constant stress from salinity intrusion, forest exploitation, sea level rise (SLR), freshwater flow reduction, siltation, oil spills, and cyclones (Harun-or-Rashid et al. 2009, Dasgupta et al. 2014.
We recorded four core functional traits: height (H, tree canopy height), wood density (WD, dry mass per unit volume of wood), specific leaf area (SLA, leaf area per unit leaf mass), and leaf succulence (LS, water content on a leaf area basis). Height, a whole-plant trait, is central to a plant's resource (carbon) acquisition and hence its ability to compete for light (Díaz et al. 2015). Height regulates valuable ecosystem services including carbon sequestration capacity (via its relationship with tree biomass) and wildlife diversity (for example, bird and mammal species prefer certain canopy layers; Moles et al. 2009). SLA, a core trait of the "leaf economic spectrum" (Wright et al. 2004), reflects the efficiency of a leaf for light capture per unit biomass invested (Poorter et al. 2009), hence regulates the relative growth rate of trees (Nicotra et al. 2010). WD is an essential trait of the "wood economic spectrum" (Chave et al. 2009) that controls hydraulics, architecture, defense, and growth potential of trees (Poorter et al. 2010). Denser wood provides mechanical stiffness to stems, hence increasing trees' ability to survive under stress (Lasky et al. 2014, Lawson et al. 2015. LS offers an accumulation of larger amounts of water and dissolved ions in leaves under hypersaline environments. Therefore, LS is a prominent indicator of species' resistance to salinity and drought (Das 1999, Naskar andPalit 2015).
Our candidate environmental drivers for the models comprised three broad categories of variables: (I) resources, (II) regulators, and (III) disturbance (Twilley and Rivera-Monroy 2005). Resources (i.e., nutrients) are assimilated by plants. Here, we selected three essential plant macronutrients, soil NH 4 , P, and K, for their importance in mangrove growth and development (Reef et al. 2010). Regulators are non-resource variables that control plant eco-physiology (Guisan and Thuiller 2005). Here, we selected soil salinity, silt, and pH for their critical roles in mangrove establishment (Sarker et al. 2016). We also selected upriver position (URP, the proportional distance of PSPs from the river-sea interface) as a regulator variable to account for the influence of the downstream-upstream gradient on species compositions and functions (Duke et al. 1998). Tropical coastal ecosystems are prone to anthropogenic disturbances such as tree harvesting (Feller et al. 2017), which offers opportunities for tree recruitment through gap creation, thus influencing forest structure and functions (Duke 2001). Therefore, we used historical harvesting (HH) as the disturbance variable to account for its possible influences on current species compositions and functions.

Trait and environmental data measurements
We measured all functional traits following standard trait measurement protocols (Pérez-Harguindeguy et al. 2013). From each PSP, we recorded trait measurements for three individual trees per species. For each tree, we collected a wood sample and at least three mature leaves from sun-exposed branches for laboratory measurements. All trait samples, as well as the environmental data, were collected in 2014 (January-June). We determined canopy height (m) using a Suunto clinometer (Suunto, Vantaa, Finland). For specific leaf area, we (1) measured green masses of the leaves immediately after collection in the field, (2) derived fresh area of the leaves using their images (captured by a digital camera, Nikon D5500, Nikon, Tokyo, Japan) in Adobe Photoshop CS6 (Adobe Inc, San Jose, California, USA), (3) quantified dry mass of the leaves using a digital balance (precisioñ 0.001 g) after keeping the green leaves in an oven at 65⁰C for 72 h, and (4) calculated SLA as fresh leaf area (cm 2 )/dry mass (g). Leaf succulence was calculated as [leaf green mass (g) − leaf dry mass (g)]/fresh leaf area (dm 2 ) (Wang et al. 2011). For WD, we (1) collected wood cores using an increment borer, (2) estimated fresh volume of the cores using the formula of a cylinder (V = πr 2 l, where r is the radius of the core and l is the length of the core), (3) quantified dry mass of the cores using a digital balance (precision~0.001 g) after keeping the fresh cores in an oven at 105°C for 72 h, and (4) calculated WD as dry mass (g)/fresh wood volume (cm 3 ).
For the environmental drivers, the data were collected in 2014 and published by Sarker et al., (2019b). NH 4 , P, K, salinity, silt, and pH were measured from the soil, for which we analyzed nine samples from each PSP (soil depth, 15 cm) and averaged them to account for withinplot variation. URP of the PSPs were measured from GIS data indicating downstream, intermediate or upstream location relative to the sea. Finally, HH levels were determined from stump count records (measuring illegal harvesting between 1986 and 2014). See Sarker et al. (2019a) for further details.

Model development
Our approach is an extension of the multivariable regression models typically used to quantify species-habitat associations (Matthiopoulos et al. 2020) on the basis of one response and multiple explanatory variables. We have augmented this approach by allowing it to model several interdependent response variables (multiple functional traits from multiple species), using a Bayesian hierarchical modeling framework. This allowed a single model, fitted simultaneously to all our data, to capture trait-environment relationship (TER) for four traits, nine tree species, and eight environmental drivers, while accounting for trait-trait (TTR) relationships across the species. Based on this inferential template, several trait-based models were developed taking a range of TER and TTR assumptions (see Table 1). These models are not an exhaustive set of the trait-based models that can be developed using our approach; instead, we chose models that address frequently asked questions about how traits interact with environmental drivers in extreme environments such as the mangrove forests (Weiher and Keddy 1995, Twilley and Rivera-Monroy 2005, Grime 2006, Minden et al. 2012, Verberk et al. 2013, Naskar and Palit 2015, Duke 2017, Xu et al. 2017, Bruelheide et al. 2018, He et al. 2019. The TER components of the models include six different features to accommodate six different assumptions about how multiple traits of multiple species respond to multiple environmental drivers by considering that (I) TER was fixed for all species, (II) TER was variable across the species, (III) the variability in TER for each species was around a common mean response for each trait-environmental driver or trait-trait pair; and the amount of TER variability for each species around those means was different for (IV) each environmental driver only, (V) each trait only or for (VI) each environmental driver and each trait, jointly.
The TTR components of the models investigated four different types of trait-trait interactions by considering that each relationship was (A) fixed for all species, (B) fixed for all species but with added flexibility in the relationship strength implied by the trait-trait covariance matrix, (C) variable across the species, and (D) variable across the species but with added flexibility in the relationship strength implied by the trait-trait covariance matrix.
We started with a baseline model (Model I) that does not account for TTR across the species, assuming that in extreme environments, traits of all constituent species can show identical responses to environmental drivers without any contributions of trait-trait trade-offs across the species. Overall, Model I and the combination of the six TER (I, II, III, IV, V, and VI) and four TTR (A, B, C, and D) model features resulted in a total of 25 traitbased models, tested using the nearly 50,000 tree records and the other trait and environmental data that we collected from the 110 PSPs in the Sundarbans mangrove ecosystem.
Conducting formal model selection among these competing versions of the model, allowed us to understand how TER and TTR can differ across species, to generate species-specific and community-wide trait predictions, and finally to translate these spatial predictions into spatial forecasts of whole-ecosystem productivity. For model selection, we used the Deviance Information Criterion (DIC), which allows several models to be compared at once (Johnson and Omland 2004). Our biologically informed model-building strategy helps us to avoid the risks of "data-dredging" (data analysis with few or no a priori questions) that often leads to overfitted models with a poor predictive ability (Burnham and Anderson 2002).

Models I and I.A-I.D: TER is fixed but TTR is fixed/variable across the species
Models I and I.A-I.D assume that traits of all constituent tree species living in stressed environments show identical responses to environmental drivers because of common, requisite adaptations to physical conditions (Weiher and Keddy 1995, Grime 2006, Duke 2017, Xu et al. 2017. The difference between Model I and the rest of the models is that Model I does not take into account the TTR across the species in defining the trait-environment linkages across the species. These simple models use multiple functional traits (without including species-specific trait identity) as response variables and environmental drivers as predictor variables. We index sites by i ¼ 1, :::, n n ¼ 110 ð Þ , each with an associated environment e i ¼ fe i,1 , :::, e i,k , :::, e i,K g, a vector of K environmental variables for site i. We index species by m ¼ 1, :::, M M ¼ 9 ð Þ: Traits are indexed by j ¼ 1, :::, J J ¼ 4 ð Þ: The vector of traits for each species at each site is denoted t mi = {t mi,1 ,. . .,t mi,j ,. . .,t mi,J }.
We model the joint likelihood of the trait set for a species at a site from a multivariate normal distribution Here, μ i = {μ i,1 ,. . ., μ i , J } is the vector of trait means, that in this model are assumed to be independent of species but functions of the local environment. The J Â J precision matrix Ω is the inverse of the variance-covariance matrix Σ between traits where the variances σ 2 j on the diagonal tell us about the marginal dispersion of each trait individually, whereas the covariances cov jJ À Á tell us about the links between traits j and J. Independent traits will have nearzero covariance, any trade-off between functional traits will be represented by a negative value of the covariance and the existence of synergistic interactions between traits will have a positive covariance. Note that we set Xxxxx 2021 THE FOURTH-CORNER PROBLEM Article e01454; page 5 In Eq. 3, α = {α 1 ,. . ., α J } are the intercepts for each trait, B is a J Â K matrix of regression coefficients, and e i represents the associated environment in the sites. That is, B jk is the coefficient for the regression of the jth trait on the kth environmental covariate. Each coefficient's intercept is generated from a normal prior with a hyperparameter τ for the precision, itself assigned a Gamma prior with a low mean to allow variability between coefficients.
τ TE ∼ Gamma 0:0001, 0:01 ð Þ : For all models after Model I, the precision matrix, Ω, was assigned a Wishart prior distribution (e.g., Horswill et al. 2019). In Model I.A, the distribution has five degrees of freedom (df) Setting the df parameter for the Wishart prior to one plus the dimensions of the variance-covariance matrix (n = 4) attains a uniform prior distribution on the individual correlation parameters; i.e., an equally likely probability between −1 and 1 (Gelman and Hill 2007). The only difference between Model I.A and Model I.B is that in Model I.B, instead of a fixed number of df, we incorporate variability in the effective df by assigning them a Gamma prior (df is constrained to be higher than On ecological time scales, species sorting may cause changes in community composition because species that lack suitable adaptations or trait values are eliminated. On evolutionary time scales, natural selection may cause changes in the traits themselves (Webb et al. 2010). However, selection pressures do not act independently on single traits, but rather, on the species as a whole, whose survival in a specific environment is controlled by multiple interacting traits. Hence, the adaptive value of a particular trait may vary across species depending on the other traits possessed by the species (Verberk et al. 2013).
Therefore, Model I.C and Model I.D incorporate variability in trait-trait relationships (TTR) for different species to capture the effects of taxonomic variability on trait plasticity. Like Models I.A and I.B, here we model the joint likelihood of the set of traits (t mi ) characterizing a species at a particular site using a multivariate normal distribution where, μ i is the vector of trait means and Ω m is the precision matrix that (unlike Models I.A and I.B) is, in this model, assumed to be species specific.
Each Ω m is the inverse of the equivalent covariance matrix Σ m (Eq. 10) that now tells us about the marginal dispersion of each trait of each species individually (i.e., intraspecific trait variability) and the covariances tell us about the interaction levels between the multiple traits of multiple species (i.e., interspecific trait variability). Ω m was assigned a Wishart prior distribution with 5 df (Eq. 11). The structure and descriptions of the linear predictors for the trait set and the priors are as in Eqs. 3, 4, and 5 In Model ID, instead of a fixed df (as was the case in Model I.D), we incorporated a variability in the df by assigning a Gamma prior (Eqs. 12 and 13).
Models II.A-II.D: TER varies for each species, but TTR is fixed/variable across the species Different mangrove species show different anatomical, morphological, and physiological adaptations to persist in the harsh intertidal areas (Naskar andPalit 2015, Mollick et al. 2021). These differential adaptations and resource partitioning may cause different coexisting species to show dissimilar functional responses to environmental drivers, leading to trait divergence (MacArthur and Levins 1967, Chesson 2000, Violle and Jiang 2009, Reef et al. 2010. Therefore, in contrast to Models I.A-I.D, which assume the traits of all constituent tree species show similar responses to environmental drivers, Models II.A-II.D formalize the idea that trait responses to environmental drivers may differ across the species although TTR is fixed in Models II.A-II.B and variable across the species in Models II.C-II.D. In Models II.A-II.B, the joint likelihood of the trait set for a species at a site (t mi ) is a multivariate normal distribution (Eq. 14) with μ mi as the vector of trait means and the precision matrix Ω. For the mth species in the ith site, the linear predictors for the trait set here are where α m is the vector of trait intercepts for the mth species and B m is a J Â K matrix of regression coefficients for that species. Hence, α m,j is the intercept for the jth trait of the mth species and B m,jk is the coefficient for the jth trait of the mth species on the kth environmental covariate. Here, we allow regression coefficients to vary across species. Each coefficient's value is generated from a normal prior with a hyperparameter τ TE (Eq. 16) for the precision and τ TE is then generated from a Gamma prior as in Eq. 5 Because species-specific variability in TTR is not accounted here, the definition of Ω in Models II.A-II.B follows Eq. 6. In Model II.B, we incorporate a variability in the df by assigning a Gamma prior as Eqs. 7 and 8.
As well as considering TER, Models II.C and II.D also incorporate variability in trait-trait relationships (TTR) for different species. Here, the joint likelihood of the trait set and the structure and descriptions of the linear predictors follow Eqs. 14 and 15, except that we incorporate species-specific trait identity to the precision matrix Like II.A-II.B, here we allow regression coefficients to vary across multiple species and the coefficients are generated from a normal prior with a hyperparameter τ TE (Eq. 16) for the precision and τ TE is then generated from a Gamma prior as Eq. 5. Ω m was generated assigning Wishart prior distributions and 5 df (Eq. 11). In Model II.D, we incorporated flexibility in the df by assigning them a Gamma prior as in Eqs. 12 and 13.
Models III.A-III.D: TER for each trait-environmental driver pair varies around its own common mean, but TTR is fixed/variable across the species In extreme environments such as mangrove forests, species can usually tolerate a certain level of stress via specialization . Hence, it is likely that all species show an average functional response to environmental filters (Naskar and Palit 2015). Nevertheless, the strength of the response around that average may vary for different species and for different traits depending on how severe the stress conditions are and how one trait facilitates or constrains another (Reef et al. 2010). We incorporated these ecological hypotheses into models III.A-III.D to VI.A-VI.D. Models III.A-III.D tested the general hypothesis that, although the mangrove species differ in their response to the environmental drivers regarding their traits (assumption of Models II.A-II.D), variability is around a common interspecific average response. However, TTR is considered fixed in Models III.A-III.B and variable across the species in Models III.C-III.D.
We model the joint likelihood of the trait set for a species at a site (t mi ) as Eq. 14 for Models III.A-III.B, and as Eq. 17 for Models III.C-III.D, while the definition and description of the linear predictors (μ mi ) for all these models follow Eq. 15. The difference from Model II.A-II.D is that here we define a common interspecific average trait response for the species by assigning a normal prior with hyperparameters B jk and τ TE (Eq. 18). B jk is then generated from a normal prior with another precision hyperparameter τ TE (Eq. 19), which is then generated using a Gamma prior (Eq. 20). τ TE is generated from a Gamma prior as Eq. 5 τ TE ∼ Gamma 0:0001, 0:01 ð Þ : To incorporate the fixed TTR, the precision matrix Ω is defined as in Eq. 6. in Model III.A, Eqs. 7 and 8 in Model III.B, and to incorporate the variable TTR across the species, Ω m is defined as Eq. 11 in Model III.C and as Eqs. 12 and 13 in Model III.D.
Models IV.A-IV.D: TER varies around the common mean for each environmental driver, but TTR is fixed/variable across the species Traits associated with extreme environmental stress may show convergence, while traits associated with resource partitioning and interactions may show divergence (Flowers et al. 2010, Weiher et al. 2011. To incorporate these fundamental assumptions, Models IV.A-IV.D include not only the possibility that the amount of variability in TER around the average trait response (Models III.A-III.D) vary for each environmental driver, but also assumes that the TTR is either fixed (IV.A-IV.B) or different (IV.C-IV.D) for each species.
Like Models III.A-III.D, here we model the joint likelihood of the trait set for a species at a site (t mi ) as Eq. 14 for Models IV.A-IV.B and as Eq. 17 for Models IV.C-IV.D, while the structure and description of the linear predictors (μ mi ) for all these models follow Eq. 15. The Ω in Model IV.A is defined as Eq. 6 while for Model IV.B it follows Eqs. 7 and 8. To incorporate speciesspecific variability in TTR, Ω m is defined as Eq. 11 in Model IV.C and as Eqs. 12 and 13 in Model IV.D.
The only difference from Models III.A-III.D is that here we account for environmental-driver-specific variability in defining the common interspecific average trait response for the species by assigning a normal prior with hyperparameters B jk and τ TE,k (Eq. 21). B jk is generated using Eqs. 19 and 20 while a Gamma prior is assigned for τ TE,k (Eq. 22) τ TE,k ∼ Gamma 0:0001, 0:01 ð Þ 8 k: Models V.A-V.D: TER varies around the common mean for each trait, but TTR is fixed/variable across the species Models V.A-V.D test whether the amount of variability around the average trait response of the species is different for each trait and also assume a fixed TTR (V.A-V.B) or a variable TTR (V.C-V.D) across the species. Like Models IV.A-IV.D, we model the joint likelihood of the trait set for a species at a site (t mi ) as Eq. 14 for Models V.A-V.B, and Eq. 17 for Models V.C-V.D. The definition and description of the linear predictors (μ mi ) for all these models follow Eq. 15. The Ω in Model V.A is defined as Eq. 6 while for Model V.B it follows Eqs. 7 and 8. Ω m is defined as Eq. 11 in Model V.C and as Eqs. 12 and 13 in Model V.D.
Compared to IV.A-IV.D, here we replace the hyperparameter τ TE,k with the hyperparameter τ TE,j (Eq. 23) to account for trait-specific variability in defining the common interspecific average trait response for the species. A Gamma prior is assigned to τ TE,j (Eq. 24) Models VI.A-VI.D: TER varies around the common mean for each trait and environmental driver, but TTR is fixed/variable across the species Models VI.A-VI.D combine the concepts of Models IV.A-IV.D and Models V.A-V.D and test whether the amount of variability around the average trait response is also different for each trait [assuming not all traits (i.e., leaf and stem traits) of a species show a similar degree of sensitivity under stress (Naskar and Palit 2015)] and each environmental driver (assuming not all environmental drivers (i.e., abiotic and biotic) modify the trait values equally ).
Like Models IV.A-IV.D and Models V.A-V.D, we model the joint likelihood of the trait set as for Eq. 14 for Models VI.A-VI.B, and Eq. 17 for Models V.C-V.D. The linear predictors for Models VI.A-VI.D follow Eq. 15. The precision matrix Ω in Model VI.A is defined as Eq. 6 while for Model VI.B it follows Eqs. 7 and 8. Ω m is defined as Eq. 11 in Model VI.C and as Eqs. 12 and 13 in Model VI.D.
The difference in Models VI.A-VI.D from Models IV.A-IV.D and Models V.A-V.D is that here we incorporate both trait-and environmental driver-specific variability in defining the common interspecific average trait response for the species by assigning a normal prior with hyperparameters B jk and τ TE,jk (25). The definition of B jk follows Eq. 19. The new hyperparameter τ TE,jk that replaces τ TE,jk (Models IV.A-IV.D) and τ TE,jk (Models V.A-V.D) is assigned a Gamma prior (26) τ TE,jk ∼ Gamma 0:0001, 0:01 ð Þ 8 j, k:

Model fitting and selection
We fitted the models using Markov Chain Monte Carlo (MCMC) in the symbolic language JAGS version 4.2.0 (Plummer 2015) with the runjags interface version 2.0.2-8 (Denwood 2016) in software R version 3.6.0 (R Core Team 2020). JAGS code for the models is included in Data S1:1.1_Models.R. To evaluate the convergence of the two MCMC chains, we inspected the trace plots and ensured that the potential scale reduction factor (Gelman-Rubin statistic) for all monitored parameters was <1.05 (Gelman et al. 2004). For reducing the chance of any relic effects of initial values, we discarded the first 5,000 samples as burn-in and then simulated two independent MCMC chains for 60,000 iterations. For assessing the goodness-of-fit of the best model, we used R 2 values from regressions of observed vs. predicted traits.
Model selection in the Bayesian framework is complex, computationally challenging and less automatic for Bayesian hierarchical models than for the simpler, likelihood-based generalized linear models (Tenan et al. 2014, Hooten and Hobbs 2015, Broms et al. 2016, Han 2017. Although there are many choices available for Bayesian model selection (cross-validation, information criteria, regularization, etc.), each has its own strengths and limitations (Hooten and Hobbs 2015). In this study, we compared the performance of the models using the DIC (Deviance Information Criterion; Spiegelhalter et al. 2002), a widely used Bayesian analog of the Akaike Information Criterion for likelihoodbased models. We centered and standardized the response and predictor variables to improve model convergence, to facilitate comparisons between the strength of different effects and to allow the intercepts to be interpreted as true baseline responses. We obtained posterior distributions and 95% credible intervals (CIs) for the model parameters.

Mapping traits and productivity
Using kriged surfaces (Sarker et al. 2016) of the environmental variables, the best Bayesian model predicted functional trait distributions for the four most prominent species (that, together, contributed 98% of the total count). Using previously derived species distribution models and maps by Sarker et al. (2016), we averaged the summed posterior trait values of each species weighted by their relative abundance in each grid cell to produce community trait maps for the entire ecosystem. In the species distribution models, the response of the species abundances to environmental drivers was modeled using a likelihood-based generalized additive modeling framework. At first, all possible candidate models with all possible combinations of predictor variables were fitted and ranked by the Akaike Information Criterion (AIC). Then a set of competing models were determined using the ΔAIC ≤ 2 criterion (Burnham and Anderson 2002) to determine a set of competing models for each species. To inspect the relative support for each model in the confidence set, Akaike weights (AICw; Burnham and Anderson 2002) were calculated. When there was only one model with ΔAIC ≤ 2, it was clear that it outperformed all possible candidate models. However, in the case of multiple competing models, AIC-weighted model averaging was performed to reduce model selection uncertainty and bias. Finally, the averaged parameter estimates were used to produce species density maps based on the interpolated surfaces of the environmental variables.
To develop whole-ecosystem biomass productivity maps, we first built species-specific trait-based productivity models and then linked these models with the species traits and density maps to generate productivity predictions. Here, we modeled the biomass production of individual trees (taking basal area as a proxy) as a linear function of the functional traits. This posthoc model was fitted using a Bayesian GLM framework. We fitted all possible candidate models with all possible combinations of predictor variables for each species using the bayesx function in R2BayesX package version 1.1-1 (Umlauf et al. 2015) in R. For each of the model, 12,000 MCMC iterations were carried out after a burn-in sample of 2,000. The convergence of the MCMC chains was checked using the potential scale reduction factor (Gelman-Rubin statistic) for all monitored parameters. The performance of the traitbased productivity models for each species was compared using the DIC (Appendix S1: Table S2). Then, using the posterior trait surfaces of the species, the best model predicted productivity for the four most dominant species. Finally, using the species density surfaces (Sarker et al. 2016), we summed the posterior productivity values of each species weighted by their relative abundance in each grid cell to produce productivity maps for the entire ecosystem.

Forecasting traits and productivity
The detrimental effects of salinity intrusion and siltation on mangrove species distributions and diversity are well established (Sarker et al. 2016(Sarker et al. , 2019b. Moreover, biotic homogenization has been underway in the Sundarbans (Sarker et al. 2019a). The SLR rate along the Bangladesh coast in the 20th century (5.93 mm/yr) was substantially higher than the global average (1.0-2.0 mm/yr). The National Adaptation Program of Action has projected 32 cm of SLR by 2050 (Karim and Mimura 2008). Sundarbans is a sea-dominated delta, where freshwater river flows help to modulate salt-water toxicity and keep the ecosystem suitable for mangrove trees. However, the freshwater supply from the transboundary rivers into the Sundarbans has considerably declined (3,700 m 3 /s to 364 m 3 /s) since the construction of the Farakka dam (1974) in India (Mirza 1998). As a result, the average soil salinity level has already increased by 60% since 1980 (Aziz and Paul 2015). Several geomorphological or hydrological prediction models (IWM 2003, Wahid et al. 2007, Mukhopadhyay et al. 2019 have projected a 5-10% decadal increase in salinity. The river network annually transports about 2.4 billion metric tons (or Mg) of sediments in the Sundarbans (Mitra and Zaman 2016) and sediment burial of aerial roots (inhibits root aeration) is a major reason for limited mangrove growth and regeneration failure (Mitra andZaman 2016, Sarker et al. 2019b). Assuming a continuation of such historical pressures that may further affect plant functions via physiological stress and shifting relationships among species, we updated the environmental and species distribution maps (Appendix S1: Fig. S2) for five future stress scenarios (E1, 10%; E2, 20%; E3, 30%; E4, 40%; and E5, 50% increase in salinity and siltation), and then used them to forecast traits and forest productivity under these novel environments. Based on the available geomorphological or hydrological model-based forecasts of salinity and siltation in the literature (IWM 2003, Wahid et al. 2007, Islam and Gnauck 2011, Ghosh et al. 2019, Mukhopadhyay et al. 2019), we considered E5 as the worst stress scenario for the ecosystem by 2050.
Plant trait values that enhance species growth and establishment, such as higher tree canopy height (indicating higher carbon sequestration capacity via its relationship with tree biomass; Rahman et al. 2015) and higher SLA (indicating increased mass-based lightsaturated photosynthetic rate; Pérez-Harguindeguy et al. 2013) contribute to higher biomass productivity in mangrove forests (Komiyama et al. 2008. In contrast, trait values that enhance salinity and drought tolerance, such as higher WD (indicating less hydraulic failure; Choat et al. 2018) and higher LS (implying an accumulation of larger amounts of water and dissolved ions in leaves under salinity stress; Naskar and Palit 2015) lead to lower biomass productivity in mangrove forests (Santini et al. 2012, Mitra 2013. Therefore, a decrease in height and SLA and an increase in WD and LS together indicate reduction in biomass productivity. Specifically, higher height and SLA values and lower WD and LS values together indicate low stress, while lower height and SLA values and higher WD and LS values together indicate high stress (Larjavaara and Muller-Landau 2010, Pérez-Harguindeguy et al. 2013, Hoeber et al. 2014. To evaluate the robustness of these forecasts, we used samples from the joint posterior distribution of all model parameters to simulate traits and productivity distributions for both current and future environmental scenarios. Using these simulations, we then mapped the posterior probability of deterioration of traits (a decrease in height and SLA values and an increase in WD and LS values) and productivity at each grid cell, averaged across all grid cells in the ecosystem, in response to each of the future stress scenarios, to uncover the uncertainty related to our forecasts.
In addition, to understand the relative contribution of the changing trait values and shifts in species densities on the projected productivity of the ecosystem, we predicted, combined, and mapped the productivity of the individual species (through the best productivity models of the species) (1) by considering the forecasted species density surfaces only (changes in trait values were not considered), (2) by considering the forecasted values of the traits only (shifts in species densities were not considered), and (3) by considering both forecasted traits and species density surfaces for 2050. We also calculated the proportion of grid cells across the whole ecosystem where the expected response is a deterioration of productivity when we consider shifts in trait values and species densities separately. All maps were generated using the raster package version 2.4-18 (Hijmans 2017) in R.  Fig. 2 shows the relative performance of the traitbased models in terms of the deviance information criterion (DIC). Model I that did not account for the TTR across the species and only focused on the TER was the weakest performer, indicating the importance of incorporating trait plasticity information in the trait-based models. Models (IC, ID, IIC, IID, IIIC, IIID, IVC, IVD, VC, VD, VIC, and VID) that incorporated variability in trait-trait relationships (TTR) for different species to capture the effects of taxonomic variability on trait plasticity performed better than the models (IA, IB, IIA,  IIB, IIIA, IIIB, IVA, IVB, VA, VB, VIA and VIB) that assumed identical TTR across the species. Models (II.A-II.D) that assumed variable TER across the species performed better than the models (I.A-I.D) that assumed identical TER across the species. However, the performances of the models substantially improved when we assumed all species living in extreme environments show an average functional trait response to environmental filters (III.A-III.D), and the strength of the response around that average was allowed to vary for each environmental driver (IV.A-IV.D), each trait (V.A-V.D), or both (IV.A-IV.D).

Model comparison summary
Model IV.D (DIC = 4954.41) was identified as the preferred model, showing evidence that tree species growing in a stressed ecosystem like a mangrove forest have a typical functional response to the environmental filters with inter-specific variability around this average, and the amount of variability is further contingent upon the nature and magnitude of the filters. We further observed interspecific variability in trait-trait relationships across the species.
The explanatory power of the best model was high (regressions of observed vs. predicted trait values yielded the following: canopy height R 2 = 0.68, specific leaf area R 2 = 0.64, wood density R 2 = 0.94, and leaf succulence R 2 = 0.67) and plots of predictions vs. observations indicate that the model made unbiased predictions (Appendix S1: Fig. S3). The potential scale reduction factor (Gelman-Rubin statistic) for all monitored regression parameters was <1.05, indicating model convergence.

Drivers of functional traits
Soil salinity and siltation are the dominant environmental drivers that limit species canopy height (Fig. 3, Appendix S1: Table S3). However, the intensity of such detrimental effects varies across the species. Salt-induced reduction in canopy height is highest in the climax tree species (H. fomes) of the ecosystem (posterior mean μ = −0.57, 95% CI −0.68 to −0.46). In contrast, increasing URP (i.e., more available freshwater) favors height increase in H. fomes (μ = 0.14, 95% CI 0.08-0.24), and also in the late-successional species X. mekongensis (μ = 0.13, 95% CI 0.05-0.21) and the generalist species E. agallocha (μ = 0.09, 95% CI 0.03-0.15). Soil pH, salinity, siltation, and K have limiting effects on specific leaf area (SLA) of many species. In contrast, we identified a strong positive response of SLA to URP and NH 4 , suggesting that nitrogen-rich upstream habitats favor plants' resource acquisition in tropical coastal environments.
While traits related to resource acquisition (height and SLA) show a strong negative response to increasing salinity, siltation, and pH, we discovered a clear positive response of the traits related to resource conservation (WD and LS) to these dominant environmental filters (Fig. 3, Appendix S1: Table S3), suggesting a substantial growth-survival trade-off across the species. Nevertheless, the strength of this positive response or positive filtering (adaptive shifts in trait values to increase resistance to stress) varies across the species, trait types, and environmental drivers. For example, salinity is the dominant positive filter for H. fomes (μ = 0.18, 95% CI 0.14-0.22), siltation for B. sexangula (μ = 0.10, 95% CI FIG. 3. Effects of environmental drivers (URP, upriver position) on the functional traits (height; canopy height; SLA, specific leaf area; WD, wood density; LS, leaf succulence). Each circle represents a species. Species are ordered (left to right, then top to bottom) based on their overall abundance (e.g., E. agallocha 1 is the most abundant and S. apetala 9 is the least abundant species). The colors of the circles indicate the type and strength of the trait response (red for negative and blue for positive responses, see scale on right). A dashed circle around all nine species indicates that the average species trait response to this environmental driver is strong (negative or positive, i.e., 95% credible intervals do not include 0). Posterior means and 95% credible interval values for the parameters are presented in Appendix S1: Table S3. Full species names are Heritiera fomes, Excoecaria agallocha, Ceriops decandra, Xylocarpus mekongensis, Amoora cucullata, Avicennia officinalis, Bruguiera sexangular, Cynometra ramiflora, and Sonneratia apetala.

Covariation among the functional traits
SLA and LS show negative associations in most species (Table 2), and for H. fomes we see broader tradeoffs between traits related to resource acquisition (height and SLA) and traits related to resource conservation (WD and LS). We also see an apparent positive relationship between the resource acquisition-related traits in A. officinalis, and between the resource conservation related traits in X. mekongensis and A. cucullata.

Trait and productivity maps
Our approach offers species-specific as well as community-wide spatial trait predictions under current and future stress scenarios regarding salinity intrusion and siltation, and links them with the trait-based productivity models to forecast ecosystem productivity. It further evaluates the robustness of such forecasts by propagating parameter uncertainty through sampling from the joint posterior distribution of all model parameters and simulating traits and productivity distributions for both current and future stress scenarios.
Our spatial maps (Fig. 4) reveal that the most productive tree communities are currently distributed in the freshwater-dominated eastern and northern Sundarbans. If historical increases in salinity and siltation are maintained, our model predicts a substantial productivity loss over the entire ecosystem under all future stress scenarios (10%, 20%, 30%, 40%, and 50% increase in salinity and siltation; Appendix S1: Fig. S8) although there is the possibility of local productivity gains in a few of the communities, particularly under the least stress future scenarios (Appendix S1: Fig. S9C). The worst stress scenario (a 50% rise in salinity and siltation) is predicted to cost the ecosystem 29.5% of its current total productivity by 2050 (Appendix S1: Fig. S8A). Our analysis on the relative contributions of the shifts in species densities and changes in species trait values to this projected productivity loss revealed that shifts in species densities alone would cause a decrease of productivity in 62% of the grid cells (Appendix S1: Fig. S10B). However, changes in trait values alone would cause a loss of productivity in all areas of the ecosystem (Appendix S1: Fig. S10C), indicating the importance of integrating trait plasticity information into the trait-based modeling frameworks.
The eastern and northern regions currently support the tallest mangrove communities. However, if the habitat degradation trend is maintained, it is highly likely that these will turn into dwarf communities, following an average of 36% height loss by 2050 (Appendix S1: Figs. S8B, S11). In turn, we predict a community-wide  Here, we defined significance as the event of the 95% credible intervals not including 0. Full species names are Heritiera fomes, Excoecaria agallocha, Ceriops decandra, Xylocarpus mekongensis, Amoora cucullata, Avicennia officinalis, Bruguiera sexangular, Cynometra ramiflora, and Sonneratia apetala.  Vol. 0,No. 0 increase in the values of LS (14%) and WD (5%), indicating highly plastic responses of the mangrove species to ensure efficient water use and mechanical support under the extreme stress scenarios, by 2050. Community SLA shows more stable patterns between time points. Our species-specific trait predictions (Appendix S1: Figs. S4-S7) further determine that, although every constituent species of the ecosystem will lose height under the future stress scenarios, the loss will be substantially higher for the climax species H. fomes (47%) by 2050 (Appendix S1: Figs. S8C and S12). The increase in WD (9%) is also highest for this species. In turn, the disturbance specialist species, C. decandra, which is expected to undergo ecosystem-wide range expansion by 2050, loses the least amount of height and SLA, compared to others.

DISCUSSION
Tropical forest ecosystems are shaped by the plastic responses and complex interactions between primary producers played out on the stage set by different environmental drivers. These are interconnected processes that collectively (through forest species composition and individual species function) determine fundamental functional aspects of the entire forest, such as its productivity. Looking at subsets of this causal network (e.g., statistical relationships between single traits and the environment) can only reveal part of the biological signal in the data. By underrepresenting the complexity of the connections between traits, species, and environment (a challenge known as the "fourth-corner problem"), any predictive model of future forest dynamics stands to suffer from increased imprecision. Conversely, by failing to propagate inherent parameter uncertainty through to the prediction stages, we risk ascribing unjustifiable confidence in any such predictions. In this paper, we have developed an analytical approach that follows a holistic (but still phenomenological) view of forest ecosystems, both by developing a model as complex as the data would support and by fully propagating parameter uncertainty into the predictions. This approach offers a solution to the fourth-corner problem, contributes to a data-driven clarification of different ecological assumptions about how complex forest ecosystems function, uncovers the environmental drivers of functional traits, and generates new, testable predictions about how traits and forest productivity will be affected by environmental change in the future.
Our integrated Bayesian approach allowed different trait-based ecological hypotheses to be tested in a single analysis simultaneously using all available data. In extreme environments, such as the tropical intertidal zones, tree communities may comprise functionally redundant species (Ricotta et al. 2016). Therefore, our studied species might have shown homogeneous functional responses to the stressors. Models I and I.A-I.D (Fig. 2), which represent this expectation, received the least support in our analyses. Rather, we found better support for Models II.A-II.D, suggesting that resource partitioning (Twilley andRivera-Monroy 2005, Reef et al. 2010) and specialization via different morphological, anatomical, and physiological adaptations (Naskar and Palit 2015, Chowdhury et al. 2016a, Rahman et al. 2020) may cause different mangrove species to show disparate functional trait responses to the environmental drivers. However, Model IV.D, which considers that all species show an average functional trait response to the environmental drivers and the amount of variability around the average functional trait response depends on the nature and magnitude of the drivers, received the most support from our data. In addition, the performances of each of the models improved when we accounted for the interactions between multiple traits across the species, thus confirming the importance of considering both intra-and interspecific trait variations in trait-based models.
Applying our approach to the dynamic Sundarbans ecosystem, we determine salinity, siltation, pH and upriver position as the key drivers governing the plants' trait responses. Further, we discovered diverging trait responses among the species and the strength of such responses to the environmental drivers were found to vary across the trait categories (Fig. 3, Appendix S1: Table S3). For example, among the traits, height is more affected by salinity than the other drivers and the magnitude of the effect is higher on the climax species H. fomes than the generalists E. agallocha and C. decandra. These species-specific responses of canopy height (and other traits) could be related to the variability in water stress, nutritional imbalance, and salt stress (Rasool et al. 2013, Gustafsson andNorkko 2019).
Simultaneous plastic decreases in traits related to resource acquisition (height and SLA) and plastic increases in traits related to resource conservation (WD and LS) along the soil salinity, siltation, and pH gradients, suggest conservative resource use (low resource acquisition and turnover rates and slow growth rates) as the dominant ecological strategy in many species to ensure their long-term persistence under unfavorable environments (Rosado et al. 2016). Traits supporting mangroves resource acquisition and traits supporting resource conservation showed opposing responses to increasing upriver position (i.e., more available freshwater), suggesting that the same species can increase resource acquisition and growth under benign environmental conditions.

Trait and productivity maps
Our model predicts an ecosystem-wide productivity loss under all future stress scenarios and forecasts a 29.5% loss of its current total productivity under the highest stress scenario by 2050 (Appendix S1: Fig. S8A). This productivity loss is mainly associated with community-wide decreases (36%) in height, by 2050 (Appendix S1: Fig. S8B). Our analysis further revealed that shifts in species distributions and densities alone would cause a decrease of productivity in 62% of the grid cells (Appendix S1: Fig. S10B), while plastic changes in trait values alone would cause a loss of productivity in all areas of the ecosystem (Appendix S1: Fig. S10C). The discrepancy between the relative contributions of shifts in species densities and shifts in species trait values to the forecasted productivity loss demonstrates the importance of incorporating trait plasticity information while making trait-based predictions. This also emphasizes the importance of showing extreme caution when using species distributions as a proxy for fitness, particularly in highly dynamic ecosystems such as the Sundarbans (Matthiopoulos et al. 2015(Matthiopoulos et al. , 2020. Trait maps (Fig. 4) reveal that community SLA and WD show more stable patterns between the stress scenarios. This could be the result of turnover in species composition rather than stasis in stress conditions (Appendix S1: Fig. S2). For example, while the climax species, H. fomes (medium SLA and high WD), is projected to lose most of its suitable habitat by 2050, both generalist, E. agallocha (high SLA and low WD), and disturbance specialist species C. decandra (low SLA and high WD), are likely to replace most of the degraded habitats and are projected to undergo ecosystem-wide range expansions.
Our species-specific predictions further uncover spatial variability in plastic responses of the traits under projected stress conditions. We find that, under increasing stress, every species will lose canopy height and SLA, but the loss will be substantially higher for the climax species H. fomes (52% height and 16% SLA loss by 2050; Appendix S1: Figs. S8C, S12), which is in agreement with the previous results that the species has narrower salt tolerance than the generalists because of its limited ability to balance water and salt uptake (Das 1999). Plastic increases in WD and LS values are also predicted to be highest for H. fomes, implying maximum stress on the species contributing the most biomass.

Challenges and prospects of improvement
While our approach incorporated multiple species, in common to other existing approaches (Dolédec et al. 1996, Legendre et al. 1997, Lavorel et al. 2007, Dray and Legendre 2008, Jamil et al. 2012, Pollock et al. 2012, Brown et al. 2014, Dray et al. 2014, Warton et al. 2015, it does not consider competitive interactions between species. Species-species interactions can lead to trait divergence or sometimes trait convergence when traits have a strong influence on the competitive ability of the coexisting species (Grime 2006), which in turn may affect the species composition and biomass productivity in forest communities. However, by allowing the flexibility to incorporate physiologically informed priors based on sound biological first principles (Kindsvater et al. 2018), our modeling framework can inform more mechanistic modeling approaches of tree dynamics (Zavala and Bravo de la Parra 2005, Purves et al. 2008, Adler et al. 2013, Chauvet et al. 2017, Lines et al. 2020) that will better capture tree interaction dynamics.
We acknowledge that our empirical modeling approach does not incorporate the physiological mechanisms giving rise to trait trade-offs and is unable to account for the connections between functional traits and the recognized coexistence mechanisms such as spatial or temporal environmental heterogeneity, resource partitioning and natural enemies (Adler et al. 2013). We must therefore stop short of speculating about the mechanisms that generate niche partitioning (MacArthur 1958, Chesson 2000 and environmental filtering (Weiher and Keddy 1995) and thus how traits mediate plant community assembly and coexistence. Nevertheless, by explicitly analyzing the hierarchical relationships between traits and environmental filters at the species, community, and ecosystem levels, and by simultaneously quantifying the intra-and interspecific trait-trait tradeoffs, our approach can better inform purely mechanistic trait-based tests of coexistence mechanisms (Adler et al. 2013, Kraft et al. 2015a,b, D'Andrea and Ostling 2016. We recorded trait measurements for only three individuals of each species at each Permanent Sample Plot (PSP). Inclusion of more individuals for each species at each site may offer better insights into within-species variations in traits. We emphasize that future local scale individual species-based studies focusing on the anatomical and physiological traits should consider a bigger sample size.
In addition, several fundamental unknowns may also interfere with our forecasts. The first is the climate-or local stress-induced trait acclimation and adaptation that could potentially affect the accuracy of the projected trait values and productivity. However, variability in acclimation potential among different tree species and different trait types have yet to be discovered in mangrove ecosystems and elsewhere (Van Bodegom et al. 2014). The second unknown is the degree to which environment-trait relationships and trait-trait interactions will remain the same under future climatic and stress conditions (i.e., whether they undergo abiotic/biotic filtering).

Practical applications
Measurable and transparent forest conservation targets require an accurate a priori determination of the environmental stressors limiting forest growth and development (Lewis 2005, Sarker et al. 2019b). Our integrated approach can also provide forest managers with a muchneeded functional basis to tailor science-driven management and conservation actions, particularly for conservation priority areas like the Sundarbans. Inadequate knowledge of how different mangrove species respond in dynamic coastal environments has previously led to unsuccessful conservation efforts in many tropical regions (Lewis 2005), including the Sundarbans (Islam et al. 2014). Therefore, determining the environmental drivers that limit mangroves' growth and establishment is crucial before designing and implementing habitat restoration and reforestation programs. For example, our finding that soil salinity, alkalinity and siltation are the most limiting environmental drivers of mangroves' primary growth and resource acquisition, especially for the climax species, can target the selection of species and sites for future reforestation and restoration initiatives. Our finding that currently the productivity hotspots are confined to the hyposaline eastern and northern regions in the Sundarbans and the projected 29.5% decline in productivity, lends strong support for immediate habitat enhancement programs (e.g., river dredging for increasing freshwater flows and for reducing silt deposition) in the eastern and northern regions. Most importantly, the productivity maps can be important tools for implementing the United Nations REDD+ (Reduced Emissions from Deforestation and Degradation) initiatives and guiding national strategy while negotiating for payments for ecosystem services (Rahman et al. 2015).
Our hierarchical Bayesian model can simultaneously estimate the parameters associated with TER and project the fitted TER to make spatial predictions for multiple traits of multiple species. By estimating TER parameters and providing trait predictions simultaneously, our integrated approach provides automatic estimates of parameter uncertainty from model fitting through to the prediction stages, thus preventing us from assigning unjustifiable confidence in any such predictions. Such an explicit representation of our prediction uncertainties will help the resource managers in taking appropriate reforestation and rehabilitation programs. In addition, this simultaneous modeling process can offer greater computational efficiency than many of the current approaches focusing on a single trait of a single species at a time and talking pair-wise trait combinations to infer TER (Wüest et al. 2018). This high computational efficiency of our integrated approach offers the flexibility to incorporate more traits, species, and environmental variables. Our current regional data set comprises data for four quantitative traits from 20 mangrove species (out of 69 mangrove species worldwide) and 110 sites. Like established global trait databases (e.g., TRY [Kattge 2020] and GLOPNET [Reich et al. 2009]) and mangrove-centric trait database (Quadros and Zimmer 2017), this data set includes multiple measurements for the same trait, species, and site, and follows the standard trait measurement protocols (Pérez-Harguindeguy 2013). Our easily measurable trait set (height, SLA, wood density, and leaf succulence) are common in all tree species whatever the scale is (Stahl et al. 2014) and also available for different forest ecosystems. For example, the Tundra Trait Team database (Bjorkman et al. 2018b) contains observations of 18 plant traits on 978 tundra species and the most frequently measured traits (>1,000 observations each) include plant height, specific leaf area, leaf fresh, and dry mass, leaf dry matter content, and wood density.
Based on the strength of simultaneous modeling and predictions, high computational efficiency and data requirement commonalities, our integrated Bayesian approach offers a generalizable foundation for powerful modeling of trait-environment linkages under changing climate and for predicting their consequences on ecosystem functions and services in other critical forest ecosystems of the world such as tundra, the most rapidly warming biome on the planet (Bjorkman et al. 2018a). The increasing contribution of collaborators to the global trait databases (Kattge et al. 2020) and fast expansion of the spatially explicit global data sets for species distributions (Hudson et al. 2014), and high-resolution environmental layers (Shangguan et al. 2014, Basher et al. 2018 provide opportunities for such applications, especially within dynamic global vegetation models (Scheiter et al. 2013, Berzaghi et al. 2020).

CONCLUSIONS
Our integrated Bayesian approach provides the foundations for trait-based predictions in plant ecology through simultaneously modeling trait-trait and traitenvironment relationships (for multiple traits, species, and environmental drivers) at organismal, community, and ecosystem levels, thus resolving many fundamental methodological problems of the existing ordination/permutation and univariate approaches, and bridges community, ecosystem ecology, and functional biogeography. These novel developments allow us to (1) integrate different ecological assumptions about how traits interact with the environment into the modelbuilding process, (2) quantify the strength of the traitenvironment relationships, (3) identify the trait-trait trade-offs across the species, (4) predict species-specific and community-level trait distributions with accurate estimates of uncertainty, and (5) forecast wholeecosystem productivity under future environmental conditions. Applying this approach to the Sundarbans, we discovered (1) tree species growing in a stressed ecosystem like a mangrove forest have a typical functional response to the environmental filters with interspecific variability around this average, and the amount of variability is further contingent upon the nature and magnitude of the filters, (2) substantial intraspecific trade-offs among the functional traits in many tree species, (3) serious detrimental effects of increasing salinity, siltation, and soil alkalinity on functional traits associated with plants' resource acquisition and growth, (4) plastic enhancement of traits related to stress tolerance, indicative of plants prioritizing persistence over growth, and (5) ecosystem-wide drop in biomass productivity under all anticipated stress scenarios with an average of 29.5% productivity loss by 2050 in the worst scenario. These findings advance our understandings of how species living in stressed ecosystems respond to environmental change and provide policymakers and managers with a much-needed functional basis for developing strategic approaches and setting targets for forest conservation, restoration, and ecosystem management. Because of its generality and high computational efficiency, our integrated approach is applicable to any ecosystem across broad global scales, to predict shifts in ecosystem functioning and services under a rapidly changing climate.