Climate-driven declines in arthropod abundance restructure a rainforest food web

Fig. 1 compares trends in mean maximum yearly temperature (MnMaxT) for the El Verde Field Station and Bisley Tower meteorological stations, both located at an altitude of 350 m. Between 1978 and 2015, MnMaxT at El Verde rose by 2.0 °C at an average rate of 0.050 °C per year. Between 1993 and 2015, the rate of temperature increase at Bisley Tower was 0.055 °C per year, not significantly different from the rate at El Verde.

As several authors have pointed out, increased exposure to extreme temperatures may have a greater impact on fitness than gradual increases in average temperatures (6, 19, 20). At El Verde the proportion of maximum daily temperatures equal to or exceeding 29.0 °C increased significantly between 1978 and 2015 (SI Appendix, Fig. S8E). At the Bisley station, the proportion of maximum daily temperatures equal to or exceeding 29.0 °C grew steadily from the 1990s until 2015, with the proportion averaging 0.03 for 1993–1994 and 0.44 for 2014–2015. While average temperatures increased 8–10% over the same period, the proportion of extreme temperatures increased 14 times. Because data from the Bisley and El Verde stations record just one maximum and minimum temperature per 24-h period, we could not calculate the duration of time that temperatures remained above the 29 °C extreme. However, there was a significant positive linear regression of proportion of days per year when the temperature equaled or exceeded 29 °C on MnMaxT at El Verde (SI Appendix, Fig. S8K), suggesting that our measure reflects significant exposure to extremes rather than ephemeral spikes in temperature (6, 20).

From 1975 to 2015, rainfall at the El Verde station averaged 3,653 mm/y with no significant positive or negative trends. However, the number of dry days (rainfall less than 1 mm) has been decreasing since 1975 (SI Appendix, Fig. S8G), while the variation in daily rainfall has been increasing (SI Appendix, Fig. S8H). Over the 36-y time span encompassed by the data analyzed here, there have been five major hurricanes and eight severe droughts with two outlier dry years, 1994 and 2015, when total rainfall was 1,404 mm and 2,035 mm, respectively (21). These disturbances are all associated with La Niña and El Niño events (SI Appendix, Fig. S9).

Fig. 2 compares arthropod dry weight biomass captured per 100 sweeps for July 1976 and January 1977 in the Luquillo forest study area (22), with our samples from July 2011 and 2012, and January 2012 and 2013. Compared with the earlier samples, arthropod biomass for July 2011–2012 and January 2012–2013 fell four and eight times, respectively. Fig. 3 shows that the declines in biomass occurred across all 10 of the major taxa captured in the sweep samples (Fig. 3 A and B vs. Fig. 3 C–E and F, respectively). A two-sample Kolmogorov–Smirnov test indicated that the distributions of biomass proportions across taxa were not significantly different between the two sample periods (Kolmogorov–Smirnov z = 0.671, P = 0.759), suggesting that all taxa declined proportionately.

Sticky-trap samples for the ground (Fig. 4A) and canopy (Fig. 4B) were indicative of a collapse in forest arthropods. The catch rate for the ground traps fell 36 times, from 473 mg per trap per day in July 1976 to 13 mg per trap per day in July 2012, and approximately 60 times, from 470 mg per trap per day to 8 mg per trap per day, between January 1977 and January 2013. In July 1976, the canopy trap catch rate was 37 mg per trap per day compared with 5 mg per trap per day in July 2012, and 21 mg per trap per day in January 1976 vs. 8 mg per trap per day in January 2013.

Analysis of Schowalter’s (23) canopy data also revealed significant decreases in invertebrate abundance, scaled with respect to foliage weight sampled between 1990 and 2010 (Fig. 5A), and a significant, nonlinear decline with increasing temperature (Fig. 5B). In keeping with our sweep net results, all 10 of the most-abundant canopy taxa exhibited negative trends in abundance between 1990 and 2010. The binomial probability of this pattern occurring by chance alone is 0.00097. Employing Willig et al.’s (24) census data, we also conducted quasi-Poisson regressions of walking stick abundance against time and against MnMaxT during the census periods (Fig. 5 C and D). Both regressions were negative and significant.

