Annual Tropical‐Rainforest Productivity Through Two Decades: Complex Responses to Climatic Factors, [CO2] and Storm Damage

Yearly changes in tropical carbon‐cycling have been a major biotic determinant of the interannual variation in the rise of atmospheric carbon dioxide ([CO2]). The environmental responses underlying these changes remain poorly understood. A 21‐year field study (1997–2018) in a Costa Rican rainforest has produced the first multi‐decade record of landscape‐scale annual aboveground productivity for a tropical ecosystem. While none of the four production components (wood production; leaf‐, twig‐, and reproductive litterfall) increased over two decades, all significantly varied among years. Multi‐factor environmental models explained a large fraction of the inter‐year changes in two production components. Annual changes in three climatic factors and [CO2] jointly accounted for 84% of the variation in annual wood production. While the strongest response of wood production was to decline in years of slightly warmer nights, it also significantly declined in years with more dry‐season hours hotter than 28°C. The yearly changes in reproductive litterfall were largely (52%) determined by the combined effects of increasing [CO2] and declining annual radiation. Twig litterfall declined over the 21 years. Annual leaf litterfall, the largest component of field‐assessed aboveground production, showed no significant environmental responses. An extreme storm in the study's final year disrupted the environmental controls on production that prevailed through the prior 20 years. While important uncertainties are introduced by some of the field methods for assessing aboveground forest production, the findings of this first long‐term field study of annual landscape‐scale productivity raise deep concerns for the future of tropical forests as global warming proceeds.

increased with higher yearly tropical mean temperatures (partial r = +0.71, after factoring out the tropical rainfall anomaly; Piao et al., 2020). A parallel analysis of the past 50 years (Anderegg et al., 2015) found that the IAV of the inferred global land C balance was most strongly correlated with annual means of tropical daily minimum temperature (T min ); in years of warmer tropical nights, inversion studies indicated greater tropical net C emissions and increased global land respiration. These findings implicate a large negative impact on the global land C sink due to the response of tropical ecosystem respiration to nighttime warming. Increases in nighttime tropical plant respiration would in turn reduce tropical NPP. As yet no comparable global analysis has indicated negative effects of high daytime temperatures on tropical NPP, the link expected when warming decreases tropical photosynthetic C uptake. Nevertheless, such effects are reasonably expected to be detectable soon (Huang et al., 2019). Leafand stand-level (eddy-flux) studies in tropical forests (Doughty & Goulden, 2008;Loescher et al., 2003;Miller et al., 2021;Pau, Detto, et al., 2018;Slot & Winter, 2017;Tribuzy, 2005) have all indicated substantial drops in C uptake by tropical forests at air temperatures (T air ) above 28-29°C, a level now often exceeded in the daytime in that biome (D.A. Clark et al., 2013;Fontes et al., 2018;Pau, Detto, et al., 2018).
Variation in annual water availability is another likely driver of the IAV of tropical C balance. While annual tropical rainfall itself was not significantly related to the tropical C-balance anomaly over recent decades (partial r = −0.14, after factoring out the tropical temperature anomaly), the impacts of annual changes in other potentially important aspects of tropical water availability have been little studied (Piao et al., 2020). The elevated Vapor Pressure Deficit (VPD) that accompanies higher temperatures is likely to increasingly constrain productivity in many biomes (Grossiord et al., 2020;Novick et al., 2016;Yuan et al., 2019). Annual variation in soil moisture and in dry-season rainfall are additional potential drivers.
Multiple additional factors could also be contributing to the IAV of tropical C cycling. A major assumption of the global land models has been that tropical NPP has been increasing due to rising [CO 2 ]: the "CO 2 fertilization hypothesis" (Walker et al., 2021). This effect, however, could be lessened or prevented by the nutrient limitation affecting large portions of the tropics, by other aspects of C-sink limitation, or by decreasing C residence times (Fleischer et al., 2019;Körner, 2009;Walker et al., 2021). While variation in light levels is a dominant control of short-term (30-min) daytime tropical-forest C balance in eddy-flux studies (e.g., Doughty & Goulden, 2008;Loescher et al., 2003;Vourlitis et al., 2005), field studies have not yet demonstrated a link between the IAV of tropical C-cycling and annual changes in radiation. A factor that could increasingly impact tropical NPP is forest damage from cyclonic storms (Negrón-Juárez et al., 2018), which are expected to increase in frequency and intensity with warming (McDowell et al., 2020). By causing both elevated tree mortality and major crown damage to surviving trees, such storms could cause repeated short-to medium-term drops in forest productivity.
To quantify tropical C-cycling and the environmental drivers underlying its interannual changes, in situ studies are critically needed in the major tropical biomes (forests, savannas, other semiarid ecosystems). The tropical-forest biome has a particularly large role in both regional and global C dynamics. Tropical forests have been estimated to contain 25% of land biosphere C and to account for ca. 33% of global terrestrial NPP (Bonan, 2008). Field studies to quantify this biome's C-cycling and to identify the environmental drivers of its IAV face two major challenges. Studies of ecosystem processes in tropical forests need to be conducted at the landscape scale because of these forests' internal edaphic mosaics (e.g., Itoh et al., 2012;Valencia et al., 2009;Yoneda et al., 2017) and their characteristically localized, sporadic disturbances (Chambers et al., 2013;Marvin et al., 2014). A single small (≤1-2 ha) plot per forest, the sampling approach that has dominated tropical-forest studies to date (e.g., Brienen et al., 2015;Hubau et al., 2020;Malhi et al., 2015;Qie et al., 2017), will not suffice for quantifying a given forest's landscape-level processes. Second, to quantify the IAV of forest C-cycling and to detect its environmental drivers requires continuous monitoring through many years. For example, to research the IAV of aboveground NPP (ANPP*; D.A. Clark et al., 2001), litterfall must be continuously collected at short (e.g., 2-week) intervals, the growth of all woody stems needs to be measured at least annually over the same period, and the study needs to span many years (D.A. Clark & D.B. Clark, 2011). A few years of coordinated field productivity measurements (e.g., Aragão et al., 2009;Doughty et al., 2014;Girardin et al., 2010;Malhi et al., 2009) can provide useful estimates of NPP components but are insufficient for assessing their IAV and the underlying environmental drivers.
To date, these requirements have been met by only one tropical-forest ANPP* study: the CARBONO Project in the old-growth tropical rainforest at the La Selva Biological Station in Costa Rica. The first 12 years of the landscape-scale CARBONO measurements, together with local meteorological monitoring, enabled a first assessment of the responses of tropical-forest ANPP* to yearly variation in climatic factors and [CO 2 ] (D.A.

10.1029/2021JG006557
3 of 18 Clark et al., 2013). The study has now produced more than two full decades  of uninterrupted field measurements of annual landscape-level ANPP*. Here we report the findings from this unique tropical-rainforest study. The 21-year observations indicate: no detectable [CO 2 ] impact on the largest productivity component; highly predictable long-term control of two of the four ANPP* components by multiple environmental factors; disruption of these controls by an extreme storm event; and no directional increase in this tropical forest's productivity.

