Solar radiation and soil moisture drive tropical forest understory responses to experimental and natural hurricanes

Tropical forest understory regeneration occurs rapidly after disturbance with compositional trajectories that depend on species availability and environmental conditions. To predict future tropical forest regeneration dynamics, we need a deeper understanding of how pulse disturbance events, like hurricanes, interact with environmental variability to affect understory demography and composition. We examined fern and sapling mortality, recruitment


INTRODUCTION
Cyclonic windstorms (named hurricanes in the Western Hemisphere) are disturbances that affect terrestrial ecosystems in space and time (Bellingham et al., 1995;Brokaw & Grear, 1991;Lugo, 2008;Tanner et al., 1991;Uriarte et al., 2004). Hurricane occurrence over millennia has shaped ecosystem structure and species composition in tropical, subtropical, and temperate coastal and inland areas, leading to selection pressure on affected biota (Griffith et al., 2008;Hogan et al., 2018;Ibanez et al., 2019;McLaren et al., 2019;Tanner & Bellingham, 2006;Zimmerman et al., 1994). Hurricanes, sudden and catastrophic forest disturbance events, are both impressive and ecologically dramatic because they produce visible ecosystem change and can induce rapid fluctuations in the demographics of biota and abiotic conditions (Hogan et al., 2018;Patrick et al., 2022;Uriarte et al., 2019). Structural and species compositional changes following hurricanes can affect emergent ecosystem function (e.g., biogeochemical cycles, carbon uptake, erosion rates, and hydrology) and are thus important to understand to predict potential long-term ecosystem responses to climate change Lin et al., 2020;Michener et al., 1997;Uriarte et al., 2004).
Forests differ in resistance (i.e., the magnitude of perturbation from historical baselines given a particular disturbance event) and resilience (i.e., time and ability to recover to a comparable or baseline forest state) to hurricanes (Attiwill, 1994;Laurance & Curran, 2008;Nikinmaa et al., 2020;Peterson, 2019;Webb, 1958). This is primarily a function of the historical disturbance regime of the forest (i.e., the prevailing frequency and severity of disturbance events), and the time since the last disturbance (i.e., the ecosystem state at the time of the hurricane) (Hogan et al., 2018;Ibanez et al., 2019;Lin et al., 2020;Peterson, 2019), although tree species diversity has been shown to modulate disturbance severity (Tanner & Bellingham, 2006;Zimmerman et al., 1994). The timing, frequency, and severity of cyclone disturbances are key to shaping the long-term forest resistance and resilience states via their effects on forest damage and subsequent forest succession, which interact with longer-term anthropogenic-driven environmental change pressures (e.g., climate warming and elevated CO 2 ) on ecosystems (Dale et al., 2001;Everham & Brokaw, 1996;Johnstone et al., 2016;Lugo, 2008). As a result, individual hurricane disturbances can be understood as distinct disturbance pulses overlaid on less dramatic, but longer-term environmental change pressures .
Hurricanes are expected to become an increasingly important force for ecological and social disturbance in a future affected by climate change. Recent work has shown that hurricanes are slowing in transit speed (Kossin, 2018), increasing in wind speed and intensity (Elsner et al., 2008;Kossin et al., 2013;Sobel et al., 2016;Velden et al., 2017), increasing in frequency (Emanuel, 2013), and disturbing previously unaffected areas (Altman et al., 2018), especially as they migrate further inland in coastal areas (Wang & Toumi, 2021). Therefore, it is likely that hurricanes and their interactions with other climatechange-driven disturbances (e.g., drought) will become increasingly important agents of change for forest dynamics, shaping forest structure and composition at local and global scales McDowell et al., 2018). In this context, it is critical to understand how an increased frequency or severity of hurricanes will affect forest ecosystems in the coming decades and how other environmental drivers exacerbate or mitigate hurricane effects .
Hurricane disturbances in forest ecosystems induce vegetation dynamism through their direct effects (i.e., damage) and by altering environmental conditions (e.g., solar radiation, soil moisture, debris deposition, and nutrient cycling) (Everham & Brokaw, 1996;Lugo, 2008). Individual hurricane disturbances represent distinct events that drastically influence successional trajectories of a forest by shaping forest structure (e.g., maximum tree heights and biomass), tree population dynamics (e.g., recruitment and mortality rates), and tree growth, over decades and centuries (Bellingham et al., 1995;Brokaw & Grear, 1991;Lin et al., 2020;McLaren et al., 2019;Zimmerman et al., 1994). For example, tropical forest structure and canopy architecture are shaped by cyclone disturbance frequency to varying degrees (de Gouvenain & Silander Jr., 2003;Ibanez et al., 2019;Peereman et al., 2022). Hurricane-induced compositional changes can lead to varying successional trajectories of forests, altering species composition, and relative abundances of forest taxa in novel and ecologically important ways, especially when combined with other disturbances (Vandermeer et al., 2000). Timing, frequency, and severity of hurricanes are key modulators of hurricane effects on vegetation compositional dynamics in space and time . Given that hurricanes are increasing in severity and frequency, it is imperative to understand forest responses to multiple hurricanes over time, and how the independent environmental drivers of hurricane disturbance and recovery (e.g., canopy opening vs. debris deposition, increased rain throughfall, and soil moisture) modulate vegetation response (Turton, 2008).
Multiple disturbances span a gradient that can be nonadditive to compound each other or can counteract each other in determining trajectories of ecosystem change over time (Buma, 2015;Micheli et al., 2016;Tye et al., 2016). Shifts in community composition postdisturbance may develop slowly over time then return to baseline prehurricane composition, or there may be key threshold responses (i.e., "tipping points") in response to disturbanceinduced environmental changes, which result in very longterm changes that may not return to predisturbance conditions in our lifetime (e.g., arrested succession; Brokaw et al., 2012;Reyer et al., 2015;Zimmerman et al., 2021). The drivers behind these different successional trajectories and the subsequent consequences for ecosystem function are not entirely clear. If the increased frequency of hurricane disturbances is leading to changes in tropical forest community composition via environmental alterations, then changes may be most readily observed in the understory community of forests (i.e., tree saplings, ferns, and other ground-dwelling plants), given their general faster response times than large trees and that understory tropical plant communities have different environmental controls than the larger adult tree community (Baraloto & Goldberg, 2004;Brokaw & Scheiner, 1989;Montgomery & Chazdon, 2001;S. Murphy et al., 2016). The Canopy Trimming Experiment (CTE) was designed (Shiels et al., , 2015; see description in Methods) to examine tropical forest community and ecosystem change in response to repeated hurricane disturbances. Here, we evaluate how multiple sequential canopy opening events via both experimental trimming and natural hurricanes affect tropical forest understory plant dynamics and species compositional change. We ask two research questions:

Study site
The CTE is located in the subtropical montane forest of the Luquillo Experimental Forest (LEF) in northeastern Puerto Rico (18.3215 N,65.8141 W) Shiels & Gonz alez, 2014). The area has an aseasonal everwet climate, where typically no month receives <100 mm of rainfall (S. F. Murphy et al., 2017). The intrusive volcaniclastic geology of the forest has led to the development of highly weathered tropical ultisol and inceptisol soils with low fertility (Mount & Lynn, 2004;Scatena, 1989;Seiders, 1971). The CTE is located within the mid-elevation (340-485 m above sea level) tabonuco forest type of the LEF, which is dominated by Dacryodes excelsa and Prestoea acuminata var. montana, has a single-tiered canopy that lacks canopy emergent trees and reaches about 20 m (Brokaw & Grear, 1991).