During a previous study (25), we sampled arthropods within a tropical dry forest in the Chamela-Cuixmala Biosphere Reserve in western Mexico using the same sticky-trap methods employed at Luquillo. Since 1981, the Chamela-Cuixmala Biosphere Reserve has experienced a 2.4 °C increase in average annual air temperature, from 24.5 °C in 1981 to 26.9 °C in 2014 (SI Appendix, Fig. S1). We returned to the Chamela forest in August 2014 to repeat our original samples within the same study area. Overall, we found an eightfold decrease in the dry-weight biomass of arthropods captured per day in ground traps in 2014 compared with the average for July 1987 and July 1988 (SI Appendix, Fig. S2).

Given these declines in arthropod resources, parallel declines in insectivore abundances would be expected. The dominant insect predators within the Luquillo rainforest are anoline lizards, eleutherodactylid frogs, and birds (26). Here we analyze changes in abundance through time for forest-dwelling anoles, the most common forest frog, Eleutherodactylus coqui, and the bird community near the El Verde Field Station.

Within our Luquillo study area, the most common anoles during 1976–1977 were Anolis gundlachi, Anolis evermanni, and Anolis stratulus (22). In 2011 and 2012, we censused anoles within the same 225-m² study area using the Schnabel method (27) employed in 1976/1977 by Lister (22). Average density of anoles in 2011–2012 had decreased ∼2.2 times compared with 1976–1977, from an estimated 12,490–5,640 anoles per hectare. Mean anole biomass within the study area dropped 32%, from 401 g to 277 g (Fig. 6 and SI Appendix, Table S1). The density of A. gundlachi, the most common species, had dropped twice, and the density of A. evermanni declined 10 times. A. stratulus, with a density in the 1970s of 200 individuals per hectare, had disappeared from the interior forest, edge habitats, and other regularly visited habitats lying within 4 km to the east and west of the study area.

We conducted quasi-Poisson regressions of E. coqui counts vs. time for censuses conducted at El Verde (Fig. 7 A and B) by Stewart (28), and at El Verde (SI Appendix, Fig. S3 A and B) and the Bisley watershed (SI Appendix, Fig. S3 C and D) by Woolbright (29, 30). Data from Woolbright’s El Verde study area were previously analyzed to investigate the impact of Hurricane Hugo on Luquillo frogs (29, 30). Here we present an analysis and interpretation. To our knowledge, no analyses of the Bisley data have been previously published.

The quasi-Poisson regressions of E. coqui numbers vs. time and MnMaxT for Stewart’s Activity Transect were both negative and significant (Fig. 7 A and B) (28). The regressions of E. coqui numbers vs. time and temperature for the Sonadura West study area were also negative and significant (SI Appendix, Fig. S3A). For the Sonadura Old population, only the regression of numbers on time was significant (P = 0.013), while the regression of abundance on temperature was marginally insignificant (P = 0.067). Because no climate data were taken at the Bisley Tower until 1993, we were unable to conduct regression analyses of E. coqui numbers vs. temperature for the Bisley Stream and Bisley forest study areas. Regressions on time were both negative and significant (SI Appendix, Fig. S3 C and D).

The temporal pattern of E. coqui abundances in Woolbright’s (29, 30) El Verde study area (SI Appendix, Fig. S3B) departs from the monotonic, negative trends found in the other E. coqui populations, and required a segmented regression consisting of three phases. In the first phase, E. coqui numbers had a significant negative trend followed by a rapid increase during phase 2 from 6 to 57 individuals between January and October 1990, most likely due to immigration from surrounding areas following Hurricane Hugo. In phase 3, the population resumed a steady decline at about the same rate as phase 1, until the census ended 6 y later. The rapid return to the preperturbation rate of decline suggests the ongoing influence of climate change on birth and death rates.

To explore trends in avian abundance, we analyzed mist netting data taken near the El Verde Field Station by Waide (31) from 1990 to 2015, an analysis that had not been previously published. As Fig. 7C shows, under the same yearly sampling effort there was a significant decline in total insectivorous birds captured per year between 1990 and 2005, with captures falling 53% from a high of 137 birds in 1990 to 64 birds in 2005. During this period, seven of the eight most common bird species at El Verde experienced reductions in numbers. The quasi-Poisson regression of insectivorous birds caught per year on MnMaxT during the period of sampling was also negative and significant (Fig. 7D).

If the declines in arthropod resources negatively impact avian abundance, then species with primarily granivorous or frugivorous diets should suffer smaller declines than those consuming a greater proportion of insects. For example, the ruddy quail dove (Geotrygon montana), which eats grains and fruits almost exclusively, showed no decline in abundance between 1990 and 2015 at El Verde, while Todus mexicanus, the Puerto Rican tody, whose diet is composed entirely of insects (26), experienced a 90% reduction in catch rate. Using data on diets of Puerto Rican birds (26) and a Theil-Sen nonparametric regression, we found that for six common species the relative decline in catch rate compared with G. montana increased significantly with an increase in the proportion of insects in their diet (SI Appendix, Fig. S4).

