Succession and seasonality drive tropical butterfly assembly after an extreme hurricane

The composition and biodiversity of insect community assemblages are mediated by a complex set of biotic and abiotic factors. Among these factors are forest structure and atmospheric variables (like temperature and humidity), which are heavily influenced by frequent hurricane activity in the Caribbean. Despite this, changes in Caribbean insect assemblages as forests recover from hurricane disturbance are poorly documented. Butterflies represent a charismatic model taxon in biodiversity and conservation, and are thus an ideal subject for exemplifying these shifts in insect abundances and diversity across ecological succession. Here, we monitored butterfly communities in two Puerto Rican forests differing in structure (i.e., canopy height, tree size) to assess butterfly diversity, abundances, and community level wing traits (size and color) over 1 year, beginning 6 months after Hurricane Maria. Monthly sampling over the course of 1 year revealed no relationships between abundances and canopy openness or humidity; instead, species abundances fluctuated seasonally and were nonlinearly correlated with temperature. In contrast, wing size and color were linearly correlated with abiotic variables. Specifically, wings were larger in cooler and more open conditions. Wing color saturation and brightness were negatively correlated with humidity. Our results suggest that, first, a functional approach may provide better insight into the factors mediating species responses to disturbances. Second, further disentangling abundance seasonality from impacts of extreme disturbances necessitates long‐term monitoring.

functional approaches has been proposed as a strategy for elucidating overall trends and underlying mechanisms.Specifically, functional traits (e.g., morphological characteristics) can indicate which environmental variables drive community assembly (Córdova-Tapia & Zambrano, 2015;Lebrija-Trejos et al., 2010).Functional diversity (the range of trait values within a community) is also considered key for understanding the relationship between diversity, community structure, ecosystem functioning (Chapin III et al., 2000;Diaz & Cabido, 2001;Naeem & Wright, 2003;Tilman et al., 1997) and responses to disturbance (Cadotte et al., 2011), in part because high functional diversity is associated with ecological resilience (Schmitt et al., 2020).
Hurricanes are one such disturbance for which a combined taxonomic and functional approach is useful for elucidating community-wide responses.Extreme hurricanes are increasing in intensity and frequency, a consequence of global climate change (Dale et al., 2001;Novais et al., 2018).These storms immediately affect the structure and dynamics of tropical forests (Heartsill Scalley et al., 2010;Lugo, 2008;Turton, 2008;Vogt et al., 1996).For example, the composition of Caribbean forests-and the insect communities that inhabit them-are profoundly influenced by shifting environmental variables during and after these extreme climatic events (Barberena-Arias & Aide, 2002;Imbert & Portecop, 2008;Luke et al., 2016).Hurricane disturbances decrease forest canopy height and density which, in turn, alters abiotic and biotic conditions (Novais et al., 2018;Shiels & González, 2014).Subsequent changes in environmental conditions therefore have both an immediate impact on communities (e.g., mortality, redistribution) and long-term impacts during succession (Schowalter et al., 2014).
Although extreme disturbances like hurricanes are increasing in tropical regions, our understanding of tropical biodiversity recovery following these events is limited.Likely because of the complexity of tropical systems compared to temperate systems, fewer empirical studies have documented community recovery after extreme disturbances in tropical regions (Devries & Walla, 2001;Willig et al., 2007).The paucity of studies in tropical latitudes is especially true for insects, which is surprising because insects are the most diverse class on Earth and are particularly vulnerable to climate change (Harvey et al., 2023).While many studies have focused on chronic drivers of global change on insect populations (e.g., Sánchez-Bayo & Wyckhuys, 2019), few studies have documented responses to acute disturbances.
Lepidoptera is a model group for understanding potential drivers of biodiversity responses to extreme disturbances (Calero-Mejía et al., 2014;Cárdenas-Lugo et al., 2015).Butterflies are key ecosystem and biodiversity indicators because they have short life cycles (Devries et al., 1999), are easy to monitor in the field (Devries et al., 1997), are relatively well studied taxonomically (Palacios & Constantino, 2006), and, as ectotherms, are sensitive to abiotic variables like temperature, humidity, and light (Brown & Freitas, 2000;Kremen et al., 1993).Yet, few studies have documented tropical butterfly responses following disturbance (Devries & Walla, 2001;Iserhard et al., 2013;Luk et al., 2019) even though Lepidoptera is among the most studied insect orders and 90% of all butterfly species occur in the tropics (Bonebrake et al., 2010).In Puerto Rico, general decline of Lepidoptera abundance and diversity was reported following hurricanes (Barberena-Arias & Aide, 2002;Schowalter & Ganio, 1999;Southwood et al., 1979).However, pioneer and gap specialists benefited from a more open forest canopy, causing population outbreaks of some species (Schowalter, 2017;Torres, 1992), and suggesting that life history strategies likely shape responses to hurricane disturbance.Indeed, both vegetation cover (Salvato & Salvato, 2007) and life history strategies defined by functional traits (Ewers & Didham, 2006) are key drivers of species responses to disturbances.
Although the application of functional traits for understanding successional trajectories has been applied broadly in plants (May et al., 2013;Vicente-Silva et al., 2016) and some insects like bees (Williams et al., 2010;Winsa et al., 2017), whether traits predict responses to disturbances has not been well studied in Lepidoptera.For butterflies and moths, functional traits like wing size and color are key aspects of dispersal (Aguirre-Gutiérrez et al., 2017;Sekar, 2012;Van Dyck & Wiklund, 2002) and thermal performance (Talloen et al., 2004).Research on functional trait responses in Lepidoptera following a disturbance is scarce, though Steffan-Dewenter and Tscharntke (1997) reported declining wing size with succession.
Declining body size in animals is emerging as a universal response to global change (Gardner et al., 2011), suggesting that smaller species may be favored in a world of frequent disturbance.
Hurricane Maria (September 2017), one of the strongest hurricanes on record in Puerto Rico, created a rare opportunity to examine underlying drivers of butterfly responses during early succession.
Understanding these complex drivers is crucial for informing conservation efforts on tropical islands, as they are often characterized by steep climatic gradients (Khalyani et al., 2016) and harbor large numbers of rare and endemic species, making them particularly vulnerable to global change (Fordham & Brook, 2010).
In the present study, we monitored butterfly populations in two nearby Puerto Rican forests after Hurricane Maria.We sought to understand trends in butterfly recovery over time during an early stage of succession.We also examined the influence of abiotic factors (temperature, humidity, and canopy cover) on butterfly abundances, diversity, and key functional traits (i.e., wing size and color).
We predicted that (1) butterfly abundances and species richness would increase unidirectionally or, at least, asymptotically over the sampling period as the forest recovered, (2) that functional diversity metrics would increase alongside taxonomic richness, and (3) that wing traits would change over time.Specifically, we predicted that the high light environment immediately following disturbance should select for smaller, less melanized species, eventually replaced by larger and more melanized species as the forest recovered (Figure 1).
In addressing these questions, we considered both taxonomic and functional diversity as an integrated approach for understanding successional trajectories and complex responses to disturbance.

