This study was conducted in the Luquillo Experimental Forest (LEF; 18°18′N, 65°50′W) in northeastern Puerto Rico. Categorized as a subtropical wet forest according to the Holdridge Life Zone System (Holdridge, 1967), this site holds the TRACE experimental plots, which were established in 2015 in a secondary forest that has regenerated naturally after agricultural development during the first half of the 20th century (Kimball et al., 2018). Located at 100 m elevation, the site contains clay-rich and acidic soils classified as Ultisols (Scatena, 1989), with >60% of root biomass concentrated in the top 10 cm soil layer (Silver and Vogt, 1993; Yaffar and Norby, 2020). In the study area, the dominant canopy trees are Syzygium jambos, Ocotea leucoxylon, and Casearia arborea, and the site is also dominated by the palm Prestoea montana.

Mean annual temperature of the site is 24°C with relatively low seasonal variation, and mean annual rainfall is 3,500 mm (García-Martinó et al., 1996). During the 2015 drought, total precipitation measured at a neighboring site (El Verde Field Station, ∼15 km away) was 48% less than the mean in the previous decade (O’Connell et al., 2018). On September 06, 2017, Hurricane Irma brought >300 mm of rain to the northeast of the island, and just 2 weeks later on September 20, 2017 Hurricane María brought ∼1,500 mm of rain in only 48 h (Hall et al., 2020).

The understory plant community was evaluated inside each of the six 12 m² hexagonal TRACE plots (i.e., three heated and three control), which range in slope from 15 to 26° and are separated by at least 10 m. All plots contain naturally occurring understory vegetation but were installed in areas with no adult trees to allow space for experimental equipment. To evaluate the effects of warming on understory plants and soil, three of the plots were warmed with infrared (IR) heaters +4°C above ambient temperatures of the three control plots. Temperatures at the heated plots were constantly controlled to maintain a 4 ± 0.1°C increase in hourly average temperatures compared to control plots, as sensed by IR thermometers (Kimball et al., 2018). The design of the warming system was adapted from Kimball et al. (2008), with six IR heaters (Model Raymax 1010, Watlow Electric Manufacturing Co., St. Louis, MO, United States) per warmed plot installed on crossbars at approximately 3.6 m from the ground and aimed downward to maintain uniform and constant temperatures for long periods of time. Control plots contained the same infrastructure (i.e., concrete footings, posts, and crossbars) but with sheet metal plates the same size as the heaters to control for potential effects of shading and equipment installation. More details on the layout of the experimental site, infrastructure of the plots, and success of the warming method can be found in Kimball et al. (2018). Experimental warming began in September 2016, was paused in September 2017 due to the passage of Hurricanes Irma and María, and was resumed in September 2018 after 1 year of monitoring hurricane disturbance effects.

All woody seedlings (≥10 cm and <20 cm height) and saplings (≥20 cm height) inside the TRACE experimental plots were initially tagged and identified to species in June 2015 (Census 0; C0), during the peak of the drought. Previously tagged seedlings were reassessed for drought damage in November 2015 (C1), after the rains resumed. The following censuses were conducted in June 2016 (C2; control year, 10 months after the drought ended), June 2017 (C3; 8 months after the warming started), May 2018 (C4; 7 months after the hurricanes, in the absence of warming), and April 2019 (C5; 7 months after warming resumed and 14 months after the hurricanes). Seedlings that surpassed 20 cm height in a census year were reclassified as saplings. Additional details about the seedling surveys can be found in Figure 1A and Bachelot et al. (2020). Raw data from these seedling and sapling surveys is available on the ESS-DIVE data archive at https://data.ess-dive.lbl.gov/datasets/doi:10.15485/1873618.

Functional trait data for the woody tree and shrub species found in our plots were obtained from the Luquillo Forest Dynamics Plot (LFDP) database, which were measured with standard procedures (Cornelissen et al., 2003; Swenson et al., 2012) within the Luquillo Experimental Forest (∼15 km away from our study site; Swenson et al., 2019). We selected functional traits that capture life history differentiation among tree species (Kraft et al., 2010; Wright et al., 2010), which included: wood density (WD; g/cm³), specific leaf area (SLA; cm^–2/g¹), and leaf succulence (LS; g H20/cm² leaf area), as well as leaf phosphorous (P; percentage P of oven-dry mass), leaf nitrogen (N; percentage N of oven-dry mass), and leaf carbon (C; percentage C of oven-dry mass) concentrations. Wood density was measured on 10 samples per species using an increment borer in trees of 10–20 cm diameter at 1 m height, or using branch material for smaller trees (Swenson et al., 2012). Leaf traits were measured on 25 samples per species using sun-exposed leaves or on leaves from the crown of large trees (Swenson et al., 2012).