To go beyond univariate regressions of population abundance on ambient temperature, we constructed multiple quasi-Poisson regressions using independent variables identified by a hierarchical partitioning (HP) analysis (32, 33). Given a set of predictors, HP averages the effect of each predictor on goodness of fit measures over all possible 2^K regression models, where “K” is the number of predictors (34, 35). This averaging alleviates the effects of multicollinearity and allows the explanatory impact of a given variable to be partitioned into an effect that depends only on the predictor itself and a joint effect shared with one or more other predictors. For each predictor, independent and joint effects are expressed as a percent of the total deviance explained. Variables with large independent effects are more likely to be causal than variables with small independent effects (33, 34). Hence, the technique is not designed to build an optimal predictive model, but to complement the model building process by revealing potential explanatory/causal relationships between a dependent variable and a set of predictors.

Using the nine climate-related variables listed in SI Appendix, Table S3 as independent variables, and abundance as the dependent variable, we conducted HP analyses for the three E. coqui populations, the canopy arthropods, walking sticks, and the El Verde birds. Results are shown in SI Appendix, Fig. S5. MnMaxT had considerably greater independent explanatory power than any of the other variables. Across all six populations, the average independent effect of MnMaxT on abundance was 25%, while the averages for four other temperature related variables ranged from 11.4 to 14.3%. The independent explanatory power of average daily rainfall was highest in the E. coqui and walking stick populations (9.8–24.4%). The Southern Oscillation Index (SOI) also emerged as an impactful, potentially causal factor with independent effects ranging from 10.2 to 17.3% in two of the three E. coqui populations, the walking stick population and the canopy arthropods.

The final multiple quasi-Poisson regressions, with supporting statistics, are given in SI Appendix, Table S5. For the El Verde birds, none of the variables in SI Appendix, Table S3 were significant once MnMaxT was entered into the model. Notably, the multivariate models for the three E. coqui populations and the walking stick population all contained MnMaxT and PropGT29 (proportion of days per year when daily temperature exceeds 29 °C) as major predictors. Multiple regressions for the El Verde E. coqui and walking stick populations also included the SOI as a significant predictor of abundance (SI Appendix, Table S5).

To further analyze the effect of El Niño/Southern Oscillation (ENSO) on arthropod populations, we conducted cross-correlation analyses between the SOI and the abundances of canopy invertebrates and walking sticks (SI Appendix, Fig. S6). For the canopy invertebrates, there is a significant (P < 0.05) positive cross-correlation with the SOI at a time lag of 0. This indicates that positive SOI values (La Niña episodes) have an immediate positive impact on the abundance of canopy invertebrates, and negative SOI values (El Niño episodes) have an immediate negative impact. There is also a significant negative cross-correlation at a negative 5-y time lag, indicating that higher abundances precede El Niño events and lower abundances precede La Niña events, as would be expected. For walking sticks, there is also a significant positive cross-correlation between abundance and the SOI, again indicating that higher and lower values of the SOI lead to higher and lower abundances respectively. However, for walking sticks, the ENSO events precede the impact on population abundance by approximately 2 y.

Regression and correlation analyses inevitably leave the issue of causality open to question (33, 34). In lieu of intractable experiments spanning large spatial and temporal scales, Vinod causality offers a statistical method that can help identify causal variables and reveal the direction of causal paths (36). Vinod’s generalized correlation and kernel causality methods (37, 38) improve and extend earlier attempts to discern causal paths between two or more variables, including those of Granger (39), Pearl (40), and Hoyer et al. (41). With a success rate of 70–75%, Vinod’s methods are the most accurate to date in terms of correctly identifying causality in the CauseEffectPairs benchmark database (42). Here we explore the evidence that temperature is the causal agent affecting abundance using Vinod’s generalCorr R package.

Given an xy dataset, generalCorr calculates a unanimity index (UI) that quantifies the likelihood that either x or y is causal. This index will always lie in the range [−100, 100]. If all three of Vinod’s causality criteria are met, then the UI will have an absolute value of 100, signifying the strongest evidence for the causality of a given variable. Three decision rules based on the value of the UI determine the direction of the causal path. If the UI lies in the interval [−100, −15], then y causes x, and if UI is in the interval [15, 100] then x causes y. If the UI lies within the range [−15, 15], the causal direction is indeterminate.

