Temperature seasonality drives taxonomic and functional homogenization of tropical butterflies

To better understand the potential impact of climate change on butterfly assemblages across a tropical island, we model the potential for taxonomic and functional homogenization and determine climate‐ and trait‐mediated shifts in projected species distributions.


| INTRODUC TI ON
Worldwide shifts in species distributions in response to contemporary climate change are altering the taxonomic and functional composition of plant and animal communities, ultimately leading to biotic homogenization (i.e. the convergence of ecological communities over time; Parmesan, 2006;Wardle et al., 2011).
Whereas taxonomic diversity describes the numbers and abundances of species, functional diversity describes the range of life history traits among species (McGill et al., 2006).The extent to which taxonomic and functional measures of diversity are linked is not entirely clear, but these connections may help understand mechanisms underlying the re-assembly of ecological communities under intensifying global change (Devictor et al., 2010;Green et al., 2022).At local scales, a growing number of field-based studies conducted mostly in temperate regions reveals shifts in species and trait distributions among diverse taxonomic groups (Green et al., 2022;MacLean & Beissinger, 2017), yet the rate of change between taxonomic and functional metrics often differs.
At regional scales, data collected from networks of permanent plots show that relationships between taxonomic and functional diversity are idiosyncratic, varying across taxa (Hevia et al., 2016) and ranging from high covariation (Arnan et al., 2018;Reymond et al., 2013) to low or none (de Castro et al., 2020;Kaltsas et al., 2018;Nunes et al., 2017).These idiosyncrasies make it challenging to predict the re-assembly of ecological communities based on taxonomic patterns alone.A complementary approach to field-based sampling involves the use of species distribution models (SDMs) to predict both taxonomic and functional diversity patterns at larger spatial and temporal scales.In tropical regions, the application of such an approach is particularly urgent given that many animals are operating near their critical thermal maximum (Deutsch et al., 2008) and that narrow thermal tolerance breadths may constrain species' responses to shifting climates (Bota-Sierra et al., 2022).Documenting distribution shifts and testing trait-based hypotheses of rapid ecological change using a modelling approach is an obvious first step for understanding the complexity and consequences of climate change responses across the tree of life.
Insects are one of the most species-rich group of organisms on Earth, with an estimated 5.5 million species (Stork, 2018).
Two traits in particular -body size and melanization -are thought to play important roles in mediating responses to climate change.First, the temperature-size rule hypothesis posits that declining body size should emerge as a universal response to warming temperatures (Atkinson, 1994;Gardner et al., 2011;Ohlberger, 2013).As a result, body size is predicted to be smaller in warmer environments.This relationship holds even when elevation or latitude is used as a proxy for environmental temperature (Angilletta & Dunham, 2003).The temperature-size rule is thought to arise because larger organisms incur greater metabolic demands, have higher equilibrium body temperatures and gain and lose heat more slowly (Nijhout, 2003).Second, the thermal melanism hypothesis predicts that declining melanization (i.e. the reduction in pigment) may be another widespread response to warming temperature (Clusella-Trullas et al., 2008), given its relationship with solar absorption (MacLean et al., 2018).Thus, less melanization in ectotherms is favoured in warmer environments because it enables slower heating rates and lower body temperatures.For insects, some studies indicate warmer environments select for reduced body size (Attiwilli et al., 2022;Brehm et al., 2019;Henriques et al., 2022;Leingärtner et al., 2014;Xing et al., 2016) and melanization (Brakefield & de Jong, 2011;Stelbrink et al., 2019;Xing et al., 2016;Zeuss et al., 2014), while others report the opposite (MacLean et al., 2018) or show no relationship (Brehm & Fiedler, 2004;Dufour et al., 2018).In short, general patterns of insect body size or melanization responses to climate change have not yet emerged.
Butterflies are an excellent group of insects to model climaticand trait-mediated responses to climate change as evidenced by efforts in temperate (e.g.Aguirre-Gutiérrez et al., 2017;Girardello et al., 2019;Mattila et al., 2011;Ockinger et al., 2010) and tropical (e.g.Meineke et al., 2018) regions.First, butterflies have many historical records and reliable identification guides.The multitude of occurrence records enables the quantification of current and future trends in taxonomic diversity and distributions at broad spatiotemporal scales.Furthermore, butterflies are well represented biotic homogenization, Caribbean, ectotherm, Lepidoptera, melanism, SDM, trait in a two-dimensional image, facilitating the measurement of trait data from museum specimens.These data may then be used to understand the mechanisms underlying changes in species distributions and the re-assembly of ecological communities.Given that 90% of butterfly species are tropical (Bonebrake et al., 2010) and tropical ectotherms may be more vulnerable to warming (Deutsch et al., 2008;García-Robledo et al., 2016), climate change in the tropics may lead to rapid taxonomic and functional homogenization of the butterfly fauna, as certain species or traits are more favoured over others.
In this study, we use thousands of digitized collection records to model whether climate change will result in the taxonomic and functional homogenization of a tropical butterfly fauna by focusing on three questions.First, we ask if forecasts point to taxonomic or functional homogenization by modelling and comparing the current and future distribution of taxonomic and functional richness and turnover.Second, we ask what are the climatic-mediated effects on modelled taxonomic and functional composition, and specifically, does evidence support temperature as the primary driver of current and future changes of projected species and trait distributions?Lastly, we ask whether wing traits mediate changes in modelled species distributions, owing to their roles in thermoregulation, and examine interactions among species distributions, traits and climate.Across a topographically complex Caribbean island, we predicted that modelled taxonomic and functional richness and turnover would decline in cooler environments.This expectation was based on the idea that environmental selection favours species that are more melanized and larger-winged in cooler environments due to thermoregulatory constraints, per the temperature-body size rule and the thermal melanism hypothesis.
As a result, under future climate scenarios, taxonomic and functional homogenization should increase in cooler environments due to upslope distributional shifts of lowland species.Simultaneously, warmer lowland environments may also exhibit reduced species and trait turnover due to biotic attrition.These expectations may be complicated by interactions among species distributions, traits and climate, which we model here, in addition to land use change, which we discuss further.Doing so can augment our understanding of potential shifts in species' distributions and the possible selection of functional traits under future climates.We can use these modelled results to develop testable hypotheses to identify species most at risk, anticipate cascading impacts on ecosystems and guide conservation priorities, which could serve as important benchmarks for other Caribbean and tropical regions.

