Significance

Arthropods, invertebrates including insects that have external skeletons, are declining at an alarming rate. While the tropics harbor the majority of arthropod species, little is known about trends in their abundance. We compared arthropod biomass in Puerto Rico’s Luquillo rainforest with data taken during the 1970s and found that biomass had fallen 10 to 60 times. Our analyses revealed synchronous declines in the lizards, frogs, and birds that eat arthropods. Over the past 30 years, forest temperatures have risen 2.0 °C, and our study indicates that climate warming is the driving force behind the collapse of the forest’s food web. If supported by further research, the impact of climate change on tropical ecosystems may be much greater than currently anticipated.

Abstract

A number of studies indicate that tropical arthropods should be particularly vulnerable to climate warming. If these predictions are realized, climate warming may have a more profound impact on the functioning and diversity of tropical forests than currently anticipated. Although arthropods comprise over two-thirds of terrestrial species, information on their abundance and extinction rates in tropical habitats is severely limited. Here we analyze data on arthropod and insectivore abundances taken between 1976 and 2012 at two midelevation habitats in Puerto Rico’s Luquillo rainforest. During this time, mean maximum temperatures have risen by 2.0 °C. Using the same study area and methods employed by Lister in the 1970s, we discovered that the dry weight biomass of arthropods captured in sweep samples had declined 4 to 8 times, and 30 to 60 times in sticky traps. Analysis of long-term data on canopy arthropods and walking sticks taken as part of the Luquillo Long-Term Ecological Research program revealed sustained declines in abundance over two decades, as well as negative regressions of abundance on mean maximum temperatures. We also document parallel decreases in Luquillo’s insectivorous lizards, frogs, and birds. While El Niño/Southern Oscillation influences the abundance of forest arthropods, climate warming is the major driver of reductions in arthropod abundance, indirectly precipitating a bottom-up trophic cascade and consequent collapse of the forest food web.

From pole to pole, climate warming is disrupting the biosphere at an accelerating pace. Despite generally lower rates of warming in tropical habitats (1), a growing body of theory and data suggests that tropical ectotherms may be particularly vulnerable to climate change (2). As Janzen (3) pointed out, tropical species that evolved in comparatively aseasonal environments should have narrower thermal niches, reduced acclimation to temperature fluctuations, and exist at or near their thermal optima. Consequently, even small increments in temperature can precipitate sharp decreases in fitness and abundance. These predictions have been verified in a variety of tropical reptiles, amphibians, and invertebrates (4⇓⇓⇓–8).

Given their abundance, diversity, and central roles as herbivores, pollinators, predators, and prey, the response of arthropods to climate change is of particular concern. Deutsch et al. (5) predicted that, for insects living at mid-to-high latitudes, rates of increase should grow as climate warms, while in the tropics insects should decline by as much as 20%. Reduction in population growth, combined with elevated metabolic rates, could potentially lower abundances and raise arthropod extinction rates. If these predictions are realized, climate warming may have an even more profound impact on the functioning and biodiversity of the Earth’s tropical forests than currently anticipated.

Although arthropods comprise more than two-thirds of all terrestrial species and are centrally important to the ecological well-being of the Earth’s ecosystems, long-term data on population abundance and extinction rates are severely limited (9). Studies documenting declines in insects have focused on temperate species (10⇓⇓–13), and have identified climate warming, along with habitat disturbance and insecticides, as primary causal mechanisms (9, 10). While demonstrated impacts of climate change on tropical forests include reductions in plant diversity (14), changes in plant species composition (15), and increases in tree growth, mortality, and biomass (16), little is known about the impact of climate warming on rainforest arthropods (17, 18). Here we analyze long-term data on climate change, arthropod abundance, and insectivores within the Luquillo rainforest in northeastern Puerto Rico, with the aim of determining if increases in ambient temperature may have driven reductions in arthropod numbers and associated decreases in consumer abundance.

Results

Climate Trends in the Luquillo Forest.

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.

Air temperature trends within the Luquillo forest. MnMaxT vs. year for the (A) El Verde (1978–2015) and (B) Bisley Tower (1994–2012) meteorological stations. The stations are ∼8.2 km apart at elevations of 350 and 352 m, respectively. The least-squares regression lines have been drawn through the data points. The difference between the slopes of the two regression lines is not significant (t = 0.281, P = 0.78).

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).