SI Appendix, Table S4 summarizes the outcome of the generalCorr analyses of causation. The unanimity indices equaled 100 in five of the six groups, indicating that the causal path Temperature → Abundance is supported by all three of Vinod’s criteria for causality. Inclusion of mean daily rainfall as a covariate improved the UI for all analyses except for walking sticks, which had a UI of 100 without the covariate, and for the E. coqui at the Sonadura Old study site, which retained a low UI and indeterminate causal path with and without the covariate.

Based on 14 y of light-trap data within a Panamanian rainforest, Wolda (43) observed that the stability of tropical insects varied considerably. He concluded that this result reflected the absence of a “dominant environmental process” affecting many species simultaneously. We suggest that climate warming has emerged as such a dominant force affecting arthropod populations in the Luquillo forest and, most likely, in other tropical forests experiencing significant increases in ambient temperature.

Several lines of evidence support this hypothesis. First, from a single species of walking stick to a community of over 120 canopy taxa, long-term declines have occurred in parallel with rising temperatures. Average ambient temperature is also a significant predictor of changes in the abundance of canopy invertebrates and walking sticks. Second, the simultaneous declines of all 10 arthropod taxa in our sweep samples, the 10 most-abundant taxa in Schowalter’s (23) canopy samples, and arthropods occupying all strata of our forest study area, all point to an overriding environmental factor that has had ubiquitous, adverse effects on forest arthropods regardless of taxonomic affiliation, stratum occupied, or type of niche exploited. Third, declines in arthropod abundance have occurred despite major decreases in their predators. This “counter cascade” conflicts with decades of research on trophic dynamics (44), as well as studies documenting the major impact of anoles, frogs, and birds on arthropod abundance (45, 46). Stewart and Woolbright (47), for example, estimated the E. coqui population alone harvested 114,000 prey items per hectare per night in the Luquillo forest, and Reagan and Waide (26) calculated that A. stratulus, A. gundlachi, and A. evermanni collectively consumed ∼443,300 invertebrates per hectare per day at El Verde. If climate warming is not driving arthropod declines, then whatever the factors may be, they must completely override the reduction in mortality from lowered predation levels. Fourth, our analyses are in keeping with predictions of Deutsch et al. (5) that tropical insects should be especially vulnerable to climate warming, adding to field and laboratory studies over the past several years that have consistently supported the predicted sensitivity (8, 48–53). Finally, our arthropod samples from the Chamela-Cuixmala reserve in Mexico add support to our hypothesis that climate warming is having a serious negative impact on invertebrate populations in tropical forests.

Other studies have also documented declines in arthropod abundance across a broad range of geography, habitats, and taxa (10–12, 54, 55, 56–58). Research on causal factors has focused on anthropogenic disturbance and pesticides (57, 58). Given its long-term protected status (59), significant human perturbations have been virtually nonexistent within the Luquillo forest since the 1930s, and thus are an unlikely source of invertebrate declines. Due to the ongoing reduction in agriculture and associated farmland, pesticides use in Puerto Rico also fell up to 80% between 1969 and 2012 (SI Appendix, Fig. S7). Most pesticides have half-lives measured in days, not decades (60), making it improbable that, despite precipitous declines in their use, remaining residues are responsible for waning arthropod abundance. Recently, Hallman et al. (61) published an extensive analysis of long-term data on the biomass of flying insects sampled at 63 protected reserves in Germany. Overall, they found a 75% decline in biomass between 1989 and 2016. Their statistical models largely eliminated land-use change as an explanatory variable but did not thoroughly analyze a number of climate-change variables. Increases in average temperature had a positive influence on insect biomass, an effect predicted by Tewksbury et al. (62) for temperate insects. However, the reasons for the alarming losses in insect biomass are currently unknown and, of course, may differ from causal agents in tropical habitats.