| Study site
The island of Puerto Rico is located within the Caribbean Biodiversity Hotspot, one of the 'hottest hotspots' (Myers, 2003).With an area of 8876 km 2 , the island encompasses varied climatic, topographic and edaphic conditions over small spatial scales, creating a landscape that harbours diverse ecological and physiographic regions (Ewel & Whitmore, 1973;Helmer et al., 2002;Miller & Lugo, 2009).
In Puerto Rico, variation in temperature is mostly driven by elevation (0-1340 m), whereas variation in precipitation is driven by a combination of prevailing winds from the north-east and interior mountains, producing a rain shadow effect where xeric conditions occur in the south and south-west and moist conditions in the north and north-east (Daly et al., 2003).Climatic variability determines the presence of lowland moist seasonal evergreen forests, submontane and lower montane wet evergreen sclerophyllous forests and lower montane wet evergreen forests (including both tall and elfin cloud forests) along the north and east sides and lowland dry semideciduous forests, lowland dry mixed evergreen drought-deciduous shrubland with succulents and lowland dry semi-deciduous forests along the south-west side (Helmer et al., 2002, and see Figure S1.1).
These plant formations occur within three major physiographic regions (Murphy, 1916): the lowland coastal and interior plains; the Cordillera Central, a mountain range that extends in an east-west direction; and the Sierra de Luquillo, an isolated mountain range in the north-east.The karst formations in the northern region (mogotes), and the Sierra de Cayey -a south-eastern extension of the Cordillera Central and the Caguas-Juncos Valley, the largest interior valley on the island -are also identified as unique physiographic regions (Smith & Pico, 1976) (Figure 1).

| Sources of butterfly data
The main source of data came from the Stuart J. Ramos Collection (SJRC; 3410 records; SJRC) which was donated by the late S. J.  (Franz & Yusseff Vanegas, 2009).The collection represents efforts of the late entomologist Stuart José Ramos Biaggi, an expert on the taxonomy, biogeography and ecology of Puerto Rican butterflies (Franz & Yusseff Vanegas, 2009;Ramos, 1996).We digitized the Stuart J.
Ramos Collection following standardized protocols (iDigBio; www.idigb io.org), which included creating digital images of each specimen positioned on a standard grey background, with a colour card standard and scale, under constant light conditions.In addition, for each specimen, we transcribed the geographical location described by its pinned data label and georeferenced location data using the batch georeferencing tool (GeoLocate) in the Symbiota Collections of Arthropods Network (Gries et al., 2014;SCAN, 2019) where images and records are available.

| Species occurrences
The total number of compiled occurrence records was 8193.We examined records for accuracy and spatial bias, which resulted in the removal of 549 duplicate records and the elimination of records of unidentified species, species with less than 10 occurrences and nocturnal species.In addition, we spatially thinned species occurrence records using spThin (Aiello-Lammens et al., 2015) to reduce the effects of spatial sampling bias.This process takes the set of occurrence records per species and identifies a new species-level subset that meets the specified minimum nearest-neighbour distance (5 km) while retaining the largest number of records.The final cleaned set of occurrence records totalled 2003 for 62 species belonging to 5 butterfly families (16 Hesperiidae, 4 Lycaenidae, 25 Nymphalidae, 1 Papilionidae and 16 Pieridae).These five families represent the most species-rich butterfly families in Puerto Rico (Ramos, 1996).