| Study sites
We monitored butterfly communities in two forests located in Puerto Rico's Cordillera Central: the old mining area in Bosque del Pueblo (Adjuntas, Barrio Vegas Arriba, 18.18° N: 66.68° W) and Bosque Escuela La Olimpia (Adjuntas, Barrio La Olimpia, 18.143° N: 66.71° W) (Figure S1).We chose these sites based on their proximity to each other and differing forest structure (i.e., canopy cover, tree size) prior to the hurricane.Both forests are classified as wet subtropical forests (Ewel & Whitmore, 1973), and can be considered secondary or novel forest (Lugo, 2009) (Daly et al., 2003).Hurricanes Irma and Maria moved across the area in September 2017, with high wind speeds and rainfall causing compounded damage to foliage, increased tree mortality, widespread flooding, and an unprecedented number of landslides, altering forest structure and plant species composition (Bessette-Kirton et al., 2019;Hall et al., 2020).

| Specimen collection and preservation
Due to limited mountain access immediately following the hurricane, we began monitoring and collecting specimens 6 months post-Hurricane Maria (April 2018).We sampled monthly at each sampling site for 1 year, except for July and August 2018, due to logistical constraints.We captured butterfly specimens using nets, fruit bait traps, and visual counts along major trails between 800 and 1500 h, a time of peak insect activity which encompassed species that prefer different environmental conditions (Villarreal et al., 2004).We used Van Someren Rydon fruit bait traps to catch fruit-feeding adults with fermented banana (Andrade et al., 2014).We deployed five traps at each site for each sampling occurrence, placed 20 m apart, hanging from trees at a height between 1 and 3 m above ground level (following Villarreal et al., 2004).They remained for at least 3 h, as suggested by previous studies (Gaviria Ortiz, & Henao., 2014;Millán, Chacón, & Giraldo, 2009;Palacios & Constantino, 2006).For visual counts, we performed standardized transect walks along a 2 km route, with approximately seven person hours dedicated to sampling per month per site.
We collected three individuals per species, usually one female, one male (when sex could be differentiated in the field) and an additional individual, to prepare as voucher specimens.We then froze, relaxed, mounted, and dried the collected specimens before depositing them at the Museo de Zoologia at the University of Puerto Rico-Rio Piedras (museum code: UPRRP).

