Hurricanes fertilize mangrove forests in the Gulf of Mexico (Florida Everglades, USA)

Significance Despite the destructive effect of hurricanes on mangrove forests in tropical and subtropical latitudes, hurricanes are major drivers controlling soil fertility gradients in the Florida Everglades mangroves, and therefore represent a positive influence in maintaining observed mangrove spatial distribution and productivity patterns. Hurricane-induced mineral inputs to near-coast mangroves in the Everglades enhance phosphorus (P) concentrations in soils, increase plant P uptake, promote soil elevation gains relative to sea level, and facilitate rapid forest recovery following disturbance. The response of mangroves to large-scale P fertilization from hurricanes may be an important adaptation of neotropical mangroves in the Gulf of Mexico and the Caribbean region to withstand the impacts of both sea-level rise and P limitation. Hurricanes are recurring high-energy disturbances in coastal regions that change community structure and function of mangrove wetlands. However, most of the studies assessing hurricane impacts on mangroves have focused on negative effects without considering the positive influence of hurricane-induced sediment deposition and associated nutrient fertilization on mangrove productivity and resilience. Here, we quantified how Hurricane Irma influenced soil nutrient pools, vertical accretion, and plant phosphorus (P) uptake after its passage across the Florida Coastal Everglades in September 2017. Vertical accretion from Irma’s deposits was 6.7 to 14.4 times greater than the long-term (100 y) annual accretion rate (0.27 ± 0.04 cm y−1). Storm deposits extended up to 10-km inland from the Gulf of Mexico. Total P (TP) inputs were highest at the mouth of estuaries, with P concentration double that of underlying surface (top 10 cm) soils (0.19 ± 0.02 mg cm−3). This P deposition contributed 49 to 98% to the soil nutrient pool. As a result, all mangrove species showed a significant increase in litter foliar TP and soil porewater inorganic P concentrations in early 2018, 3 mo after Irma’s impact, thus underscoring the interspecies differences in nutrient uptake. Mean TP loading rates were five times greater in southwestern (94 ± 13 kg ha−1 d−1) mangrove-dominated estuaries compared to the southeastern region, highlighting the positive role of hurricanes as a natural fertilization mechanism influencing forest productivity. P-rich, mineral sediments deposited by hurricanes create legacies that facilitate rapid forest recovery, stimulation of peat soil development, and resilience to sea-level rise.