| Species distribution models (SDMs)
The 2003 records were used to generate SDMs which were created following a three-step procedure.First, we selected bioclimatic variables and checked for collinearity, eliminating highly correlated variables; second, we used the remaining variables to generate species distribution models using an ensemble approach; and lastly, we tested model performance.We describe each of these steps in more detail below.Current climate conditions were based on WorldClim (v 2.0; Hijmans et al., 2005) which includes elevation and 19 interpolated climatic variables.The WorldClim bioclimatic dataset has a spatial resolution of ~1 km 2 and represent annual trends (e.g.mean annual temperature and annual precipitation), seasonality and extreme conditions (e.g.temperature of the coldest and warmest month and precipitation of the wet and dry quarters) calculated as an average across a 30-year period .These 1-km 2 pixels were used to define our present-day butterfly communities (see Zurell et al., 2020).Future climate conditions were based on the Representative Concentration Pathway RCP 4.5, which is a more conservative scenario compared to other RCPs, and the global climate model HadGEM2-ES, which performs well at predicting bimodal seasonality for precipitation and temperature trends in Puerto Rico (Hayhoe, 2013), whereas many other GCMs do not.To reduce collinearity among predictor variables, we calculated variance inflation factors (VIF) with a threshold of less than two (Queen et al., 2002).As a result, six variables were included in species distribution models representing seasonality (isothermality, temperature and precipitation seasonality and temperature annual range); extreme conditions (precipitation of driest month and precipitation of the warmest quarter); and elevation.Together, these variables describe the bimodal precipitation and temperature seasonality characteristic of the island (Hayhoe, 2013).We used the biomod2 library in R to generate SDMs (Thuiller et al., 2016) using an ensemble approach under current and future climate conditions.The algorithms we chose to model species distributions were maximum entropy (MAXENT), random forest (RF) and a generalized linear model (GLM).Maximum entropy performs well with few occurrence records and was developed to model presence-background data (Phillips et al., 2006).Random forests are machine learning algorithms which use bootstrap-based classification and regression trees.They are widely used in SDMs and are robust to overfitting (Lawler et al., 2006).GLM is an extension of classic regression modelling, which uses linear or higher-degree polynomials to model the relationship between response and predictor variables.Spatially thinned occurrence points were treated as presences, and 10,000 background points were randomly selected as pseudo-absences.We used 80% of presences and pseudo-absences for calibration; the remaining 20% of data were used to test model performance.We then calculated a consensus of predictions based on the weighted area under the curve (AUC) average of the three algorithms, selecting model outputs where the AUC values were greater than 0.6, based on 45 SDMs per species (2790 SDMs total).The number of models meeting or exceeding this threshold ranged from 2 to 43 models per species, which were averaged to create final species distribution maps.
Lastly, we used the true skill statistic (TSS, Allouche et al., 2006) to evaluate models, discarding models with TSS < 0.60.TSS varies from −1 to 1, with values less than 0 representing no discrimination and 1 representing perfect discrimination (Araújo & New, 2007).We used default settings in biomod2 (as suggested by Thuiller et al., 2009) to transform the ensemble model predictions (both current and future) into binary maps of presence/absence using species-specific thresholds that maximize specificity and sensitivity (Liu et al., 2011).Given the small area of Puerto Rico, its biogeographic history (most butterfly species dispersed to the island from elsewhere, Miller & Miller, 2001) and high forest cover (>50%, Helmer et al., 2002), we assumed perfect climate tracking (i.e.no dispersal limitations).

| Functional traits
While all specimen collections were sources of occurrence records, we only used the Stuart J. Ramos Collection to measure wing traits because it was (and is) the only digitized collection.
We used FIJI (Abràmoff et al., 2004) to measure dorsal forewing length as the length between the forewing tip and the site of wing attachment on the thorax, and forewing width as the widest section of the forewing (Figure S1.3).To facilitate generalizations, we categorized species based on their wing length.Species with a forewing length less than the mean (3 cm) were classified as 'smallwinged', while those with a forewing length greater than the mean were classified as 'large-winged'.We also measured hue, saturation and brightness (HSB) on a 1 cm 2 region of the dorsal forewing near the thorax, avoiding major veins and visible scale damage, as this region is important for thermoregulation (Kingsolver, 1987;Wasserthal, 1975;Watt, 1969).The HSB colour space provides a standardized metric for comparing colour values and is correlated with melanin content (Foster et al., 2023;McGraw et al., 2005).Hue represents the colour (i.e.red and yellow), saturation represents the intensity and brightness represents the lightness or darkness of the colour.Because melanin absorbs light, organisms with high melanin content will have lower HSB values (i.e.darker colours) and vice versa.Organisms with low melanin content will have higher HSB values (i.e.lighter colours).This method, however, emphasizes pigment colouration and ignores structural colouration or iridescence.The impact of structural colouration on butterfly thermoregulation is not well known (Krishna et al., 2020).We focused on females because their parental investment and, thus, energetic requirements are often greater than those of males, and male colouration can be biased by sexual selection (Clutton-Brock & Parker, 1992).Of the 62 species, only 4 (all Nymphalidae) had visually noticeable size differences, with females being larger than males.In a few instances in which images of females were not available, or the sex could not be determined from visual characteristics such as markings or size differences, we selected individuals randomly within species for trait measurements (Female:Unknown = 48:14).We aimed for a maximum of 10 individuals per species when possible (median = 10; mean = 8.1 individuals per species).