For each treatment (i.e., disturbance type), we calculated the community-weighted means (CWM) of each functional trait by averaging the trait value across all species present in each plot weighted by the number of stems in that same plot. Measuring the CWM trait value is one of the most common approaches for quantifying functional composition at the community level and analyzing trait-environment relationships (Funk et al., 2017). These values were calculated separately for seedlings and saplings, to evaluate changes in CWMs observed along ontogeny.

To assess how plant community responses to climate warming and disturbance differed across ontogeny, we conducted all the analyses described below separately for seedlings and saplings. First, we investigated how climatic events influenced broad biodiversity patterns by running linear mixed models of species richness (i.e., number of species) and diversity (i.e., Shannon’s Diversity Index) as a function of treatment (i.e., control, drought, hurricanes, and warming). Since each of the six plots was censused six times from 2015 to 2019, we used plot identity as a random effect to account for repeated measurements. All linear mixed models were run in R (R Core Team, 2020) with the lmerTest package (Kuznetsova et al., 2017), using the Satterthwaite method for calculating degrees of freedom and corresponding p-values from t-tests.

We then investigated how each event was associated with changes in species composition. As heterogeneity in dispersions across groups can alter the results of composition analyses (Anderson and Walsh, 2013), we tested for homogeneity of dispersions in species composition and found no significant differences in dispersions across censuses (seedling: F_5,30 = 0.54 and P = 0.77, sapling: F_5,30 = 0.11 and P = 0.97) or across treatments (seedling: F_3,32 = 0.62 and P = 0.57, sapling: F_3,32 = 0.08 and P = 0.94). To detect significant changes in species composition at the beginning and the end of the study period, we used permutational multivariate analysis of variance (PERMANOVA) using data from the first (C0) and last (C5) censuses. Then, to evaluate how variation in plant community composition was explained by time (six census dates) and treatment (i.e., control, drought, hurricanes, and warming), we used distance-based redundancy analysis (RDA; Rao, 1964). We also used PERMANOVAs to investigate changes in species composition following each individual treatment. This enabled us to assess specifically how species composition changed following drought (by comparing composition in all plots in C0 and C2), warming (by comparing composition in all plots in C2 and C3), hurricanes (by comparing composition in all plots in C2 and C4), and warming before and after hurricanes (by comparing composition in the heated plots in C3 and C5). All analyses were run in R (R Core Team, 2020) with the vegan package (Oksanen et al., 2020), using 999 permutations and based on standardized Euclidean dissimilarity matrices of the raw abundances (count of plant individuals) aggregated per plot in a given census (Clarke and Warwick, 2001). We visualized all the results using non-metric multidimensional scaling (NMDS; Kruskal, 1964).

To evaluate how individual species responded to treatments and to predict the consequences for the functional composition of the forest, we used generalized joint attribute modeling (GJAM; Clark et al., 2017) in order to implement a Bayesian predictive trait model (Clark, 2016). This approach models species weights as composition data, using multivariate normal distribution, before using them to predict traits. Specifically, we modeled species weights (w_i, a vector of relative species abundances in plot i) of all species found inside the plots as a function of treatment and plot as random effect (contained in the predictor vector x_i) using a multivariate normal distribution:

where β is the matrix of coefficients and ∑ the species covariance matrix. Once this first step is achieved, a variable change allows the fitted model to be transformed to traits from which trait predictions can be made. Specific details about the model implementation can be found in Clark et al. (2017). This method allows evaluating changes in species abundances not only in response to the variables of interest (with treatment and plot as random effect to account for repeated measurements) but also to the other species. Since species influence each other, using a multivariate normal distribution is more appropriate than using a linear regression on individual species. Furthermore, this method has been shown to perform better than directly modeling traits (Clark, 2016). To fit the model, we used three chains with 2,000 iterations and 500 burn-ins. Model fit was assessed with predictive checks of both species counts and CWMs of functional traits. All analyses were conducted in R (R Core Team, 2020), using the GJAM package for predictive trait modeling (Clark et al., 2017).