D isturbances are large-scale episodic events that create abrupt changes in community structure, regulate ecological processes in ecosystems (1), and generate biological legacies that interact with environmental conditions, thus defining trajectories of ecosystem recovery (2,3). Understanding spatiotemporal ecosystem responses to disturbances and how resilience maintains specific ecosystem states and regime shifts is challenging and still one of the knowledge gaps in ecological theory (4,5). Since resilience is the "capacity of a system to experience shock while retaining essentially the same function, structure, feedbacks, and therefore identity" (5)(6)(7), it is paramount to determine the effect of natural disturbances on ecosystem response, particularly in coastal regions. Given their location at the boundary between terrestrial and marine ecosystems along tropical and subtropical latitudes, mangrove wetlands are particularly prone to potential regime shifts as a result of increasing human impacts (8) and landscape level interactions with natural, high-energy "pulsing" disturbances, such as hurricanes (9)(10)(11)(12). However, how those shifts occur is unknown in part because of the lack of data and information to discern the magnitude and direction of ecosystem responses (13). Mangrove species are well adapted to recover quickly from disturbance due to their resilient traits (e.g., resprouting from epicormic shoots, high rates of water-use, and nutrient-use efficiency) (10,14,15). Hurricane-force winds change forest structure through defoliation, tree snapping, and uprooting, which in turn influence species composition, successional patterns, nutrient cycling, tree mortality, and potential loss in soil elevation (9,(15)(16)(17)(18)(19). The degree of impact depends on physical conditions of the storm (i.e., intensity, wind velocity, size) and the proximity of the forest to the hurricane's path (20). In addition to the immediate physical damage of hurricanes, storms can also deliver phosphorus (P)-rich mineral sediments that build elevation and provide nutrients that fertilize soils, promoting vegetation regrowth and fast recovery postdisturbance (13,(21)(22)(23)(24). This stimulating effect of storms in northern neotropical mangrove forests (i.e., the Florida Coastal Everglades, FCE) was first hypothesized by statistically relating tree height to soil P concentrations along a fertility gradient in the Shark River estuary, southwestern FCE, where the tallest trees grow in soils with the highest soil total P (TP) concentrations in the lower estuary (25,26). Concentrations of Ca-bound P in this Significance Despite the destructive effect of hurricanes on mangrove forests in tropical and subtropical latitudes, hurricanes are major drivers controlling soil fertility gradients in the Florida Everglades mangroves, and therefore represent a positive influence in maintaining observed mangrove spatial distribution and productivity patterns. Hurricane-induced mineral inputs to nearcoast mangroves in the Everglades enhance phosphorus (P) concentrations in soils, increase plant P uptake, promote soil elevation gains relative to sea level, and facilitate rapid forest recovery following disturbance. The response of mangroves to large-scale P fertilization from hurricanes may be an important adaptation of neotropical mangroves in the Gulf of Mexico and the Caribbean region to withstand the impacts of both sea-level rise and P limitation. estuary region were 40 times higher than in the upper estuary, suggesting a shift from organic to mineral P down the estuary gradient. Thus, Ca-bound P inputs during storm events to the mouth of the Shark River from the Gulf of Mexico (GoM), rather than upland inputs of nutrients, are thought to influence mangrove development and productivity patterns in the southwestern Everglades (25,26).
Mangroves in karstic settings are generally deprived of large volumes of terrestrial sediment supply and nutrients (e.g., P) compared to minerogenic environments (e.g., river-dominated mangroves) that receive an allochthonous sediment supply (27). Therefore, plant-soil feedbacks are driven by the calcium carbonate soil matrix, where P supply is often limited due to carbonate chemisorption reactions and the sequestering of P into recalcitrant organic pools, which reduces the availability of exchangeable P to plants (28). The southwestern FCE in the GoM is a karstic environmental setting where P concentration in wetland soils and benthic sediments is highly limiting, as observed in other locations throughout south Florida, the Caribbean, and other GoM regions (29,30). Six years after this correlative observation between forest development and soil fertility was documented, Hurricane Wilma (October 2005), a category 3 storm, landed on the southwest coast of Florida, creating large-scale physical damage to mangrove forests along the Shark River and other estuaries (22,31). Wilma's storm surge also deposited large amounts of mineral sediment equivalent to twice the average surface (top 10 cm) soil P density (0.19 mg cm −3 ) before impact. Furthermore, vertical accretion from this storm was 8 to 17 times greater than the annual vertical accretion rate (0.30 ± 0.03 cm y −1 ) averaged over the last 50 y (22). Mangrove forests along the estuary were heavily impacted close to the coastline, as evidenced by >90% defoliation and cumulative mortality over 10 y (24). Yet, mangrove forest canopy and net primary productivity returned to pre-Wilma conditions 5 y after impact (2010), indicating a significant forest resilience capacity (i.e., potential ability of the mangrove region to recover or bounce back to prestorm conditions) (13).
On September 10, 2017, 7 y after mangrove recovery from Wilma, mangroves along the southwest Florida coast, including the Shark River, were impacted by Hurricane Irma (hereafter Irma), a category 3 storm. Irma made landfall on the southwest Florida coast near Marco Island with sustained, hurricane-force winds of 180 to 193 km h −1 (112 to 120 mph). The storm surge produced maximum water levels up to 3 m above ground level along the southwest coast of Florida (32). Irma's winds resulted in large-scale physical damage to mangrove forest structure (i.e., defoliation, tree snapping, and uprooting) while the storm surge deposited a thick layer of allochthonous mineral sediment across the FCE. Indeed, this landscape mineral deposition was similar to that caused by Wilma in 2005 after its passage across the FCE (22). This recurrent storm pattern allowed a systematic landscape-level evaluation of the interaction of storm surge impacts with wind damage since our long-term (19 y) mangrove experimental plots and transects were already positioned across the impacted area.
Hurricanes are quite common in south Florida, where the tropical storm recurrence frequency is 6 to 8 y (33). Because the frequency of strongest hurricanes (categories 4/5) is expected to increase as a result of climate change (34)(35)(36)(37), Irma's impacts provided a unique opportunity to further quantify the positive and negative effects of hurricanes on mangrove forest productivity. In fact, we provide strong evidence of the positive role of hurricane-induced P fertilization (i.e., Ca-bound P) in maintaining landscape gradients of soil fertility that control mangrove vegetation patterns in the Everglades. We hypothesized that the distribution and thickness of new storm deposits on mangrove soils would be similar to that observed after Wilma, with a significant reduction in deposition with distance inland from the mouth of estuaries. Since we assumed that the mineral material deposited on mangrove soils after the storm surge was resuspended sediment from the coastal shelf, we also expected similar total P and N density as those measured after Wilma in 2005. Hydrographs indicated that Irma's storm surge amplitude was significantly ameliorated by mangrove vegetation along the estuaries. Therefore, we expected maximal deposition in mangrove areas adjacent to the mouth of estuaries, and lack of storm deposits 14 to 18 km from the GoM. We also anticipated a decreasing P gradient with higher deposition at the fringe mangrove zone next to the channels/estuary relative to the interior forest along the estuary. We also used long-term species-specific leaf litter TP content (2008,2014,2018) after hurricane impact to determine potential plant P uptake and its relationship with soil TP density (2007 to 2017). We expected higher leaf enrichment in sites with high soil TP concentration, especially after Irma's impact. Moreover, given the magnitude of the estimated P loading rates, we predicted that these rates would be at least similar to rates measured in treatment wetlands receiving agriculture runoff high in TP upstream from the Everglades. Our results provide a quantitative assessment of TP loading rates in mangrove forests resulting from hurricanes. Quantifying the relative magnitude of this allochthonous mineral input will help advance our understanding of P cycling in oligotrophic wetlands and determine the relative balance and interaction of long-term effects, both positive (P fertilization, soil accretion) and negative (tree mortality, defoliation) (13), on mangrove forest structure and productivity. Because hurricane disturbance in coastal regions is frequent and increasing, a consideration of the interaction of these effects is critical to evaluate the role of pulsing high-energy disturbances (hurricanes) and pressing changes (sea-level rise) in maintaining forest species composition, productivity, and soil elevation in neotropical mangrove wetlands under climate change.  above the forest soil surface and ∼0.75 m at mid-and upstream locations (SRS-5 and -4). Tidal gauges in adjacent Harney and Broad Rivers, close to the Shark River, showed maximum surge of 1.1 m and water levels from 0.25 to 0.90 m above average. Water levels in SRS-6 and SRS-5 peaked at about 5:00 PM and 9:00 PM on September 10th, respectively, whereas a more gradual increase was observed upstream at SRS-4, peaking at 9:00 AM the following day. At the Taylor River sites, water levels peaked early morning on September 10th at the downstream site (Taylor River Slough [TS/Ph]-7) and during early afternoon at the upstream location (TS/Ph-6). Overall, the storm surge lasted ∼6 to 9 h (Fig. 1).

Results
Storm Sediment Deposition and Nutrient Inputs. Irma's storm surge delivered a large-scale mineral sediment deposition that was evident up to 10-km upstream from the GoM in southwestern FCE, with higher deposition in near-coast mangroves relative to upstream mangrove areas. There was also a lateral gradient in deposition at all transect sites, with storm sediments declining in thickness from the fringe zone to the forest interior ( Fig. 2 and SI Appendix, Table S1). Storm sediments at SRS-5 (midstream estuary) extended up to 100 m into the mangrove forest, whereas at SRS-6 storm deposition was more widespread (3.8 ± 0.3-cm deposits within the first 50 m, decreasing to 2 cm at 350 m into the forest). In the Harney River, deposition was higher downstream (WSC-10; Water, Sustainability, and Climate Project; http://sfwsc.fiu.edu/) compared to midstream (WSC-9), but consistently decreased with distance into the forest at both sites (Fig.  2). Along the Broad River transect (WSC-13, downstream), sediment deposition peaked at 50 m in the forest interior (4.5 cm) and decreased to 1.6 cm at 350 m. In Taylor Ridge, sediment deposition was highly variable, with higher deposition in the creek side (2.6 cm) and middle section of the transect (2.5 cm at 80 m) compared to the Bay side (0.7 cm) (Fig. 2). Sediment deposition in discrete mangrove locations (i.e., sites sampled only at 30 to 40 m from edge) was also different among basins. Along the Broad River, deposition was lower upstream (WSC-11: 1.9 cm) relative to the midstream section of the estuary (WSC-12: 3.8 cm). A similar trend was observed at the Harney River, with no deposition at the most upstream site (WSC-7). At the Shark River mouth (SRS-7), deposition was 3.7 ± 0.5 cm.
Bulk density (BD) of hurricane deposits (0.6 ± 0.1 g cm −3 ) was significantly higher in all transect sites compared to underlying surface (top 10 cm) mangrove soils (0.3 ± 0.1 g cm −3 ), except in Taylor Ridge, where both storm sediments and surface soils showed high BD values (SI Appendix, Tables S1 and S2). BD data from discrete sites showed the same trend between storm sediments and mangrove soils compared to transect sites (SI Appendix, Table S3). Organic matter (OM) content showed an opposite trend relative to BD, with lower values in storm sediments (17.4 ± 1.5%) compared to surface soils (42.5 ± 6.5%) across all transect and discrete sites (SI Appendix, Tables S1-S3).
Total organic carbon (TOC) concentrations were lower in storm sediments (31.9 ± 2.3 mg cm −3 ) than in mangrove soils (40.0 ± 1.2 mg cm −3 ) across transect and discrete sites (SI Appendix, Tables S1-S3). In contrast, total inorganic carbon (TIC) was double in storm deposits (40.6 ± 6.8 mg cm −3 ) compared to surface soils (22.1 ± 12.1 mg cm −3 ) in all transect sites, except at Taylor Ridge (SI Appendix, Table S2). Discrete sites showed a similar pattern (SI Appendix, Table S3). Overall, TIC was the largest fraction in storm sediments compared to soils (except at Taylor Ridge) and accounted for 39 to 65% of the total soil C (TC) pool across all mangrove sites.
Mean total N (TN) concentrations were significantly different among transect sites, with higher values at Taylor Ridge, SRS-5, and WSC-9 relative to the rest of the sites ( Fig. 3A and SI Appendix, Table S1). However, there was no significant difference in either TN between layers or the interaction between sites and layers. Discrete sites showed similar magnitudes in TN between both layers (Fig. 3A). TP density varied significantly among sites and layers ( Fig. 3B and SI Appendix, Table S1). In contrast to TN, TP was on average 1.5 times higher in storm sediments (0.29 ± 0.04 mg cm −3 ) than in surface soils (0.19 ± 0.02 mg cm −3 ). All transect and discrete sites had higher TP in storm sediments compared to soil surface (Fig. 3B), except at Taylor Ridge, which showed an opposite trend. Overall, TP density in surface soils consistently decreased with distance inland from the mouth of all estuaries (Fig. 3B). When comparing estuaries, there were significant (F 3, 31 = 33.3, P < 0.001) differences in TP density between Taylor Ridge and the three estuaries in southwestern FCE, with the lowest mean value in Taylor Ridge (0.10 ± 0.02 mg cm −3 ) relative to other estuaries (0.34 ± 0.02 mg cm −3 ) (SI Appendix, Fig. S2). Although we detected a spatial trend in these three estuaries, the differences were not significant. TP loading rates varied among estuaries, with the highest mean rate in the Broad River (118.4 kg ha −1 d −1 ), followed by Harney (90.5 kg ha −1 d −1 ), Shark River (72.8 kg ha −1 d −1 ), and Taylor Ridge (18.5 kg ha −1 d −1 ).
The Ca-bound inorganic P fraction was significantly different among transect sites and layers (SI Appendix, Tables S1 and S2). Mean Ca-bound P concentrations were on average 2.3 times higher in storm sediments (0.18 ± 0.03 mg cm −3 ) than in surface soils (0.08 ± 0.01 mg cm −3 ) across all sites. Discrete sites followed a similar trend compared to transect sites. Overall, Ca-bound P in surface soils was 2.4 (WSC-10, Harney), 4.3 (SRS-7, Shark), and 5.5 (WSC-13, Broad) times higher at downstream sites of each estuary compared to upstream sites. In contrast, concentrations in storm sediments did not differ considerably from downstream to upstream locations across all estuaries (SI Appendix, Tables S2 and S3).
Long-Term (2007 to 2017) Soil Properties. Soil BD was consistently lower upstream at SRS-4, intermediate at SRS-5, and higher at SRS-6 over the 11 y of sampling at Shark River mangrove sites. Soil OM content showed an opposite trend to that of BD (Table  1 and SI Appendix, Fig. S3 A and B). Soil TC was more variable from year to year and was significantly higher at SRS-6 relative to SRS-4 and SRS-5 (Table 1 and SI Appendix, Fig. S3C). Mean TN was significantly higher at SRS-4 compared to other sites, whereas mean TP was significantly different among all three  sites, with consistently higher concentrations at SRS-6 (0.23 ± 0.01 mg cm −3 ), followed by SRS-5 (0.15 ± 0.01 mg cm −3 ), and lower values at SRS-4 (0.12 ± 0.01 mg cm −3 ) (Table 1 and SI Appendix, Fig. S4 A and B). Soil N:P atomic ratios showed a significant increase with distance inland from the mouth of the estuary in all years, with overall means ranging from 33 (SRS-6) to 51 (SRS-5) to 82 (SRS-4) across years (Table 1 and SI Appendix, Fig. S4C).

Discussion
Landscape Sediment Deposition. In 2005, Hurricane Wilma's storm surge delivered sediment deposits across FCE mangroves (22,31). High surface-water levels (∼3 m) were observed within mangrove areas closer to the mouth of Shark River estuary in contrast to lower values (0.5 m) at upstream sites (22). This pattern was also registered by tidal gauges in the Harney, Broad, and Shark River estuaries and similar to storm surges after the passage of Hurricanes Andrew (1992) and Irma (2017; present study) across the southwestern FCE, including the Shark River estuary (22,31,38). However, the magnitude of Irma's surge was lower inside the forest (water levels ∼1 m) at mangrove areas closer to the mouth of the estuaries, where the storm reached maximum winds of 51 to 55 m s −1 during landfall (SI Appendix, Fig. S5). Irma's storm surge caused a large-scale mineral sediment deposition (<10 cm), with maximal values at the mouth of all estuaries relative to mid-and upstream locations. This deposition in mangrove forests was evi-dent up to 10 km from the GoM in southwestern FCE, with no deposition in mangrove areas 14-to 18-km upstream from the mouth of estuaries (e.g., SRS-4) (Fig. 4). This deposition is comparable to impacts observed after Andrew and Wilma (22,38). Greater deposition is observed in mangrove areas adjacent to the GoM (2 to 5 km) compared to upstream regions, and less deposition in the interior forest compared to the fringe mangrove zone (Fig. 4). Our results underscore the unique effect of these recurrent storms that create legacies of sediment and nutrient (mainly P) distribution in mangrove soils, despite different physical characteristics (e.g., direction, angle of approach, intensity) during landfall.
Notably, Irma's sediment deposition will have a potential contribution to the long-term net gain in soil elevation in mangroves of the southwestern FCE as has been reported by the longterm effects of Hurricanes Wilma and Irma on soil elevation change in the study area (39). Moreover, soil vertical accretion from Irma ranged from 1.8 to 3.9 cm, resulting in 6.7 (SRS-5) to 14.4 times (WSC-10) greater accretion during a single storm event compared to the long-term (100 y) annual accretion rate (e.g., SRS-6: 0.27 ± 0.04 cm y −1 ) estimated using radioisotopes (40). These mineral sediments increased the overall elevation capital at mangrove sites, ultimately contributing to the long-term net gain in elevation of the full soil profile, and potentially reducing their vulnerability to sea-level rise (39). Whereas our results showed a relative gain in elevation across FCE mangroves as a result of Irma, other soil processes, including erosion and compaction There was also a lateral gradient in deposition at each mangrove site, with higher deposition in the fringe mangrove zone compared to the interior forest (C and D). (E) In southwestern Everglades estuaries (e.g., Taylor River), storm-derived sediments were not observed in scrub mangrove sites. These storm deposition events contribute to phosphorus fertility gradients that control mangrove forest structure (e.g., tree height) and productivity patterns along estuaries and between the southwestern and southeastern Everglades (26,44). (39,41), soil shrinking and swelling, and subsidence (42) need to be considered in the long-term when assessing net changes in soil elevation. Indeed, soil vertical accretion in SRS-6 (4.3 cm) as a result of Wilma was reduced 1.3 times after 1 y due to erosion and compaction from shallow subsidence (21). Nevertheless, our findings show that changes in sediment volume as a result of Irma's storm surge have a significant overall positive effect in soil elevation in the long-term, contributing to the stability of neotropical mangroves.
Landscape Total Phosphorus Loading Rates. Our loading rates reflect not only the relative difference in TP density among estuaries, but also an expected spatial pattern where the lowest sediment deposition was observed in upstream regions of estuaries in southwestern FCE (e.g., WSC-11, Broad River) and eastern Florida Bay (i.e., Taylor Ridge). These differences in deposition are caused by a reduction in storm wind fields and energy attenuation due to the drag force induced by the presence of mangrove trees within the flow field, and other physical factors that regulate hydroperiod along the estuary (e.g., forest density, tree architecture, along-channel distance upstream, channel dimensions) (Fig. 4). Using data from the Advanced Circulation Model simulation (Coastal Emergency Risks Assessment; https://cera.coastalrisk.live/) of Irma's surge, we determined that maximum inundation depths as a result of storm surge at the vicinity of Broad, Harney, and Shark River mouths were 3.2, 3.0, and 3.0 m, respectively. This surge inundation pattern follows the same gradient in loading rates estimated for these estuaries in the southwestern region. In contrast, storm surge in Taylor River (southeastern FCE) and nearby mangrove areas (e.g., TS/Ph-7) was <1 m. This gradient in maximum inundation depths from north to south and west to east Everglades coincides with mean water levels registered by sensors in each estuary ( Fig. 1 and SI Appendix, Fig. S6). In fact, we found a negative linear relationship (R 2 = 0.62, P = 0.0024) between storm surge anomaly and sediment deposition (SI Appendix, Fig.  S7), suggesting the regulatory effect of the surge magnitude on forest maximum sediment deposition across our sites.
The spatial pattern of total P density and loading rates, both along and across estuaries, indirectly trace the storm impact across the FCE landscape. In areas where deposition was lower, mangrove tree height and extension were also low (e.g., WSC-8 and 11; upstream Harney and Broad Rivers, respectively) ( Fig. 4A and SI Appendix, Fig. S8). Moreover, when overlapping mangrove tree height distribution and Irma's TP deposition along those estuarine regions, a distinct positive association emerges between the highest TP concentrations and highest forest canopy (Fig. 4). This strong spatially explicit association demonstrates the positive role of hurricane-induced TP fertilization in maintaining forest productivity patterns. This is further supported by the low loading rate, associated with a lesser impact, estimated for the Taylor Ridge area (18.5 kg ha −1 d −1 ) located ∼70 to 80 km from Irma's direct path. Indeed, this area experienced weaker winds (30 to 40 m s −1 ) compared to the southwest Florida coast (SI Appendix, Fig. S5). The geophysical configuration of Florida and Madeira Bays enabled the storm surge to transport sediment over great distances into the lower Taylor River Slough, with maximal deposition in the shore fringe of Taylor Ridge (SI Appendix, Fig.  S8D). Mangrove tree height is predominantly low (<3 m) overall in Taylor River (22,43,44), while trees on the Ridge can reach heights up to 5 m (45). This difference in sediment deposition has been previously reported after the passage of Hurricane Wilma in 2005 (22). Similarly, Davis et al. (46) documented a similar sediment deposition pattern across the Taylor Ridge area as a result of Hurricane Irene in 1999.
Our TP loading rate estimates represent robust first-order values underscoring the influence of hurricanes in regulating P cycling in this oligotrophic karstic landscape (29,47). This allochthonous P input may also enhance belowground root growth and accumulation, thus increasing soil elevation and offsetting changes in sea-level rise, as shown by field P fertilization studies in karstic oceanic mangrove islands of Belize (48). To determine the potential impact of this "positive" storm-induced "pulsed fertilization" mechanism on mangrove forest development and longterm productivity, we compared our estimates with loading rates reported for wetlands in the Everglades Water Conservation Areas (WCAs) and other agriculture fields. The WCAs located south of the Everglades Agricultural Area are designed as surface water reservoirs that receive P-enriched waters directly from the Everglades Agricultural Area and from urban runoff. Increased nutrient loading above normal levels in the northern Everglades has led to the dramatic expansion of Typha species (49,50). We assumed that Irma's loading rate occurred in 1 d, given that the pulsing effect of the storm surge delivering a massive amount of TP-enriched sediment lasted between 6 and 9 h (Fig. 1). Since the Everglades is an oligotrophic ecosystem (29,30), any increase in TP could have significant impacts on vegetation structure (e.g., species shift, abundance, biomass) and function (e.g., productivity, carbon storage) (13,22,29,44,51). We found that Irma's TP mean loading rates in FCE mangrove wetlands (75 ± 21 kg ha −1 d −1 ) were 40 times higher than those used in "conventional" treatment wetland studies in the WCAs (1.9 ± 0.7 kg ha −1 d −1 ) in the northern Everglades ( Fig. 5 and SI Appendix, Table S4).
Moreover, using P loading rates (phosphate fertilizer) from agriculture fields for comparison (100 ± 7 kg ha −1 days-cycle −1 ), we show the massive TP input caused by Irma in 1 d (Fig. 5). This comparison is staggering given that fertilization of agriculture fields, expressed on per area basis, occurs over a growing season/ harvest cycle (corn: 105 to 120 d; sugar cane: 365 to 455 d) to enhance and maintain annual crop production ( Fig. 5 and SI Appendix, Table S4). Irma's P deposition made a significant contribution to the nutrient pool of mangrove soils across an extensive area. Indeed, the average contribution of P deposition from Irma's event to the soil nutrient pool ranged from 49.3% (Harney sites) and 49.6% (Shark) to 91.8% (Broad) across southwestern FCE estuaries. In contrast, TP deposited in the Taylor Ridge area only accounted for 11.2% of the total soil P pool. This TP input during hurricane landings is thought to maintain high levels of biomass and productivity in mangrove areas closer to the mouth of estuaries in southwestern FCE, in contrast to upstream areas dominated by scrub mangroves (43,44). Our Fig. 5. Comparison of total phosphorus (P) loading rates across different habitats, including (A) Everglades mangrove wetlands and treatment wetlands, and (B) agricultural crops (corn and sugar cane). For mangrove wetlands, rates were calculated using field data collected along transects and discrete sites to estimate total P loading rates as a result of Hurricane Irma's storm surge. The total area sampled was estimated using Arc View and the QGIS Geographic Information System. These rates represent total P deposition along a mangrove band of variable width adjacent to the estuary where actual field data were obtained. long-term data show that mangrove forests along the Shark River are indeed characterized by a distinct gradient in soil P fertility and corresponding shifts in N:P ratios from downstream to upstream locations (26,44,52,53), as confirmed after Irma (SI Appendix, Fig. S4 B and C). These P fertility gradients are associated with decreased mineral deposition (Ca-bound P) upstream of the estuary during hurricanes (22). Comparison of surface soil N:P ratios (Fig. 3C) across our mangrove sites further illustrate this pattern, where lower ratios (range: 20 to 27) are observed in mangrove areas downstream of estuaries and higher ratios (range: 30 to 60), indicative of P-limited conditions, in mid-and upstream (e.g., SRS-5, WSC-11) regions.
In addition to the conspicuous soil P fertility gradient contributing to differences in forest productivity in our study area (24,25,44), the spatiotemporal association between TP concentrations in green and brown-senescent leaves (leaf litter) and soil porewater soluble reactive P (SRP) concentrations (see below) also reflect a direct tree response to P loading by hurricanes, as illustrated by both Wilma's (2005) and Irma's impact (2017). Previous studies in the mid 1990s (54) and early 2000s (55) showed differences between TP concentrations in green leaves and fresh litter, including a significant interaction between nutrient retranslocation efficiency and species composition in the same sites along the Shark River estuary. For example, TP content of green leaves for the three mangrove species present in all Shark River sites were higher in forest stands at the mouth of the estuary (green leaves: 800 to 1,500 μg g −1 dry mass (dm); litter: 400 to 800 μg gdm −1 ) compared to lower concentrations (green leaves: 600 to 1,000 μg gdm −1 ; litter: 400 to 500 μg gdm −1 ) 9.9 km upstream (i.e., site SRS-5) from the coastline (54). These results are similar to mean TP concentrations we measured in leaf litter of all mangrove species in the same sites 3 (2008), 9 (2014), and 13 (2018) y after Wilma's impact (Fig. 6), when mangrove defoliation was >90% close to the estuary mouth. In the case of Wilma, forest canopy in these sites recovered by 2010 (24). Wilma's impact (22,24) was similar to Irma's impact (September 2017) by causing widespread defoliation and tree mortality (<10%). Our TP concentrations in leaf litter from litterfall productivity studies in 2008, 2014, and 2018 were similar to values reported in the 1990s, including the highest concentration in Avicennia germinans-senescent leaves (910 ± 36 μg gdm −1 ) and in the case of Rhizophora mangle, the lowest concentrations upstream (SRS-4) versus downstream (SRS-6) of the Shark River ( Figs. 6 and 7). The lowest mean concentration was consistently registered in R. mangle leaves across all sites and years (Fig. 7). Similarly, TP concentration in leaf litter of all species increased from SRS-4 to SRS-6 ( Fig. 6). This species-specific difference in foliar nutrient content indicates that mangrove species respond physiologically different to TP availability, although in a similar spatial pattern. For example, Laguncularia racemosa, a shadeintolerant species, reaches the highest growth rate (13,56) during canopy openings caused by lightning, natural tree mortality, or hurricane impact, thus triggering a competitive advantage that increases its dominance as observed in SRS-6 (13,44).
This TP difference among species is apparent when anayzing leaf litter concentrations at the SRS-6 site where the highest hurricane deposition occurred after Irma's impact (Fig. 2). Indeed, L. racemosa leaf litter TP concentration reached a mean value of 1,565 ± 331 μg gdm −1 in the period from January to May and then from October to November (1,509 ± 252 μg gdm −1 ) in 2018 (Fig. 7); these estimates contrast with the lower mean annual values estimated for this species in 2008 (635 ± 55 μg gdm −1 ) and 2014 (627 ± 37 μg gdm −1 ), which are 50% lower compared to the annual mean measured in 2018 (1,186 ± 145 μg gdm −1 ) (Fig. 7). Both A. germinans (1,339 ± 111 μg gdm −1 ) and R. mangle (655 ± 62 μg gdm −1 ) also show a signficant increase in litter foliar TP concentrations in early 2018, 3 mo after Irma's impact, thus underscoring the interspecies differences in nutrient uptake as indicated by mangrove isotopic and biomarker studies in this estuary (57,58). Furthermore, the increase in foliar TP in all species after Irma's impact could be explained by the significant increase in soil porewater SRP concentrations (Fig. 8), especially in the case of SRS-6, where values were almost seven times higher in the 2017 wet season (14.5 ± 0.5 μM) than the average value measured in the period 2008 to 2016 (2.0 ± 0.2 μM) (Fig. 8); this value is significantly higher compared to the SRP values reported for eutrophic coastal environments like the Mississippi River coastal deltaic plain, where concentrations in the water column can reach up to 5 μM (59). This SRP "spike" indicates that the Ca-bound P in sediments deposited by Irma became readily available for plant uptake as shown in late 2017 and early 2018 across all species in both SRS-5 (1.7 μM) and SRS-6 (4.5 μM, dry season) (Fig. 8).
These site differences in SRP concentrations are consistently observed from 2008 to 2018 (Fig. 8). Although we do not have mangrove foliar litter TP data for the year after Wilma's impact (i.e., 2006), our 2018 results post-Irma indicate that the forest canopy P stochiometry was back to predisturbance levels after ∼2 y. This recovery is indirectly reflected on the conspicuous "U" shape of the monthly TP values for all species where the lowest value occurs in June-July at the peak of litterfall productivity (24,44) (Fig. 7); this temporal pattern in TP is associated with retranslocation efficiency at different levels (SI Appendix, Fig.  S9). For example, in 2008 and 2014, L. racemosa mean efficiency (65%) was similar to R. mangle (65.9%), despite its lower foliar TP content over the year, whereas A. germinans showed the lowest mean retranslocation value (46%), possibly as a result of its highest concentration throughout the year (SI Appendix, Fig.  S9). This seasonal pattern in foliar TP and retranslocation efficiency has been observed consistently in other years (e.g., 2001, 2010, 2012, 2014) before and after hurricane impacts like Wilma and Irma (55).* Therefore, mangrove wetlands with soil porewater SRP limitation and lack of hurricane P inputs as those observed upstream of the estuary (SRS-4) (Fig. 8), have higher P retranslocation efficiency and foliar N:P ratios compared to nutrient-rich, near-coast mangroves (e.g., SRS-6). These tradeoffs are indicative of the strong link between soil P availability and the phenotypic plasticity of mangrove leaves in response to nutrient gradients across the FCE landscape as a mechanism to enhance nutrient conservation.
Based on these findings and previous work, we propose that mangrove areas that do not receive this Ca-bound P pulse have less soil P accumulation triggering P limitation conditions that constrain mangrove forest development (e.g., scrub mangroves in Taylor River). This pattern is apparent in the polygons used to estimate P loading rates along the estuaries where low-stature trees were generally present at the end of the transects or upstream regions of estuaries (SI Appendix, Fig. S8); this is the case in SRS-4, where maximum tree height is <6 m (43,44) and storm deposits were not observed after Wilma (22) or Irma's impact. The relative influence of the hurricane-induced P fertilization mechanism to near-coast mangroves in the FCE is comparable to other natural P subsidies (i.e., guano) delivered by seabirds to mangrove islands within karstic oligotrophic environments in the Yucatan Peninsula, Mexico (60). This P fertilization by seabirds increases soil P concentrations (up to eight times the regional mean during nesting season) in mangrove islands, thus reducing nutrient limitation and inducing changes in mangrove tissue stoichiometry and physiological processes, as observed in Shark River mangroves.
Despite the well-documented destructive effect of hurricanes on FCE mangrove ecological attributes (16,24,31,61), our results show that hurricanes are major drivers controlling fertility gradients in Everglades mangroves (22), and therefore create a positive effect in maintaining observed mangrove spatial distribution and productivity patterns (43,44). Hurricanes can setback/reset the succession stages in mangroves, but help to self-maintain the system in a relatively steady state at the landscape scale (14), where initial conditions of forest structure (e.g., basal area, canopy gaps) and rates of sapling recruitment, along with P storm deposition, will largely control the rate of recovery and change of ecological attributes following disturbance (13,56). Furthermore, the interaction between hurricane physical properties (intensity, duration, wind velocity) (62) and forest gaps and patches create a highly dynamic forest matrix postdisturbance. In fact, mangrove structural development in the FCE does not mature to a typical "old growth forest" (13,63) as those observed in other regions where mangrove forests are not impacted by hurricanes, and as a result, the forest canopy can reach heights >50 m (64). Additionally, the landscape variability in sediment deposition (inland and laterally) observed across FCE mangroves could also be the result of the interaction between the geomorphology of the coast, local microtopographic changes, and storm physical properties (46) (SI Appendix, Fig. S10). Thus, tropical storms along GoM coastlines may impart both positive and negative effects on mangrove forests, which need to be considered when assessing and modeling interactions between biogeochemical and intranutrient cycling processes (e.g., C storage and sequestration, P resorption efficiency) (23,24), especially over decadal scales. Hurricanes can enhance soil elevation relative to sea level and maintain high mangrove biomass and productivity rates, particularly when tree mortality and defoliation could be quickly offset in the long term by P fertilization by influencing reproductive output and tree and seedling growth rates (23, 24) (SI Appendix, Fig. S10). Increases in tree basal area poststorm due to increases in recruitment and density as a result of P fertilization will increase mangrove structural complexity enhancing sediment trapping in future storms. Indeed, sedimentation rates have been correlated with mangrove roots and pneumatophore density in mangrove forest stands (27). Mangrove functional roles, such as sinks or source of atmospheric C, can also be impacted by storms, ultimately modifying the C balance with the atmosphere and the capacity of mangrove forests to mitigate C emissions (65). This feedback mechanism should be explored in other coastal areas where mangroves are frequently impacted by hurricanes controlling tree height as a result of structure damage (64,66,67). As mangroves migrate inland with sea-level rise, this hurricaneinduced P deposition may offset nutrient limitation and enhance mangrove development and productivity (SI Appendix, Fig. S10). Moreover, climate-induced changes in environmental conditions have the potential to change the distribution, abundance, and composition of mangrove wetlands worldwide (68)(69)(70).
We propose that the mangrove response to large-scale P loading from hurricanes may assist neotropical mangroves in the GoM and the Caribbean region over long time-scales to withstand the impacts of both sea-level rise and nutrient limitation. Indeed, this persistence and mangrove colonization and transgression in the GoM dates back to the mid-Holocene (71), indicating the complex interaction among geomorphology, within-site environmental gradients, and regional disturbance regime (SI Appendix, Fig. S10). This is particularly significant to south Florida due to the high recurrence of tropical storms, where complex interactions between hurricane legacies, climate change, and the rapid increase in sea-level rise during the last decade (72) will potentially affect mangrove structural and functional attributes and recovery trajectories after disturbance, thus affecting the provision of ecosystem services, such as carbon sequestration. Further assessment of this complex mechanism is needed to quantify the long-term interaction between natural positive (fertilization) and negative (forest mortality/structural damage) effects in neotropical mangrove wetlands.

Methods
Study Sites. This study was conducted in mangrove forests of Everglades National Park, south Florida ( Fig. 4A and SI Appendix, Fig. S1). This area is part of the Florida Coastal Everglades Long-Term Ecological Research (FCE-LTER) program (73) (https://fcelter.fiu.edu/). Mangrove forests in the southwestern FCE form a continuous band of riverine mangroves that extends along a freshwater-estuarine gradient, from the lower reaches of freshwater marshes in SRS to the GoM, a distance of about 15 to 20 km that sets the limits of periodic saltwater and tidal influence (44). The FCE-LTER sites are located along two major drainage basins, the Shark River Slough (southwestern FCE) and Taylor River Slough (southeastern FCE). In December 2000, three mangrove sites were established each along the SRS (SRS-4, SRS-5, SRS-6) and TS/Ph (TS/Ph-6, TS/Ph-7) basins to monitor structural and functional attributes and soil properties (44,53). The SRS-6 site is located 4.1 km from the mouth of the estuary, while SRS-5 and SRS-4 are 9.9 and 18.2 km from the river mouth, respectively. A new permanent site (SRS-7) was established at the mouth of the estuary after Hurricane Wilma to assess mangrove resilience since this area was heavily impacted. Mangrove sites along the Taylor River are located about 1.5 km (TS/Ph-7) and 4 km (TS/Ph-6) inland from Florida Bay. Another site at the Buttonwood Ridge was established ∼1 km east from the mouth of Taylor River. The Ridge is a high elevation (∼0.5 m in height) depositional berm that extends ∼60 km across the southern tip of Florida and isolates these scrub forests from the direct influence of Florida Bay (22,46). Similar to Shark River sites, four and three riverine mangrove sites were established along the freshwater-estuarine gradient in Harney and Broad Rivers (southwestern FCE), respectively. Mangrove sites along the Harney River were located ∼14 km (WSC-7), 10 km (WSC-8), 6 km (WSC-9), and 2 km (WSC-10) from the mouth of the estuary. In Broad River, sites were established at 10 km (WSC-11), 7 km (WSC-12), and 3 km (WSC-13) from the river mouth. Detailed information on sites is provided in SI Appendix, SI Text.
Hydrology and Storm Surge. Changes in water levels associated with Irma's storm surge were measured at the FCE-LTER mangrove sites. Continuous water levels relative to soil surface have been measured at 1-h intervals at all sites since December 2000 using ultrasonic water level recorders (Infinities USA). Instruments were installed in the interior of each mangrove site about 50-to 80-m inland from the edge (22). At Harney and Broad Rivers, longterm real-time water-level data were obtained from permanent gauges installed along the main river channels using the Everglades Depth Estimation Network (https://sofia.usgs.gov/eden).
Sediment Deposition and Nutrient Inputs. We used a comprehensive spatial sampling design (i.e., transect sites and discrete sites) to compare gradients in storm deposits and nutrient inputs across all mangrove sites (Fig. 4A). A similar approach was previously used in the study area after the passage of Hurricane Wilma in October 2005 to characterize storm-related sediments and nutrient inputs in mangroves (22). We measured the thickness, distribution, and physico-chemical properties (BD, OM content, TC, TN, TP, TOC, TIC, and Ca-bound P) of storm sediments using duplicate cores collected at different distances along transect sites and discrete sites across the FCE. Sampling was conducted in January (Shark and Harney River sites), early March (Taylor Ridge), and early April 2018 (Broad River sites). All soilsediment cores were collected with a piston core (2.5-cm diameter × 15-cm length) and sectioned into two layers, storm sediments (variable depths) and surface (top 10 cm) mangrove soils, and the depth of each layer was registered. Each layer was stored separately in prelabeled 50-mL centrifuge tubes, placed on ice, and transported to the laboratory for further analyses. The storm layer was easily distinguished from the mangrove layer due to its gray color, fine sandy texture, and organic-free content. For the transect data, we used a randomized-block ANOVA design to test for differences in physico-chemical variables among sites, distance along transects, and layers (soil vs. sediment). Data collected within the discrete sites were analyzed separately with a two-way ANOVA, with sites and layers as main factors. Unless otherwise stated, data presented are the means (±1 SE) of untransformed data. Statistical analyses were performed with PROC MIXED (SAS Institute). Specific details about the sampling design and chemical and statistical analyses can be found in SI Appendix, SI Text.
Long-Term Soil Properties. Soil TC and nutrient (TN and TP) concentrations, BD, and OM content, and porewater SRP were measured at the Shark River mangrove sites from 2007 to 2017 (74). In July or August each year, small, 3-cm diameter cores (top 10-cm surface soils) were collected in triplicate from SRS-4, SRS-5, and SRS-6 using 60-mL syringe barrels that were pushed into the soil while holding the plunger at the soil surface to minimize compaction. Soil cores were processed in the laboratory using standard analytical protocols. See SI Appendix, SI Text for specific chemical and statistical analyses information. Porewater SRP concentrations have been measured during the dry and wet seasons from 2008 to 2018 at four permanent sampling stations within each of the two plots following procedures described by Castañeda-Moya et al. (44).
Litterfall Collection and Nutrient Content. Detailed information on litterfall collection, processing, and P analyses is provided in SI Appendix, SI Text.
Landscape TP Loading Rates. We used data collected along transects and discrete sites to estimate TP loading rates as a result of Irma's storm surge and sediment deposition in the Broad, Harney, and Shark River estuaries and the Taylor Ridge (Little Madeira Bay). Based on the actual sediment deposition data measured at different distances along transects and the position of transects along each estuary, we delimited a number of polygons to estimate total deposition per estuary and assess spatial variability (SI Appendix, Fig. S8). The polygons were delineated following the contour of the estuary and the length of the transect; this length per estuarine region was considered the width of the polygons. The total area per polygon was estimated using ArcGIS (ArcGIS Desktop 2010, ESRI Inc.). We only estimated the total deposition along the side of the estuary where actual field data were obtained. We considered these values conservative, first-rate estimates of this "new" TP input (opposite to regenerated in situ) into the mangrove wetlands at the landscape level.