Arthropods.

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.

Mean dry-weight arthropod biomass per 100 sweeps taken in the same sample area in the Luquillo rainforest during July 1976, January 1977, July 2011, and January 2013. One SE around the mean biomass is shown for each bar. Total sweeps taken in each period was 800, except for July 1976, when 700 sweeps were taken. Data for 1976 and 1977 are from Lister (22).

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.

Comparison of the average dry-weight biomass of arthropods caught per 12-h day in 10 ground (A) and canopy (B) traps within the same sampling area in the Luquillo rainforest. Numbers above the bars give the mean daily catch rate in dry weight of arthropods per day for the respective dates. Data for 1976 and 1977 are from Lister (22).

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.

Trends in the abundance of canopy arthropods and walking sticks in the Luquillo forest El Verde study area. (A) Linear regression of the total number of canopy arthropods captured per foliage weight sampled at El Verde against the period when the samples were taken. (B) Cubic regression for the total number of canopy arthropods captured per foliage weight sampled against the MnMaxT during the period when the samples were taken. (C) Quasi-Poisson regression of total number of walking sticks vs. the period when the population was sampled. (D) Quasi-Poisson regression of total number of walking sticks vs. the MnMaxT during the period when the population was sampled. The 95% confidence intervals are shown around the best-fit regression lines. For the Poisson regressions, Pr(χ) is the result of a likelihood-ratio χ2 test of whether the independent variable improves the Poisson model beyond an intercept-only model. P < 0.05 indicates a statistically significant regression.

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).

Insectivores.

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.

Anoles.

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-m2 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.

Comparison of anoles censuses conducted in the same Luquillo rainforest study area during July 1976 and January 1977, with censuses conducted during July 2011 and January 2012. The total number of anoles in each census, estimated by the Schnabel multiple capture–recapture method (27), is given above each histobar. Brown represents male and female A. gundlachi. Green represents A. evermanni. Gray represents A. stratulus. SI Appendix, Table S1 gives the 95% confidence intervals around estimated abundances for A. gundlachi and A. evermanni.

E. coqui.

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.

Population trends for E. coqui and birds near the El Verde Field station. (A) Quasi-Poisson regression of estimated total number of E. coqui individuals against time from censuses conducted by Stewart (28) in the Activity Transect. (B) Quasi-Poisson regression of the estimated number of E. coqui individuals against MnMaxT during the time periods when Stewart’s censuses were conducted. (C) Quasi-Poisson regression for the total number of birds captured during equal length, 4-d sessions of mist netting (31) near the El Verde Field Station against the period when the mist netting was conducted. (D) Quasi-Poisson regression of the total number of birds captured during Waide’s (31) 4-d sessions vs. MnMaxT during the year of mist netting. The 95% confidence intervals are shown around the best-fit regression lines. Pr(χ) is the result of a likelihood-ratio χ2 test of whether the independent variable improves the model beyond an intercept-only model. P < 0.05 indicates a statistically significant regression.

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.

Birds.

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 2K 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.

Vinod Causality.

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.

Discussion

Climate Change and Arthropods.

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.

Anolis Population Declines.

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.

Eleutherodactylus Declines.

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).

Avian Declines.

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).

Multiple Regression and Cross-Correlation.

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).

Conclusion

In general, a more integrated understanding of how climate warming and global cycles affect population and food web dynamics awaits more extensive time series that encompass all levels of tropical forest ecosystems, from producers to top predators. The impacts of major climatic perturbations on the forest food web, such as the hurricanes and prolonged droughts driven by ENSO (SI Appendix, Fig. S9), also need further study. While the initial damage from Hurricane Hugo was devastating (129), the Luquillo forest vegetation recovered remarkably quickly, with 70% of trees at El Verde producing new leaves after just 7 wk (130), and after Hurricane Georges insect populations at all sites studied by Barberena-Arias and Aide (131) recovered in less than 1 y. The resilience of the forest was also manifested after Hurricane Maria when researchers found many locations well on their way to recovery after 2 mo (132). However, as climate warming continues, the frequency and intensity of hurricanes in Puerto Rico are expected to increase (133), along with the severity of droughts and an additional 2.6–7 °C temperature increase by 2099 (134), conditions that collectively may exceed the resilience of the rainforest ecosystem.