Natural hurricanes affecting the LEF and CTE
The LEF experiences frequent hurricane disturbances with a historic return interval estimated at 60 years (Booze et al., 2004;Scatena et al., 2012). Two strong and damaging hurricanes (San Ciprian and San Felipe) affected the LEF in the 1930s (Scatena & Larsen, 1991). Several storms have occurred since the 1930s; however, none severely damaged the LEF until Hurricane Hugo in 1989 (Zimmerman et al., 1994). Thus, the LEF was in a relatively mature successional state in 1989 with tall trees and high levels of tree biomass before Hurricane Hugo severely damaged the forest (Brokaw & Grear, 1991). Another strong (category 3) hurricane, Georges, affected the LEF in 1998, and although being similar in strength to Hurricane Hugo, the pathway of Hurricane Georges across Puerto Rico resulted in less damage to the LEF (Hogan et al., 2016). The CTE was established at the end of 2002, and treatments began 2 years after plot establishment Shiels & Gonz alez, 2014; see more details below). Of the seven hurricanes that have occurred since the establishment of the CTE (Knapp et al., 2010), Hurricanes Irene (22 August 2011), Irma (7 September 2017), and Maria (20 September 2017) passed nearest to the LEF and CTE, with minimum distances from the eye of the hurricane to the CTE of 16.6, 37.4, and 74.7 km, respectively. Hurricane Irene did not lead to noticeable forest damage at the CTE. However, Hurricanes Irma and Maria collectively caused significant damage to the entire LEF including the CTE. Irma partially defoliated the forest canopy and led to some canopy branch breakage, then-2 weeks later-Maria caused widespread canopy and tree stem damage, with complete canopy loss at the landscape level and unprecedented levels of tree mortality (Y. Feng, Negr on Ju arez, et al., 2018;Uriarte et al., 2019). Hurricane Maria likely caused more severe damage to the LEF than any storm since Hurricanes San Ciprian and San Felipe in the 1930s (Scatena et al., 2012;Scatena & Larsen, 1991;Uriarte et al., 2019). As Hurricanes Irma and Maria both occurred in September 2017, and Maria caused much more damage, we consider them to have created one combined impact on the forest. Henceforth, we will refer to the combined impact of both Hurricanes Irma and Maria as the response of the understory to Hurricane Maria.

The CTE: An experimental simulation of hurricane effects
The CTE is an ongoing experiment that simulates the impacts of hurricanes on the forest canopy by removing (i.e., trimming) tree branches to mimic hurricaneinduced canopy loss (Richardson et al., 2010;Shiels et al., 2010;Shiels et al., 2015). The experiment was designed with three replicates (blocks A, B, and C), each comprised of four plots (30 Â 30 m, which included a 5-m edge buffer around a 20 Â 20 m data collection area). The 20 Â 20 m data collection area was divided into 16 subplots (subplot 4.7 Â 4.7 m with 0.4-m-wide walking paths between the subplots) (see figure 3 in Shiels & Gonz alez, 2014;Zimmerman, 2017). Before the implementation of experimental treatments, plots were assessed, vegetation structure was recorded, and abiotic environmental equipment was installed and monitored to collect approximately 1 year of pretreatment data. Then, from October 2004 to June 2005, four treatments were applied, each to one plot of the three replicate blocks, in a factorial block design: (1) Trim + Debris, the hurricane simulation, where trees were trimmed to remove all branches, and the trimmed material (debris) was deposited on the forest floor within the experimental plot area; (2) Trim + No debris, where trees were trimmed and debris collected and added to the forest floor of plot; (3) No trim + Debris where trees were not trimmed but debris was added; and finally, (4) No trim + No debris, an untreated control where trees were not trimmed and no debris was added. Treatments were implemented (in order: blocks B, C, and A, with each block taking an average of 75 days to trim the canopy and deposit the debris).
For trimmed plots, all trees ≥15 cm diameter at breast height (dbh; measured at 1.3 m above the ground surface) had branches <10 cm diameter trimmed off at the main stem, while trees 10-15 cm dbh were cut off at 3-m height and palms ≥3-m height had all leaves removed (see Shiels & Gonz alez, 2014, for further details). In 2004, trimmed debris was weighed and uniformly distributed upon the Debris addition plots Shiels & Gonz alez, 2014). A suite of measurements within the plots occurred at frequent intervals (see more details below) until 2014 when trees in the plots were retrimmed. After the retrim, debris manipulation was discontinued. Thus, only the Trim + Debris (Hurricane simulation) treatment was re-implemented, while the No Trim + No Debris (Control) plots remained as before. For our statistical analyses, we focused on these two treatments-Trim + Debris (Hurricane simulation) and No Trim + No Debris (Control)-to examine the understory plant community demographic and compositional responses to repeated hurricane disturbances.

Vegetation census methods: Saplings
For each of the 12 CTE (20 Â 20 m) plots (3 replicate blocks), tree censuses were carried out approximately annually (Zimmerman, 2020b). Tree censuses of the CTE commenced on 1 March 2003, 15 October 2004, 17 September 2009, 17 October 2008, 12 November 2009, 1 October 2014, 26 October 2015, 24 October 2016 December 2017, and 1 November 2018. Plots were censused systematically in the same order each year to reduce differences in the length of inter census periods, with each census completed within 6 months. Freestanding woody vegetation census protocols follow the Center for Tropical Forest Science (now Smithsonian Forest-GEO) methodology as outlined by Condit (1998). Saplings are classified as stems ≥1 cm but <10 cm in dbh; these were identified to species, assessed for living/ dead status, assessed for damage, and tagged with a unique number. New stems were included in a census (as recruits) when they reached a size of ≥1 cm dbh. Heliconia caribaea, an abundant and iconic understory monocot of the LEF, was censused in the CTE, recording alive/dead status, and the number of culms for each alive individual >50-cm tall, the minimum size at which H. caribaea individuals are included in tree censuses. In all analyses included here, we define saplings as freestanding (woody and understory shrub) stems ≥1 cm but <10 cm dbh plus Heliconia individuals >50-cm tall.

Vegetation census methods: Ferns
Terrestrial ferns are long-lived perennials that can dominate the herbaceous layer of the LEF understory (Sharpe & Mehltreter, 2010;Sharpe & Shiels, 2014). Ferns represent understory herbaceous vegetation responses in our analyses. Fern abundance was censused within each subplot (16 subplots 4.7 Â 4.7 m) of every CTE plot (20 Â 20 m excluding the 5-m border buffer), as described in Sharpe and Shiels (2014). This gave a total effective fern census area of 353 m 2 /plot. All ferns rooted in the soil that were ≥10 cm in height at the longest leaf length were identified and counted. Ferns were counted early in January at the end of the calendar years 2003 through 2009, except for 2004 when the first canopy trimming treatment was applied. After the 2009 count, fern censuses were discontinued but were resumed in 2014 with an October census just before the second trim that began in November 2014 and continued through 2018 (Sharpe, 2021).
In addition to the count of ferns across the whole plot, individual terrestrial ferns ≥10 cm (longest leaf length) were tagged, identified, counted, and measured at the longest leaf length in 5 designated seedling monitoring quadrats (1.0 Â 3.5 m), located within 5 of the 16 subplots of each plot (Zimmerman, 2020a). Monitoring took place in October each year from 2003 to 2019, except for 2004 when it occurred in June (for more details, see Sharpe & Shiels, 2014;Shiels et al., 2010). Plant mortality was noted, and new plants were tagged and recorded as recruits when they reached 10 cm in leaf length.