| Abiotic variables
To evaluate forest recovery, we measured canopy cover using a fisheye lens (ECO-FUSED, U.S.) attached to a smartphone.We took digital images at two locations each month at each site and analyzed them using the Gap Light Analyzer (Frazer et al., 1999) which calculates the percentage of occupied pixels.We then averaged the two values to obtain one canopy openness value per month per site.
Additionally, we measured temperature and humidity using data loggers (HOBO Prov2;Onset Computer Corporation, 2010).We replaced one data logger in September 2018 in Bosque del Pueblo because the deployed device was stolen.

| Functional traits and diversity
We photographed all voucher specimens using standardized protocols (Seltmann et al., 2017) and used the resulting images to calculate wing size using image analysis software (ImageJ/FIJI; Ferreira & Rasband, 2012).We measured two forewing traits on all specimens: wing length, calculated as the length in centimeters of the longest axis from the point of wing attachment at the thorax (disc cell) to the distal tip (R4 vein), and wing width, calculated as the width in centimeters of the longest line (R2 vein) that could be drawn perpendicular to the A2 anal vein (Altizer & Davis, 2010; Figure S2).We then calculated length-to-width ratio as an indicator of wing shape.We also measured wing color using image analysis software on a 1 cm 2 square area near the thorax, avoiding veins and regions with damaged or missing scales.We converted the raw RGB (red, green, blue) digital images to HSB images (hue, saturation, brightness, and intensity color space; Schindelin et al., 2012).
Next, to analyze functional community structure through time, we calculated community weighted means for wing traits (CWMs) and functional diversity metrics (FDMs) based on species composition and abundances at each site for each sampling occurrence (Garnier et al., 2004), using the vegan (Oksansen et al., 2020) and FD (Laliberté et al., 2014) packages in R version 4.2.2 (R Core Team, 2022).CWMs can be used to detect environmental selection for certain traits (Perović et al., 2015) and, we argued, should vary during succession.FDMs describe various dimensions of trait space.
Functional richness and evenness are related to the distribution of traits in space, though functional richness tends to be independent of abundances and is thus strongly related to observed species richness (Mason et al., 2005).Functional divergence is a measure of similarity among dominant species.High functional divergence indicates greater functional differentiation and more diverse resource use among species (Córdova-Tapia & Zambrano, 2015).Low functional divergence reflects a high degree of functional similarity, indicating homogenization of the lepidopteran community.In contrast with other components, functional dispersion considers relative species abundances and is therefore less sensitive to the presence of rare species (Laliberté & Legendre, 2010).Further, functional dispersion can be used to estimate the degree of functional specialization among species; species that are functional specialists will be farther from the center of gravity of the community (Córdova-Tapia & Zambrano, 2015;Villéger et al., 2010).Thus, low functional dispersion suggests an assemblage consisting of mostly generalists.

| Statistical analysis
First, to analyze taxonomic community structure through time, we calculated species richness for each sampling event at each site using true Hill numbers (Chao et al., 2014).We then calculated extrapolation curves according to the iNEXT algorithm, using the iNEXT package in R (Hsieh et al., 2016) to estimate total species diversity over time.We also used an Analysis of Similarity (ANOSIM) to determine differences in butterfly taxonomic community composition and abundances (Calero-Mejía et al., 2014) between seasons (wet/dry), forests (Bosque La Olimpia/Bosque del Pueblo), and over time (first/second half of sampling).Finally, we performed Non-metric Multidimensional Scaling (NMDS) analysis using vegan to illustrate differences in species composition between the two study sites.The NMDS stress value measures how well the multidimensional data is represented in reduced ordination space, with a value of <0.2 indicating a reliable fit.Bray-Curtis was selected as the distance metric for both the ANOSIM and NMDS based on the degree of coverage (Figure S3).Second, to test for abiotic and trait-mediated drivers of butterfly recovery, we performed Generalized Additive Models (GAMs) using the mgcv package (Wood, 2011) in R. All variables were scaled to a mean of zero and a standard deviation of one.We looked for relationships between abiotic variables (temperature, humidity, and canopy cover) and butterfly response variables (taxonomic diversity, abundances, CWM trait values, and FDMs).We chose to use GAMs because we aimed to simultaneously capture linear and nonlinear relationships, a necessity for disentangling complex environmental drivers of community composition.To test for changes in response variables over time, and for direct correlations between the abiotic variables of interest, we used linear regression analyses.

