Trait Response to Nitrogen and Salinity in Rhizophora mangle Propagules and Variation by Maternal Family and Population of Origin

Many coastal foundation plant species thrive across a range of environmental conditions, often displaying dramatic phenotypic variation in response to environmental variation. We characterized the response of propagules from six populations of the foundation species Rhizophora mangle L. to full factorial combinations of two levels of salinity (15 ppt and 45 ppt) reflecting the range of salinity measured in the field populations, and two levels of nitrogen (N; no addition and amended at approximately 3 mg N per pot each week) equivalent to comparing ambient N to a rate of addition of 75 kg per hectare per year. The response to increasing salinity included significant changes, i.e., phenotypic plasticity, in succulence and root to shoot biomass allocation. Propagules also showed plasticity in maximum photosynthetic rate and root to shoot allocation in response to N amendment, but the responses depended on the level of salinity and varied by population of origin. In addition, propagules from different populations and maternal families within populations differed in survival and all traits measured except photosynthesis. Variation in phenotypes, phenotypic plasticity and propagule survival within and among R. mangle populations may contribute to adaptation to a complex mosaic of environmental conditions and response to climate change.


INTRODUCTION
Many plant species thrive across an extensive range of environmental conditions, often displaying dramatic phenotypic variation (McKee, 1995;Smith and Snedaker, 1995;Richards et al., 2005;Feller et al., 2010). This is particularly true in coastal ecosystems that are characterized by temporal cycles and spatial variation in tidal inundation, temperature, nutrient availability, and salinity (Pennings and Bertness, 2001;Krauss et al., 2008;IPCC, 2014;Proffitt and Travis, 2014;Wuebbles et al., 2014). In addition to the naturally dynamic nature of coastal habitats, anthropogenic activities can increase the input of nutrients and alter watersheds, further contributing to environmental variation (Bertness et al., 2002;Barbier et al., 2008;Gedan et al., 2009Gedan et al., , 2011Crotty et al., 2020). Within these dynamic conditions, foundation plant species such as mangroves provide ecosystem services. These services include providing habitat for many juvenile fish species, biotic filtering of pollutants, and buffering of storms (Ellison et al., 2005;Zedler and Kercher, 2005;IUCN, 2007;Alongi, 2008Alongi, , 2013Costanza et al., 2008;Gedan et al., 2011;Bertness, 2020). Foundation species are defined not only as those that dominate a community assemblage numerically or in biomass, but they also determine diversity of associated taxa through a variety of interactions (Ellison, 2019). Further, foundation species modulate fluxes of nutrients and energy in their ecosystem (Ellison, 2019). Hence, these species disproportionately contribute to maintaining habitat integrity and ecosystem resilience (Bertness and Callaway, 1994;Keith et al., 2017;Ellison, 2019;Bertness, 2020;Qiao et al., 2021). Understanding how these species cope with challenges from anthropogenic impacts is key to preserving the ecosystems they create and define (Gedan et al., 2009(Gedan et al., , 2011Guo et al., 2021).
Understanding the mechanisms of response in coastal foundation species has become increasingly important for conservation and management strategies as these species must cope with rising sea levels, increased warming, and anthropogenic disturbances (Gedan et al., 2011;Kirwan and Megonigal, 2013;Osland et al., 2013Osland et al., , 2017. The Food and Agriculture Organization of the United Nations (FAO) estimates that as much as 35% of global mangrove forest habitat has been destroyed since the early 1980's for the development of human settlements, agriculture and aquaculture, and industrial shipping harbors, although the rate of loss appears to have slowed in the last decade (Food and Agriculture Organization of the United Nations, 2007; Polidoro et al., 2010;Ellison et al., 2015;FAO, 2020). In some regions, mangrove trees are also harvested for wood and charcoal (Ellison et al., 2015), resulting in habitat fragmentation and isolation of existing remnant fragments (Friess et al., 2012;Haddad et al., 2015). The resultant loss of diversity could pose risks for these coastal foundation species in the future, particularly as sea levels are projected to rise between 0.2 and 2 m over the next century due to anthropogenic climate change (Melillo et al., 2014).
The vulnerability of coastal foundation plant communities to global change has been debated. Several authors have suggested that the combination of eutrophication and sea-level rise may have a synergistic effect that results in enhanced losses of coastal habitats, and requires further research (Deegan et al., 2012;Kirwan and Megonigal, 2013;Kirwan et al., 2016;Crosby et al., 2017;Schuerch et al., 2018). While coastal eutrophication may enhance growth of foundation species, nutrient enrichment studies have reported a range of impacts on coastal systems depending on the local conditions (Anisfeld and Hill, 2012;Kirwan and Megonigal, 2013). For example, nutrient cycling and marsh stability were affected by local sediment characteristics, soil nutrients, microbial processes, and shifts in allocation of the plant species (McKee et al., 2007;Turner, 2011;Deegan et al., 2012;Lewis et al., 2021).
Predicting species level responses to environmental challenges requires an understanding of the amount of phenotypic variation within and among populations, which may reflect both phenotypic plasticity and heritable differences in phenotype (Richards et al., 2006;Nicotra et al., 2010;Banta and Richards, 2018). Several studies have shown that plant species harbor heritable differences in eco-physiological traits (Arntz and Delph, 2001;Geber and Griffen, 2003;Caruso et al., 2005), and variation in plasticity of traits (Sultan, 2001;Matesanz and Sultan, 2013;Nicotra et al., 2015;Matesanz et al., 2021). However, the amount of variation in natural populations for traits that are important for response to future climates is not well known (Davis and Shaw, 2001;Parmesan, 2006;Lovelock et al., 2016). Plants that inhabit coastal and intertidal zones have putative physiological adaptations that enable them to grow and reproduce in the anoxic and saline conditions that characterize these habitats. These adaptations include adjustment of water required by the plant (Antlfinger andDunn, 1979, 1983;Glenn and O'Leary, 1984;Donovan et al., 1996;Ball et al., 1997), adjustment of carbon uptake and nutrient absorption (Donovan et al., 1997;Lovelock et al., 2006Lovelock et al., , 2016Flowers and Colmer, 2008), and changes in resource allocation (Cavalieri and Huang, 1979;Glenn and O'Leary, 1984;Donovan et al., 1996Donovan et al., , 1997Flowers and Colmer, 2008;Lovelock et al., 2016). In addition, mangrove species can moderate anoxia that results from flooding via root growth, and altered peat formation has allowed mangrove communities to keep pace with sea level rise (McKee et al., 2007). Mangroves have been shown to respond to changes in nitrogen (N) availability by altering relative growth rate, photosynthetic rate, and resource allocation (Feller, 1995;McKee, 1995;Feller et al., 2003), which could be an important response to anthropogenic activities, such as runoff from agriculture and other types of land use change (Feller et al., 2003;Alongi, 2013).
Given spatial differences in salinity, anoxia, and N in the intertidal habitat, plasticity of traits that allow for tolerating such conditions may be adaptive. We expect intertidal plants like mangroves to show plasticity in response to salinity and N conditions (Antlfinger, 1981;Pennings and Richards, 1998;Richards et al., 2010b;Vovides et al., 2014;Lovelock et al., 2016). Proffitt and Travis (2010) found plasticity in growth rate and reproductive output within and among natural Rhizophora mangle mangrove populations in the Tampa Bay region of Florida in the United States. However, they also found both site of origin and maternal tree of origin affected R. mangle growth and survival, and that these effects varied by intertidal position (significant maternal family by elevation interaction; Proffitt and Travis, 2010). On the other hand, in a study of other Tampa Bay populations, we found low genetic diversity, and little population differentiation (Mounger et al., 2021). Although genetic diversity was low, we discovered high epigenetic diversity based on DNA methylation polymorphisms (Mounger et al., 2021). This type of epigenetic diversity has been associated with phenotypic and functional diversity, and could be a mechanism underlying phenotypic plasticity (Zhang et al., 2013;Nicotra et al., 2015;Herrera et al., 2017;Jueterbock et al., 2020). Several studies have suggested that epigenetic diversity could be particularly important in genetically depauperate species, providing a nongenetic source of response to the diverse conditions experienced by natural populations (Gao et al., 2010;Verhoeven et al., 2010;Richards et al., 2012;Jueterbock et al., 2020). Further, some epigenetic differences have been shown to be heritable. In fact, we discovered that the differences in DNA methylation in R. mangle propagules were predicted by maternal tree indicating a high degree of inheritance of differences in DNA methylation (Mounger et al., 2021).
In this study, we used plants from the same Tampa Bay populations to characterize within and among population level variation in putative adaptive traits in response to combinations of salinity and N in a full factorial design. Given the dynamic environment inhabited by R. mangle and the evidence of heritable non-genetic differences among populations, we hypothesized that propagules collected from different populations would respond differently to salinity and N amendment treatments. Our study was designed to test three predictions. First, R. mangle seedlings will be plastic in response to salinity and N amendment in putative adaptive traits that conserve water and adjust allocation of N. Second, response to salinity and N amendment will co-vary as plants shift resources to maintain osmotic balance. Finally, populations will vary in putative adaptive traits, and in plasticity of these traits, due to population differentiation.

Study Species
The red mangrove, Rhizophora mangle L. 1753 (Malpighiales, Rhizophoraceae), is an evergreen shrub or tree found along tropical and subtropical coastlines across the Americas, East Africa, Bermuda, and on a number of outlying islands across the South Pacific (Proffitt and Travis, 2014;Tomlinson, 2016;DeYoe et al., 2020) that can grow to heights of 24 m (Bowman, 1917). Poleward expansion of the species is limited by freezing events (its current northern range limit is roughly 29 • N latitude; Proffitt and Travis, 2014;Kennedy et al., 2017). Rhizophora mangle is considered a self-compatible species, with selfing rates in Tampa Bay estimated to be as high as 80-100% (Proffitt and Travis, 2005;Nadia and Machado, 2014). However, colder temperatures and contaminants from anthropogenic sources correspond with increased flowering and outcrossing, potentially resulting in higher genetic diversity particularly in the smaller populations at range limits from 28 to 30 • N Travis, 2005, 2014). Rhizophora mangle stands in our study area have a mean number of about 600 reproducing trees per kilometer of estuary (Proffitt and Travis, 2014). Pollinated R. mangle flowers mature in approximately 95 days, producing the buoyant hypocotyl also known as a propagule (Raju Aluri, 2013). The viviparous propagule germinates and matures on the maternal tree before it drops off, is dispersed pelagically, and becomes established as a seedling (McKee, 1995). Rhizophora mangle excludes salt in the root system through selective uptake of potassium (K + ) to sodium (Na + ) ions (Wise and Juncosa, 1989;Flowers and Colmer, 2008;Krauss et al., 2008;Medina et al., 2015), and allocates resources to manage osmotic potential (Bowman, 1917;Flowers and Colmer, 2008).

Sampling Design
We collected propagules from six populations of R. mangle  (Figure 1). The sites varied in salinity, mean tidal range, and neighboring species. We measured salinity at each site with a refractometer finding that salinity ranged from 20 to 40 parts per thousand (ppt) across the sites at the time of collection. This area has a humid subtropical climate with mean monthly temperatures ranging from 15.6 • C in January to 28.5 • C in August (1991-2020 monthly normals, U.S. NOAA National Centers for Environmental Information, station GHCND:USC00088824, Tarpon Springs, FL, USA), and annual precipitation of 1379 mm (annual mean 1991-2020, Tarpon Springs station). Precipitation falls as rain, with 60% falling during June through September, and 40% evenly distributed among other months. Monthly mean relative humidity ranges from 67% in April and May to 76% in August and September (1948-2018, U.S. NOAA Comparative Climate Data for the USA through 2018, station GHCND:USW00012842, Tampa International Airport). Tides are semi-diurnal, with 0.57 m median amplitudes (U.S. NOAA National Ocean Service, Clearwater Beach, Florida, station 8726724). Sea-level rise is 4.0 ± 0.6 mm per year (1973-2020 trend, mean ± 95% confidence interval, NOAA NOS Clearwater Beach station).
Honeymoon Island had a near monoculture of R. mangle while the remaining sites contained mixtures of two other mangrove species that are common in Florida: Laguncularia racemosa L. and Avicennia germinans L. We refer to plants from different sites as members of different populations based on our previous work which found differences among sites based on molecular markers (Mounger et al., 2021). At each population, we collected 20 propagules directly from each of 10 maternal trees separated by at least 10 m from each other to maximize the range of genetic variation sampled within each population (Albrecht et al., 2013). Propagules from each maternal tree were at least half-siblings but they could be more closely related due to the high selfing rate in the study area (Proffitt and Travis, 2005).
We refrigerated the propagules at 4 • C for up to 14 days until we planted them in the greenhouse at the University of South Florida Botanical Gardens. The greenhouse temperature was maintained between 18 and 29 • C. In the greenhouse, propagules from four of the maternal trees at AC and nine of the maternal trees at FD failed to establish, so we returned to sample propagules from eight new maternal trees at FD on August 12 and 29, and from the same original maternal trees at AC on October 17. Hence, while most of the propagules were in the greenhouse from the end of June until mid-October before they were exposed to treatments, some propagules from these FD and AC maternal families had less acclimation time before treatments began.

Experimental Treatments
We measured the length of each propagule and planted each in a 0.5-L pot with a 50:50 mixture of sand and peat soil. We watered the plants each day with tap water until we started applying the salinity and N amendment in mid-October. We set up the experiment in five spatial blocks. Within each block we randomized the position of plants such that each block had one replicate of each family for each treatment combination [i.e., a full factorial randomized complete block design with N = 6 populations × 10 maternal families × 4 environmental treatment combinations (2 salinity × 2 N fertilization) × 5 blocks × 1 replicate/block = 1200 plants]. The treatments were a factorial combination of low and high salinity with either no addition or addition of N. The salinity treatments were made with Morton solar salt (NaCl). The low salinity treatment was 15 ppt and the high salinity was 45 ppt, reflecting a slightly wider range of salinity than we measured in the field sites considering that we only measured salinity in the field at one time point and we expect the range to be slightly greater. The N treatments were made from equal amounts (in moles N) of urea and ammonium nitrate in tap water approximating 3 mg N per pot each week. This level is similar to a rate of 75 kg N per hectare per year, based on an estimated rate of loss by soil erosion and water runoff from corn crop residue in the USA (Pimentel et al., 1989). The source water for our irrigation solution meets drinking water standards, and thus has low mineral N and phosphorous (P) concentrations, with mineral N (nitrate-N plus ammonium-N) less than 1 mg N/L and P (phosphate) less than 0.1 mg P/L (Penuela Useche, 2015). Furthermore, the aquifer that provides our water had a median sulfate concentration of 8.2 mg/L per liter (interquartile range 2.3-17.8) (Berndt et al., 2015). This sulfur (S) would not be targeted for removal via water treatment, as the Environmental Protection Agency recommended S level is <250 mg S (as SO 4 )/L. This level of S would avoid plant S limitation since plants need about 1 mol S per 15-20 mol N.
At the start of treatments, we recorded seedling initial height from the soil to end of any growth. To avoid osmotic shock, the salinity treatment was applied twice a week and gradually increased by five ppt each treatment. The low salinity level (15 ppt) was reached in 2 weeks and the high salinity level (45 ppt) in 6 weeks. We started N treatments after the first week (when salinity treatments were 10 ppt), and applied N once per week from October 15, 2015 to May 1, 2016. We watered on nontreatment days with enough water to saturate the soil, but not flow through. Once per week, we watered with sufficient water to flow through the soil to prevent salt buildup. To determine if the N amendment was lost between treatments, we collected the flow-through leachate for a subset of eight plants, two of each combination of salinity and N treatment. We analyzed leachate for total dissolved N via combustion and luminescent detection (Skalar Formacs TN analyzer, Breda, Netherlands).

Traits Measured
We characterized each plant as alive, dormant, or dead at the end of the experimental treatments. We assigned plants that showed no growth and no desiccation to the dormant category. We measured five traits related to salt tolerance and overall performance for alive plants: change in height from beginning to end of treatments (cm) (hereafter, height growth), leaf mass per area or LMA (dry leaf mass in g/total leaf area cm 2 ), succulence (grams of water in all leaves/total leaf area cm 2 ), root to shoot ratio based on dry biomass, and total dry biomass (g) at harvest. We only used healthy leaves for succulence and LMA. We defined healthy leaves as attached, a minimum of 50% green, and fully developed. For dry above ground biomass, we included leaves that were attached, but not 50% green or fully developed. We measured the biomass of above and below ground tissues of all harvested plants after the tissues were dried at 60 • C until they maintained constant mass. Finally, we measured the total dry mass of leaves after drying in silica desiccant beads for a minimum of 7 days to constant mass.
In addition, we used a LI-COR 6400 to measure maximum photosynthetic rate (micromoles CO 2 /m 2 s) for a subset of the plants just prior to harvest. We determined that the appropriate photosynthetically active radiation (PAR) for saturation in these plants was 1000 micromoles PAR/m 2 s based on light curves generated from six data points from each of two plants (one low salinity -no N and one high salinity -high N). We then measured maximum photosynthetic rate on one plant with at least two healthy leaves for each surviving maternal line for each treatment (n = 29 low salinity -no N; n = 31 high salinity -no N; n = 26 low salinity -high N; and n = 32 high salinity -high N, for total N = 118 total plants). We took maximum photosynthetic rate measurements at a CO 2 rate of 400 micromoles/(m 2 s) and a flow rate of 500 micromoles/s. We took the measurements on a healthy leaf from the second node of each plant after the leaf had been clamped in the LI-COR for 1 min to ensure conditions had stabilized. We measured maximum photosynthetic rate for the 118 plants in random order by block over six consecutive days from April 23 to 28, 2016, between 8:30 and 11:30 in the morning.

Statistical Analysis
We performed all statistical analyses in R, version 4.0.3 (R Core Team, 2020). All analyses reported here used either the General Linear Model (GLM) or Generalized Linear-Mixed-Model (GLMM) frameworks. GLMs are a generalization of the familiar analysis of variance framework, which allow one to use error distributions other than the restrictive Gaussian used in ANOVA. Moreover, the GLM link function allows one to fit a linear model to a dependent variable that has an essentially nonlinear relationship with the predictors. For example, using a log link, one can model the log of the expected value of the data as a linear function; this is distinct from ANOVA on log-transformed data, which models the expected value of the log. GLMMs are a further generalization that allow one to model the variances, rather than expected values, of some predictors. These predictors are sometimes (incorrectly) assumed to be sampled randomly; conceptually the real distinction between these predictors and the others is whether one is interested in the mean values or variances. We fit all models as GLMMs because we were interested in the variance among the maternal families, rather than in their particular means (Bolker, 2015). However, we also fit GLMs that excluded the family terms, to determine whether maternal family contributed to improved model fits. We also treated Block as a "random" term where possible, but when such models failed to converge numerically (due to the small number of blocks) we fit the analogous model with Block as a fixed term (Bolker, 2015).
We checked the residuals to assess normality on traits as appropriate; we did not transform height growth or photosynthesis (lmer and glm), but we used the log link function (glmer) for analysis of succulence, LMA, root to shoot ratio, and total biomass. We used the function lmer or glmer implemented within the lme4 package (Bates et al., 2015) to fit a series of models and identify the best fit model for survival and for each trait ( Table 1).
For survival and each trait, we began with a saturated model that included as fixed factors salinity, N, and their interaction as well as the random effects of block, population and maternal family nested within population. We ran subsequent models removing individual terms from the model (see details in the Supplementary Material file "Rhizophora mangle data analysis code and annotated results"). In several cases, models with random terms for block or population failed to converge ( Table 1), most likely because there were relatively few blocks and populations. In these cases, we proceeded by treating these as fixed effects. We used AIC to evaluate the fit among all of the models that we ran for survival and for each trait. The model with the lowest AIC was considered to have the best fit. We used ANOVA (type III) in the package car (Fox and Weisberg, 2019) to evaluate the significance of fixed effects when they were included in the best model. When the random terms were found to be significant, we provided graphics of the relative differences among the maternal families, populations or blocks in the Supplementary Material file "Rhizophora mangle data analysis code and annotated results" to visualize the relative differences among these groups.
For survival, we assessed three states that were coded as 0 for live plants, 1 for dormant plants, and 2 for plants that died during the experiment. We modeled survival using random effects logistic regression. In one set of models, we included dormant plants as alive, and in another we excluded them. The results were qualitatively similar, so we report only the case where dormant plants were treated as alive (see results for both in Supplementary Material file "Rhizophora mangle data analysis code and annotated results"). For change in height, we included the covariate of height at the beginning of treatments in October.
To gain insight about variance explained by models, we calculated the R 2 approximations proposed by Nakagawa and Schielzeth (2013). As these authors explained, in the context of GLMMs this leads to two different sorts of R 2 , a marginal R 2 that reflects variance explained by fixed factors only, and a conditional R 2 that reflects variance explained by both fixed and random factors. We reported each of these as appropriate, e.g., where the best model included only fixed factors we reported only the marginal R 2 .

RESULTS
We found that some combination of salinity and N treatments were included in the best models for three of the six traits we measured: succulence, root to shoot biomass ratio, and maximum photosynthetic rate (Table 1). In addition, maternal families and populations were significantly different for survival and all traits We started with a saturated model that included salinity, N, and their interaction as fixed factors as well as the random effects of block, population and maternal family nested within population. We ran subsequent models after removing individual terms from the model (see details in Supplementary Material file "Rhizophora mangle data analysis code and annotated results"). Fam, maternal family; ht, height; LMA, leaf mass per area (dry mass/area); Pop, population; Pop/Fam, maternal families nested within source populations; RTS, root to shoot biomass ratio; SNG, number of survivors in April. Terms in parentheses are random terms; (1| term) indicates that a random intercept is estimated for each term. For some models, Block or Pop were fit as fixed effects if estimation as random effects failed (generally due to the small numbers of blocks and source populations). AIC is the difference between the saturated model (or closest model to it when saturated model was singular) and the AIC for the given model. measured with the exception of photosynthesis (Table 1). Below we report the results of model selection and assess significance of effects retained in the best model.

Treatment Validation
To ensure that our N amendment treatments were not flushed out during the once weekly flow through watering we measured the total dissolved N in leachate from a subsample of the seedlings. We found that N was not significantly different between the low salinity -no N and the high salinity -high N amended plants, and therefore, confirmed that we did not lose the N due to watering between treatments [Mean Square = 0.11, F ndf 3/ddf 48 = 0.23, Pr(>F) = 0.88].

Survival
Of 1,149 plants, 76 were unequivocally dead in April, 907 unequivocally alive, and 166 dormant. The number of plants that showed active growth ranged from 63% in HI to 92% in FD (Figure 2). The number of plants that did not grow, but also did not appear to be dead ranged from 3% in WB and WI to 10% in HI. We modeled survival by including these dormant plants alternatively as either alive or dead. In both cases, the best-supported model was one including a fixed effect for population and a random effect for maternal family ( Table 1

Trait Responses to Treatment
The only useful predictor among the fixed effects for the change in height was the height at the start of treatments in October (Table 1; marginal R 2 = 0.82; ANOVA chisq = 4432.67; p < 0.0001). The conditional R 2 was 0.92, and the marginal R 2 was 0.82, indicating that 82% of the variance in change in height was explained by a positive relationship with the fixed effect of height at the beginning of the experimental treatments (estimate = 0.9009). The variance components (scaled to sum to 1) for population and maternal family within population were 23 and 37% of the random variance, respectively. A plot of the random effects suggested that the populations were largely similar, with the exception of FD which was very different (Supplementary Figure 2). It is not obvious that any population had more among-family variability, but our design was limited to determine this.
For LMA (log link), the best model included only the random effects of block, population and maternal family within population. The conditional R 2 was 0.99. In this model, maternal family explained 51% of the variance, population explained 30% of the variance and block explained 17% (Supplementary Figure 3).
For succulence (log link) the best model included the fixed effect of salinity and the random effects of population and maternal family nested within population. The marginal R 2 = 0.39, while the conditional R 2 = 0.99. The scaled variance components for maternal family within population and population were, respectively, 0.91 and 0; in other words there was very large variance among maternal families within populations, but not among populations (Supplementary  Figure 3). While increased salinity resulted in a statistically meaningful reduction in succulence, it was only when conditioned on family that it added to the explanatory power of the model.
Root to shoot biomass ratio (log link) was the only variable for which the data supported the saturated model as the best model (Figure 3). The conditional R 2 = 0.206, while the marginal R 2 = 0.017. An ANOVA to evaluate the fixed effects revealed that the main effect of salinity was significant (chisq = 7.38; p = 0.007) indicating increased root to shoot allocation in response to salinity, but not the main effect of N. In addition, the interaction of salinity × N was significant (chisq = 6.19; p = 0.01). At ambient N levels (no N added), root to shoot ratio was not affected by increasing salt concentration. With the addition of N fertilizer, root to shoot ratio increased in response to high salinity ( Figure 3B). However, the small size of the marginal R 2 suggested that these effects were mainly meaningful when conditioned on the random terms. The random terms population, maternal family nested within population, and block account for 13, 24, and 6% of the random variance, respectively (Supplementary Figure 5).
The best model for total biomass (log link) included only the random effects of population and maternal family within population. The conditional R 2 = 0.12. The terms for population and family nested within population accounted for 13 and 21% of the random variance, respectively (Supplementary Figure 6).
The photosynthesis data were the most limited in sample size since we were only able to assess one individual with at least two healthy leaves for each surviving maternal line for each treatment (N = 118). The best model was one including the fixed effects salinity, N, and their interaction, but no random effects ( Table 1). An ANOVA to evaluate these fixed effects showed the main effects of salinity and N were not significant but the interaction was (chisq = 4.46; p = 0.035). At ambient N levels (no N added), maximum photosynthetic rate declined with increasing salt concentration, but this negative impact of salinity was absent or reversed upon the addition of N fertilizer ( Figure 3C).
In summary, we found plasticity to our treatments for three of the six traits we measured. Succulence decreased in response to higher salinity ( Figure 3A), but this response varied largely by family (Supplementary Figure 3). For root to shoot biomass ratio and maximum photosynthetic rate the responses to experimental treatments depended on changes in both salinity and N (Figures 3B,C). Root to shoot ratio also varied overall in response to salinity as well as by family and population (Supplementary Figure 4).

DISCUSSION
We assessed the growth and survival of R. mangle propagules to full factorial combinations of two salt and two N levels because these are two important abiotic properties that potentially have important impacts on mangrove survival and traits, and by extension, on the biodiversity and ecosystem function of coastal wetlands. In addition to natural variation in these conditions, anthropogenic activities may result in more extreme levels of these conditions from runoff and flooding (Ellison et al., 2005;Krauss et al., 2006;Gedan et al., 2009;Kirwan and Megonigal, 2013;Lewis et al., 2014). We predicted that R. mangle seedlings would be plastic in response to salinity and N amendment in putative adaptive traits that conserve water and adjust allocation of N. Our study showed that only succulence and root to shoot ratio were plastic in response to salt, regardless of N treatment. Allocation to root and shoot biomass and maximum photosynthetic rate were also plastic, but response of both traits to N amendment depended on the level of salt. This supported our second prediction that response to salinity and N amendment will co-vary as plants shift resources to maintain osmotic balance. We also found support for our third prediction that populations would vary in putative adaptive traits, and in plasticity of these traits, due to population differentiation. Importantly, every trait except for photosynthesis varied among population and maternal families within populations. This was also true of survival. In fact, maternal family and population were the most consistent predictors for variation in traits and survival.

Phenotypic Plasticity in Response to Treatments
We expected that R. mangle propagule survival and growth would decrease in response to high salinity and increase in response to N fertilization. We also expected that N could alleviate some of the effects of salinity as indicated by an interaction of the two conditions. However, we found no response to treatments in survival, height growth, LMA or total dry biomass. This may be because the propagules were supported by resources provided by the maternal tree, which in R. mangle can support growth for at least a year (Ball, 2002;Proffitt and Travis, 2010). If the seedlings were supported by these maternal reserves, height growth would likely be more correlated to propagule length at collection which would be corrected for in the start-of-experiment (time zero) height measurements that we included as a covariate. Because our treatment duration was only 6 months, the lack of growth response to treatments was consistent with dependence on maternal reserves. However, we did see seedling response to treatment in succulence, allocation to root and shoot biomass and maximum photosynthetic rate.
Increased succulence is a common response to water deficiency under high salinity conditions (Rosenthal et al., 2002;Vendramini et al., 2002;Ottow et al., 2005;Karrenberg et al., 2006;Richards et al., 2008Richards et al., , 2010a, but in our experiment we detected reduced succulence in response to high salt. However, this could be because R. mangle excludes salt. Many halophytes absorb high concentrations of salt for osmotic adjustment that would normally be toxic. One strategy is to compartmentalize these ions along with increased succulence (Flowers and Colmer, 2008). On the other hand, many plants avoid this toxicity by excluding salt at the roots (Cavalieri and Huang, 1979;Donovan et al., 1996;Tester and Davenport, 2003). For example, several other halophytes that are salt excluders, including the succulent plant Salicornia europaea L., and another member of Rhizophoraceae, Kandelia candel (L.) Druce, do not increase succulence or leaf thickness in response to high salinity (Glenn and O'Leary, 1984;Kao et al., 2001). In fact, with N fertilization K. candel decreased leaf thickness when salinity was increased (Kao et al., 2001). Thus, one possible explanation for our results is that the N fertilized seedlings were able to reallocate resources (e.g., increased N) to salt tolerance mechanisms like compatible solutes, and still maintain positive water balance in the high salt condition without increased succulence.
Although we saw a statistically significant response to our treatments in three of the six traits, the interaction between salt and N fertilization only effected R:S and maximum photosynthetic rate. We expected maximum photosynthetic rate to increase in response to N fertilization because the enzyme RuBisCO, which catalyzes the dark reactions in photosynthesis, requires a large amount of N (Sage et al., 1987;Andersson, 2008). Further, a meta-analysis across 356 diverse species representing most biomes and growth forms showed increased maximum photosynthetic rate with increased N is a general phenomenon (Walker et al., 2014). This could lead to increased biomass allocation to shoots relative to roots. Despite this expectation, there was no overall response to high N level independent of salt treatment. One reason might be that photosynthesis overall was limited by other nutrients, not just N, and thus increasing N alone might not have been enough to elicit a response (see e.g., Lovelock et al., 2006). In a field study, dwarf R. mangle did not respond to N alone, presumably because other nutrients were limiting (Feller, 1995;Lovelock et al., 2006). Dwarf mangroves did increase biomass in response to fertilization with all three nutrients: N, phosphorus (P), and potassium ions (K + ) (Feller, 1995), and in response to P addition alone (Lovelock et al., 2006). We also expected that response to N amendment would depend on salinity. We found that in plants treated with high salt, maximum photosynthetic rate was slightly enhanced by high N. Possibly, the additional N enabled the plants to synthesize N-rich compatible solutes for osmotic regulation and continued photosynthetic gain of carbon (Flowers and Colmer, 2008). Plants in high salt also responded differently in allocation of biomass with increased N; instead of increasing shoot biomass, they increased root biomass on average.

Variation Within and Among Sites
Phenotypic variation that is maintained in common garden from within and among populations would indicate R. mangle has heritable trait diversity which could allow for response to changing environmental conditions. Seedling survival depended on population and varied among maternal families within the six populations. We found variation in height growth, succulence, LMA, root to shoot allocation, and total dry biomass was largely determined by maternal families within populations. Proffitt and Travis (2010) also found seedling survival varied among maternal families, as well as by location in the intertidal zone. But in their study after 3 years, growth and survival did not reflect initial propagule size (Proffitt and Travis, 2010). Our results support these previous findings that propagule length is positively correlated to short term performance. This findings suggests that maternal reserves in the R. mangle propagule can help the seedling survive, and larger propagules contain more maternal reserves than smaller propagules (Ball, 2002;Proffitt and Travis, 2010). Because our study was a short-term, controlled greenhouse study, maximum photosynthetic rate might be the best indicator as an immediate response. Variation in maximum photosynthetic rate can ultimately manifest as variation in growth and allocation of resources, particularly once the seedling has depleted maternal reserves. The seedlings did not show significant differences in height growth or total dry biomass in response to treatments. But given the plasticity we saw in maximum photosynthetic rate, and the 1-year and 3-year growth results found in a previous study of nearby populations (Proffitt and Travis, 2010), we might have found a more dramatic response to these treatments with more time (Ball and Farquhar, 1984). However, previous work in several different systems have argued that greenhouse studies are often unable to recreate relevant field conditions, so such responses may only be registered in the field (Schittko et al., 2016;Rinella and Reinhart, 2017;Forero et al., 2019;Dostál et al., 2020).
The amount of heritable phenotypic diversity and differentiation we discovered in this study is an important indicator of the potential for this species to respond to changing conditions. The high level of diversity in response among populations and maternal families within populations may be surprising since we previously reported low genetic diversity among these plants based on molecular markers. Low genetic diversity is expected to limit the potential for different responses among individuals. On the other hand, we discovered high epigenetic diversity (Mounger et al., 2021), which could contribute to phenotypic and functional diversity, and could be a mechanism underlying the type of phenotypic differences and plasticity we found here (Zhang et al., 2013;Nicotra et al., 2015;Herrera et al., 2017;Jueterbock et al., 2020). In addition, we know very little about the interactions with the microbiome in the species, but microbes have been highlighted as important symbionts in these and other challenging environments (Bowen et al., 2017;Angermeyer et al., 2018;Jung et al., 2021). Soil microbial activity could have dramatic impacts on the future nutrient availability and stability of these coastal sediments (Deegan et al., 2012;Bowen et al., 2017;Hughes et al., 2020;Lewis et al., 2021). In fact, a recent study suggested that bacterial community composition differed among R. mangle maternal genotypes but not with genetic diversity (Craig et al., 2020).

CONCLUSION
Previous work suggested that R. mangle growth and survival depended on an interaction of intertidal elevation and maternal genotype, suggesting variation in response to flooding conditions (Proffitt and Travis, 2010). However, in addition to changes in flooding, anthropogenic activities are causing changes in salinity and N level in R. mangle ecosystems (McKee et al., 2007). The interacting effects of salinity, N level, and elevation are complex and potentially non-additive (McKee et al., 2007). Our experimental findings suggest that the important traits of succulence, root to shoot biomass allocation and photosynthetic rate respond to salinity, N level, or the combination of these conditions. We also discovered that seedling survival and the magnitude of almost all responses varied among populations and even maternal families within populations.
This variation in survival and important traits among families and among populations is particularly interesting considering the importance of this foundation species for the functioning of the coastal ecosystem. Our previous work showed a lack of genetic diversity, which might be alarming considering the expected limitations of low genetic variation (Allendorf et al., 2010;Estoup et al., 2016). However, accumulating studies provide important evidence that genetic variation must be interpreted with caution (Hufford and Mazer, 2003;Estoup et al., 2016) and that the emphasis on only variation in DNA sequence can be misguided (Keller, 2002(Keller, , 2014Sultan, 2015;Bonduriansky and Day, 2018). Non-genetic sources of response may contribute to the phenotypic diversity we report here that is particularly relevant under different salinity and N conditions. This may provide another source of resilience for R. mangle and other critical species to changing environmental conditions and contribute to future adaptation to a complex mosaic of environmental conditions.

DATA AVAILABILITY STATEMENT
All of the data collected and analyzed for this study along with R code are included in the Supplementary Material, Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CR, KL, GF, and DL conceived the study and designed the experiments. KL and JM collected plants and maintained the experiments. KL and GF analyzed the data. CR and KL wrote the first draft of the manuscript. All authors provided input and revisions to the manuscript.

FUNDING
This work was supported by funding from the National Science Foundation (United States ) IOS-1556820 (to CR), DUE-1930451 (to DL), and the Federal Ministry of Education and Research (BMBF, Germany; MOPGA Project ID 306055 to CR).