Environmental measurements: Soil moisture and solar radiation
Data on environmental parameters including temperature, rainfall, canopy throughfall, canopy openness and understory solar radiation, soil biogeochemistry, and soil moisture have been extensively monitored in the CTE since 2002 Shiels et al., 2015;Shiels & Gonz alez, 2014;Van Beusekom et al., 2020). Environmental data were collected on the 20 Â 20 m plot interiors, excluding the 5-m border zones to minimize edge effects. Solar radiation and soil moisture are among the most important factors that can elicit a response from the understory plant community. Accordingly, our analysis focused on the effects of soil volumetric water content (VWC) and incoming solar radiation on the understory vegetation (Gonz alez & Van Beusekom, 2021).
From 2003 to 2006 at 3-month intervals, gravimetric soil water content (GWC) was measured manually from collected soil samples (Richardson et al., 2010). Then in 2015, soil volumetric moisture sensors were installed at 5 cm depth (Decagon Devices 5TM, METER Environment, Pullman, WA, USA), and Campbell Scientific SC616 soil water content reflectometers were installed at 0-15 cm depth (Campbell Scientific, Logan, UT, USA) to continuously monitor soil VWC. In 2015, soil bulk density and GWC were repeatedly measured near installed sensors, for comparison of concurrent GWC and VWC measurements for 1 year. To create a single data series, GWC measurements were converted to VWC using soil bulk density and smoothed VWC data from the soil reflectometers. Soil VWC data were checked against satellite data from the Advanced Microwave Scanning Radiometer 2 (AMSR2), confirming that increases and decreases in soil VWC tracked satellite observations over time (Van Beusekom et al., 2020).
For incoming solar radiation, hemispherical canopy photographs from within the CTE plots (before 2015) were combined with measurements from solar pyranometers (L1200X, Campbell Scientific, Logan, UT, USA). Beginning in December 2003 (i.e., 1 year before treatment implementation), hemispherical photographs were taken at low light (i.e., dawn) using a fisheye lens at 1-m height above the center of five random subplots and the center and corners of the 20 Â 20 measurement area of each CTE plot (n = 10 per plot per sampling period; Shiels et al., 2010;Van Beusekom et al., 2020). Additional photographs were taken at roughly 6-month intervals (i.e., 2-3 times per year) between 2005 and 2012. The hemispherical canopy photographs were converted from color to black and white images using threshold iterative separation of background and foreground with the Ridler and Calvard method (Bachelot, 2016). Converted images were analyzed using the Hemiphot algorithm using R (ter Steege, 1997Steege, , 2018, which estimates canopy openness and leaf area index and calculates photosynthetically active radiation (PAR) for one or several days using the sun's position and track over the canopy. A complete time series of PAR was calculated by linearly interpolating PAR for each day of the year as a fraction of the previous and next photograph's PAR on that day, which roughly accounts for canopy cover recovery from trimming and seasonal dynamics (Van Beusekom et al., 2020). Photosynthetically active radiation was translated to solar radiation by multiplying by a factor of 2 (Escobedo et al., 2009) and calibrating values to solar radiation measurements from above the forest canopy. A complete time series of solar radiation data was constructed by validating and interpolating values derived from canopy photographs with measurements from the pyranometers (post-2015). For further details, see Van Beusekom et al. (2020).

Data analyses
All data analyses were done in R 4.0.3 (R Core Team, 2020). For our first research question, which concerns how repeated simulated and natural hurricanes have affected understory plant demographic rates, we calculated the demographic rates of ferns and saplings for each plot in the CTE using the following equations: where N t is the number of plants or stems alive at a given second of two censuses, N dead is the number of dead plants or stems for that census, N recruits is the number of newly recruited plants or stems for that census, and t is the census interval in months (Condit et al., 1995). We chose a monthly census interval because of the variation in calendar dates for tree censuses among years. Fern demographic rates were calculated using a total of five plant census subplots, within each of the CTE treatment plots. Sapling demographic rates were calculated using tree census surveys of the entire area of each 20 Â 20 m CTE plot. Mortality and recruitment rates are difficult to measure, variable, and not distributed normally; therefore, the assumptions of standard statistical tests (e.g., ANOVA) are not typically met (Condit et al., 1995;Sheil & May, 1996). The Scheirer-Ray-Hare test, which is a two-factor nonparametric statistical test that can test for differences among factors and their interactions when designs are balanced (Sokal & Rohlf, 1995), was used to statistically test for demographic rate differences among treatments over time and their interaction. The test was run separately for mortality and recruitment rates for both saplings and ferns. This is just like doing a multiple factor Kruskal-Wallis test, which we explored and gave congruent results. Fern and sapling densities were calculated at the hectare scale for all plots and then for each species by treatment. We searched for common species (>100 individuals in the CTE) whose demographic responses to the treatments (i.e., stem density changes) typified the demography for specific plant life-history strategies. Based on those selection criteria, three fern species (Thelypteris deltoidea, Cyathea borinquena, and Cyathea arborea), six woody species (Cecropia schreberiana, P. acuminata var. montana, Psychotria berteroana, D. excelsa, Casearia arborea, and Tetragastris balsamifera), and the clonal understory monocot H. caribaea were selected as species which exemplified responses. Their densities were plotted by treatment over time.
For our second research question, which examined changes in species composition over time for the CTE hurricane simulation treatments and natural hurricanes versus control plots, we used principal coordinate analysis (PCoA), also referred to as metric multidimensional scaling (MDS) ordination (Gower, 1966). Community (i.e., plot Â species) matrices for hurricane simulation and control plots were constructed by combining fern and sapling communities for censuses that took place in the same calendar year. Combining fern and sapling censuses led to understory plant community counts for 2003, 2004, 2007, 2008, 2009, 2011, 2012, 2013, 2014, 2015, 2016, 2017, and 2018. Species abundances per plot were relativized, and the Bray-Curtis distance (Bray & Curtis, 1957) was used to compute compositional dissimilarity among plots over time. The PCoA was done using the "capscale" in the vegan R package (Oksanen et al., 2020) before plotting the plot scores in the MDS ordination space. We also plotted the species scores, noting where tree and fern taxa sorted out concerning the first two ordination axes.
Multidimensional scaling ordination site scores were regressed against environmental predictors (soil moisture and solar irradiance) using linear mixed-effects models (LMMs). For each of the first two MDS axes, the bestfitting LMMs that included fixed predictors of the environmental variable (either soil VWC or solar radiation), year, treatment (hurricane simulation or control CTE plot), and interactions between year and the environmental variable, and between year and treatment were selected. Random intercept terms for year and block were considered. The "step" function from lmerTest package (Kuznetsova et al., 2017) was used to find the best-fitting models based on the Akaike information criterion, by doing backward selection of random effects, then fixed effects. Models themselves were fit using the lme4 package (Bates et al., 2014). We tried to fit models that included both soil VWC and solar irradiance and their interaction in the same LMM; however, we were limited by the degrees of freedom in the dataset (n = 60 for the solar irradiance models and n = 42 for the soil moisture models), and could not fit these models satisfactorily (i.e., the models resulted in rank-deficient model matrices). As a result, we used separate LMMs for solar radiation and soil moisture.
The best-fitting LMMs for each MDS axis using the environmental predictors all included fixed effects for treatment, the environmental predictor, and an interaction term between treatment and the environmental predictor. The fixed effect for year dropped out during the model selection procedure because it did not improve model fit. A random intercept term for block was included in each model. The random intercept for replicate block allows for residual variation in MDS site scores among experimental units in the LMM; however, it removes such variation when estimating the effects of treatment, the environmental predictors, and their interaction. We were particularly interested in the interaction between treatment (i.e., the effect of canopy trimming) and the environmental predictors (i.e., solar radiation and soil moisture) on understory plant community compositional change in the CTE.
Lastly, we conducted permutational multivariate analysis of variance (PERMANOVA) to evaluate the effects of treatment, time, and environmental variables (i.e., solar radiation and soil moisture) including interactions. This was done to verify the results from the LMMs of MDS axes using the entire community composition dataset. The PERMANOVAs were done using the "adonis" function in the vegan R package (Oksanen et al., 2020). PERMANOVAs were done separately for soil moisture and solar radiation to maintain congruence with the LMMs. Permutations were constrained within experimental blocks (A, B, and C), because of high β-diversity among experimental replicates. A total of 9999 permutations were used.

Treatment effects on the environment over time
Trimming of the canopy in 2004 and 2014 in the experimental plots increased solar radiation and soil moisture in the forest understory (Appendix S1: Figures S1, S2 and Table S1). In 2003 before the first canopy trim, solar radiation levels were generally low, averaging 25.1 and 21.4 W/m 2 in control and hurricane simulation plots, respectively. At this time, the mean difference in solar radiation was 3.7 W/m 2 greater in hurricane simulation than in control plots. Similarly, soil moisture averaged 0.45 and 0.51 cm 3 /cm 3 in the control and hurricane plots, respectively, being on average 0.14 cm 3 /cm 3 greater in hurricane simulation than in control plots (Appendix S1: Figure S1). Then, the mean difference among treatments maximized in 2005 after the first trim, increasing to be 7.8 W/m 2 and 0.25 cm 3 /cm 3 greater in simulated hurricane plots. In 2005, incoming solar radiation and soil moisture measured 36.9 W/m 2 and 0.12 cm 3 /cm 3 for hurricane plots, compared to 29.1 W/m 2 and 0.08 cm 3 /cm 3 for the control plots (Appendix S1: Figure S1).
Differences in solar radiation and soil moisture among treatments attenuated over time between trimming implementations, but then increased again in 2015 following the second trim, where the average differences in solar radiation and soil moisture were 22.9 W/m 2 and 0.07 cm 3 /cm 3 greater in simulated hurricane plots. At this time, average incoming solar radiation levels were 51.0 W/m 2 in hurricane simulation plots and 28.1 W/m 2 in control plots, and average soil moisture measured 0.41 cm 3 /cm 3 in hurricane simulation plots and 0.35 cm 3 /cm 3 in control plots (Appendix S1: Figures S1 and S2). Maximum solar radiation levels in the hurricane simulation plots reached about 200 W/m 2 following canopy trimming (Appendix S1: Figure S2). In 2017, Hurricane Maria induced increases in solar radiation and soil moisture in both treatments, wherein differences diminished among treatments (Appendix S1: Figures S1 and S2). In 2018, 1 year post-Maria, average incoming solar radiation levels were 72.0 W/m 2 in the hurricane simulation plots and 63.9 W/m 2 in control plots, and average soil moisture measured 0.44 cm 3 /cm 3 in both hurricane simulation and control plots (Appendix S1: Figure S2). In the years following Hurricane Maria, minimal differences have emerged among treatments and there has been a slightly decreasing trend in both solar radiation and soil moisture in all CTE plots, because of forest canopy recovery and recruitment dynamics.

Sapling and fern demographic responses
We compare the dynamics of understory saplings and ferns in simulated hurricane plots, which were trimmed in 2004 and 2014, to control plots having no canopy trim. All plots were affected by Hurricane Maria in 2017. In the prolonged absence of simulated or natural hurricane disturbances (i.e., not immediately following canopy trimming or Hurricane Maria), understory plant mortality and recruitment rates in both treatments of the CTE were generally low, fluctuating near zero (Figure 1). Treatment implementation led to a significant difference in sapling mortality and recruitment rates (Table 1). Sapling mortality and recruitment rates also differed significantly over time (Table 1). Baseline mortality of saplings from 2003 to 2016 (i.e., not in intervals following simulated or natural hurricanes) in the control plots was 0.0060 stems/month, whereas it was 0.0115 stems/month in the hurricane simulation plots (Figure 1). Thus, baseline mortality for saplings was approximately twice as high in hurricane simulation plots relative to control plots. Mortality rates of ferns between treatments were nearly equal, with mortality rates averaging 0.0051 and 0.0052 plants/month from 2003 to 2016 in control and hurricane simulation plots, respectively. Fern mortality rates were not statistically different among treatments or over time; however, fern recruitment differences between treatments were marginally significant and fern recruitment differed significantly over time (Table 1). Following hurricane simulations, mortality rates of both saplings and ferns tended to be higher in the hurricane treatment than in the control treatment ( Figure 1). Baseline recruitment rates for saplings averaged 0.0308 stems/month in the control plots and 0.0195 stems/month in the hurricane simulation plots. For ferns, baseline recruitment was 0.0019 plants/month in the control plots and 0.044 plants/month in the hurricane simulation plots (Figure 1).
For both saplings and ferns, demographic rates increased with the experimental canopy trimming and Hurricane Maria (Figure 1). For saplings, recruitment rates approximately doubled in the hurricane simulation plots compared to control plots for census intervals following the canopy trimming treatment (2007 and 2015; Figure 1). For understory ferns, recruitment in the control plots was nonexistent for census intervals following experimental manipulations (2006 and 2015) but measured 0.0037 and 0.0071 plants/month in the hurricane simulation plots for the first and second canopy trimming treatments, respectively. However, understory recruitment related to Hurricane Maria was generally greater in the control plots than it was in the hurricane simulation plots (Figure 1). For example, sapling recruitment rates for the two census intervals following Hurricane Maria averaged 0.0308 plants/month in the control plots compared to 0.0195 plants/month in the hurricane simulation plots. For ferns, the difference in demographic rates after Hurricane Maria was more evident, with recruitment averaging 0.0465 plants/month in the control plots, which was approximately 15 times higher than the recruitment rate of 0.0034 plants/month in the hurricane simulation plots (Figure 1).

Key species driving demographic responses and compositional turnover in the CTE
For saplings, the overall recruitment response to canopy trimming treatments was pronounced. Stem densities increased fourfold in the trimmed hurricane simulation plots (Figure 2a). The increase in stem densities in the hurricane simulation plots (Figure 2a) is primarily driven by the recruitment of pioneer species that responded strongly to the increased solar radiation because of canopy opening after trimming, such as the understory shrub, P. berteroana (Figure 2d), and the pioneer tree species C. schreberiana (Figure 2b). Both species recruited heavily; then, stem densities thinned over time as the canopy of the trimmed plots closed, and some of the pioneer saplings grew to larger, canopy-dwelling trees. In the hurricane simulation plots, the recruitment pulse of P. berteroana following the first canopy trim persisted for about four times as T A B L E 1 Scheirer-Ray-Hare test statistics for differences in recruitment and mortality rate for saplings and ferns in the Canopy Trimming Experiment (CTE) by treatment and time (and their interaction).  Stem densities for both Ca. arborea (Figure 2f) and T. balsamifera (Figure 2g) were similar in all plots before canopy trimming, and both showed an initial pulse in recruitment in the hurricane simulation plots after trimming by increasing four times for Ca. arborea and doubling for T. balsamifera, but these saplings did not suffer such high mortality as the trimmed canopy gradually closed, and their densities have remained elevated since then. In contrast to the pioneer and intermediate lightdemanding tree species mentioned above, many species show almost no change in stem densities in response to the effects of canopy trimming to simulated hurricane disturbance. Two of the most dominant species that are characteristic of this LEF forest type, the sierra palm (P. acuminata var. montana; Figure 2c) and the tabonuco (D. excelsa; Figure 2e), did not show any changes in stem densities in response to either simulated or natural hurricanes. Additionally, H. caribaea (Figure 2h), a heliophilic dominant understory monocot (Musaceae family) measured as part of the tree census, recruited into hurricane simulation plots following canopy trimming, where it had been absent and shows several clear increases over time in response to hurricane effects. In control plots, it increased in abundance following Hurricane Maria (Figure 2h).
Of the 20 fern species observed during the study, 10 were resident at the beginning, while only 6 of these were present at the end (Appendix S1: Table S2). The most common resident fern species were T. deltoidea and the trunkless tree fern Cy. borinquena. Densities of T. deltoidea steadily declined in the control plots throughout the study (Figure 3b). By contrast, T. deltoidea density initially declined in hurricane simulations plots but increased after the first trim in 2004, declining again by the second trim to which it showed little response; however, the magnitude of increase in the trimmed plots after hurricane Maria in 2017 was of about the same as to the first trim (Figure 3b). There were very few individuals of Cy. borinquena in the control plots, and this species showed little response to the two trimming treatments or Hurricane Maria (Figure 3c). All other fern species occurred in lower densities throughout the study (Appendix S1: Table S2) except the pioneer fern species Cy. arborea, a tall tree which recruited briefly into the hurricane plots following the first and second canopy trims, and also appeared in very large numbers in the control plots 1 year after Hurricane Maria (Figure 3d).
The species composition of the fern community also changed over time (Appendix S1: Table S2). Before CTE treatments (i.e., [2003][2004], there were more resident fern species (seven) present in the hurricane simulation plots than in the control plots (five), but by 2018 (after Hurricane Maria), only the core three resident species remained in the hurricane plots. Both the control and the hurricane plots maintained core populations of understory resident ferns (T. deltoidea, Cy. borinquena, and Maria. While a small number of four low-growing pioneer species had appeared by the end of the study, it was the tree fern species Cy. arborea that dominated the fern F I G U R E 4 Metric multidimensional scaling (MDS) ordination for the understory community of the Canopy Trimming Experiment, including tree saplings, ferns, and H. caribaea over time. (a) Ordinations by treatment (control-No Trim + No Debris, and hurricane simulations-Trim + Debris). Note the separation among block and the general trend for plots to generally move up in the ordination space over time in the absence of experimental of natural hurricane disturbance, but down in the ordination space following canopy trimming and Hurricane Maria. Block replicates are labeled (blocks A, B, and C). Points are colored by and connected with arrows (and numerals) between consecutive vegetation censuses. (b) Species scores in the ordination space. Axes are square-root transformed to help visually distinguish small values. Species are colored based on their life form (see Appendix S1: Table S2 for a complete table of taxonomic information corresponding to six-letter species codes). Notable species codes clockwise around the periphery of the two-dimensional species space are as follows: NEPBRO-Nephrolepis brownii, CASARB-Casearia arborea, GUTCAR-Guatteria caribaea, SLOBER-Sloanea berteroana, SCHMOR-Schefflera morototoni, NEPRIV-Nephrolepis rivularis, CYAARB-Cyathea arborea, PSYBER-Psychotria berteroana, CECSCH-Cecropia schreberiana, DACEXC-Dacryodes excelsa, DANGEN-Danaea geniculata, THEDEL-Thelypteris deltoidea, PREMON-Presotea acuminata var. montana, ADIOBL-Adiantum obliquum, MYRLEP-Myrcia amazonica, and MANBID-Manilkara bidentata. See Appendix S1: Figures S3 and S4 for species loadings along each MDS axis plotted separately. community in the control plots after Hurricane Maria. Some of Cy. arborea that recruited into the hurricane plots after trim 1 had disappeared by 2009 but reappeared in very small numbers after the second trim. Thus, by 2018, the control plots had more species (nine) because of the influx of pioneer species than the hurricane plots (six), which had both lost resident species and gained only two other species.

Changes in understory plant community composition in CTE treatments over time
Community changes in species composition visualized using MDS ordination showed that composition varied among the three blocks (i.e., spatially) and over time (i.e., in plots from 2003 to 2018) in response to hurricane disturbance ( Figure 4a). For the hurricane simulation treatment, MDS site scores shifted most dramatically in 2004 and 2015 in response to the canopy trimming. The movement of the hurricane simulation plots was uncoordinated in the MDS ordination space over time. Following the 2004 experimental trimming, blocks A and C shifted positive on MDS axis 2 (i.e., upwards), but block B moved negative; additionally, blocks A and C moved left on MDS axis 1, while block B shifted right (Figure 4a). After the first experimental trim, hurricane treatment plots were more or less consistent in moving down and to the left in the MDS ordination space. For the control plots, all experimental blocks shifted positively on MDS axis 2 over time. Some variation in movement among blocks existed about MDS axis 1 of the ordination space, as block B moved mainly in a positive direction (i.e., right), block A moved negative (i.e., left), and block C remained stationary. For the control plots, MDS site scores shifted most noticeably in response to Hurricane Maria in September 2017 (Figure 4a), showing drastic shifts down on MDS2 and to the right for blocks A and B. Movement of hurricane simulation plots in the ordination space due to Hurricane Maria in 2017 was much less, being similar in magnitude to compositional change due to previous experimental implementations. Furthermore, there were considerable differences in community composition among blocks A, B, and C, as the plots from each block separate in the MDS ordination space.
Responses of pioneer species, like C. schreberiana and S. morototoni, and the understory shrub (P. berteroana), which are highly responsive to hurricane disturbance, influence the positive side of axis 1 and the negative side of axis 2 in the MDS species space (Appendix S1: Figures S3 and S4). Seven (three ferns, three trees, and a shrub) of the 11 species that influence the positive side of MDS1 also influence the negative side of MDS2 (C. schreberiana, Cy. arborea, Cy. borinquena, N. rivularis, P. berteroana, S. berteroana, and S. morototoni). Three (one tree and two ferns) of 11 species that influence the negative side of MDS2 also influence the negative side of MDS1 (D. excelsa, D. geniculata, and T. deltoidea). Species exerting more influence in the MDS species space tended to be those with higher densities in the CTE plots, those that show a clear preference for the treatments, or those that respond strongly over time (Appendix S1: Table S2).
Thus, the MDS species space shows a forest understory species community composition trade-off, primarily capturing sapling and fern pioneer species responses, but also certain responses of late-and secondary-successional species (e.g., M. bidentata, Ca. arborea, G. caribaea, N. brownii, and A. obliquum). The secondary-successional and late-successional species tend to be located on the top and right of the species MDS ordination space (Figure 4b), whereas the pioneer taxa are found lower and to the left (Appendix S1: Figures S3 and S4).

Solar radiation and soil moisture effects and their interaction with the CTE treatments
For the LMMs that included effects for solar radiation, the fixed effects for treatment were the strongest significant effects in the LMMs (see Appendix S1: Figure S5 and Table S3). The hurricane simulations resulted in a positive shift for MDS site scores on both axis 1 and axis 2 ( Figure 5). The effect of variation in solar radiation among plots alone, after accounting for the main treatment effect, was minimal (variable not included for MDS1 and β = À0.01, p = 0.19 for MDS2). The interaction between solar radiation and the hurricane simulation treatment was statistically significant in affecting MDS2 axis scores ( Figure 5; Appendix S1: Table S3). Thus, the majority of variation in understory composition and its variability over time can be attributed to the simulated hurricane treatments of the CTE (i.e., the categorical fixed effect of treatment, β = 0.37, p < 0.01, for MDS1 and β = 104.91, p < 0.05, for MDS2). The interaction of solar radiation with experimental treatment was significant for compositional variation along MDS2 (β = 0.02, p = 0.02), but not for compositional variation along MDS1, where the interaction was not included in the best-fitting LMM ( Figure 5).
Similar to the LMMs for soil radiation, the fixed effects for experimental treatment were among the strongest significant effects in the LMMs for soil moisture (MDS1: β = À2.26, p < 0.05; MDS2: β = À2.92, p < 0.01; see Appendix S1: Figure S6 and Table S4). The effects of soil moisture alone on compositional variation along MDS axes were marginally significant (MDS1: β = À2.08, p = 0.18; MDS2: β = À2.04, p = 0.16). Yet, the interaction between soil moisture and treatment was significant for MDS variation along both axes (MDS1: β = 6.14, p = 0.02; MDS2: β = 5.60, p = 0.03), wherein understory composition shifted toward the negative side of both MDS axes in control plots with increasing soil moisture but had the opposite trend in experimentally trimmed F I G U R E 5 Interaction plots for linear mixed-effects models. Points show plot mean solar radiation or soil volumetric water content (VWC) values and multidimensional scaling (MDS) sites scores for the corresponding census interval (n = 42 for soil moisture and n = 60 for solar radiation). Models were fit to explain MDS loadings for Canopy Trimming Experiment (CTE) hurricane simulation and control plots over time using annual mean solar radiation (upper panels) or soil VWC (lower panels) Â Treatment as the explanatory variables with a random intercept for experimental block. Probabilities in the lower right of the figure panels show the statistical significance of the interaction between treatment and environment (either mean soil VWC or solar radiation) on the MDS site scores. See Appendix S1: Tables S6 and S7 for complete model summary tables and effect plots. plots ( Figure 5). Thus, increases in soil moisture because of hurricane-driven canopy opening led to divergent trends in the species composition of the regenerating understory. Regarding compositional variation among replicated blocks, along MDS1 block B was different from blocks A and C, and along MDS2 block A was distinct from blocks B and C (Appendix S1: Figure S7).
The PERMANOVA results largely confirmed the results from the LMMs (see Appendix S1: Table S5 for solar radiation and Appendix S1: Table S6 for soil moisture). Treatment was the most significant factor in driving compositional dissimilarity (F 1,59 = 8.13, p = 0.001 for solar radiation, and F 1,41 = 5.39, p = 0.001 for soil moisture). Additional effects of solar radiation alone, time, or interactions between treatment and time and solar radiation and time were all marginally significant (Appendix S1: Table S5). Soil moisture alone had a nonsignificant effect; however, it interacted significantly with treatment to affect compositional dissimilarity (F 1,41 = 3.18, p = 0.003). Effects for time and its interaction with treatment were also significant factors affecting community composition in the soil moisture PER-MANOVA (F 6,41 = 1.58, p = 0.003 and F 6,41 = 1.04, p = 0.006, respectively; Appendix S1: Table S6).

DISCUSSION
The simulated hurricane and natural hurricanes affecting the CTE plots have created a novel testing ground for the compounding effects of multiple hurricane disturbances on tropical forest plant communities. Using a 17-year dataset (2003-2019) on understory forest dynamics of saplings and ferns, our first research question addressed how demographic rates of understory saplings and ferns vary in response to hurricane disturbance. We hypothesized that after repeated disturbance from experimental canopy trimmings followed by a naturally occurring hurricane, the understory would show less variation in mortality and recruitment rates and less change in species composition than in the untrimmed control plots. This hypothesis rests on the assumption that the species composition of trimmed plots would maintain a higher proportion of pioneer species than in control plots, and the ability of pioneer species to recruit to natural hurricane disturbance would be lower in the previously trimmed plots than in the controls (Hogan et al., 2016;Uriarte et al., 2009). We found some support for this hypothesis, in that mortality and recruitment rates differed among treatments and over time for saplings, and recruitment rates of ferns differed over time (Table 1). Demographic rates were more evenly spread over the 17 years of monitoring for the trimmed, hurricane simulation plots than for control plots. Additionally, trimmed plots had a muted demographic (i.e., mortality or recruitment) response to the natural hurricanes in 2017 compared to the untrimmed control plots (Figures 1-3).
Our second research question asked to what degree differences in solar radiation and soil moisture that result from canopy loss drive compositional change in the forest understory. We hypothesized that solar radiation would be the primary driver of understory plant community compositional change, as it may change more dramatically in the understory after a hurricane than soil moisture, and because soil moisture is not believed to be limiting in this wet tropical forest. Contrary to the hypothesis, we found that although solar radiation had significant effects, soil moisture interacted with hurricane disturbance to also affect the composition of recruits. Both soil moisture and solar radiation to the understory increased in the hurricane simulation treatments relative to control plots (Appendix S1: Figures S1, S2 and Table S1). Although measured solar radiation levels alone had no additional effect on altering community composition after the direct effects of canopy opening (i.e., treatment) were accounted for, the changes in solar radiation following canopy opening significantly shifted species composition of the regenerating understory toward pioneer saplings (e.g., P. berteroana, C. schreberiana, and S. morototoni) and ferns (e.g., O. aculeata and Cy. arborea). Soil moisture, however, interacted with the hurricane simulation treatment to move the understory community composition toward a more pioneer species composition ( Figure 5), pointing to a secondary control on understory species recruitment, which may be independent of solar radiation responses in this historically and evolutionary-adapted hurricane-disturbed tropical wet forest.

Demographic responses of the understory to repeated hurricane disturbances
For saplings, the recruitment of understory shrub species (e.g., P. berteroana) and larger pioneer tree species (e.g., C. schreberiana) is an iconic successional signature of hurricane disturbances in Caribbean tropical forests (Bellingham et al., 1995;Brokaw et al., 2012;Hogan et al., 2016;Zimmerman et al., 2021). The recruitment response of trees in the CTE understory to canopy openings from simulated and natural hurricanes damages more than tripled stem densities (Figure 2) (Chevalier et al., 2022;Zimmerman et al., 2014). This rapid colonization by pioneer species results in a community shift toward species with low wood densities and faster lifehistory strategies in the short term Hogan et al., 2018). These pioneer species grow fast and rapidly replace and infill the forest canopy. The closure of the forest canopy recreates the shaded environment that late-successional species (e.g., M. bidentata and D. excelsa) with higher wood density require to recruit and grow , which occurred within 4 years following the first canopy trim . Mid-successional species such as Ca. arborea show a curious response to the effects of simulated hurricane disturbance (Figure 2) in stem densities, potentially benefiting long-term under cyclic hurricane recurrence. Casearia arborea and T. balsamifera recruited to higher stem densities after the first canopy trim, and these high densities remained elevated for the duration of the CTE monitoring ( Figure 2).
Simulations of the forest dynamics in the LEF using the SORTIE model parameterized using demographic data from the nearby Luquillo Forest Dynamics Plot  showed that these mid-successional species are preferentially selected for under cyclic, repeated hurricane disturbance versus infrequent, noncyclic disturbance. Thus, these findings from the CTE support previous simulations of forest demographics, at least over the short term (i.e., one to two decades, given several canopy opening events). Such preferential shifts in the recruitment of the sapling community can change the successional trajectories of the tropical tree community (Chazdon, 2008;Hogan et al., 2016;Lai et al., 2021), which in turn may have lasting effects on the forest structure and carbon dynamics of the forest (Chevalier et al., 2022;X. Feng, Uriarte, et al., 2018;Uriarte et al., 2009), especially if hurricane disturbance is coupled with other environmental stressors (e.g., drought). The ecosystem demography model predicts a decline in the carbon storage capacity of LEF over the next several decades, which is driven by an increase in early and midsuccessional tree species at the expense of latesuccessional tree species (X. . These declines are predicted to be accelerated by a 30% increase in drought frequency and continued warming (of 1.2 C from 2016 mean temperatures) (X. . Downscaled climate models do predict a decrease in total annual rainfall, which is driven primarily by an increase in the frequency of dry days (those with <5 mm of rainfall) during the onset of the wet season (May-June) (Ramseyer et al., 2019). Therefore, drought will likely interact with hurricane disturbances to affect understory species composition via preferential recruitment and mortality dynamics in the LEF in novel ways into the coming decades (Henareh Khalyani et al., 2019).
High fern mortality rates after the first canopy trim in 2005 and the even higher rates after the second trim in 2014 ( Figure 1) could have been due to the unavoidable damage to ground-level plants from arborist trimming activity (Halleck et al., 2004) but could also have been due to the rapid increase in grass cover . Fern mortality was highest immediately following debris deposition of Hurricane Maria (Figure 1) just as adult (i.e., stems >10 cm dbh) tree mortality doubled after Hurricane Maria , suggesting that wind, a factor not included in the hurricane simulation experiment, may be an indirect driver of fern mortality rates (Royo et al., 2011).
Thelypteris deltoidea, Cy. borinquena, and Cy. arborea are the three most abundant fern species observed in our study and showed different responses to hurricane disturbance ( Figure 2). The resident understory fern species (T. deltoidea and Cy. borinquena) maintained relatively high densities throughout the study (compared to lesscommon fern species) and have been recorded as common in other studies in the LEF where fern observations have been included and Cy. arborea, though never common, appears in these studies as well (Halleck et al., 2004;Portugal Loayza, 2005;Royo et al., 2011;Smith, 1970). Unlike most of the other fern species present in the CTE plots, these three common species have very limited ranges. The distribution of T. deltoidea is limited to Puerto Rico and a few islands in the Lesser Antilles. Cyathea borinquena is endemic to Puerto Rico, and Cy. arborea is limited to the Caribbean islands and northern Colombia. In most fern hotspots (Suissa et al., 2021), high levels of habitat variation (e.g., topography and elevation gradients) have provided a wide variety of niches for potential fern speciation. However, it is possible that for these three common fern species in our study that hurricane-induced habitat heterogeneity has instead led to various fern adaptations that have allowed their persistence and understory dominance in wet forests with hurricane disturbance regimes of the Greater Antilles. Indeed, this has been reported for many Caribbean tree species (Bellingham et al., 1995;Griffith et al., 2008;Hogan et al., 2018;Tanner et al., 1991;Vandermeer et al., 2000).
It appears that resident understory fern species have limited opportunities for recruitment from spores under a closed canopy. However, Sharpe and Shiels (2014) reported that spores are produced in great numbers in response to the increase in solar radiation from canopy openings. The relatively high densities of the resident fern species, T. deltoidea and Cy. borinquena, at the start of our study, could be the result of spore release after Hurricane Georges (1998), which occurred only 5 years before the start of the CTE. Although the canopy damage was not as severe during Hurricane Georges as in Hurricanes Hugo (1989) and Maria (2017) (Hogan et al., 2016;Uriarte et al., 2019), it likely still created spore-driven dispersal and recruitment for these fern species. The steady decline of both T. deltoidea and Cy. borinquena in the control plots until Hurricane Maria arrived in 2017 resulted in fern densities returning to only half of the initial numbers recorded in the CTE. An undisturbed canopy for almost 20 years would result in ever lower light levels and a deeper shade habitat, leading to lower numbers but more reproductively mature ferns. Relatively little soil surface leaf litter cover in an undisturbed forest would also facilitate spore germination (although treatments that created initial differences in leaf litter cover and litterfall rates subsided over time; Lodge et al., 2014;Silver et al., 2014). If sufficient spore production and dispersal are followed by successful germination and sporeling establishment after Hurricane Maria, the populations of resident ferns will likely again increase throughout the LEF.
In contrast to these resident fern recruitment limitations, the tree fern Cy. arborea is a successful pioneer that mostly inhabits open areas such as roadsides and landslides throughout the LEF (Walker et al., 2013). The role of the tree fern Cy. arborea as a pioneer species in colonizing disturbed areas has been well documented; it has prolific and rapid regeneration on disturbed hillslopes following landslides (Walker et al., 2010) and also a preference for colonizing tipped-up mounds pits following hurricane disturbance (Bellingham et al., 1995;Walker, 2000). The height of the tree fern facilitates wind dispersal of huge numbers of spores (Conant, 1976), and although almost all leaves of Cy. arborea individuals are lost in a natural hurricane, they are rapidly replaced at a rate of one frond per month (Conant, 1976). These leaves also rapidly mature, release spores and the fronds fall all within about 6 months (Walker et al., 2013), which facilitates rapid colonization after a disturbance. Tree ferns do not thrive in the typical low light conditions of the understory when the canopy closes, however, and may, therefore, only facilitate tree seedling and other fern sporeling establishment in open areas or tree-fall gaps (Walker et al., 2010). Spores of the smaller species known to colonize landslides such as Pityrogramma calomelanos var. calomelanos also appeared and were able to germinate, albeit in small numbers after Hurricane Maria, exemplifying the positive effect of natural hurricane winds on fern recruitment and new species establishment.
Simulated versus natural hurricane responses: The effects of solar radiation and soil moisture on understory dynamics in the LEF Hurricane Maria in 2017 provided more of an opportunity for the recruitment of pioneer species than did either of the experimental trim events. The pioneer species Cy. arborea appeared in small numbers in the hurricane plots following the second trim, but ferns over 10-cm height recruited more in the control plots in the 2 years after Hurricane Maria (Figure 2b). A similar trend was observed for C. schreberiana (Figure 2b) and several lesscommon pioneer species, such as Simarouba amara, which did not respond strongly to the experimental trimming but increased copiously following the natural hurricane disturbance of 2017. The experimental trims opened a relatively small area of the forest canopy, so only small areas of the forest understory were subjected to changes in microclimate (i.e., light and temperature), soil moisture, and any interaction between the two. Thus, there were likely relatively few available microsites for the spores of the common tree ferns to colonize, for seedlings from the seedling bank to flourish, or for seeds in the soil seed bank to germinate and establish. However, when propelled by natural hurricane-force winds, many more spores from further dispersal distances could increase recruitment of ferns into the small CTE plots after Hurricane Maria. Similarly, bigger changes in the understory microclimate driven by larger areas of canopy opening, with more diffuse lateral light entry into the forest understory, could stimulate the germination of seeds and establishment and growth of tree seedlings (Comita, Uriarte, Thompson, et al., 2009;Comita, Uriarte, Forero-Montaña, et al., 2018).
A big question is why recruitment of the tree fern-Cy. arborea (Figure 3b), and the pioneer trees such as C. schreberiana (Figure 2b), in 2018 after Hurricane Maria was more common in the CTE control plots than in previously trimmed hurricane simulation plots. The control plots are located in the same three blocks of forest as the trimmed, hurricane simulation plots Shiels & Gonz alez, 2014), so we were surprised at these differences. A dense grass cover appears in hurricane simulation plots after trimming Zimmerman, 2020a) perhaps creating a barrier to leaf expansion in newly germinated spores of Cy. arborea. Or perhaps the comparatively low amount of surface litter on the forest floor of the control was a wetter and more shaded substrate that was more suitable for spore germination, than the previously disturbed, and likely drier, substrate of the hurricane plots. There may also have been a greater accumulation of spores in the control plot soils that had not germinated as they did after the two canopy trims in the simulated hurricane plots. This could mean that there was a larger soil spore bank available in the control plots ready to germinate when the natural hurricane occured.
Similar processes as described above for ferns have likely affected sapling demographics in the control relative to canopy trimmed plots (Luke et al., 2016;Royo et al., 2011). In 2017 at the time of Hurricane Maria, there could have been more seeds in the soil seed bank and seedlings in the seedling bank of the control plots, which had not already grown to sapling stage than in trimmed, hurricane simulation plots. After Hurricane Maria, these suppressed seedlings in the control plots would be able to grow rapidly. Similarly, accumulated seeds of pioneer species would germinate and grow quickly to reach sapling size, being recorded as recruits in the post-Maria tree census. In all of the plots-the hurricane simulation plots with two canopy trims and a hurricane, and the control plots with only one hurricane impact-sapling densities will decline as the canopy closes and reduces light, causing shade-induced mortality of pioneer species, and density-dependent mortality for abundant shade-tolerant species (Comita, Uriarte, Thompson, et al., 2009;Comita, Uriarte, Forero-Montaña, et al., 2018;Luke et al., 2016;Muscarella et al., 2013;Royo et al., 2011;Uriarte et al., 2018). Compositional shifts toward more pioneer and secondary-successional species in the CTE may be partly due to increased throughfall and wetter soils in the trimmed plots, with soil VWC being >0.1 cm 3 water/cm 3 soil wetter in trimmed plots (Appendix S1: Figure S2; Van Beusekom et al., 2020). Soil moisture also has effects on seedling survival in the LEF, with seedlings growing in dry sites having slightly increased survival with increasing soil moisture (i.e., in wetter years), and vice versa. For example, dry years preferentially benefit the survival and potentially the recruitment of liana seedlings in the LEF, especially in wetter areas (Umaña et al., 2019); however, the predominant trend is for high soil moisture to decrease seedling survival, especially in the wettest microsites, especially for late-successional shade-tolerant species .
Data from our analyses of compositional changes in the CTE plots show that higher soil moisture levels interacted with the experimental canopy trimming light increases to shift community composition toward pioneer-type tree species associated with secondary succession and past land-use legacies in the LEF (see the lower left quadrant of MDS space, Figure 4b, or negative sides of MDS axes in Appendix S1: Figures S3 and S4; Hogan et al., 2016;Uriarte et al., 2009), and away from old-growth, more shade-resistant species (see upper right quadrant of the MDS ordination space, Figure 4b, or positive sides of NDMS axes in Appendix S1: Figures S3 and  S4). Our results are similar to those of Arasa-Gisbert et al. (2021), who reported canopy openings caused by tree removal in a Mexican tropical forest caused a shift toward pioneer and secondary-successional generalist species because of the lower recruitment of old-growth, late-successional species. Thus, in addition to light availability, soil moisture content is a distinct, yet critical factor that governs the composition of the regenerating understory plant community following hurricane disturbance. Drier conditions tend to support increased relative abundances of many late-successional species in the LEF (e.g., M. bidentata and Sloanea berteroana) (Umaña et al., 2019;Uriarte et al., 2018), therefore hurricanedriven increases in throughfall likely also facilitate pioneer species recruitment after canopy disturbance.
In conclusion, the demography of the understory plant community-tree saplings, ferns, and the grounddwelling H. caribaea-responds dynamically to repeated hurricane disturbances, indicating that an increased frequency of hurricanes will reduce understory recruitment after each disturbance event, if there is not sufficient time to allow for stand thinning and development, or regeneration of the seed and seedling banks in between hurricanes Brokaw & Grear, 1991). This reduction in recruitment to repeated hurricanes may be greatest for late-successional species (e.g., M. bidentata and S. berteroana), when increased solar radiation and soil moisture in the understory create environmental conditions that favor the recruitment of pioneer species at the expense of late-successional species. Hurricaneinduced increases in canopy openness and solar radiation and soil moisture both drive community compositional changes in the understory plant community. Increased solar radiation is the main driver, however, increased throughfall and lower transpiration from canopy trimmed trees increase soil moisture availilibity to shape the species composition of the regenerating understory. A significant future drying trend is projected for Caribbean forests, including wet forests at Luquillo (Ramseyer et al., 2019), so the relative strength of climate drying versus increased throughfall and lower transpiration because of hurricane disturbance will infleunce shortterm understory responses to hurricanes (e.g., Maria) and longer-term forest dynamics.
AUTHOR CONTRIBUTIONS Aaron B. Shiels, Grizelle Gonz alez, and Jess K. Zimmerman designed the CTE. Grizelle Gonz alez, Joanne M. Sharpe, Sarah Stankavich, Samuel Matta Carmona, John Bithorn, and Jamarys Torres-Díaz collected the data in the field. Samuel Matta Carmona, John Bithorn, and Sarah Stankavich entered, and quality checked the data. J. Aaron Hogan, Joanne M. Sharpe, and Ashley Van Beusekom analyzed the data for this manuscript. Aaron B. Shiels, Grizelle Gonz alez, and Jess K. Zimmerman secured and provided funding and provided project oversight/management. J. Aaron Hogan and Joanne M. Sharpe wrote the manuscript. All authors read the manuscript, provided comments, and contributed intellectually to this work.

ACKNOWLEDGMENTS
We recognize and thank the many volunteers and LTER technicians who helped make the CTE possible because of their hard work in the field. We thank Jill Thompson and two anonymous reviewers for comments, which have improved this work. We acknowledge the USDA Forest Service, especially the International Institute of Tropical Forestry, for their dedicated support of decades of research in the Luquillo Mountains, including the CTE. Carlos Estrada and Carlos Torrens provided invaluable assistance with sensor monitoring in the field. The Luquillo LTER is supported by the following research grants from NSF to the former Institute for Tropical Ecosystem Studies, within the Department of Environmental Sciences at the University of Puerto Rico: BSR-8811902, DEB 9411973, DEB 0080538, DEB 0218039, DEB 0620910, DEB 0963447, and DEB-1546686. Additional support was provided by the University of Puerto Rico.

CONFLICT OF INTEREST
The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT
All data used in this research are previously published, publicly available, and cited properly herein as follows: The CTE fern data are available from the Environmental Data Initiative (EDI) (Sharpe, 2021): https://doi.org/ 10.6073/pasta/0ea1a4b97112cc66e5566129cf88ec1c. The CTE plant and seedling measurement dataset (which was used for fern recruitment and mortality rates) is available from EDI (Zimmerman, 2020a): https://doi.org/10.6073/ pasta/e4eb7f6d0e8287a36c2ec12be3f781ea. The CTE tree dataset (from which sapling abundances and sizes are derived) is available from EDI (Zimmerman, 2020b)