| Taxonomic and functional richness and turnover
To model the current and future distribution of taxonomic and functional composition across the island, we calculated taxonomic and functional richness and turnover for each grid cell (1 km 2 resolution) based on model outputs.Species richness was expressed as the number of species projected to occur in each grid cell.We obtained estimates of species richness by stacking and summing distribution maps.Taxonomic beta diversity accounts for differences in the numbers of species between sites (Baselga, 2010;Leprieur et al., 2012) and is partitioned into nestedness and turnover.Nestedness represents the difference in species richness between assemblages (Leprieur et al., 2012;McKnight et al., 2007), whereas turnover represents true species replacement between assemblages.The latter is insensitive to species richness gradients, such as those we expected to occur across the island's climatic gradients, and contributes more to beta diversity than nestedness (Baselga et al., 2012).We thus focused our analyses on turnover rather than nestedness and, towards this end, we used the Sørensen pairwise dissimilarity index among grid cells.A value of 0 indicates neighbouring pixels or communities had similar species composition, whereas a value of 1 indicates communities with distinct species composition.
Functional richness and turnover were calculated following a similar approach.Functional richness was calculated as the total branch length of a tree linking all projected species occurring within a grid cell, which was calculated using functions in the betapart and BAT packages (Baselga, 2010;Baselga & Orme, 2012;Cardoso et al., 2015).The tree was calculated using a hierarchical cluster analysis using trait values for species within each grid cell.Functional turnover or beta diversity was then calculated as the turnover of traits between a focal grid cell and its eight neighbouring cells using the Sørensen pairwise dissimilarity index (Sorensen, 1948;following Petchey & Gaston, 2006;Safi et al., 2011)

| Climatic-mediated effects on species and trait distributions
Our next set of analyses was aimed at examining environmental drivers, in particular temperature, on present and relative per cent change in projected butterfly taxonomic and functional richness and turnover, as well as wing trait variation across Puerto Rico.We used generalized additive models (GAMs; Hastie & Tibshirani, 1986) in which bioclimatic variables from WorldClim were used as predictors and layers for taxonomic and functional richness and turnover were used as responses.Similar to the species distribution models, multicollinearity among bioclimatic variables was removed using a VIF of less than 2 (Queen et al., 2002), resulting in four predictor variables (temperature seasonality, temperature annual range, precipitation of the warmest quarter and elevation).Because this was a pixel-scale analysis, these four variables differed from the six variables used for species distribution models (derived based on occurrence records rather than pixels).The four variables were scaled and centred to a mean of 0 and a standard deviation of 1. GAMs were fitted using the mgcv package (Wood & Wood, 2015) with a thin-plate spline smoother (which tends to have better RMSE performance) and default dimension (k) of the basis for the smoothing term.Estimation of smoothing penalties was performed using restricted maximum likelihood (method = 'REML') which has been shown to perform better than other methods and treats the model as a mixed effects model with the non-linear parts of the spline treated as special random effects.We generated all possible additive models using MuMIn (Barton & Barton, 2015) and used model selection according to the Akaike information criterion (AIC, Burnham & Anderson, 2004).
Because GAMs use smooth functions that can be complex and nonlinear, interpreting model coefficients can be complicated (Hastie & Tibshirani, 1986), making it difficult to understand how the predictor variables affect the response variables.Because of this, we also fit an ordinary least squares (OLS) regression to aid in interpreting the relative effect size and direction (i.e.negative and positive) of relationships between response and predictor variables.This dual approach can help avoid misinterpreting the importance of predictor variables.

| Trait-mediated effects on species distributions
Our third and final set of analyses was aimed at testing whether traits mediate projected distribution changes including upslope migration, leading to taxonomic or functional homogenization.To this end, we first ran a principal component analysis (PCA) on all traits (i.e.wing length, hue, saturation and brightness) to interpret a reduced number of axes of variation more easily.All variables were scaled to a mean of 0 and a standard deviation of 1 prior to analyses (z-score transformation) due to differences in magnitudes and units among variables.We calculated Spearman correlation coefficients among PCA scores and area, elevational range, maximum elevation, mean elevation and minimum elevation for projected present and future distributions (and relative per cent change) based on binary transformed ensemble model predictions.