| Variation across months and sites
Species richness and the number of individuals observed peaked in each forest during September 2018 (wet season) and again in January 2019 (dry season) (Figure 2a,b; see also: Figure S4).Temperatures remained higher throughout the sampling period in Bosque del Pueblo when compared to Bosque La Olimpia, while humidity was generally lower (Figure 2c,d).Percent canopy openness was higher in Bosque La Olimpia (53% at the start) compared to Bosque del Pueblo (29%) but declined in both sites during the study period (Figure 2e).
In general, abundance fluctuations (lows and highs) were detected during wet and dry seasons in both forests.The ANOSIM showed a difference in species composition between seasons (R = 0.38, p = .004).While the ANOSIM did not reveal a difference in composition between forests (R = 0.11, p = .071),the NMDS illustrated segregation between the two sites, with a stress value of 0.1587 indicating a reliable fit in reduced ordination space (Figure 3).This disparity is likely because the NMDS is a more sensitive test, able to detect subtle differences in species composition between sites.In the NMDS, the first half of the sampling months were clustered in negative space of axis 2 while the second half of sampling was clustered in positive space of axis 2. We further confirmed this shift over time in the ANOSIM (R = 0.66, p = .001).These results indicate that species composition differed between seasons, between locations, and through time.While simple species richness (q = 0) was similar between sites (32 species in Bosque del Pueblo; 33 species in Bosque La Olimpia), Hill's diversity was higher in Bosque del Pueblo at higher orders (q = 1 and q = 2; Figure S3), indicating a more even assemblage.

| Environmental variables
To determine whether shifts in environmental conditions during succession mediated butterfly species richness or the number of individuals observed (abundance), we measured temperature, humidity, and canopy openness in each site.We found that there were no relationships between species richness or numbers of individuals with humidity or canopy openness, however, in the GAM analysis, temperature was nonlinearly correlated with abundance (p = .0087;Table 1; see also: Figure S5F).In the linear regression analysis, temperature and canopy openness were uncorrelated (R 2 = 0.048, p = .40),as were temperature and humidity (R 2 = 0.0013, p = .88),but canopy openness was positively correlated to humidity (R 2 = 0.31, p = .020).

| Functional traits
We found only two linear relationships between trait CWMs (community weighted means) or FDMs (functional diversity metrics) over time.Wing light intensity decreased (Figure 4c; R 2 = 0.26, p = .014)and functional divergence increased (Figure 4g; R 2 = 0.17, p = .057)over the sampling period.Instead of unidirectional or asymptotic trends, seasonal peaks and valleys were more prevalent (Figure 4).

| Taxonomic diversity and composition
We examined whether abiotic factors (temperature, humidity, canopy cover) or wing traits (size, color) mediate the recovery of butterfly diversity post-disturbance.We began with a taxonomic approach to determine whether species abundances and richness increased over time during succession, and whether abiotic variables drove those changes.
Contrary to expectations, neither abundances nor richness exhibited unidirectional or asymptotic changes over time, rather, seasonal peaks and valleys occurred.In addition, abundances were not linearly related to any abiotic variable (temperature, humidity, canopy openness), but were nonlinearly related to temperature, indicating relatively low abundances at mid-range temperatures and perhaps further reflecting seasonal peaks (Figure S5F).Richness (q = 0) 1.868 a 1.000 1.000 0.246 42.9% Note: for edf values of 1.000, the direction (positive or negative) of the correlation is provided in parentheses.For edf values greater than 1.000 (indicating nonlinear correlations), the direction is provided only if there were visually obvious trends in the partial dependence plots.