The central question addressed by our research is why simultaneous, long-term declines in arthropods, lizards, frogs, and birds have occurred over the past four decades in the relatively undisturbed rainforests of northeastern Puerto Rico. Our analyses provide strong support for the hypothesis that climate warming has been a major factor driving reductions in arthropod abundance, and that these declines have in turn precipitated decreases in forest insectivores in a classic bottom-up cascade. This hypothesis also provides a parsimonious explanation for why similar cross-taxa, concordant decreases in reptiles, anurans, and birds have occurred in Costa Rican rainforests (3, 63), and are likely occurring across a broad range of tropical ecosystems. Overall, there is an urgent need for more widespread monitoring of arthropods and insectivores throughout the tropics (135, 136). As the sixth mass extinction continues to decimate the world’s biota (137, 138), these data will be crucial to understanding the impact of climate change on terrestrial food webs (139), ecosystem dynamics (140), and biodiversity (8), and to formulating conservation strategies aimed at mitigating the effects of future climate forcing.

Materials and Methods

Study Sites.

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−1 and 26 °C at lower limits, to 450 cm y−1 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).

Arthropod Samples.

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).

Anolis Abundance.

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.

Climate Data.

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.

Canopy arthropods.

Data were collected by Schowalter (23) near the El Verde field station between February 1991 and June 2009. Several articles have analyzed these samples with respect to invertebrate diversity, functional groups, arthropod composition in gap and intact forest, and recovery from disturbance (145), but none have looked for trends in overall abundance. Here we summed all arthropods sampled each year across taxa, forest type, and tree genera.

Walking sticks.

We analyzed data from a census of walking sticks (Lamponius portoricensis) carried out by Willig et al. (24) between 1991 and 2014 in the 16-ha Luquillo Forest Dynamics Plot (LFDP) near the El Verde Field Station. Sampling was conducted during the wet and dry seasons and captured individuals were classified as adults or juveniles. To analyze walking stick abundance through time, we summed all juveniles and adults across seasons and land classes.

E. coqui abundance.

We analyze census data for the Puerto Rican frog E. coqui taken by Woolbright (29, 30) between 1987 and 1997 at study areas near the El Verde Field Station and in the Bisley watershed, as well as E. coqui census data taken by Stewart (28) in another study area at El Verde. For Woolbrights’s (29, 30) census data at the Sonadura Old study area, we removed one outlier for the 1995 census data when E. coqui numbers increased from 29 to 236 and returned to 45 individuals in the next census.

Mist netting.

We analyze mist netting data taken by Waide (31) in the LFDP from June 1990 until May 2005. Nets were run from dawn to dusk for a total of 4 consecutive days each year during late May and early June. Only the post-Hurricane Hugo capture rates for June 1990 have previously been published (31). We calculated the total number of individuals caught per year across all species across all 15 y of data. Given that the ruddy quail dove (G. montana) is an exclusive granivore (SI Appendix, Fig. S4) and could confound patterns in insectivore abundance, we excluded this species from counts of captured birds. Recaptures of the same individuals within a given sampling period were also eliminated from the analysis.

Statistical Analyses.

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 R2 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 χ2 test statistic P values.

Acknowledgments

We thank Ray Huey for his advice throughout this study; Patricia Burrowes, Luisa Otero, and Hector Alvarez for their support in Puerto Rico; Sarahi Torobio, Eduardo Avila, Helda Hernandez, and Katherine Lister for assistance in Mexico; Felipe Cano (US Department of Agriculture Forest Service) for facilitating research at Luquillo; Wayne Arendt (US Department of Agriculture Forest Service) for sharing his knowledge of Luquillo avian populations; Prof. H. D. Vinod for guidance on using his generalCorr R Package; and our reviewers for their comments, which broadened our understanding of climate warming effects and strengthened our analyses. We received support from the Department of Biological Sciences at Rensselaer Polytechnic Institute, the graduate program in the biological sciences at the Instituto de Biologia, Universidad Nacional Autonoma de Mexico, and the Estacion de Biologıa Chamela of the Instituto de Biologia- Universidad Nacional Autonoma de Mexico. This research was supported by National Science Foundation Grant 1038013.

Blood-sucking sand flies from disparate global regions have a predilection for feeding on the marijuana plant (Cannabis sativa), and the findings hint at a potential avenue for controlling sand flies, which can transmit leishmaniasis.