| Taxonomic and functional richness and turnover
Modelled taxonomic and functional richness and turnover were variable across the island, both under current and future climate scenarios.Under current conditions, species richness was projected to be the highest in the southern and western lowlands, whereas species turnover was projected to be the highest in the Cordillera Central and the Sierra de Cayey (Figure 2a).In contrast, functional richness was projected to be the highest in the north-west karst region, and functional turnover was projected to be the highest in the large Caguas-Juncos eastern valley (Figure 2a).Wing length was projected to be the highest along the northern karst region where hue values were also projected to be highest in the yellow range (Figure 2b), whereas saturation and brightness were projected to be the highest in the eastern portions of the Cordillera Central.
The patterns described above exhibited important changes under future conditions.In general, our models projected an overall increase in species richness at the site level.The largest increases in species richness occurred in the northern karst region and the Cordillera Central.In contrast, our models projected a decrease in species turnover, with the greatest decreases in species turnover projected to occur along the western and southern coasts.Likewise, our models projected an increase in functional richness across the island, although the greatest increases were projected to occur in the eastern Caguas-Juncos Valley (Figure 2a).Functional turnover, like species turnover, was projected to decrease at the site level across the island, yet the greatest decreases were projected to occur in the northern karst region and the eastern Caguas-Juncos Valley (Figure 2a).Although functional turnover was projected to decrease at the site level across the island, it was projected to increase in the Sierra de Luquillo (Figure 2a).Wing traits also exhibited variable patterns across the island; however, all traits were most notably projected to change in the northern karst region or the eastern Caguas-Juncos Valley (Figure 2b).
These patterns are best summarized by examining the relative per cent change between modelled future and current taxonomic and functional richness and turnover at the island level.
The extent to which taxonomic and functional measures of diversity are linked is not well studied.Here, we found that projected taxonomic richness was positively correlated with projected functional richness (R 2 = .41,p < .001),taxonomic turnover was positively correlated with functional turnover (R 2 = .51,p < .001),but richness was negatively correlated with turnover for both taxonomic (R 2 = .38,p < .001)and functional indices (R 2 = .33,p < .001)(i.e.taxonomic richness was negatively correlated with taxonomic turnover and functional richness was negatively correlated with functional turnover; Figures S1.4 and S1.5).

| Climatic-mediated effects on species and trait distributions
Based on results from both GAMs and OLS regressions, temperature seasonality (bio4) was the strongest predictor of projected current taxonomic and functional richness and turnover, although with opposing effects (Table 1).Temperature seasonality was negatively associated with projected measures of richness (based on OLS regressions; species: β = −5.61,p < .001;functional: β = −0.29,p < .001)and positively associated with those of turnover (species: β = −0.06,p < .001;functional: β = 0.04, p < .001).Similarly, temperature seasonality was the strongest predictor of projected relative per cent change in species richness (β = 3.49, p < .001),functional richness (β = 0.45, p < .001)and functional turnover (β = −0.37,p < .001),but not projected changes in species turnover (Table 1).Temperature seasonality was positively associated with projected change in species (β = 3.49, p < .001)and functional richness (β = 0.45, p < .001),and negatively associated with projected change in functional turnover (β = −0.37,p < .001).Projected change in species turnover was more strongly (and positively) , future (2061-2080) and per cent change (%) of (a) species richness (number of species), functional richness (sum of the trait-distance matrix branch lengths), species and functional turnover (Sørensen pairwise dissimilarity index) and (b) wing length (cm), hue, saturation and brightness (colour metrics all unitless, range 0-255) for 62 butterfly species using stacked ensemble distribution models.Note the different units and scales for each panel.associated with elevation (β = 0.03, p < .001).Under current conditions, predictors of wing traits varied by trait (Table 1).In contrast, temperature seasonality was the strongest predictor of projected change in all wing trait values, but the direction of the relationship varied from negative (length and saturation) to positive (hue and brightness, Table 1).

| Trait-mediated effects on species distributions
Species were primarily separated by wing hue, saturation and brightness along Axis 1 of the PCA and secondarily by wing length along Axis 2 (Figure S1.6, Table S1.2).Spearman correlation coefficients among PC scores indicated significant relationships among the second PC axis (strongly associated with wing length), current distribution area (rho: .26,p < .05)and modelled future elevation including modelled future maximum elevation and change in elevational range (Tables S1.3 and S1.4).Wing length was negatively correlated with projected changes in maximum elevation (R 2 = −.34,p < .01)and elevational range (R 2 = −.34,p < .01).Wing saturation was positively correlated with wing length (R 2 = .28,p < .01)and negatively correlated with hue (R 2 = −.49,p < .0001;Table S1.4).Notably, largerwinged species were projected to exhibit decreased distribution area with declines in maximum elevation, whereas smaller-winged species were projected to exhibit increased distribution area with upslope expansion (Figure 4).

| DISCUSS ION
We used digitized museum collections and species distribution models to model the extent to which climate change could lead to taxonomic and functional homogenization of tropical butterfly communities.The modelling approach used here pointed to both taxonomic and functional homogenization, as well as a negative relationship between wing length and temperature seasonality.Finally, we found that changes in modelled distribution area were negatively correlated with wing length.Altogether, our work indicates that temperature seasonality is a prominent driver of biotic homogenization TA B L E 1 Results of generalized additive models (GAM) and ordinary least squares regression (OLS) of projected present and projected change in response variables: taxonomic and functional richness and turnover, and wing length, hue, saturation and brightness and explanatory variables: elevation, precipitation of the warmest quarter (bio18), temperature seasonality (bio4) and temperature annual range (bio7).Note: For brevity, degrees of freedom (df) and, as measures of overall fit, the adjusted R 2 and the deviance explained are presented for GAMs (see Table S1.1 for F-values for each variable).Estimates of the regression coefficients are presented for OLS regressions.All coefficients were significant (p < .001)unless otherwise noted ('ns').Bold values indicate non-significance (p > .05).
by favouring smaller, less melanized species that expand their distributions upslope.