To assess changes in broad biodiversity patterns, we evaluated variations in species richness and diversity throughout the study period using linear mixed models. As expected, responses to drought, hurricanes, and experimental warming depended on plant ontogenetic stage (Figure 2). Hurricane disturbance was marginally associated with a decrease in seedling richness (t₂₇ = −2.01, P = 0.055), while seedling diversity showed no significant response to any disturbance type. In contrast, sapling richness (t₂₈ = 3.40, P = 0.002) and diversity (t₂₈ = 3.40, P < 0.001) were positively associated with experimental warming, and sapling diversity was marginally associated with drought (t₂₇ = 1.87, P = 0.072). Full results of linear mixed models can be found in Supplementary Table 1.

As expected, seedling and sapling species composition significantly changed between the first and last census (Table 1, History model). However, the effects of disturbance type and time on species composition were not as pronounced (Figure 3). Redundancy analyses showed that 27% of the variation in seedling composition was explained by disturbance and time (RDA, pseudo-F_6,29 = 1.82, P = 0.027, Table 1). Specifically, disturbance type, although marginally significant, explained 11% of the variation in seedling composition (pseudo-F_3,29 = 1.51, P = 0.103) and time captured 16% of the variation in seedling composition (pseudo-F_3,29 = 2.13, P = 0.023). Meanwhile, sapling composition did not show a strong response to disturbance and time (RDA, pseudo-F_6,29 = 1.82, P = 0.188, Table 1). When evaluating the individual effects of each disturbance type on species composition, we detected no significant effects of drought, hurricanes, or warming on either seedling or sapling composition (Table 1).

Several species showed significant changes in abundance associated with disturbance at both the seedling and the sapling ontogenetic stages (Figure 4 and Supplementary Table 2). In general, drought and warming were more often associated with higher abundance in both seedlings and saplings, while the effects of hurricanes were more varied (Figure 4). Across the three disturbance types, warming was significantly associated with about twice as many species as drought and hurricanes. Warming resulted in high abundance of six species in the seedling stage and seven species in the sapling stage, and with low abundance of two species in the seedling stage and three species in the sapling stage. In contrast, drought was significantly associated with high abundance of four species and low abundance of one species in the seedling stage and with high abundance of two species and low abundance of one species at the sapling stage. Finally, hurricanes were associated with high abundance of two species in the seedling stage and three species in the sapling stage, and with low abundance of three species in the seedling stage and one species in the sapling stage. These species covered a wide range of life forms and successional status (Table 2).

Some of these species demonstrated a significant response to more than one disturbance type. For example, seedlings of one of the most common species in the study area, the early successional shrub Psychotria brachiata (PSYBRA), showed significantly lower abundance during drought but higher abundance under warming conditions. Meanwhile, seedlings of the early successional tree Tabebuia heterophylla (TABHET) showed increased abundance during drought but lower abundance after hurricane disturbance. Observed responses of plant species to the different disturbances also showed variation across ontogenetic stages. In particular, the previously mentioned P. brachiata shrub had higher abundance at the seedling stage but lower abundance at the sapling stage under warming conditions. Lastly, a few other species showed a consistent response to the disturbances across ontogenetic stages, such as the early successional shrub Piper glabrescens (PIPGLA), which showed higher abundance during drought and lower abundance during warming at both the seedling and sapling stages (Figure 4).

We found no disturbance effect on any of the six sapling CWM functional traits for predicted changes in the functional structure of the forest (Figure 5 and Supplementary Table 3). Specifically, drought had no effect on the functional structure of either seedling or sapling communities in our study site. In contrast, hurricane and experimental warming were significantly associated with changes in CWM functional traits of seedling communities (Figure 5). Both the hurricanes and the experimental warming led to an increase in CWM leaf P concentrations and leaf succulence, and a decrease in CWM leaf C concentrations in seedlings. The seedling CWMs for wood density, specific leaf area and leaf N concentrations showed no association to any of the disturbance types included in our study.

Despite the observed shift in species composition between 2015 and 2019, we found few significant changes in seedling and sapling communities specifically associated with drought, hurricanes, and experimental warming. This result suggests that the woody understory might be more resistant to climate warming and disturbance than expected. As mentioned in previous studies, plant communities in this forest may have already developed pre-existing tolerance to a wide range of stressors over a long history of repeated disturbance events, increasing their resistance and resilience through time (Ostertag et al., 2005; Zimmerman et al., 2020, 2021). Alternatively, the lack of marked shifts in seedling and sapling communities associated with these disturbances might suggest lagged responses that take longer to emerge (Xie et al., 2020).