TA B L E 1 Summary of generalized additive model (GAM) analyses correlating community weighted mean (CWM)
wing traits, functional diversity (FD) metrics, abundance, and taxonomic species richness with abiotic variables (temperature, humidity, and canopy openness).Analyses were conducted using the mgcv package in R version 4.2.2.Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 ' a ' 0.1.An edf (effective degrees of freedom) value close to 1.000 indicates a linear relationship; the greater the edf, the more nonlinear the relationship is between response and predictor variables.
variables and butterfly abundances is the duration of the study period.Butterfly communities oscillate (in abundance and species richness) seasonally and long-term (Bonebrake et al., 2010;Cárdenas-Lugo et al., 2015;Davila & Asencio, 1992;Devries et al., 1997;Devries & Walla, 2001;Young, 1982).For example, after Hurricane Patricia in Mexico, Novais et al. (2018) found that catastrophic hurricanes temporarily increased insect species diversity and abundance in tropical dry forests.However, insect communities continued to oscillate seasonally, as we also show here.Other research on the Lepidoptera of Mexico also emphasizes the interacting effects of succession and seasonality (Essens et al., 2010).While these processes are infrequently investigated in tandem using long-term data, complex interactions between the two have also been observed in tropical insect groups outside of Lepidoptera.For example, in dung beetles (De Siqueira Neves et al., 2010) and ants (Marques et al., 2017), species richness generally increased with succession while varying drastically between seasons.Thus, more long-term monitoring is needed to explore potential relationships between insect abundances and richness following disturbances and concomitant changes in abiotic factors within the context of tropical seasonality.
We were also interested in whether shifts in species composition pointed to differing responses to disturbance.We found that species composition varied between seasons and forests.This may be due to differences in forest structure and plant species composition prior to the hurricane.As a result, the two forests differed in subsequent changes of percent canopy openness, thus supporting different butterfly assemblages.Additionally, while butterfly species richness was similar between the two study sites, evenness was higher in Bosque del Pueblo when compared to Bosque La Olimpia (Figure S3).This could be due to greater habitat heterogeneity in Pueblo-a forest which contains both areas of dense vegetation and more open, fire-disturbed areas-thus supporting a more diverse plant assemblage with a more robust range of ecological niches.
In both forests, some species were highly seasonal, most appearing in January (the early cool, dry season).The appearance of these species may be partially explained by the recovery of vegetation after Hurricane Maria.Alternatively, more favorable environmental conditions could occur during the dry season, suggesting that butterflies may occur at low abundances during hurricane season, perhaps in diapause.
If this is true, then this suggests that butterfly assemblages may have adapted to hurricane disturbances by shifting their peak abundances to a time when hurricanes are unlikely to occur, thereby reducing mortality.Yet, other species were observed in relatively high abundances throughout the year.These species may be generalists, able to feed on a wide range of host plants and persist in a variety of habitat and environmental conditions.Indeed, high abundances may be attributed to generalist feeding habits, abundant resources such as nectar, or abundant host plants (Ríos-Málaver, 2007).Though beyond the scope of the present study, including plant phenology and plant species composition in addition to abiotic factors may better disentangle seasonal changes in herbivorous insect species composition (Bonebrake et al., 2010;Cárdenas-Lugo et al., 2015;Olarte-Quiñonez et al., 2016;Ospina-López & Reinoso-Flórez, 2015;Palacios & Constantino, 2006) from changes due to disturbance.
It should be noted that the lycaenids are possibly undersampled in the present study.In a database of over 8000 museum records from Puerto Rico, lycaenids represented 4.7% of all specimens, whereas they represent 1.4% in our field observations.In addition, there are 20 species of Lycaenidae in the museum record, but only  two species appeared in our field observations.Neither of these species belong to the diverse Theclinae tribe, which contains six species in the Puerto Rico record.This disparity could reflect several artifacts, such as our sampling being insufficient for these smaller winged, often elusive species, collector bias in the collection record, or these species specializing in habitats and/or regions not encompassed by our two sampling sites.