Long-term abundance data taken by Pounds et al. (63), Whitfield et al. (4), and Stapley et al. (20) indicate that Anolis populations in other tropical forests have declined in parallel with anoles at Luquillo. At Monteverde, Costa Rica, populations of Anolis tropidolepis and Anolis altae were diminishing by 4.4% and 3.1% per year, respectively, for 11 and 13 y before both species became locally extinct. Over a 35-y period at La Selva, Costa Rica, Anolis capito, Anolis humilis, and Anolis apletophallus decreased by 23%, 78%, and 71%, respectively (3), and in the rainforest on Barro Colorado Island, Panama, A. apletophallus abundance declined by 71% between 1971 and 2011 (19). Proposed drivers of these trends were reduced mist frequency at Monteverde, diminished leaf litter due to climate change at La Selva, and the effects of the ENSO combined with a rise in minimum temperatures on Barro Colorado.

Estimated rates of decline for Anolis populations at Luquillo, Monteverde, La Selva, and Barro Colorado are given in SI Appendix, Table S2. All of these species are nonbasking, thermal-conforming forest inhabitants (6). In a macrophysiological and focal species analysis of Puerto Rican anoles and other lizards from 12 neotropical locations, Huey et al. (6) proposed that thermal-conforming lizards inhabiting shaded tropical forests should be especially vulnerable to climate warming. This vulnerability results from lower critical thermal maxima than basking species, limited acclimation via thermoregulation, and reduced thermal safety margins. As ambient temperatures increase, thermal stress can lead to extinction by lowering performance, restricting activity time, and lessening food intake. In an encompassing study spanning local extinctions of lizard populations on five continents, Sinervo et al. (7) reached broadly similar conclusions for baskers and thermal-conformers in temperate and tropical habitats.

In the Huey et al. (6) and Sinervo et al. (7) hypotheses, rising ambient temperatures indirectly lower food intake via changes in thermoregulatory behavior, while we propose that reduced consumption and declines in abundance result directly from climate-driven declines in food abundance. Although the mechanisms differ, the results are the same. As climate warms, lizards should increasingly become nutritionally challenged, resulting in reduced body condition, growth rates, egg production, fat body mass, physical performance, and immune system response. Consequently, death rates from starvation, predation, and disease should increase along with concomitant decreases in reproduction.

Decreases in Puerto Rico’s Eleutherodactylus populations began in the 1970s, with three species classified as threatened by 1985 (64). Stewart’s data (28), documenting the collapse of E. coqui, was also an early warning sign of negative trends in Puerto Rican frogs. Currently, 3 of the original 16 species of Puerto Rican Eleutherodactylus are considered extinct, and 13 species are listed as critically endangered, endangered, vulnerable, or near threatened (65). The monolithic pattern of rapid declines of E. coqui populations shown in Fig. 7A and SI Appendix, Fig. S3, when extrapolated, suggest their extinctions may be imminent by 2030–2040.

The ongoing losses of anuran biodiversity on Puerto Rico and other West Indian islands (66, 67) are part of a global amphibian decline that has continued largely unabated since at least the 1970s (68–72). To our knowledge, no studies have explored reductions in arthropod prey as a possible causal mechanism. Most research has focused on the fungus Batrachochytrium dendrobatidis, which has precipitated mass-mortality events and catastrophic reductions in over 200 frog species world-wide (63, 73, 74). While B. dendrobatidis has been detected in E. coqui populations at higher elevations in Puerto Rico (66, 75), the fungus is restricted to areas above 600 m (76) and, at present, there is no evidence that it has impacted the lower elevation E. coqui populations analyzed in this study.

Pounds and coworkers (63, 73, 77–79) have presented a wide-ranging body of evidence indicating that climate warming is a central factor promoting outbreaks of Batrachochytrium and the mass die-offs that follow. The relationship between chytrid outbreaks and climate is complex, but in general epidemics are often associated with repeated warm years and drier weather. Given ongoing climate change, the elevational range of Batrachochytrium will almost certainly change, as well as its impact on Puerto Rican frogs, especially juveniles (80). The key to understanding future impacts lies at the intersection of a changing climate with the specific temperature- and moisture-related physiological limitations of the fungus, and how these constraints in turn affect its reproduction (81, 82).

Based on our literature searches, no studies have explored reductions in arthropod prey as a possible causal mechanism behind amphibian declines, nor the possibility that malnutrition could interact with rising temperature and drought, compounding the impacts of disease through depression of immunocompetence (83).

Bird populations in many tropical and temperate ecosystems have been declining in parallel with those in Puerto Rico. In Canada, 44% of 460 bird species have had negative trends since 1970, and 42 species common to the United States, Canada, and Mexico have lost 50% or more of their populations in 40 y (84). Over a 20-y period, 54% of bird species native to Britain have also decreased (55). In the rainforest at La Selva, Costa Rica, 10–20% of resident birds declined between 1960 and 1999 (85), and in the Guanica dry forest in southwestern Puerto Rico mist-net data indicate that dry season bird populations dropped by two-thirds over two decades (86).