| Taxonomic and functional richness and turnover
The first goal of this study was to model current and projected trends in taxonomic and functional richness and turnover across the topographically complex island of Puerto Rico.Under current climate conditions, butterfly species richness was projected to be the highest along the south-west coast, mirroring patterns of bird diversity (Acevedo & Restrepo, 2008;Campos-Cerqueira et al., 2021;Kepler & Kepler, 1970;Smith et al., 2017).A combination of high precipitation seasonality and reduced abundance of amphibians competing with birds for invertebrate prey is suggested as an explanation for high bird richness in this region (Kepler & Kepler, 1970).In comparison, butterfly species turnover was projected to be the highest across two main mountain ranges, namely the Cordillera Central and the Sierra de Cayey.High species turnover in mountainous terrain is not uncommon (e.g.Dewan et al., 2022) and can be attributed to environmental heterogeneity (Wilson & Fox, 2021).Functional richness was projected to be the greatest in the humid, north-west region of the island, where taxonomic richness was also projected to be high.This region is characterized by limestone hills or haystacks (mogotes), consisting of narrow valleys surrounded by hills with steep slopes (Lugo et al., 2001).In general, high biodiversity in this region (including endemic birds, amphibians and reptiles) is attributed to protection from strong winds and reduced agricultural impact on steep terrain (Rivera & Aide, 1998).These conditions have led to persistent forest cover which may insulate populations from disturbances and this stability may, over long periods of time, result in greater diversity (e.g.Roubik et al., 2021;Willig et al., 2019).Finally, functional turnover was projected to be the greatest in the humid and fertile Caguas-Juncos valley, an alluvial valley once used for cattle and sugarcane.As in other tropical regions (e.g.Bush et al., 2015), land-use legacies influence contemporary patterns of plant species composition (Pascarella et al., 2000) and thus suitable host plant distributions, which likely play a role in the high butterfly functional turnover found in the Caguas-Juncos valley.In Puerto Rico, the legacy of deforestation and land use continues to impact plant species richness and distributions (Grau et al., 2003) with cascading effects on insects that depend or specialize on different plants (Acevedo & Restrepo, 2008;Barberena-Arias & Aide, 2003).For example, one study found that although secondary forests in Puerto Rico became structurally similar to old-growth forests, they differed in plant species composition, often dominated by non-indigenous species (Aide et al., 2000;Zimmerman et al., 2000).(Brandeis & Turner, 2013;Rojas-Sandoval et al., 2022).It is possible that changes in forest cover and plant abundances during this time led to shifts in butterfly abundances from early successional species specializing in agricultural crops to later successional species, supporting a more generalist butterfly fauna given the increase in non-indigenous plant species composition.Furthermore, forest regeneration combined with frequent hurricane disturbance leads to complex successional processes (Flynn et al., 2010) -2000).A related caveat is that the butterfly species included in our modelling approach encompasses a bit more than half of the total estimated number of butterfly species on the island (see Figure S1.7).The included butterfly species reflect those with at least 10 occurrence records and tend to be more common.Thus, understanding whether rarity in the collection record translates to rarity in abundance could better determine whether our results underestimate patterns of richness and turnover.
Compared to diversity patterns of current modelled butterfly communities, future butterfly communities were projected to become more taxonomically and functionally similar or homogeneous.Biotic homogenization was projected to occur throughout the island but especially in the coastal lowlands and interior valleys.These results are in stark contrast with those for the flora of Puerto Rico where a greater floristic homogenization was projected to occur at higher elevations (Henareh Khalyani et al., 2019), which may be explained by differences in the bioclimatic variables used in the SDMs.Henareh Khalyani et al. (2019) used annual means of temperature and precipitation, whereas we used metrics describing their variation (seasonality, range), extremes (water availability during the warmest quarter) and elevation (a major determinant of temperature on the island; Daly et al., 2003).Shifting climatic variables have interacting effects on insect growth and mortality (Pureswaran et al., 2018), so more detailed studies on how climate means and variability interact is a pressing challenge.Another possibility is that plants and insects have different abiotic limitations, and higher elevation habitats are currently too cold or climatically variable for the butterfly species studied here, pointing to potential future spatial mismatches between insects and their host plants under future climate scenarios (e.g.Gérard et al., 2020;Schweiger et al., 2008).Alternatively, these contrasting findings may be driven by differences in plant and butterfly diversity patterns or differences in the degree of host plant specialization.Lowland habitats are warmer, supporting a diverse butterfly community but reduced plant diversity (Weaver & Murphy, 1990), perhaps because the butterfly assemblage there consists of more generalists.In comparison, upland habitats harbour fewer butterfly species even though these areas are floristically more diverse (Gould et al., 2006;Weaver & Murphy, 1990), pointing to a more specialized butterfly assemblage.
Given that species distributions reflect a range of variables in addition to climate, actual species re-distributions will depend on life histories, behaviour, dispersal and concomitant shifts in other species including host plants, predators, parasitoids and pathogens as well as changes in land use.For example, behavioural thermoregulation can diminish the direct impacts of warming and facilitate species persistence in unfavourable thermal environments (Bonebrake et al., 2014;Wenda et al., 2021).Additionally, in this study, we estimated adult butterfly responses but caterpillar responses may have equally important consequences for population persistence in warming climates (e.g.Agosta et al., 2017;Fenberg et al., 2016;Van Dyck et al., 2015;Wilson et al., 2019).Finally, the results shown here may be influenced by evolutionary relationships with certain families exhibiting greater sensitivity to changes in climatic variability or, conversely, exhibiting greater sensitivity to the presence of particular host plants (which differ across environments), confounding the observed linkages between traits and climate.