| Functional traits
While a taxonomic approach revealed seasonal variation in butterfly assemblages, a functional approach may better explain the underlying drivers of these shifts.We asked whether wing size or color shifted, reflecting changes in environmental conditions following disturbance.We predicted that smaller, less melanized species should be favored by the hot, high light environments of early succession due to thermoregulatory demands, with assemblages becoming larger and more melanized as succession progressed.We found few unidirectional trends in wing size or color over time.However, we showed that CWM traits were correlated with temperature and humidity.Butterfly assemblages were indeed smaller in warmer environments, though, contrary to our predictions, they were also smaller in the more shaded conditions typically associated with later stages of succession.
Wings were also brighter and more saturated when humidity was lower.These results highlight the role of seasonal environmental conditions in selecting for functional traits and suggest that a longer sampling period is necessary for elucidating morphological trends across succession.
Wing size and saturation declined in the hottest months (July to September), then increased during the coolest months (October to April), indicating seasonal trends during butterfly recovery.In line with these dual processes, wing size was positively correlated with canopy openness (larger species in more open canopies) suggesting that trait-mediated recovery occurred within the background of seasonal shifts.Alanen et al. (2011) found that small, less mobile lepidopteran species colonized habitat fragments slower than large and more mobile species during succession.Indeed, in our study, the smallest winged species (wing length <2 cm) were only observed 1 year after the hurricane (October 2018) and only in one forest (Pueblo).This suggests that larger-winged species with a greater ability to disperse long distances may have an advantage immediately following a major disturbance, in contrast with our prediction that smaller species would be favored.Other studies have shown that butterfly body size decreased during succession (e.g., Steffan-Dewenter & Tscharntke, 1997) and, in Coleoptera, wing length was greater in disturbed habitats (Venn, 2007), further supporting the idea that small-winged species may be more impacted by disturbances or require longer recovery times.Perhaps some species disappeared, or declined so dramatically that a much longer recovery period will be required, as reported in the Florida Keys after various hurricanes (Salvato & Salvato, 2007).Whether this occurred in the present study system, however, is difficult to determine given the lack of pre-hurricane data.
Interestingly, wings were more elongated in cooler, drier conditions, and rounder and stubbier with increasing temperature and humidity.The negative linear relationship between length-to-width (or aspect) ratio and humidity was particularly strong.This relationship is surprising as several existing analyses have, in contrast, found positive correlations between wing aspect ratio and various proxies of humidity.For example, one study which exclusively evaluated the tropical species Bicyclus anynana (Nymphalidae) found increased wing aspect ratio under wetter plant conditions (Kuczyk et al., 2021).
Other studies found lower aspect ratio in one temperate species (Pararge aegeria; Nymphalidae) when raised on drought-stressed plants (Gibbs et al., 2012) and at higher latitudes (and thus, lower average humidity and temperature in the study region ;Vandewoestijne & Van Dyck, 2011).However, there are no existing community-wide studies which have directly related shifts in humidity to variation in aspect ratio among Lepidoptera.The present study therefore offers a novel perspective by demonstrating a link between humidity and wing morphology across multiple species.Wing aspect ratio is related to important aspects of dispersal, including maneuverability, flight speed, and gliding efficiency (Betts & Wootton, 1988;Le Roy et al., 2019).Here, rounder wings may have been advantageous for dispersal capacity in dense, high moisture air conditions.Contrary to our predictions, wing light intensity decreased, becoming duller over time.According to thermoregulatory theory (Kingsolver, 1995(Kingsolver, , 1996) ) there is a disadvantage to intense coloration in early stages of succession due to high thermal energy and increased visibility to predators (Vane-Wright & Boppré, 1993).Yet, it is unclear to what degree the decline in intensity reflects seasonal variation in temperature and humidity rather than successional processes.Wing color saturation and brightness were negatively correlated to humidity, providing further evidence that this particular wing trait may be more indicative of seasonal changes.In the tropics, lower humidity may be related to higher solar radiation, resulting in brighter individuals due to changes in wing melanization (e.g., due to thermoregulatory effects).
Increased melanin can increase body temperatures during basking, thus increasing flight activity in cool temperatures (Kingsolver, 1995), though excess melanization may result in overheating when ambient temperatures and solar radiation are high (Kingsolver, 1996).In addition to environmental factors, sexual selection (Allen et al., 2011) and larval development (Hoffmann, 1973;Koi & Daniels, 2017) may also influence the results shown here as these factors are associated with variation in wing pigmentation and size.