Common drivers of plummeting bird populations mentioned in the literature are habitat loss, effects of temperature on reproduction (87–89), water stress (90), increased predation (91, 92), a disconnect between rates of climate change and climate tracking (93), mismatches between breeding period and prey abundance (94), and atmospheric pollutants and insecticides (95). For farmland birds in Europe, reduction in prey due to herbicides, plowing, and loss of marginal habitats (96) appear to be major factors. We found only two articles that ascribed the collapse of avian populations to the impact of climate warming on prey abundance. Pearce-Higgins et al. (97) analyzed the demographics of the golden plover in Europe and showed that trends in population dynamics were closely associated with temperature-based reductions in prey abundance. Pearce-Higgins et al. (97) also analyzed the diets of 17 upland, insectivorous bird species in the United Kingdom, calculated a climate-sensitivity index for their prey, and showed that the index was correlated with recent negative population trends. In the forests of Victoria, Australia, Mac Nally et al. (98) documented major declines in the resident avifauna regardless of spatial scale, disturbance level, species conservation status, or foraging guild. They attributed the ongoing declines to the impact of increasing temperatures and decreasing rainfall on arthropods and other food resources.

We also note that in many studies of passerine declines, such as Mac Nally et al. (98), insectivores are often the group exhibiting the largest decreases in abundance, as appears to be the case in the Luquillo forest. In Canada and the northeastern United States, the total reduction in aerial insectivores between 1970 and 2010 was greater than 60%, the largest decline of all passerines analyzed (84, 99). Nebel et al. (100) concluded that these declines resulted from impoverishment of flying insect populations, and suggested acid rain as a possible mechanism. Between 1923 and 1998 in Singapore, insectivores experienced the highest level of extinctions of any foraging guild (101), and in fragmented forests in Costa Rica and Brazil, insectivores suffered greater declines and local extinctions than other avian groups (102–104). In a lowland tropical rainforest at La Selva, Costa Rica, Sigel et al. (85) also found that at least 50% of declining bird species were insectivorous.

Given decades of research linking avian survivorship, reproduction, and abundance to food resources (105–107), the paucity of studies on climate-related reductions in arthropod prey is surprising. Galapagos finches provide a classic example of how climatic events associated with El Niño can lead to major fluctuations in food supplies, which in turn drive reductions in abundance of up to 40% per year in Geospiza populations (108). Anomalous rises in seawater temperature associated with ENSO have also caused severe reductions in food supplies, leading to starvation and mass mortality of marine birds (109–111). More recently, ocean warming associated with climate change has precipitated bottom-up trophic cascades, analogous to those proposed here for terrestrial systems, which have restructured both benthic and pelagic food webs and driven declines in shearwaters and auklets (112–114).

Our multiple regression results speak to the impact of both rising and falling temperatures on vertebrates and invertebrates in the Luquillo forest. While the effect of MnMaxT on abundance was negative in all models, the proportion of temperatures greater than 29 °C had both negative (El Verde E. coqui) and positive (walking sticks, E. coqui at Sonadura West and Sonadura Old) effects (SI Appendix, Table S5). This was perplexing, as we expected the frequency of extreme temperatures to track the rise in mean temperature. While this is true over the long term (SI Appendix, Fig. S8E), it was not the case during the Sonadura West and Sonadura Old censuses when the proportion of temperatures greater than 29 °C fell from 0.1 to 0.0, or during the walking stick census when the proportion dropped from 0.3 to 0.0. Only during Stewart’s earlier E. coqui census (28) did the proportion of higher temperature rise steadily from 0.0 in 1982 to 0.1 in 1992.

The time line in SI Appendix, Fig. S9 places the dates of the population studies analyzed here into the context of El Niño and La Niña years as measured by the SOI. Positive SOI values (SI Appendix, Fig. S9, blue areas) correspond to La Niña years and negative values (SI Appendix, Fig. S9, red areas) to El Niño years. As SI Appendix, Fig. S9 illustrates, the canopy and walking stick censuses spanned three La Niña and three El Niño episodes. As SI Appendix, Fig. S9 also shows, the eruption of Mount Pinatubo, which lowered global temperature by as much as 0.6 °C over the next 2 y, occurred during the walking stick, Sonadura West, and Sonadura Old censuses. The effects of Pinatubo most likely lowered the frequency of extreme temperatures at El Verde despite a simultaneous strong El Niño episode.