Tropical seedlings have been shown to be sensitive to water availability (Bunker and Carson, 2005; Slot and Poorter, 2007; Álvarez-Yépiz et al., 2018; Uriarte et al., 2018). However, the effects of drought on seedlings have mostly been assessed with experimental manipulations, which may simulate more extreme conditions than those observed in nature and thus exacerbate seedling responses to water stress. In addition, controlled experiments might obscure the balancing effects of other natural occurring mechanisms, such as indirect effects on soil pathogens and mycorrhizal associations that influence seedling survival (Liang et al., 2015; Bachelot et al., 2020). Observational studies provide more realistic data for understanding these cross-trophic dynamics, but extreme climatic events are unpredictable and often with long return intervals (Scatena and Larsen, 1991; Boose et al., 2004). Further, the results may be confounded by other biotic or abiotic components of the ecosystem (Soong et al., 2020).

Our results show a slight increase in sapling diversity in response to drought. In the TRACE site, Bachelot et al. (2020) found higher seedling growth and survival during the 2015 drought when compared to other censuses, which could in turn influence sapling diversity. However, given that we have no data prior to the drought, we were unable to assess the drought response with respect to the baseline. Thus, the survey in the drought year recorded late-stage effects of the drought and might represent an assessment of seedlings and saplings that already possessed high fitness and drought tolerance (Zimmerman et al., 2021).

Tropical forests in Puerto Rico have shown a high degree of resistance and resilience to hurricanes (Beard et al., 2005; Zimmerman et al., 2021). Yet Hurricane María’s combined extreme precipitation and strong winds caused significantly higher stem breakage and tree mortality than Hurricanes Hugo (1989) and Georges (1998) in the Luquillo Experimental Forest (Uriarte et al., 2019). Despite the unprecedented strength of this storm, we only found a mild reduction in seedling richness following the hurricanes, and no effects on diversity or species composition. Reductions in seedling richness might have been caused by the significant increase in herbaceous cover observed after the hurricanes within the TRACE plots (Kennard et al., 2020), potentially leading to an increase in intra and interspecific competition for resources. In the initial years following a hurricane, tropical forest understory herbs, shrubs, and newly recruited tree saplings tend to increase in abundance due to higher light availability (Zimmerman et al., 2021). However, as canopy cover increases through successional time, we can expect a reduction in herbaceous cover, leaving space for shade-tolerant woody species to grow and recover (Scatena et al., 1996; Chinea, 1999; Walker and Sharpe, 2010).

Species richness and diversity of saplings, but not seedlings, increased in response to warmer temperatures. This finding contrasts with our expectation that temperature would negatively affect tropical plant diversity (Stork et al., 2009; Wright et al., 2009). Warmer temperatures might be selecting for pioneer species that thrive in hotter conditions, leaving the forest in a transition phase comprised of a mix of species where species that are not doing well have yet to die out (Connell, 1978; Corlett, 2011). The TRACE site is also a secondary forest with a long history of climatic and anthropogenic disturbance (Kimball et al., 2018). Secondary forests have been shown to retain high levels of late successional species, possibly in response to environmental filtering that has selected for species that do well under a range of conditions (e.g., Chazdon et al., 2009). Species composition in our study did not change in response to warming, which supports the idea that, while there is a shift in diversity of the saplings, the community is relatively resistant to change. Yet recent evidence of a decrease in root production and slower root recovery after hurricane disturbance in the TRACE heated plots points to complex interactive effects of warming and disturbance that could affect biodiversity patterns in the long-term (Yaffar et al., 2021). As such, whether this resistance continues under persistent and longer-term warming is yet to be determined.

Seedling and sapling communities differed in their responses to the various climatic stressors assessed in our study. This was expected, as ontogeny alters many ecological processes that influence plant demographics (Lasky et al., 2015) and responses to environmental change (Comita et al., 2009; Fortunel et al., 2020). In the case of drought, we only found a marginal increase in sapling richness and a significant response from a few species at the seedling and sapling stage. The abundance of P. glabrescens, an early successional shrub, was significantly higher during the drought at both the seedling and sapling stages. This result is in line with the theory that shrubs in Puerto Rico might have adapted to a wide range of soil moisture conditions (Zimmerman et al., 2021). However, the abundance of another early successional shrub species, P. brachiata, was reduced at the seedling stage during the drought. This range of results emphasizes the need to understand the functional characteristics of each plant species.