Study Site and Plot Network
The CARBONO Project was carried out in old-growth lowland tropical rainforest at the Organization for Tropical Studies' La Selva Biological Station. This research site in Costa Rica (10°26' N, 83°59' W, elevation range 37-150 m) adjoins ca.100,000 ha of continuous protected rainforest. Over the 21-year study period, the mean annual temperature was 25.1°C and annual rainfall averaged 4,441 mm. Lower rainfall typically occurs in January-April, a period we refer to as the "dry season"; however, for each of these 4 months the long-term average rainfall exceeds the 100-mm threshold generally used (Dewalt et al., 2010) to characterize tropical dry-season months. La Selva's volcanically derived soils place this forest at the high-fertility end of the spectrum of lowland neotropical forests (Powers et al., 2005).
The ANPP* and soil-moisture measurements were taken in the CARBONO Project's network of 18 0.5-ha plots. These plots were sited across ca. 500 ha of old-growth forest to sample, with replication (6 plots each), the three dominant edaphic conditions of the upland old-growth forest landscape: younger oxisol (alluvial) terrace; older oxisol plateau; and older oxisol slope. The plots were located using a random-block design based on coordinates drawn from the digital elevation model and the soil coverage of La Selva's Geographical Information System. The plot-network map and details of the plot installation are provided in D.A. Clark & D.B. Clark (2021). Each measurement-year was the period: 1 October, Year 1-30 September, Year 2.

Productivity Measurements
ANPP* is a field-based estimate of forest aboveground NPP (D.A. Clark et al., 2001;Huston & Wolverton, 2009), obtained through measurements of aboveground wood production and the three categories of tissues in fine litterfall (leaf, reproductive, twig). It addresses a large but unknown fraction of total NPP. Due to intractable challenges for yearly field-monitoring of other productivity components (see D.A. Clark, Asao, et al., 2017), ANPP* omits multiple aboveground components and all of belowground NPP, which is potentially large and time-varying. Through the 21 years ending in late 2018, we measured the four components of ANPP* continuously in each of the 18 CARBONO plots. The field ANPP* data from all plots are provided in two data depositions (fine litterfall: D.A. Clark, 2021; wood production: D.A. Clark & D.B. Clark, 2021). Table S2 in D.A.  provides the annual landscape-level ANPP* data.
In each plot, we collected fine litterfall (leaves, twigs <1 cm dia., reproductive materials) biweekly from nine standing basket traps (0.25 m 2 area, 0.8 m basket-frame height) and nine ground-level "traps" (0.25 m 2 vertically projected squares demarcated on the ground) for collecting leaf litterfall items >50 cm length (e.g., 6-m long palm leaves), which are not reliably captured by standing basket traps (Villela & Proctor, 1999). When a trap was damaged or blocked from above, the contents were discarded. Plot-level litterfall was quantified based on the summed area of the collected traps. At each collection, the contents from all collected litter traps in a plot were combined by trap type. The plot-level combined collection was sorted by category and oven-dried at 65°C to constant mass. For each CARBONO measurement-year, each component of annual fine litterfall was quantified as the area-corrected dry mass summed across the 26 bi-weekly collections, annualized to 365 days.
At installation of the 18 plots in 1997, all included live woody stems of ≥10 cm diameter were mapped and identified (>250 species of trees, palms, lianas). Beginning in September 1997, all individuals were then censused annually for growth, mortality, and recruits (stems newly reaching 10 cm diameter). The plots were measured in the same sequence each year, during September-November/December. At each annual census, the diameter of all live stems ≥10 cm diameter was measured with a fabric diameter tape to the nearest 1 mm (rounding down) at a permanently marked point of measurement on the stem (POM), at 130 cm height above the ground or above any basal stem irregularities or buttresses. For stems with higher POM's, the diameter was measured by a technician(s) on one to four 3-m ladders. Through the study the same two field technicians made all the diameter measurements, with consistently high repeatability. At the outset of each annual census, measurement precision was assessed by re-measuring the diameter of ca. 80 stems after an interval of 1-3 days; in each of the 22 years, 70%-89% of the re-measurements matched the first and 97%-100% were within 1 mm. When a stem's POM was found to be threatened by encroaching stem irregularities or by upward "buttress-creep", the diameter was measured both at the current POM and at a new one higher on the stem, above the irregularity, to maintain continuity of that stem's record. The diameter-measurement protocols and the qa/qc of the data are documented in more detail in the data deposition (D.A. Clark & D.B. Clark, 2021) Following each annual census, the estimated aboveground biomass (EAB, kg) of each live individual ≥10 cm diameter in the plots was calculated from the current diameter using the Tropical Wet Forest biomass allometric equation of Brown (1997). That allometry was based on 175 harvested trees, 100 of which were from forest within 30 km of La Selva. The annual estimated aboveground biomass increment (EABI, kg) for each live stem was calculated as the difference between the stem's successive annual biomass estimates, with the difference then annualized based on the actual days between the two measurements. For the few multiple-stemmed individuals (23 of 6,705), each individual's diameter and EAB were set to those of a tree with the summed basal area of its live stems. For stems with a POM change(s) during the study, yearly EAB and EABI were based on the stem's diameter measurement at the original POM, augmented by the subsequent annual diameter increments; this procedure prevents progressive biomass underestimation when the POM is moved upward, where diameter may be reduced by bole taper. Annual estimates of plot-level wood production (plot EABI) were calculated as the sum of the EABI's of all surviving individuals ≥10 cm diameter in the plot, plus the EABI's of the new recruits (stems that reached ≥10 cm diameter at that census). Each recruit's EABI for the year ending in its first measurement year was calculated as their EAB at that first census, minus the EAB of a 10-cm-diameter stem (D.A. Clark et al., 2001).

