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 (Proffitt and 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).

We collected propagules from six populations of R. mangle between June 9 and June 26 of 2015, on the west coast of central Florida (USA) within the following county and state parks: Anclote Key Preserve State Park (AC), Fort De Soto Park (FD), Honeymoon Island State Park (HI), Upper Tampa Bay Conservation Park (UTB), Weedon Island Preserve (WI), and Werner-Boyce Salt Springs State Park (WB) (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.

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₄)/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 non-treatment 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).

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²), succulence (grams of water in all leaves/total leaf area cm²), 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₂/m² 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² 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₂ rate of 400 micromoles/(m² 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.

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 non-linear 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² approximations proposed by Nakagawa and Schielzeth (2013). As these authors explained, in the context of GLMMs this leads to two different sorts of R², a marginal R² that reflects variance explained by fixed factors only, and a conditional R² 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².

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].

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; Supplementary Figure 1). The model treating dormant plants as alive had conditional R² = 0.11 and marginal R² = 0.08. The model with dormant plants as dead had conditional R² = 0.28 and marginal R² = 0.13. Because there is no residual in the equation defining logistic regression, no variance component for a random residual is calculated, and thus it is not possible to calculate a meaningful scaled variance component for family here.

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² = 0.82; ANOVA chisq = 4432.67; p < 0.0001). The conditional R² was 0.92, and the marginal R² 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² 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² = 0.39, while the conditional R² = 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² = 0.206, while the marginal R² = 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² 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² = 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).

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., 2008, 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.

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).