In a recent analysis of the growing frequency of extreme weather, Diffenbaugh et al. (115) predicted a two- to threefold increase in the probability of historically unprecedented maximum daily temperatures between 2018 and 2036, and a three- to fivefold increase between 2036 and 2055. If these predictions are realized, the impact on food webs and ecosystems will be even more rapid than expected from the steady, linear increases in mean temperatures that are currently ravaging the natural world.

Numerous studies have documented the impact of the ENSO on all of the major taxa considered here, including invertebrates (116), frogs (63, 73, 77, 117), lizards (1, 5), and birds (118). In keeping with our HP results, the multiple regressions for walking sticks and the El Verde E. coqui populations included the SOI as an important predictor of abundance, and the cross-correlations between the SOI and abundances of canopy invertebrates and walking sticks (SI Appendix, Fig. S6) underscore the fluctuating effects of La Niña and El Niño. The reason for the 2-y time lag in the effects of ENSO on walking stick abundance is not clear. Changes in temperature often lag phase shifts in the SOI by up to 9 mo, and this effect may be partly responsible.

From data gathered over 55 y and 16 locations in Puerto Rico, the National Weather Service found that temperatures during La Niña events averaged ∼2.0 °C lower during the dry season (December to August), and 1.5 °C lower during the wet season (May to November) compared with El Niño episodes. ENSO also impacts rainfall across Puerto Rico, increasing dry season precipitation by 13% during El Niño years and causing a 14% increase in wet season precipitation during La Niña years. We interpret the positive effect of La Niña on invertebrate abundance as the result of cooler, wetter conditions conducive to increased survivorship and population growth compared with the hotter and seasonally drier conditions of El Niño years. Stapley et al. (20) reached a similar conclusion concerning the effects of ENSO on Anolis populations in Panama.

While we have focused on population collapse driven by climate change, shifts in altitudinal distributions of arthropods, insectivores, producers, and top predators may also play an important role in the observed declines in abundance. In montane habitats upward migrations into cooler refugia can spell the difference between survival and extinction as temperatures rise at lower elevations (119). However, altitudinal shifts can also cause rapid changes in community structure, disrupting coevolved interactions, creating “no-analogue communities,” and contributing to community disassembly via accelerated extinctions (120, 121). The reductions in habitat area as species move up mountains also increase extinction rates (87), and can erode genetic diversity, further reducing fitness and response to ongoing changes in climate (122, 123). At present, the extent to which these systemic factors may have influenced the midelevation population declines documented here is unknown. However, Luquillo forest birds and frogs are known to exhibit elevational range shifts similar to those encountered on other tropical mountains. Since 1998, 8 of 33 Luquillo forest bird species (124) and 10 of 12 resident frog species (125) have undergone elevational shifts.

An alternative to the climate-warming hypothesis for population declines in Luquillo and other tropical montane habitats is a mechanism proposed by Lawton et al. (126). In theory, clear cutting of lowland forests raises surface temperature, which in turn increases conductive heat flux from exposed ground to the atmosphere. This reduces both evapotranspiration and evaporative cooling (126) and lowers humidity. Pounds et al. (63) have presented a detailed rebuttal of the Lawton scenario with respect to montane forests in Costa Rica. With regard to the Luquillo rain forest, the assumptions of the Lawton hypothesis simply do not hold, given a major regeneration of the lowland forest since the 1950s after a transition from an agrarian to a manufacturing economy (127). Consequently, the degree of shade, cooling, and humidification of the habitats surrounding the Luquillo Mountains have all increased (128).

Our research was carried out in the Luquillo Experimental Forest (LEF), which encompasses 10,871 acres in the Luquillo Mountains of northeastern Puerto Rico. Elevations extend from <100 m to 1,078 m. Rainfall and temperature vary with altitude, ranging from 150 cm y⁻¹ and 26 °C at lower limits, to 450 cm y⁻¹ and 18.9 °C on the highest peaks. We conducted our field work in the same study area utilized by Lister (22) in the 1970s. This location (18°19′584″N longitude, 65°228″W latitude) is situated on a northeast facing slope at 440 m. The vegetation is midelevation Tabonuco (Dacryodes excelsa) forest (141), characterized by an 8- to 10-m canopy, low ground light levels, and high plant diversity (171 tree species). We also sampled arthropods in a tropical dry forest in the Chamela-Cuixmala Biosphere Reserve in western Mexico (19°22003″–19°35011″N, 104°56013″ to 105°03025″W), to compare recent samples with those we obtained in the same study area in 1986 and 1987 (25).