| Climatic-mediated effects on species and trait distributions
The second goal of our study was to understand possible climaticmediated effects on taxonomic and functional richness and turnover and wing traits.Overall, we found that current and future-modelled measures of richness and turnover were driven primarily by temperature seasonality followed by elevation.In addition, we found that under current conditions, wing size and melanization increased with elevation.In other words, butterflies were larger and more melanized at higher elevations.However, future climatic conditions favoured smaller, less melanized species, particularly in the high Central Cordillera.In other words, butterfly assemblages became smaller and lighter at higher elevations.Our results are in line with three sets of studies.The first reported larger wing size and greater wing melanization at higher elevations (Brehm et al., 2019;Henriques et al., 2022;Leingärtner et al., 2014;Xing et al., 2016), both among temperate and tropical Lepidoptera.The second reported declining body size across temperate (Bowden et al., 2015;Wonglersak et al., 2020;Zeuss et al., 2014) and tropical (Wu et al., 2019) insects as a universal response to climate warming (Gardner et al., 2011;Ohlberger, 2013).
Although the latter study set primarily focuses on within-species changes in body size using museum specimens and historical temperature records, we demonstrate these changes can also occur at an assemblage level.Taken together, the emerging importance of thermoregulatory constraints on wing size and colour (e.g.Clusella Trullas et al., 2007;Kingsolver, 1983) suggests that future warming is likely to favour smaller, less melanized species.The third set of studies demonstrates the prominent, and often overlooked, role of climatic variability on species and functional diversity.In our study, current modelled taxonomic and functional richness and turnover, as well as projected changes in these and wing traits, were primarily driven by temperature seasonality: modelled species and functional richness increased with decreasing temperature seasonality.
Thus, thermal stability promotes species and functional diversity.
In line with other studies, thermal stability was strongly correlated with species richness in our study, possibly reflecting reduced extinction risks (Huntley et al., 2016).However, our results also showed greater trait variation among species in thermally stable environments, which contradicts empirical results across latitude and taxa (e.g.Deutsch et al., 2008) and predictions from the climatic variability hypothesis (e.g.Janzen, 1967) which argues that trait variation should be greater in thermally variable environments.One explanation for this result is that climatic stability favours specialization (Vázquez & Stevens, 2004) and, consequently, greater specialization results in finer niche partitioning, allowing more species to co-exist (MacArthur, 1972).Another possibility is that thermal stability reduces selection for traits that cope with unstable temperature regimes and these traits then vary simply by genetic drift.The impact of temperature seasonality on species turnover is much less understood, and on functional turnover even less so (Dewan et al., 2022).Seasonality may impact species and functional turnover through its impact on resource availability (Gaston & Chown, 1999).Indeed, the seasonality of resource availability was a major influence on weevil body size variation (Chown & Klok, 2003) and may be especially important for the Puerto Rican Lepidoptera assemblage given that their abundances vary between wet and dry seasons (Davila & Asencio, 1992).