Hurricane disturbance had varied effects on seedling and sapling communities. Only seedling richness was negatively affected by hurricane disturbance, which could suggest that seedlings are more sensitive than saplings due to their lower belowground biomass or higher desiccation potential (Uriarte et al., 2005). Previous studies have highlighted life history differentiation amongst tropical woody species that result in varied susceptibility to hurricanes (Zimmerman et al., 1994; Walker et al., 2003; Canham et al., 2010). These studies support our results of both positive and negative effects of hurricane disturbance on different species. For example, the early-mid successional shrubs P. brachiata and Piper hispidum showed contrasting responses to hurricanes as well as drought, which points to niche differentiation (Connell, 1978; Uriarte et al., 2004). Heterogeneity in the severity and timing of disturbance might enable these two competing shrubs to coexist in the understory.

To date, TRACE is the only field-scale in situ warming experiment of understory plants in a tropical forest ecosystem (Kimball et al., 2018). As such, knowledge of how increased temperatures might affect seedling and sapling dynamics in tropical forests is still in its infancy, which is critical to predict future global carbon balance (Clark, 2004; Reed et al., 2012; Zuidema et al., 2013). Previous studies have shown high thermal tolerance during juvenile growth (Cheesman and Winter, 2013; Slot and Winter, 2018). Similarly, we found stronger effects of warming on sapling rather than seedling biodiversity, with species showing both positive and negative responses to warming. In particular, P. brachiata responded positively to warming at the seedling stage but negatively at the sapling stage, while P. glabrescens responded negatively to warming at both the seedling and sapling stages. After 1 year of warming prior to the hurricanes, P. brachiata demonstrated greater photosynthetic thermal acclimation capacity than P. glabrescens (Carter et al., 2020). The differing capacity of these understory shrubs to physiologically acclimate to warming could explain why the abundance of P. brachiata saplings increased with warmer temperatures, while the abundance of P. glabrescens did not. Finally, this variability in species responses to climate warming and disturbance can act as a mechanism to maintain species biodiversity through non-linearity (Chesson, 2000; Bachelot and Lee, 2020).

Overall, our finding of species-specific effects of disturbance suggests that the functional composition of the forest might eventually change. In other words, disturbance-driven changes in species abundance were predicted to significantly influence the functional composition of seedling and sapling communities. We predicted that drought would lead to higher leaf succulence (i.e., higher water storage capacity), while warming and hurricanes would lead to the prevalence of pioneer and early successional species with fast-growing strategies. Contrary to our expectations, however, we found no change in leaf succulence, or any other functional trait, associated with drought. This result could reflect the lack of pre-drought data on the woody composition, such that the community already reflects traits that are best adapted to drought.

Previous studies have shown that tropical forest seedlings and saplings experience drastically different environments following a hurricane, including increased water input more similar in chemistry to rainfall when compared with water that has passed through a closed canopy forest (Heartsill-Scalley et al., 2007). However, despite the increased water input, forests recovering from hurricane disturbance tend to have drier litter layers with slower decomposition rates due to increased temperature and solar radiation (González et al., 2014), which could in turn increase water stress and select for seedling communities with higher leaf succulence, as we observed. Similar to prior studies that have followed recruitment after hurricane disturbance (Zimmerman et al., 1994), we also found a shift toward fast-growing strategy traits such as low leaf carbon and high leaf phosphorus concentrations at the seedling stage. Overall, our findings support our hypothesis that hurricanes would shift the understory plant community toward fast-growing or drought-tolerant species, at least at the seedling ontogenetic stage. Specifically, our results show two distinct shifts through time (Reich, 2014) – water stress immediately after the hurricanes selects for drought-resistant species, while subsequent changes in the environment favor the recruitment of fast-growing species (Zimmerman et al., 1994).

Surprisingly, many more species were significantly influenced by warming when compared with drought and hurricanes, and most responded with an increase in abundance. As expected, most of these favored species were pioneer trees and shrubs such as Cecropia schreberiana, Miconia racemosa, and Psychotria berteriana. As a result, we found a shift toward low leaf carbon and high leaf phosphorus, which is consistent with a community shifting toward fast-growing strategies (Strauss-Deberiedetti and Bazzaz, 1996; Bonal et al., 2007). Similar to the functional shifts observed after the hurricane, warming also favored drought-resistant species with high leaf succulence, which tend to be slow-growing species (Reich, 2014). The selection of contrasting strategies might be the result of the short length of our study – warming favors fast-growing strategies in the short-term, but drought-tolerant species with conservative strategies might be favored under prolonged warming stress. While we did observe these shifts for the seedlings, there was no significant change in CWM of functional traits for the saplings in response to warming. Together, this might suggest that significant environmental filtering is occurring at the seedling stage (Lebrija-Trejos et al., 2010), while prevailing shrubs already possess traits that are responsive to changing environmental conditions.