Lister (22) sampled arthropods within the Luquillo forest during July 1976 and January 1977. Following the same procedures and using the same study area, arthropod abundances were again estimated during July 2011 and January 2012 using both sticky-traps and sweep netting. Our 10 traps were the same size (34 × 24 cm) as Lister’s (22), and also utilized Tanglefoot as the sticky substance. Traps were laid out on the ground in the same-sized grid (30 × 24 m), and also left uncovered for 12 h between dawn and dusk before all captured insects were removed and stored in alcohol. Hoop sizes of our sweep nets (30-cm diameter) matched those used by Lister (22). Body lengths of all captured arthropods were measured to the nearest 0.5 mm using a dissecting scope and ocular micrometer. Regression equations were used to estimate individual dry weights from body lengths (142, 143).

To compare Anolis densities with Lister’s (22) estimates from July 1976 and January 1977, we sampled anoles within the same 15 × 15-m quadrat during July 2011 and January 2012. Following Lister (22), we used the Schnabel multiple recapture method (27) to estimate densities. However, instead of marking captured lizards by toe clipping, we used Testor’s enamel paint to create small (∼2 mm) spots with different color combinations directly above the dorsal base of the tail.

We analyzed climate data taken at two locations in the Luquillo forest: the United States Forest Service El Verde Field Station and the Bisley Lower meteorological tower, which is part of the Luquillo Critical Zone Observatory. The El Verde station lies 5 km southwest of our study area (18.3211° N, 65.8200° W), at an elevation of 350 m. The upper Bisley Tower is located 3.2 km southeast of our study area (18.3164 N, 65.7453 W) at 352 m in elevation. Temperature data for the El Verde station span 37 y, from 1978 to 2015 (Fig. 1A), and for the Bisley station 21 y from 1993 to 2014 (Fig. 1B). Given that the highest ambient temperatures for a given area should have the greatest impact on fitness, especially for ectotherms (144), daily maximum temperatures were utilized in our analyses. Climate data for the Estacion de Biologia Chamela were obtained from www.ibiologia.unam.mx/ebchamela/www/clima.html.

Data sets from the Luquillo long-term ecological research (LTER) online Data Center were downloaded and analyzed for trends in population abundances over time. Detailed methods employed in the various studies can be found at the LTER Data Center website (https://luq.lter.network/luquillo-information-management-system-luq-ims).

All analyses were conducted with R 3.4.2 using RStudio v1.0.153. Poisson regressions and the nonparametric Theil-Sen regression were carried out with the glm2 (146) and mblm (147) R packages. Poisson regressions were tested for overdispersion using the AER package (148). Because all dispersion tests rejected the null hypothesis, only quasi-Poisson regressions were used. Guided by Pounds et al. (63, 73), Aguilar et al. (149), Stapley et al. (20), and García-Robledo et al. (53), we constructed 11 climate-related indices from the Luquillo LTER temperature and rainfall data for use in HP and regression analyses. SOI data were downloaded from the National Center for Atmospheric Research website (www.cgd.ucar.edu). Olea et al. (150) have shown that for more than nine explanatory variables, the order in which variables are entered in an HP analysis can affect the amount of independent variance explained by a given predictor. Hence, we eliminated 2 of the 11 climate variables that had the highest correlations with the other climate variables. The final nine variables are defined in SI Appendix, Table S3. If a given variable did not make at least a 5% independent contribution to the explanation of the deviance, it was not included in the regression analyses. For the HP analyses we employed the hier.part R package (151), and for the Vinod Causality analyses the generalCorr package (37).

We did not utilize automated regression techniques, such as stepwise selection and backward elimination, due to a variety of issues including erroneous type I error rates and inflated R² measures (33). Climate variables are often highly correlated and can generate spurious results when used in regression (34, 134). To guard against such effects of multicollinearity, we constructed correlation matrices of the HP predictors and eliminated variables having correlations of 0.7 or higher. All multivariate quasi-Poisson regressions were then checked for multicollinearity using the Variance Inflation Factor and Bartlett’s Test of Sphericity. Because Akaike’s Information Criterion cannot be calculated for quasi-Poisson models, we selected models having the smallest deviance residuals and the lowest χ² test statistic P values.