| Trait-mediated effects on species distributions
The final question we asked was whether traits mediate projected changes in species distributions.Our comparisons of wing colour and size with species' projected changes in distributions showed that shifts in species distributions were associated with wing size and, to a lesser extent, wing brightness.Specifically, smallerwinged species were projected to exhibit the greatest increases in distribution area, that is, expansion of their ranges, primarily through upslope migration.In contrast, larger-winged species were projected to exhibit the greatest declines in distribution area, that is, contraction of their ranges, via downslope shifts in distribution.Examining the degree of specialization among small and larger-winged species could provide additional insight into these differences, for example, if smaller-winged species tend to be wide-ranging generalists.A phylogenetic approach could further determine whether these results are driven by specific families (e.g.small-winged Lycaenids).Interestingly, the widest-ranging species exhibited the poorest model performance.Many of these were skippers (Hesperiidae), pointing to mechanisms that may promote broad ranges or determine population-level responses to changing climate (see Figure S1.8).Furthermore, our results are in sharp contrast with a recent meta-analysis across latitudes that showed that traits, including body size, are weak predictors of distribution shifts across taxa (Beissinger & Riddell, 2021).Perhaps it is easier to detect the role of traits at smaller scales rather than the global scale of meta-analyses.The present study was restricted to a small topographically complex island, and the modelled distribution shifts occurred across elevation rather than latitude.The meta-analysis by Beissinger and Riddell (2021) also contradicts other studies on Lepidoptera that report body size as an important predictor of distribution shifts.For example, large body size predisposed species to distribution declines among butterflies in Finland (Mattila et al., 2011) and population declines among moths in England (Coulthard et al., 2019).Thus, for mobile taxa such as Lepidoptera, traits such as wing size may be good predictors of climate change-induced shifts in species distributions, as we have modelled here.Even though large body size is positively associated with dispersal ability (MacLean & Beissinger, 2017), and thus might be expected to promote distribution expansions, selection against larger wing size may be driven by thermoregulation.Warming temperatures are expected to favour smaller wing sizes, reducing thermal loads and metabolic demands.Yet the downslope shifts of larger-winged species from cooler uplands to warmer lowlands, projected here, seemingly challenge a thermoregulatory explanation.These results may point to interactive effects of wing size on temperature tolerance and desiccation resistance (Chown et al., 2011), with larger wings conferring greater drought tolerance (Hooper et al., 1999;Kellermann & van Heerwaarden, 2019;Parkash et al., 2009).Modelled distribution shifts were also explained by wing brightness, although to a lesser extent than wing size.Colour may translate into differences in heat loads and exposure to climatic conditions, which further underscores the idea that thermoregulatory constraints drive predicted species distributions in the butterfly fauna studied here.Here, we focused on species-level traits but intraspecific variation in wing size or colouration (e.g.Ellers & Boggs, 2004;Kingsolver, 1995) would be worth examining further (and the images are available for doing so).Local adaptation or plasticity could be important mechanisms for responding to climate change in situ.

| CON CLUS ION
Using digitized museum specimen records, we generated species and trait distribution models to examine the potential for taxonomic and functional homogenization among a tropical butterfly assemblage.Our modelling approach points to biotic homogenization resulting from shifts in species and trait distributions.Results also consistently showed that island-wide shifts in modelled species distributions were primarily driven by temperature seasonality and mediated by wing size.These projected relationships are tempered by the recognition that distribution shifts do not occur in isolation , single sites (e.g.Davila & F I G U R E 1 Elevation (m) in Puerto Rico and major physiographic regions including the northern karst formations (mogotes), Sierra de Luquillo, Cajuas-Juncos (CJ) Valley, Cordillera Central and the Sierra de Cayey.
using the speciesRaster package (Title 2017).As with species turnover, a value of 0 indicates that neighbouring pixels or communities have similar functional trait composition, and a value of 1 indicates that communities have distinct (i.e.no overlap in their) functional composition.The relative per cent change (RPC) between projected current and future taxonomic and functional richness and turnover, distribution area, elevational range and maximum elevation was calculated for each pixel according to Equation 1: where PF i and PP i are the projected future and projected present values of pixel i respectively.Pixel values were plotted as histograms to visualize pixel-level gains and losses in taxonomic and functional richness and turnover.We used a two-sided z-test (H 0 : μ = 0) or a Mann-Whitney U-test (for non-normal distributions) to identify statistically significant gains or losses in pixel-level taxonomic and functional richness and turnover and wing trait values.

F
I G U R E 3 Density histogram plots showing changes in modelled (a-d) species and functional richness and turnover; (e-h) wing length, hue, saturation and brightness over time.Values represent pixel-level values calculated based on mean traits for all species with at least 65% probability of occurring in each cell.Pixels are ~1 km 2 spatial resolution (883 × 923 m), and (i) total area (km 2 ) exhibiting gains or losses in each response variable.
Although occurrence records span a range of collection years from 1830 to 2019, most were collected between 1980 and 2000.During this 20-year F I G U R E 4 Local regression between modelled change in area and modelled change in maximum elevation for 62 butterfly species based on current (1970-2000) and future (2061-2080) environmental conditions, showing that predicted declines in maximum elevation are not associated with changes in distribution area.In contrast, predicted increases in maximum elevation are associated with predicted increases in distribution area (R 2 = .35,p < .001).Histogram and density distributions for predicted change in distribution area (km 2 ) and predicted change in maximum elevation (m) are shown along the corresponding axes.Most species with small wings were predicted to expand their habitat area and shift to higher elevations.Size and colour of data points indicate wing length (cm).Notably, species which were predicted to exhibit declines in maximum elevation were predominantly larger winged.period, widespread agricultural abandonment led to forest recovery, and by the early 2000s, forest cover in Puerto Rico reached 55%, even though these novel forests are dominated by non-indigenous plants at the expense of indigenous ones