| Functional diversity
Finally, some indices of functional diversity appeared to be related to reorganization and subsequent succession following the hurricane, while other indices appeared more related to seasonal shifts in temperature and humidity.Functional divergence was low for both forests shortly after the hurricane, increased over time, and had a strong nonlinear correlation with humidity.Low functional divergence reflects high functional similarity, indicating homogenization of the lepidopteran community shortly after the hurricane.Biotic homogenization, defined as an increase in similarity among biota, can occur when an environmental disturbance promotes the expansion of some (often widespread, invasive) species and reduction of other (often rare, indigenous) species (Dar & Reshi, 2014).This phenomenon is predicted to intensify and become more widespread with climate change (e.g., Hidasi-Neto et al., 2019;Savage & Vellend, 2015;Zwiener et al., 2018) and has already been reported as a result of anthropogenic land use change (e.g., Devictor et al., 2008;Gossner et al., 2016) and wildfires (e.g., da Silva et al., 2018) in some plant and animal communities.However, the role of hurricane disturbance in homogenization is not well understood.A few studies have reported that hurricanes speed up succession of indigenous plant species in Puerto Rico, leading to functional heterogenization-not homogenization (Abelleira Martínez, 2010;Aide et al., 2000;Flynn et al., 2010).It is not understood whether this effect cascades to insects.One study suggests future climate scenarios may lead to homogenization of the butterfly fauna in Puerto Rico (Hulshof et al., 2024).How this plays out on smaller scales or following more acute disturbances is unclear.
In this study, it appeared that seasonality played a more important role in determining functional dispersion even after this extreme disturbance.Functional dispersion decreased during the warmer months and increased during the cooler months.This finding suggests that functional specialists (i.e., species with functionally rare wing size or color values) are more prevalent in the cooler, dry season, a time when hurricanes are less likely to occur.Other functional indices (evenness and richness), showed no obvious trends over time, suggesting that the evaluation of multiple dimensions of functional diversity can be useful for understanding shifts in species composition (Figure 1; Villéger et al., 2008).

| CON CLUS ION
Here, we sought to examine environmental drivers of butterfly community assemblages following an extreme disturbance.We found that species abundance and community composition varied seasonally over the sampling period.Although abundances and species richness were not linearly correlated with abiotic variables, there was evidence that wing traits (size and color) were impacted by changes in environmental conditions.Wing size was positively correlated with canopy openness and negatively correlated with temperature, reflecting two interacting processes -response to disturbance and seasonal shifts.We also provided evidence that wing shape was influenced by abiotic conditions (i.e., wings were more elongated in cooler and dryer environments).We did not observe unidirectional or asymptotic trends in wing traits over time.These results further point to complex interactions between functional traits, succession, and seasonality.Taken together, the results suggest that the hurricane differentially impacted smaller sized and functionally specialized species, contributing to a growing understanding of functional homogenization in response to ecological disturbance.Our results also highlight that seasonality is an important (and often overlooked) component of succession.

| RECOMMENDATIONS AND NEXT STEPS
Given the complexity of global climate change impacts on butterfly diversity, long-term monitoring in the tropics is critical (Luk et al., 2019).Extreme disturbances, like hurricanes, will continue to increase in intensity (Novais et al., 2018), which will likely have cascading consequences on insect communities worldwide (Schowalter, 2012).This study is a first step toward better understanding how butterfly assemblages are impacted by environmental variables depending on their functional traits.Future monitoring should seek to record and integrate host plant data for a more complete understanding of fluctuations in herbivorous insect diversity and abundances.Additionally, although the butterfly fauna of Puerto Rico has been well studied from a taxonomic perspective (e.g., Davila & Asencio, 1992;Ramos, 1982;Torres & Medina-Gaud, 1998), the ecology and life history of most species are less understood (Terry et al., 2023).Future studies should seek to fill this important gap, both locally and regionally, by continuing to measure trait variation for butterfly assemblages in the tropics, a region for which there is scarce trait data (Bonebrake et al., 2010).Robust monitoring of Lepidoptera communities will require longer-term research due to high seasonality and variability within and between years and sites (Kocsis & Hufnagel, 2011).

F
I G U R E 1 Summary of expected and observed results for trends in (a) species richness and abundance; (b) functional trait values; and (c) indices of functional diversity.

F
Butterfly abundances and abiotic variables over time across the two sampling sites, Bosque La Olimpia (turquoise) and Bosque del Pueblo (blue) for (a) species richness; (b) number of individuals observed; (c) temperature (°C); (d) humidity (%); and (e) canopy openness (%).Interruption of line is due to missing data (no sampling).
Due to the predominance of wet-dry seasonality in the tropics, one possibility for the lack of relationships between environmental F I G U R E 3 Non-metric Multidimensional Scaling (NMDS) of butterfly abundances and species composition in Bosque La Olimpia (turquoise) and Bosque del Pueblo (blue).Wet season months are represented by filled circles, while dry season months are represented by filled triangles.