Meteorological and [CO 2 ] Data and the Derived Annual Metrics
Local meteorological data from the Organization of Tropical Studies' long-term weather station at La Selva were used for this study. The automated system recorded 30-min (or hourly) and daily summary data from sensors for air temperature, relative humidity, total radiation (pyranometer), and rainfall (tipping-bucket), with parallel manual data for daily rainfall (pluviometer) and for daily maximum and minimum temperatures (max and min thermometers). We screened for and gap-filled missing or anomalous data using cross-record (e.g., automated vs. manual instruments) checks and regressions. The daily climatic data produced with this qa/qc are in Table  S3 in D.A. . We monitored volumetric soil-moisture at 0-30 cm depth with one Campbell sensor buried vertically in the center of each plot (CS615 FDR sensor or CS616 TDR sensor, with cross-regression between types). These sensors were read biweekly (methods in O'Brien & Oberbauer, 2001) from March 1998 to October 2018, with limited data gaps due to instrument failure. A data deposition (D.A.  provides the documented biweekly soil-moisture data. The annual soil-moisture metrics were based on the biweekly data, averaged across the three edaphic conditions.  Table S1 in D.A. .

Data Analysis
Statistical tests were run using the program JMP 13.2.0 (SAS Institute). We used Repeated Measures Analysis, Mixed Model, to assess whether ANPP* and its four components significantly varied among years, testing for a year effect while controlling for the repeated annual data from each plot. We then used simple linear regression on year to test the five landscape-scale productivity metrics (the yearly means across the 18 plots) for directional change.
The climatic and [CO 2 ] sensitivities of ANPP* and its four components were assessed by testing the yearly values of each landscape-scale productivity metric in a standard set (Table S4 in D.A. Clark, Clark, and Oberbauer (2021)) of 66 linear models based on annual climatic factors and annual [CO 2 ]. Five models included landscape-scale EAB at the beginning of each measurement-year (EAB yr1 ), to test for an effect on productivity from temporal changes in landscape-level EAB. The set of 17 environmental factors tested in these models (Table 1) was based on expectations from the literature of those climatic drivers most likely to impact tropical-forest productivity and on the significant relationships found in the first 12 years of this study (D.A. Clark et al., 2013). Four temperature metrics involved air temperatures above 28°C because of the accumulating evidence of substantial drops in (net) photosynthesis above this threshold (see Introduction). The annual means of daily T min and daily T max have distinct relationships with C cycling (Anderegg et al., 2015;D.A.;Clark et al., 2013;Welch et al., 2010) and were both included in test models. Annual mean temperature was also included because of its widespread use in global carbon models. The VPD metrics were based on the proportion of the daytime period with the high-VPD (>1 kPa) conditions that have been linked to decreased daytime C uptake in tropical forests (Fu et al., 2018;Tang et al., 2020).
The standard model set included single-factor models based on the 17 environmental factors and models based on combinations of up to four of those factors. For [CO 2 ], the a priori expectation was a positive effect, and all tests were 1-tailed. For EAB yr1, all tests were 2-tailed. For climatic factors, a priori expectations differed among productivity components. For wood production and reproductive litterfall, negative effects were expected from higher temperatures and from greater VPD, while positive effects were expected from greater radiation and greater water availability (except, for annual rainfall: tests were 2-tailed for wood production and a negative relation was expected for reproductive litterfall due to leaching). For leaf litterfall and twig litterfall, because a given climatic factor could have simultaneous opposing effects on litter mass (see Results and Discussion), all climatic tests were 2-tailed. Given the large number of tests, we set the threshold significance level for a model at P ≤ 0.01, both for the overall model and for any included factor. This level was chosen to facilitate detection of the important environmental associations while lowering the risk of spurious relationships.
We assessed the climatic and [CO 2 ] responses three ways. (a) In a significant model, the slope of a productivity metric's association with an environmental factor indicated the rate and direction of productivity change per unit of that factor. The (partialized) slope from a multiple-factor model indicated the association with a given factor after isolating its effects from those of the other factors in the model. (b) We followed the formulation of Walker et al. (2021) to calculate for each significantly associated productivity metric a relativized response (β for [CO 2 ]; Ɣ for climatic factors) to each factor's range during the study, as: where productivity metric factormax and productivity metric factormin are the model-predicted values for that productivity metric at the maximum and minimum observed annual values of the environmental factor during the 20-year period 1997-2017, and factor max and factor min are the maximum and minimum observed annual values for the environmental factor. When the significant association was from a multiple-factor model, the relativized response was based on the partial (leverage) relationship of the productivity metric with the focal factor. (c) Productivity metrics often showed a bias across the tested models in the direction of their response to an environmental factor, regardless of significance levels. To "fingerprint" each productivity metric's set of climatic and [CO 2 ] associations, we used binomial tests of the sign of their response direction across models.

ANPP* Magnitude and Composition at Contrasting Spatial Scales
Over the full 21-year study period, landscape-scale ANPP* of this old-growth tropical rainforest (annual means across the 18 plots; Figure 1) averaged 13.25 Mg ha −1 yr −1 , with an among-year Coefficient of Variation of only 5%. In contrast, annual ANPP* varied more than 2-fold (8.68-19.93 Mg ha −1 yr −1 ) at the scale of individual 0.5ha plots.
The short-lived constituents of ANPP* (leaf-, reproductive-and twig litterfall) together dominated landscape-scale ANPP* in all years ( Figure 1). The long-term average contributions of the four components to total ANPP* were: leaf litterfall, 52%; reproductive litterfall, 10%, twig litterfall, 6%; and wood production, 33%. These contributions varied far more at the plot level than at that of the landscape. For example, while leaf litterfall accounted for 48%-58% of landscape-scale ANPP* among the 21 years, it constituted as little as 36% and as much as 67% of annual ANPP* in individual plots. We therefore assess the environmental responses and directional trends of ANPP* at the landscape scale.

Interannual Variation and Directional Trends in Climate and [CO 2 ]
The years of the study period spanned a wide range of environmental conditions at La Selva (Table 1; also  (Table 1) lacked a significant directional trend (linear regression on year). The first two-thirds of our study period (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) was identified by Ballantyne et al. (2017) as a global-scale warming hiatus. The annual averages of daily T min , T max , and T mean , however, all varied strongly among the study years at La Selva. Metrics related to water-limitation also showed significant IAV during the 21 years. Dry-season rainfall and annual soil-moisture deficit both varied more than 3.5-fold. The annual VPD metrics closely followed a sine-wave pattern (third degree polynomial fit with year, r 2 = 0.85); they were lowest in 1999-2005, highest in 2006-2015, and declined to the 1997 level over the last three study years. Only two environmental metrics had a significant (P < 0.01) directional 21-year trend. [CO 2 ] showed a tight linear increase with year, from 366 to 408 ppmv (Pearson's r 2 = 0.997, P < 0.0001). Annual averages of daily radiation significantly declined through the study, but with considerable among-year variation (Pearson's r 2 = 0.32; P = 0.008). The lowest annual mean of daily radiation occurred in 2014-2015, one of the two study years when La Selva received more than 6 m of rainfall (Table S1 in D.A. .

Interannual Variation and Directional Trends in Productivity
Landscape-level annual ANPP* and its four components (Figure 1) all varied highly significantly (P < 0.0001) among the 21 years. The among-year range in landscape-scale ANPP* over the 21-year study period was 3.46 Mg ha −1 yr −1 , a range 2.5 times greater than that observed in the first 12 years (D.A. Clark et al., 2013). Both the 21year maximum (14.76 Mg ha −1 yr −1 ) and the 21-year minimum (11.29 Mg ha −1 yr −1 ) for landscape-scale ANPP* occurred during the last 4 years of the study.

Climatic and [CO 2 ] Responses of Annual ANPP*
Until an extreme windstorm strongly affected the forest in 2018 (see Section 3.7), the IAV of annual wood production, reproductive litterfall, and twig litterfall was highly significantly linked to yearly environmental variation through the 20 years 1997-2017 (Table 2). Notably, however, the sets of drivers differed among these three productivity components. No models significantly predicted the IAV of total ANPP* or of its largest constituent, leaf litterfall.

Environmental Responses of Wood Production
A multiple-factor model (Model 1, Table 2) closely predicted the large interannual changes of this important ANPP* component over 20 years (Figure 2). This model indicated highly significant independent responses of annual wood production at La Selva to the concurrent inter-year changes in four environmental drivers. Yearly wood production declined with greater nighttime temperatures (T min ) and with more prevalent high-VPD conditions in the dry season; it increased with greater dry-season rainfall and with increasing [CO 2 ] (Figure 3). Models 2 and 3 (Table 2), based on subsets of these drivers, were also highly significant but predicted substantially less of the IAV of wood production (61%-63%) than Model 1 (84%). Single-factor models 4 and 5 (Table 2) indicated weaker but still significant associations of annual wood production with two other factors: the annual soil-moisture minimum (positive response), and the dry-season daytime fraction with average temperatures above 28°C (negative response).
When relativized to the proportional 20-year range in each of these six factors, by far the strongest environmental response of annual wood production during these two decades was to decrease with increased nighttime temperatures (Table 3). This negative response of annual wood production was more than double its positive response to the 20-year increase in [CO 2 ] (Table 3). Across the tested models, irrespective of relationship significance levels, wood-production 8 of 18 showed a significantly consistent response direction (Table 4) to five of these six environmental factors, in alignment with each factor's likely effect on NPP.
The ability of Model 1 (Table 2) to capture the 20-year IAV of wood production (Figure 2) was scale-dependent. While this model accounted for 84% of the IAV of landscape-scale wood production (Table 2), it essentially failed at the level of individual 0.5-ha plots. For 17 of the 18 plots, Model 1 was not significant when run on the individual plot's yearly data. For the one exception plot, the model was highly significant (Pearson's r 2 = 0.88, P < 0.0001). For nearly all plots (16 of the 18), however, the response direction indicated in the plot-level model results was in the predicted direction for all or all but one of the four included environmental factors.

Environmental Responses of Reproductive Litterfall
Half of this minor ANPP* component's IAV was accounted for by Model 6 (Table 2), which indicated highly significant positive responses of reproductive litterfall to greater annual averages of both radiation and [CO 2 ]. Across the 20-year ranges of both factors, the relativized increases in this ANPP* component were notably large (Table 3). Reproductive litterfall's relative response to the 20-year range in [CO 2 ] was considerably stronger than that seen with wood production (Table 3). Factors are defined in Table 1. b Probabilities are 1-tail unless otherwise noted. c This model is included here for information purposes but is not in fact significant, given the a priori expectation of a positive effect. Some of this component's significant biases in response direction (declines with higher nighttime temperatures, increases with greater [CO 2 ] and radiation; Table 4) conformed to the expected influences of these drivers on reproductive tissue production. In contrast, the tendencies to decrease with greater dry-season rainfall and greater dry-season soil-moisture could indicate artifactual effects (i.e., pre-collection losses of reproductive litterfall mass through decomposition or leaching).

Environmental Responses of Leaf Litterfall
None of the 66 tested models (Table S4 in D.A.  significantly predicted the interannual changes ( Figure 1) in annual leaf litterfall. This largest ANPP* component did, however, have a distinctive "fingerprint" of biases in response direction to environmental factors (Table 4). Four of its five significant biases (to increase with higher temperatures and to decrease with increased dry-season water) were inconsistent with the promotion of leaf production. Instead these biases could plausibly reflect the factors' promotion of either  (Table 2). Each partial regression plot shows the linear response of annual wood production to the indicated factor, after removing the effects of the other three factors.
stress-related leaf shedding or pre-collection mass losses to leaching or decomposition. Notably, this dominant ANPP* component showed no consistency of response direction to increasing [CO 2 ].

Environmental Responses of Twig Litterfall
Given the strong linear decline of this minor ANPP* component through the 20 years (Section 3.3), annual twig litterfall was necessarily negatively correlated with [CO 2 ] (Model 7, Table 2). Models 8 and 9 (Table 2) indicated negative associations of annual twig litterfall with both annual and dry-season daytime high VPD (Pearson's r 2 = 0.40-0.41). As with [CO 2 ], both of these factors had their lowest values in the study's first years and higher values later in the series (see Section 3.3). All three models, therefore, could reflect coincidental correlation rather than a causal environmental response of twig litterfall. Two of its significant response-direction biases (Table 4)-the tendencies to increase with greater nighttime temperatures and with greater radiation-aligned with those of leaf litterfall, as might be expected due to the physical connection of twigs and leaves. Regarding the other environmental factors, however, the two ANPP* components had contrasting response-direction biases (   Table 1. b Probability levels, P: * ≤ 0.01; ** ≤ 0.001; *** ≤ 0.0001. c +/N = (number of models with positive association/number of models that include the factor).

Environmental Responses of ANPP*
Given the differing environmental responses of the four components (Table 2), the lack of explanatory environmental models of total annual ANPP* is unsurprising. Its "fingerprint" of response-direction biases almost perfectly matched that of its largest component, leaf litterfall (Table 4). Total ANPP* also showed no response-direction bias to the annual changes in [CO 2 ]. The tendencies to increase with greater radiation and to decrease with more prevalent high-VPD conditions were consistent with those factors' effects on production. The significant biases toward decreasing with greater water availability and increasing with higher temperatures were instead consonant with environmental promotion of tissue-shedding or of litterfall mass loss through leaching or decomposition.

Data-Series Length and Environmental-Response Detection
The relationships indicated by Models 1 and 6 (Table 2) increased in statistical strength with more years of measurements (Table 5). While none of the relationships was significant (P < 0.01) at 8 years, all were highly significant (P = 0.005-0.0002) with the full 20-year data series. The sign of the responses to these drivers, however, was nearly perfectly consistent through all record lengths; the one exception was from the 8-year relationship of wood production and dry-season high VPD.

Lagged Environmental Responses?
When the significant models for wood production and for reproductive litterfall (Table 2) were based on prior-year environmental data (1-year lag) instead of on current-year data, none indicated significant environmental associations. For total ANPP*, wood production, leaf litterfall, and reproductive litterfall, there were no significant (P < 0.01) single-factor productivity models based on 1-year-lagged data for any of the 17 tested environmental factors (Table 1). The one marginal case was a weak positive association (Pearson's r 2 = 0.28) of annual wood production with the prior-year averages of daily radiation (P 1-tail = 0.0104). Given wood production's lack of association with current-year daily radiation (single-factor model: Pearson's r 2 = 0.02, P 1-tail = 0.55; also no response-direction bias, Table 4), its weak relationship with the prior-year data could be spurious. The two significant models for twig litterfall (Table 2) were weaker (reduced r 2 , larger P) when run on the prior-year environmental data.

Environmental Responses Disrupted by an Extreme Storm
An exceptional windstorm in May of the study's final year (Rader et al., 2020) caused record tree mortality and associated potential loss of wood production in the CARBONO plots ( The one case of inconsistent sign of effect through the time series; 2-tail P for this (non-significant) positive relationship at 8 years.  [1997][1998]. The aggregate potential lost wood production (the summed prior-year EABI's of the dead trees) was more than double that in any prior year. In spite of the year's favorable conditions (wetter dry season, cooler), wood production in that final year was the second-lowest of the entire record ( Figure 1).
All the highly significant climatic and [CO 2 ] associations of annual productivity through the preceding 20 years (Table 2) weakened or disappeared when this storm year was added to the record. For annual wood production, four of the five highly significant 20-year environmental models became non-significant. The exception, its positive association with dry-season rainfall (Model 2), was much weaker at 21 years (Pearson's r 2 = 0.36, vs. 0.63 at 20 years). Similarly, the highly significant 20-year model of annual reproductive litterfall (Model 6) was non-significant after adding the final year's data. For annual twig litterfall, both significant negative 20-year relationships (Models 8, 9) became slightly weaker with the 21st year (Pearson's r 2 = 0.38, 0.37, respectively). As with the 20-year records, there were no significant 21-year environmental models of annual leaf litterfall.

Confounding Effect of Biomass Recovery?
Landscape-scale aboveground biomass stocks in the La Selva old-growth forest were not static through the study.
With the Very Strong Niño of 1997-8, tree mortality was high in Years 1 and 2 (Table S5 in D.A.  and landscape-scale EAB dropped ca.12%; after EAB recovery over the next 11 years, biomass stocks were roughly stable for six years before the 2018 storm (D.B. D.B. Clark, Clark, &Kellner, 2021, in press). In the standard model set applied to all productivity metrics (Table S4 in D.A. , five models included landscape-scale EAB at the beginning of each measurement-year (EAB yr1 ), to test for an effect on productivity from yearly forest EAB. For annual wood production, EAB yr1 was non-significant and inconsistent in effect direction in four of the five models (positive in two models, negative in three); when EAB yr1 was added to Model 3 (Table 2), its effect was negative and marginally significant (P 2-tail = 0.0096). For both reproductive litterfall and leaf litterfall, none of the five models that included EAB yr1 was significant. For twig litterfall, EAB yr1 was a non-significant factor in the four multi-factor models; in the single-factor model, it was significantly negatively associated with annual twig litterfall, which declined through the study.

Discussion
To improve the global C models and their projections of future climate, a fundamental need is to quantify the environmental responses of tropical C-cycling (Ballantyne et al., 2021;D.A. Clark, Asao, et al., 2017;Piao et al., 2020). This multi-decade study of annual landscape-scale ANPP* in La Selva's old-growth forest has enabled the first such analysis for a tropical biome. While total ANPP* and its four components were either stable or declined through the study, they all varied highly significantly among the 21 years. Only half of La Selva's ANPP* showed significant environmental associations, and the responsive components were affected by different sets of drivers. One important production component (wood production) strongly declined with hotter years, two components (wood production, reproductive litterfall) increased with greater [CO 2 ], and one (reproductive litterfall) increased with greater radiation. Longer data series strengthened these relationships. The highly significant responses through 20 years were disrupted when an extreme windstorm caused record tree mortality.

Drivers of the IAV of Annual Wood Production
Wood production, the only aboveground productivity component involving long-lived materials, averaged 33% of La Selva's ANPP* (Figure 1). Through 20 years, 84% of the IAV of this important production component was explained by the independent effects of four environmental drivers (Figures 2 and 3).
By far the strongest relative environmental response of annual wood production (Ɣ = −2.83; Table 3) was to higher nighttime temperatures. Enhanced respiration C losses on even slightly hotter nights (Figure 3) are thus a dominant control of wood production in this tropical rainforest. The CARBONO record constitutes the first longterm field evidence from a tropical ecosystem of reduced annual aboveground productivity with warmer nights. It aligns with the strong link found by Anderegg et al. (2015) between tropical nighttime temperatures and the IAV of estimated global respiration through recent decades.
The other significant environmental responses of annual wood production at La Selva appear related to impacts on the other big plant C flux, photosynthesis. The second-strongest relative response of annual wood production (β = 1.17; Table 3) was its positive response to the 20-year rise in [CO 2 ]. This β from La Selva should not be compared to tropical-forest β values (0.27, 0.69, 1.2) calculated by Walker et al. (2021: their Table 2). Those values were derived from wood-production data compiled from mostly small (≤1-2 ha) single plots scattered around the tropics (Brienen et al., 2015;Hubau et al., 2020;Yu et al., 2019). Many of the plots were active only at the beginning or end of each study or were limited to one short inter-census period (e.g., Extended Data Figure 2 of Hubau et al., 2020). It is unclear if area-relevant responses to [CO 2 ] can be inferred from these changing sets of widely dispersed plots. Second, the La Selva β was based on a model that controlled for the concurrent responses of annual wood production to three other drivers (Figure 3). With no directional change in aboveground wood production over the 21 years, La Selva's record would have logically produced a β of 0 if simply based on the net long-term change in wood production, the approach underlying the Walker et al. (2021) values. A similar lack of long-term change in wood production has been found for two other tropical forests that have been monitored for long periods at the landscape scale (French Guiana: 37.5 ha over 16 years, Rutishauser et al., 2010;Panama: 47 ha over 30 years, Rutishauser et al., 2020).
A significant, but much weaker, relative response of annual wood production at La Selva (Ɣ = −0.31; Table 3) was to decline with years when more dry-season hours were hotter than 28°C. This is the first long-term field evidence of a daytime high-temperature constraint on annual tropical-rainforest productivity. This finding closely aligns, however, with those from numerous short-term studies (see Introduction) of a significant drop in leaf-level net photosynthesis or in forest-level daytime CO 2 uptake at T air above 28°C (or above the corresponding T leaf : ca. 30-31°C). Through the CARBONO study, T air averaged above 28°C for 31%-51% of daytime hours in La Selva's annual dry seasons (Table 1; yearly data in Table S1 in D.A. . The 20year ANPP* records indicated a reduction in annual wood production of 0.34 Mg ha −1 yr −1 with each additional 10% of the dry-season daytime period at these high temperatures (Table 3). As global warming proceeds, this effect will likely intensify. Reductions of photosynthesis were also indicated by the responses of annual wood production to three aspects of water limitation (minimum annual soil moisture, dry-season rainfall, dry-season prevalence of high-VPD periods; Tables 2 and 3). Notably, however, this important ANPP* component showed no significant response to the yearly changes in radiation (Tables 2 and 4).

Environmental Controls of the Annual Production of Short-Lived Tissues
In field ANPP* studies, the annual production of leaves, twigs and reproductive tissues is estimated from materials collected in litterfall traps. This introduces important ambiguities for assessing these production components' environmental responses, especially in tropical forests. An environmental factor can at the same time: (a) affect production through its effects on photosynthesis and/or respiration; (b) affect the physiological stress that leads to tissue-shedding; and (c) reduce the collected mass by promoting pre-collection leaching or decomposition, processes intensified by hot, wet conditions. Further, consumption by animals causes important pre-collection losses of leaves and reproductive materials in tropical forests (e.g., Metcalfe et al., 2013;Parrado-Rosselli et al., 2006).

Drivers of the IAV of Leaf-Litterfall
None of the 66 test models (Table S4 in D.A.  significantly explained the IAV (Figure 1) of leaf litterfall, the largest ANPP* component at La Selva and in most tropical forests (Huston & Wolverton, 2009). The apparent lack of environmental responses is likely at least partly due to the ambiguities from using leaf litterfall as a production surrogate. For example, increased dry-season water availability could have directly promoted leaf production (thus increasing litterfall mass), while also reducing stress-driven leaf shedding and increasing pre-collection leaching and decomposition (thus decreasing litterfall mass). Higher nighttime temperatures could have had similarly opposing effects. Annual leaf litterfall showed significant biases toward decreasing with wetter dry seasons and increasing with higher temperatures (Table 4). These biases suggest predominantly artifactual effects from environmental IAV (impacts on leaf-litterfall mass rather than on leaf production). Any inter-year differences in pre-collection mass losses from herbivory would have been another artifactual component of the leaf-litterfall IAV.
If leaf production directly responded to environmental IAV through the study years, such effects were not detectable from the annual litterfall. Alternatively, approximately constant leaf production across years might be expected in this old-growth tropical rainforest. With La Selva's Leaf Area Index of 6 (D.B. Clark et al., 2008), the costs of adding even more leaf area (materials, respiration) might not be offset by gains in light interception. The lack of directional change in leaf litterfall through 21 years of rising [CO 2 ] is consistent with the possibility that "CO 2 fertilization" benefits do not extend to leaf production in such forests.

Drivers of the Reproductive-Litterfall IAV
Reproductive litterfall averaged 10% of La Selva's ANPP*. In tropical rainforests, however, a large proportion of fruit and seed production (e.g., 93%-Lugo & Frangi, 1993%;51%-Parrado-Rosselli et al., 2006) does not reach litterfall traps. Most tropical-forest fruits and seeds are animal-dispersed, thus evolved to be eaten. In addition, pre-collection mass losses from leaching and decomposition of high-nutrient reproductive tissues are likely to be non-trivial in the hot, humid conditions. Thus, annual reproductive litterfall is likely only a small portion of annual reproductive production ("crumbs under the dinner table").
Remarkably, however, the inter-year changes in reproductive litterfall appear to reflect yearly changes in reproductive production at La Selva. More than half of the IAV of annual reproductive litterfall through 20 years was explained by a two-factor model (Model 6, Table 2) consistent with environmental effects on photosynthesis.
One highly significant response was to increasing [CO 2 ]. More years of observations brought greater cumulative change in [CO 2 ], increasing the statistical strength of the relationship (Table 5). The [CO 2 ] impact on reproductive litterfall indicated by this model (+0.01 Mg ha −1 yr −1 per ppmv) matched that found separately for wood production (Table 3). The second response of reproductive litterfall-to annual changes in radiation-was highly significant only with the full 20-year record (Table 5), which spanned a long-term radiation decline. When relativized, both responses were notably large (Table 3), reflecting their proportionately large effects on a small ANPP* component (see Walker et al., 2021). Reproductive litterfall showed no directional increase over the two decades of this study, with the [CO 2 ]-related increase being countered by the radiation decline.
Prior studies in Panamanian forests have also linked reproductive production to both radiation and [CO 2 ]. In a 2-year light-augmentation experiment in a tropical dry forest, the light-augmented branches had more reproductive buds (Graham et al., 2003). In a 28-year study in a tropical moist forest (Pau, Okamoto, et al., 2018), annual [CO 2 ] was the environmental factor most strongly related to species-level flower presences in litter traps. However, that indirect reproductive metric initially increased and then leveled out for that study's final decade, a pattern more consonant with a reproductive response to forest recovery after a disturbance.

Drivers of the Twig-Litterfall IAV
Twig litterfall, which averaged 6% of La Selva's ANPP*, linearly decreased through the study. Changing rates of pre-collection mass loss could have played a role. In a New Guinea montane forest, Edwards (1977) found that twigs that fell from branches into litter traps were already highly decomposed, having lost 36%-40% of their live mass. Also, the twigs' content of non-structural carbohydrates could have changed over time, but was not assessed. The significant environmental relationships shown by annual twig litterfall (Table 2) are not readily interpretable. The negative relationship with [CO 2 ] is counter-intuitive, and the negative correlations with both [CO 2 ] and high-VPD conditions could simply reflect their opposite temporal patterns. The long-term decline in this minor ANPP* component remains unexplained.

Storm Damage Overwhelming the Environmental Responses of ANPP*
The microburst that hit La Selva in 2018 was an extreme wind-storm event (Frank et al., 2015). In the CARBONO plots, the storm caused levels of tree mortality and associated loss of potential wood production not seen in the prior two decades (Table S5 in D.A. . The environmental responses of annual wood production and reproductive litterfall maintained through the prior period (Table 2) were greatly disrupted in the storm year. Highly significant 20-year models became non-significant for the full 21 years. Some impacts of this storm on forest C-cycling (e.g., increased stocks of woody debris, tree growth reductions from major crown damage) likely continued through at least an additional year. Storms are expected to increase in frequency and intensity as warming intensifies (McDowell et al., 2020). Through the coming decades, such events could bring repeated re-sets of forest structure, C storage and other ecosystem processes at La Selva and across the biome. It will be important to factor such disturbances into models of future tropical-forest C-cycling.

Methods for Field-Assessing Tropical-Forest Climatic Sensitivities
The CARBONO Project provides useful guidance for future studies of the environmental responses of productivity in this biome. The field measurements should be designed to be at the landscape scale. The typical tropical-forest regime of sporadic, localized disturbances produced much more variable ANPP* histories for individual 0.5-ha CARBONO plots than seen with the forest-level averages. The strongest environmental model for landscape-scale forest productivity was non-significant for all but one of the 18 plots (Section 3.4.1). Long series of field observations (Table 5) were required to detect the forest responses with high statistical confidence. Annual censuses are thus the tool of choice. The multi-year census intervals typically used in tropical-forest studies (e.g., Brienen et al., 2015;Hubau et al., 2020;Qie et al., 2017;Rutishauser et al., 2020) significantly increase the time needed to accumulate a long series of observations. In addition, averaging over multi-year intervals greatly reduces the observed climatic variation and the range of the productivity metrics, thus hampering detection of climatic responses. Multi-factor modeling can be required in order to isolate each environmental factor's influences from the effects of uncorrelated IAV in other important drivers. As we found in this study (Tables 2 and 3), the strongest environmental responses can be undetectable in single-factor analyses. Finally, although much can be learned from long-term records of litterfall, using it as a surrogate for aboveground production introduces large uncertainties for research on the environmental drivers of productivity IAV.

Implications of Our Findings for Current Global C-Cycle Models
NPP should not be modeled as a homogeneous process. In the La Selva rainforest, the four components of ANPP* had distinct environmental responses (Tables 2-4). This complexity is to be expected, given the different types of tissues being produced. Leaf litterfall, the largest ANPP* component, stands out. Its lack of significant environmental models or of a long-term directional trend raises the possibility that annual leaf production in an old-growth evergreen forest is largely unresponsive to environmental IAV.
The productivity benefit from rising [CO 2 ] was notably limited for this tropical forest. While wood production and reproductive litterfall each increased 0.01 Mg ha −1 yr −1 per added ppmv of [CO 2 ] ( Table 3), most of La Selva's ANPP* showed no [CO 2 ] response. This muted benefit is particularly notable because La Selva is relatively nutrient-rich (Powers et al., 2005) and is thus expected to show greater "CO 2 -fertilization" benefits than in the more nutrient-limited majority of tropical forests (Fisher et al., 2012;Fleischer et al., 2019;Goll et al., 2012). No productivity component directionally increased over the 21 years. Negative impacts on La Selva's ANPP* from multiple climatic stressors (Tables 2 and 3) and from the 2018 storm more than offset the limited benefits from rising [CO 2 ]. These findings from a high-fertility tropical rainforest are reason to strenuously question model projections of sustained productivity increases for tropical forests.
Finally, our findings argue for refining the climatic metrics used in the global models. While the IAV of certain specific aspects of temperature and water availability produced highly significant impacts on forest productivity (Table 2), no ANPP* component responded to the IAV of either annual mean temperature or annual rainfall.

An Imminent Decline in Tropical-Forest Productivity?
This first multi-decade assessment of the effects of environmental factors on landscape-level tropical-forest productivity raises deep concerns for the biome's future. In the La Selva rainforest, the limited productivity benefits from rising [CO 2 ] are being overwhelmed by climatic stressors. Annual wood production is already being reduced by the current prevalence of critically hot daytime periods during the dry season. As noted by Huang et al. (2019), that tropical forests are already operating near their photosynthetically optimal temperature suggests "a limited safe operating space for these ecosystems under future warming." At the same time, the strongest relative response of La Selva's annual wood production was its decline with very small increases in yearly averages of nightime temperatures. Evidence that this negative temperature impact on productivity extends across the land tropics has been inferred from atmospheric inversion analyses (Anderegg et al., 2015). Annual wood production at La Selva is also already being reduced by current high-VPD periods, a finding echoed by eddy-covariance studies in other tropical forests (Fu et al., 2018;Tang et al., 2020;Vourlitis et al., 2005Vourlitis et al., , 2011. All these climatic stressors are intensifying with global warming. Temperatures exceeding the twentieth-century seasonal maxima are expected to occur across all tropical regions in 70%-90% of the years during this century (Diffenbaugh & Scherer, 2011). All these indications point to declining productivity across tropical forests in coming decades. This would constitute both a positive feedback to global warming and a major threat to these unique biodiverse ecosystems.