Continent-wide tree fecundity driven by indirect climate effects

Indirect climate effects on tree fecundity that come through variation in size and growth (climate-condition interactions) are not currently part of models used to predict future forests. Trends in species abundances predicted from meta-analyses and species distribution models will be misleading if they depend on the conditions of individuals. Here we find from a synthesis of tree species in North America that climate-condition interactions dominate responses through two pathways, i) effects of growth that depend on climate, and ii) effects of climate that depend on tree size. Because tree fecundity first increases and then declines with size, climate change that stimulates growth promotes a shift of small trees to more fecund sizes, but the opposite can be true for large sizes. Change the depresses growth also affects fecundity. We find a biogeographic divide, with these interactions reducing fecundity in the West and increasing it in the East. Continental-scale responses of these forests are thus driven largely by indirect effects, recommending management for climate change that considers multiple demographic rates.

T he composition and structure of twenty-first century forests will depend on the seed production needed for tree populations to keep pace with climate change. North America is warming and drying out in much of the West. The dramatic impacts include large-scale die-backs 1-3 that are transforming size-species structure 4,5 . But the decade-scale trends will depend on the regeneration that follows tree death. Fecundity will determine the capacity of trees to disperse seed to the shifting habitats where they can survive in the future [6][7][8] ; risks to each species depend not only on the current distribution of fecundity but also on its trajectory 4,[9][10][11][12][13] . As with many ecological processes [14][15][16] , noisy, spatially variable fecundity trends are hard to quantify 8,17 , but this is only the first problem. Attributing trends to environmental variables is complicated by individual size, growth, and resource access [18][19][20] . Conservation efforts must anticipate not just the direct climate effects on this trajectory but also the indirect effects as growth and changing size structure also affect fecundity. Because it has thus far been impossible to estimate at continental scales, fecundity is the only major demographic process that lacks field-based estimates in models of vegetation change 5,21,22 . To address these challenges, we built the continental Masting Inference and Forecasting (MASTIF) network of primary data (Fig. 1), and we developed trend attribution (TA) to quantify climate impacts, as modulated by the condition of the organisms themselves. Application to the MASTIF network shows that indirect effects dominate, operating through stand structure and growth.
Although a substantial climate-impacts literature has focused on growth responses to short-term (interannual to a decade) climate fluctuations 3,23 , this focus is not based on evidence that fecundity effects are of secondary importance. The emphasis on tree growth comes in part from the facts that (i) data are widely available from inventory plots and tree-ring records, (ii) where absent they can often be obtained from tree rings for long periods in the past, and (iii) growth anomalies can often be at least partly explained by climate anomalies 3,23 . By contrast, fecundity is not directly observed for most species and habitats, data accumulate slowly and with substantial investment, and the effects of climate anomalies can be overwhelmed by nonlinear, internal feedbacks on reproductive effort 18,22,[24][25][26][27] . Although critical for population dynamics, a limited role for fecundity could be interpreted from stand simulators that stabilize dynamics by assuming an external seed pool [28][29][30] . This is done both for lack of estimates, but also because the contribution of fecundity to dynamics is too poorly understood to construct models that allow species to coexist; even models used to explore effects on species diversity depend on the assumption that seeds remain available even when adults are not 30 . Foundational understanding of population growth makes clear that fecundity contributes directly to fitness [31][32][33] , while the change in size can do so only indirectly. While evidence points to the direct importance of fecundity for future forests, we show here that it must be coupled with the indirect role of tree growth.
We identify and quantify the effects of climate change on fecundity at the continental scale, including the climate-condition interaction (CCI) that require individual-scale observations (Fig. 2). We hypothesized that climate change will be experienced by organisms, each in its own way. We use the term CCIs to include, for example, moisture effects that differ for deep-rooted adults and small saplings 26,34 and temperature effects that depend on light availability 35 . If CCIs are important, then they must be quantified at the scale of individuals. We find biogeographic differences in the indirect effects of climate change, slowing fecundity change in the West and increasing it in the East.

Results
TA was developed to quantify change that emerges from both direct and indirect effects and that are not amenable to traditional time-series methods. Effects of climate change are increasingly apparent, including shifts in phenology 36 and species range limits 37,38 . By contrast, the time series of species-abundance data typically lack a clear signal 16,39,40 . No ecological process suffers more from the signal-to-noise problem than seed production, where quasiperiodic, order-of-magnitude variation from year to year and tree to tree 8,18,19,25,41,42 can bury long-term trends. Autoregression models assume a fixed periodicity, but mast intervals are not fixed, not even within an individual 25,35,43 . There are as many time series as there are trees (>10 5 in this case), but they must be modeled together because there is dependence. Data are non-Gaussian (including zeros for immature trees and failed crops). Trends estimated by meta-analysis may not be comparable across studies due to divergent methods and transformations intended to force non-Gaussian data into traditional time-series models 20 . Efforts to determine whether a species is increasing or decreasing are further challenged by the uneven distribution of publications. A standard trend analysis of our sites (Fig. 3a) shows not one trend but rather a broad range, with most species (bars in Fig. 3b) increasing in some habitats while decreasing in others. Estimates are readily biased 16 due to haphazard habitat coverage (Fig. 3a).
The MASTIF network includes the primary tree-year data (a given tree in a given year) that are needed to estimate change and the contributions of CCI. Data include the canopy environment (fully exposed to deep shade) and tree size, recognizing that accelerated growth can speed reproductive maturity 44 . Ecological studies have assumed that fecundity continues to increase with stem diameter (Fig. 2b, solid line) [45][46][47][48][49] . However, horticultural practice suggests declines in large trees (Fig. 2b, dashed line), but the literature is limited 50 . Several ecological studies also suggest declines 20,51-53 , but inference suffers from few observations of large trees. The MASTIF network offers a broad range of sizes combined here with weighted regression methods that allowed us to quantify the effects of both maturation and eventual fecundity declines (see "Methods" section). Fecundity data include seed traps (STs) and crop counts (CCs) from longitudinal studies (Fig. 3a) and opportunistically through the iNaturalist MASTIF [https://www.inaturalist.org/projects/mastif] project, including 2,566,594 tree-years from 123 species. The dynamic model accommodates non-Gaussian data and serial and intertree dependence, with full uncertainty for data, model miss-specification, and parameters 20 . Continental prediction used 7,723,671 trees from inventory plots (see "Methods" section).
TA was developed to decompose direct and indirect effects on change. For transparency, three climate variables in Fig. 2a are represented here by a single state variable C. To evaluate community-wide effects, we report on log (proportionate) fecundity change, where F is seeds per tree per year. Proportionate change df/dt is analyzed because we are interested in effects on species of both high and low fecundity. Analysis of absolute change dF/dt would be dominated by the few species that produce the most seeds. To obtain community change, we average these proportionate changes over all trees on a plot. TA entails (i) model fitting and (ii) trend decomposition. Model fitting estimates responses as fitted coefficients. These responses are main effects of climate, ∂f/∂C, and size (diameter G), ∂f/∂G, and their interaction ∂f/∂(GC). Trend decomposition combines these responses with dense information on the environment, including individual states (G, C) and their rates of (1). a Trends in climate variables (since 1990) include minimum T in spring, mean summer T, and moisture deficit (D = cumulative monthly PET-P). The brown contour separates positive and negative trends. Shaded contours are green (decreasing) to brown (increasing). b Indirect effects have two elements. An arrow from b to F includes a growth effect (dG/dt) and a climategrowth interaction (C × dG/dt). An arrow from b to c depends on the uncertain relationship between tree diameter G and fecundity F shown in panel (b). If fecundity continues to increase with tree size (solid line in b), then accelerated growth (orange arrows are dG/dt) moves trees into more productive size classes, but not if fecundity eventually declines (dashed line). c Average diameter of trees (restricted to trees >20 cm) is high in the West, meaning that the effects of tree growth depend on [if] fecundity continues to increase or declines with size in panel (b). The effect of size on fecundity (arrow from c to F) comes through an interaction with climate (G × dC/dt in Eq. (1)).
change (dG/dt, dC/dt) (Fig. 2), summarized with three terms, The first two terms are both main effects that depend respectively on rates of change in climate (dC/dt) and tree growth (dG/dt). The direct effect (first term in Eq. (1)) combines the climate response with the rate of climate change (Fig. 2a). This direct effect is followed by three terms that contribute the indirect effects of size and growth. The third term holds their interaction (GC) as products of rates (dG/dt, dC/dt) and states (G, C). [Again, C is a placeholder for multiple environmental variables (see "Methods" section)]. These are (i) the size-dependent effects of climate change and (ii) the climate-dependent effects of growth.
The residual γ allows effects that are not attributed to other terms. The full effect of a climate variable C combines the direct term 1 with its indirect effects in the second and third terms.
Indirect terms are CCI, incorporating climate effects that are modulated by tree size, G × dC/dt, as when large, deep-rooted trees experience drought differently from saplings. Conversely, a change in growth rate has effects that can vary with climate, also a CCI, C × dG/dt. The indirect effects through growth g(C) = dG/dt are not shown in Eq. (1), but are given in the "Methods" section. Depending on how fecundity changes with size ( Fig. 2b), climate stimulation of growth can move small trees into more fecund size classes. If fecundity eventually declines with size, growth has the opposite effect on large trees. TA in Eq. (1) starts from a notion similar to "climate velocity" 9,54 , which replaces fecundity in Eq. (1) with distance x as dx/dt = dx/dC × dC/dt. Rather than distance-over-time in climate velocity, TA decomposes the climate and size contributions to fecundity trends over time. It relies on extracting the smooth trends from volatile seed production data. The terms in Eq. (1) are available each as a predictive distribution for each tree. There is an average over trees for each inventory plot. The climate effects on fecundity trends differ from sensitivity to interannual variation 8,[17][18][19]25,35 . A species that reduces seed production in dry years (negative response to moisture deficit D j,t in Table S2.1) may not suffer from dry climates in general-indeed, the capacity to reallocate under fluctuating conditions can be adaptive. A negative effect in TA means that species and size classes of the current forests produce less seed under the decade-scale trends occurring now, based on responses across climate and habitat variation.
Indirect effects dominate response. TA shows that continentwide trends are dominated by indirect effects. Maps of these effect in Fig. 4 have different scales, which is necessary to show the geographic patterns within maps; the scale differences must be considered when comparing maps. The direct responses in Fig. 4a are transformed by the heterogeneity of climate trends (Fig. 4b) and then, indirectly, through tree growth (Fig. 4c). The responses in Fig It is important to recognize that a positive effect of climate change occurs wherever the response to climate and the direction of climate change have the same sign. For example, the negative direct effects of spring T in the Northeast (NE) and Southwest (SW) (Fig. 4b, top) result from opposing forces: in the NE, mostly positive responses (Fig. 4a, top) combine with a negative spring T trend (Fig. 2a), that is, (+) × (−) = (−). In the SW, negative responses combine with a positive spring T trend, that is, Between is a swath of positive effects stretching from the NW toward the SE (Fig. 4b, top), where positive responses overlap with rising spring T (Fig. 2a). The direct effects of other climate variables are near zero or negative for summer T (Fig. 4b, middle). The limited direct effects of moisture deficit D is apparent from the scale differences for maps in Fig. 4b.
The foregoing direct effects are overwhelmed by the indirect effects (contrast scales in Fig. 4b, c), both as main effects on growth in the second term and interactions in the third term of Eq. (1). Whereas the full effects contribute to a positive east/ negative west divide in the effects of spring T (Fig. 4c, top) and moisture deficit D (Fig. 4c, bottom), the contribution of summer a Longitudinal studies in black (opacity proportional to numbers of sites) and opportunistic in white. Shaded ecoregions are desert/shrub/grass (browns), montane (blues), and mixed forest (greens). b Trends in mean log (proportionate) fecundity by species from sites in a range from negative (declining) to positive. As would be the case for any meta-analysis, the time scales for which trends are evaluated vary (see "Methods" section). Species belong to color-coded families below that are listed in Supplementary Data 2. There is no relationship with phylogeny (i.e., no trend in box color from left to right). Summaries in b include mean (crosshairs), 95% of site means (bold line), and range of site means (whiskers). The number of sites (n) contributing to b is shown below.
T is primarily negative in the East (Fig. 4c, center). Positive effects in the East come predominantly through spring T (Fig. 4c, top), which is transparent because both responses and climate trends tend to share the same sign (both positive: Figs. 2a and 4a, top). The full effects could not have been anticipated from the direct responses because they require consideration of how growth responds to climate and the effect of size and growth on fecundity.
To understand continental responses and the large differences between maps in Fig. 4b, c requires the decomposition of effects, which is available through the individual terms in Eq. (1). An important contributor to these differences is the pervasive fecundity declines that our analysis found for large trees (Fig. 5a). Due to management and species traits, growth stimulation in the East speeds the transition of small trees to larger, more fecund size classes (Fig. 5b). Conversely, much of the West supports trees that have passed this size. The East-West differences are amplified by maturation, which is increasing the probability of seed production in the East, but not the West ( Supplementary  Fig. S2.3b). Declining fecundity in large trees (which are also older in the West, Supplementary Fig. S2.3a) does not necessarily come from physiological decline ("senescence") because [declines can result from] crown architectural changes also occur.
Discussion. Continent-wide impacts of climate change are being governed by the indirect effects that come through the condition of individuals. Global change science has steadily improved understanding of direct responses, including photosynthetic rates,  water use, and demography 25,55,56 and the abundances of species 40 . Lagging behind is the understanding of interactions and indirect effects 18 ; individual responses to climate do not predict responses of canopies 57 , just as species responses do not predict outcomes of competition 40 . However, climate variation operates on individual organisms that see that change differently. From the multitude of individual traits that affect climate response, we show that size and growth differences provide important insights on continent-wide change. By combining detailed habitat change with the individual condition (Fig. 4), TA shows the dominance of indirect effects.
The geographically coherent picture of change that emerges from TA contrasts with inconclusive results offered by metaanalysis of trends. For example, the conflicting interpretations from recent studies of insect abundances 16,58 reflect oversensitivity to precisely which sites and species were included in each meta-analysis 58 . Indeed, simple trend analysis of our sites (Fig. 3b) is no more informative than that of ref. 58 (their Fig. 2), both showing that nearly every species is increasing and decreasing somewhere. Only long time series can provide reliable estimates of trends for erratic data, but the duration is not enough -the unrepresentative geographic distribution of sites precludes an interpretation of overall trends. TA does not attempt to extract signal or extrapolate from noisy data, [instead of] exploiting instead relationships between varying climate and individual condition. It benefits from dense geographic coverage of sites, but can provide insight without it, relying primarily on adequate coverage in covariate space rather than geographic space.
TA adds value to existing efforts because climate change is heterogeneous not only in rate but also in sign (Fig. 2). By exposing the trends masked by interannual and intertree volatility of seed production while incorporating effects of additional variables TA provides much-needed perspectives on patterns and processes that affect individuals. Because climate variables interact with one another and with the individual condition (CCI), models need to not only find coefficients for their effects (Fig. 4a) but also combine them with the changes in climate that are happening now (Fig. 4c). Current understanding suggests that fecundity is also responding to variables that could not be incorporated into this analysis, including changing CO 2 44 , irradiance and clouds 59 , and soils, depending only on data availability and distribution across sites.
TA can complement efforts based on stand simulators and species distribution models by quantifying contemporary change and extracting the reasons for it. For instance, because stand simulators rely on immigration to achieve species coexistence, fecundity estimates are mostly absent and without direct or indirect effects of climate that are based on field data. TA combines growth with fecundity estimates at the tree-year scale for understanding biogeographic consequences, thus offering an alternative perspective to stand simulators and map-based predictions of future biodiversity 60 (Fig. 4c). The full effect differs from direct effects (Fig. 4b, c) due to indirect effects of climate on tree growth. Growth changes have limited impact on fecundity trends in the West because few trees are nearing maturity (Supplementary Fig. S2.2b), and fecundity has plateaued or is decreasing (Fig. 2c). By contrast, climate changes are accelerating change toward fecund size classes in the East.
The finding that fecundity can decline in large trees, with biogeographic consequences, does not diminish their contribution to biological diversity through microhabitats for wildlife 61 . Selective removals that promote uneven-aged structures can preserve microhabitats and promote canopy heterogeneity and light penetration that stimulates growth and fecundity 35 . Growth is not currently making a strong contribution to average trends in the West; however, management priorities can be guided by disaggregating these mean trends to understand their distribution across species at risk and/or valued for their ecosystem services.
The determination that indirect effects through individual condition can dominate biogeographic responses has immediate application in forestry and conservation. As an example, scientists and managers increasingly recognize that the challenges posed by continuing trends in climate cannot be addressed with traditional nursery practice or silvicultural treatments 62 . Managing for longterm trends (as opposed to the volatile interannual variation) must consider both the direct effect of climate-induced changes on growth and the indirect effects of these changes on fecundity. This knowledge is critical because size-species structure is often under the direct control of forest managers and conservation planners, especially in eastern North America 63 , whereas climate is not 4 . TA offers concrete estimates of how fast these changes are happening now and which variables are responsible. While climate is not controllable by managers, the opportunity to influence indirect effects through stand structure can foster stronger connections between conservation planning and global change science.

Methods
Elements of TA. Identifying biogeographic trends within volatile data required several innovations in the MASTIF model 20 , building from multivariate state-space methods in previous applications 41,52 . Standard modeling options, such as generalized linear models and their derivatives, do not accommodate key features of the masting processes. First, multiple data types are not independent. Maturation status is binary with detection error, CCs are non-negative integers, also with detection error, and STs require a transport model (dispersal) linking traps to trees, and identification error in seed identification. Of course, a tree observed to bear seed, now or in the past, is known to be mature now. However, failure to observe seed does not mean that an individual is immature because there are detection errors and failed crop years 41,64 .
Second, seed production is quasiperiodic within an individual (serial dependence), quasi-synchronous between individuals ("mast years"), [and] there is dependence on environmental variation, and massive variation within and between trees 41,53,65 . Autoregressive error structures (AR(p) for p lag terms) impose a rigid assumption of dependence that is not consistent with quasiperiodic variation that can drift between dominant cycles within the same individual over time 43 . It does not allow for individual differences in mast periodicity.
Third, climate variables that affect fecundity operate both through interannual anomalies over time and as [a] geographic variation. The masting literature deals almost exclusively with the former, but our application must identify the latter: the potentially smooth variation of climate effects across regions must be extracted from the many individual time series, each dominated by local "noise." Finally, model fitting is controlled by the size classes that dominate a given site and thus is insensitive to size classes that are poorly represented. Large trees are relatively rare in eastern forests, making it hard to identify potential declines in large, old individuals 41,53 . Conversely, the shade-intolerant species that dominate second-growth forests often lack the smaller size classes needed to estimate maturation and early stages where fecundity may be increasing rapidly.
Several of the foregoing challenges are resolved in the MASTIF model by introducing latent states for individual maturation status and tree-year seed production. The dependent data types (maturation status, CCs, STs) become conditionally independent in the hierarchical MASTIF model (e.g., ref. 66 ). The serial dependence is handled as a conditional hidden Markov process for maturation that combines with CCs and STs by way of stochastic (latent) conditional fecundity. Maturation status and conditional fecundity must be estimated jointly, that is, not with separate models. The latent maturation/fecundity treatment avoids imposing a specific AR(p) structure. In the MASTIF model there is a posterior covariance in maturation/fecundity across all tree-year estimates that need not adhere to any specific assumption 20 . The dependence across individuals and years is automatic and available from the posterior distribution.
Separating the spatial from temporal components of climate effects is possible here, not only because the entire network is analyzed together but also because predictors in the model include both climate norms for the individual sites and interannual anomalies across sites 35,52 . TA depends on both of these components.
Extracting the trends from volatile data further benefits from random individual effects for each tree and the combination of climate anomalies and year effects over time. A substantial literature focuses on specific combinations of climate variables that best explain year-to-year fecundity variation, including combinations of temperature, moisture, and water balance during specific seasons over current and previous years 19,25,41 . Results vary for each study, presumably due to the differences in sites, species, size classes, duration, data type, and modeling assumptions. For TA, the goal is to accommodate the local interannual variation to optimize identification of trends in space and time. Thus, we include a small selection of important climate anomalies (spring minimum T of the current year, summer T of the current and previous year, and moisture D of the current and previous year). The climate anomalies considered here do not include every variable combination that could be important for all size classes of every species on every site. For this reason, we combine climate anomalies with year effects. Year effects in the model are fixed effects within an ecoregion and random between ecoregions (ecoregions are shown in Fig. 2 and listed in Supplementary Data 2). They are fixed within an ecoregion because they are not interpreted as exchangeable and drawn at random from a large population of possible years. They are random between ecoregions due to the uneven distribution of sites (Supplementary Data 1) 20 .
To optimize inference on size effects, the sampling of coefficients in posterior simulation is implemented as a weighted regression. This means that the contribution of tree diameter to fecundity is inversely proportional to the abundance of that size class in the data. This approach has the effect of balancing the contributions of abundant and rare sizes. Identifying size effects further benefits from the introduction of opportunistic field sampling, which can target the large individuals that are typically absent from field study plots.
MASTIF data network. Data included in the analysis come from published and unpublished sources and offer one or both of the two data types, CCs and STs (Supplementary Data 1). Both data types inform tree-year fecundity; they are plotted by year in Fig. 6.
CCs in the MASTIF network are obtained by one of three methods. Most common are counts with binoculars that are recorded with an estimate of the fraction of the crop that was observed. A second CC method makes use of seeds collected per ground surface area relative to the crown area. This method is used where conspecific crowns are isolated and wind dispersal is limited. The crop fraction is the ratio of ground area for traps relative to the projected crown area. Examples include HNHR 67 and BCEF 68 .
A third CC method is based on evidence for past cone production that is preserved on trees. This has been used for Abies balsamea at western Quebec sites 69 , Pinus ponderosa in the Rocky Mountains 70 , and for Pinus edulis at SW sites 27 .
ST data include observations on individual trees that combine with seed counts from traps. Because individual studies can report different subcategories of seeds, and few conduct rigorous tests of viability, we had to combine them using the closest description to the concept of "viable". For example, we do not include empty conifer seeds. A dispersion model provides estimates of seeds derived from each tree. ST and CC studies are listed in Supplementary Data 1. The likelihoods for CCs and STs are detailed in ref. 20 . Individually and in combination, the two data types provide estimates, with full uncertainty, for fecundity across all treeyears.
Fitted species had multiple years of observations from multiple sites, which included 211,146 trees and 2,566,594 tree-years from 123 species. Sites are shown in Fig. 2 of the main text by ecoregion, they are named in Fig. 1 and summarized in Supplementary Data 1. For TA the fits were applied to 7,723,671 trees on inventory plots. Mean estimates for the genus were used for inventory trees belonging to species for which there were not confident fits in the MASTIF model, which amounted to 7.2% of inventory trees. Detailed site information is available at the website MASTIF.
Covariates. Covariates in the model include as main effects tree diameter, tree canopy class (shading), and the climate variables in Fig. 1 of the main text and described in Table 1. A quadratic diameter term in the MASTIF model allows for changes in diameter response with size 52 . Shade classes follow the USDA Forest Inventory and Analysis (FIA)/National Ecological Observation Network (NEON) scheme that ranges from a fully exposed canopy that does not interact with canopies of other trees to fully shaded in the understory. Shading provides information on competition that has proved highly significant for fecundity in previous analyses 41,52 .
To distinguish between the effects of spatial variation versus interannual variability, spring T and moisture D are included in the model as site means and site anomalies 35 . Spring minimum T affect phenology and frost risk during flowering and early fruit initiation. Summer mean T (June-August) is included both as a linear and quadratic term. Mean summer T is linked to thermal energy availability during the growing season, with the quadratic term allowing for potential suppression due to extreme heat. Moisture D (cumulative monthly PET-P (potential evapotranspiration[-] minus precipitation) for January-August) is included as a site mean and an annual anomaly. Moisture D is important for carbon assimilation and fruit development during summer in the eastern continent and, additionally, from the preceding winter in the western continent. For species that develop over spring and summer, anomalies incorporate the current and previous year. We did not include longer lags in covariates. For species that disperse seed in spring (Ulmus spp. and some members of Acer), only the previous year was used. Temperature anomalies were included for spring, but not summer, simply to reduce the number of times that temperature variables enter the model, and these two variables tended to be correlated at many sites.
Climate covariates were derived from gridded climate products and combined with local climate monitoring where it is available. Terraclimate 71 provides monthly resolution, but it is spatially coarse. For both norms and trends, we used the period from 1990 to 2019 because global temperatures have been increasing consistently since the 1980s, and this span broadly overlaps with fecundity data (Fig. 6). CHELSA 72 data are downscaled to a 1 km grid, but it does not extend to 2019. Our three-component climate scaling used regression to project CHELSA forward using Terraclimate, followed by downscaling to 1 km with CHELSA, with further downscaling to local climate data. Even where local climate data exist, they often do not span the full duration of field studies, making the link to gridded climate data important. Local climate data were especially important for mountainous sites in the Appalachians, Rockies, Sierra Nevada, and Cascades.
Of the full list of variables, a subset was retained, depending on species (some have narrow geographic ranges) and deviance information criteria of the fitted model (Supplementary Data 2). Year effects in the model allow for the interannual variation that is not absorbed by anomalies 20 .
Model fitting and TA. As mentioned above, model fitting applied the hierarchical Bayes model of ref. 20  The climate variable referenced as C in Eq. (1) of the main text is, in fact, a vector of climate variables described in the previous section, spring minimum T, summer mean T, and moisture D ( Table 1). The anomalies and year effects in the   fitted model contribute to the trends not explained by biogeographic variation as γ in Eq. (1). For main effects in the model, the partial derivatives are fitted coefficients, an example being the response to spring minimum temperature ∂f =∂T sp ¼ β T sp . For predictors involved in interactions, the partial derivatives are combinations of fitted coefficients and variables. For example, the response to moisture D, which interacts with tree size, is ∂½F f =∂D ¼ β D þ β GD G. The response to diameter G, which is quadratic and interacts with D, Trend decomposition applied the fitted model to every tree in forest inventories from the USDA FIA program, the Canada's National Forest Inventory, the NEON, and our MASTIF collaboration. Each tree in these inventories has a species and diameter. For trees that lack a canopy class, regression was used to predict it from distances and tree diameters based on inventories that include both location and canopy class, including NEON, FIA, and the MASTIF network. Although inventories differ in the minimum diameter they record, few trees are reproductive at diameters below the lower diameter limits in these surveys, so the effect on fecundity estimates is negligible. For the indirect effects of climate coming through tree growth rates, the same covariates were fitted to growth as previously defined for fecundity, using the change in diameter observed over multiple inventories. A Tobit model was used to accommodate the fact that a second measurement can be smaller than an earlier measurement. The Tobit thus treats negative growth as censored at zero. TA to inventory plots used 7,717,677 trees. Because not all species in the inventory data are included in the MASTIF network, mean fecundity parameters for the genus were used for unfitted species. Species fitted in the MASTIF network accounted for >90% of trees in inventory plots (Supplementary Data 2).
From the predictive distributions for every tree in the inventory data, we evaluated predictive mean trends aggregated to species and plot in Fig. 2b. We extracted specific terms that comprise the components in Fig. 4 and aggregated them too to the plot averages.
General form for TA. Equation 1 simplifies the model to highlight direct and indirect effects. Again, climate variables and tree size represent only a subset of the predictors in the model that are collected in a design vector x t ¼ ½x 1;t ; ; x Q;t 0 , where the q = 1, …, Q predictors include shading from local competition, individual size, and climate and habitat variables (Table 1). On the proportionate scale, Eq. (1) can be written in terms of all predictors, including main effects and interactions, as where I q are variables that interact with x q . In this application, interactions include tree diameter with moisture deficit and diameter squared. Each term in the summation consists of a main effect of x q and interactions that are multiplied by the rate of change in variable x q . For the simple case of only two predictors, Eq. (2) is recognizable as Eq. (1) of the main text, where x 1 , x 2 have been substituted for variables G and C. In our application, predictors include additional climate and shading ( Table 1).
Recognizing that environmental variables affect not only fecundity but also growth rate, we extract the size effect, that is, the x q that is G, and incorporate these indirect effects (through growth) by expanding g = dG/dt in Eq. (1) of the main text as where ν is the component of growth that is not accommodated by other terms. This expression allows us to evaluate the full effect of climate variables, including those coming indirectly through growth.
Connecting fitted coefficients in MASTIF to TA. This section connects the continuous, deterministic Eq. (1) to the MASTIF model 20 with the interpretation of responses, direct effects, and full effects of Fig. 5. To summarize key elements of the fitted model 20 , consider a tree i at site j that grows to reproductive maturity and then produces seed depending on its size, local competitive environment, and climate. We wish to estimate the effects of its changing environment and condition on fecundity using a model that includes spatial variation in predictors that are tracked longitudinally over years t. Fecundity changes through maturation probability ρ ij (t), which increases as trees increase in size, and through conditional fecundity ψ ij (t), the annual seed production of a mature tree. Let z ij (t) = 1 be the event that a randomly selected tree i is mature in year t. Then, ρ ij (t) is the corresponding probability that the tree is mature, E[z ij (t)] = ρ ij (t)(ρ is not to be confused with the probability that a tree that is now immature will make the transition to the mature state in an interval dt = 1. That is a different quantity detailed in the Supplement to ref. 41 ). Fecundity has expected value F ij (t) = ρ ij (t) ψ ij (t). On a proportionate (log) scale, The corresponding rate equation is The discretized and stochasticized version of Eq. (1) is where dt = 1 and ϵ ij,t is the integration error. When applied to a dynamic process model, this term further absorbs process error (see above), which is critical here to allow for conditional independence where observations are serially dependent. In simplest terms, ϵ is model miss-specification that allows for dependence in data. The MASTIF model that provides estimates for TA is detailed in ref. 20 . Elements of central interest for TA are F ij;t ¼ z ij;t ψ ij;t z ij;t ¼ 1 where μ ij,t = α 0 + α G G ij,t describes the increase in maturation probability with size, Φ(⋅) is the standard normal distribution function (a probit), ϵ ij,t~N (0, σ 2 ), and h t (T) can include year effects, h(T) = κ t , or lagged effects, hðTÞ ¼ P p k¼1 κ k ψ ij;tÀk , that contribute to γ in Eq. (1) of the main text. If year effects are used, then γ includes the trend in year effects. (The generative version of this model writes individual states at t conditional on t − 1 and is given in the Supplement to ref. 20 .). If an AR(p) model is used, then γ = κ 1 (provided data are not detrended). Random individual effects in the fitted model are marginalized for prediction of trees that were not fitted, meaning that σ 2 is the sum of model residual and random-effects variance. Again, the length-Q design vector x ij,t includes individual attributes (e.g., diameter G ij,t ), local competitive environment, and climate (Table 1). There is a corresponding coefficient vector β.
Elements of basic theory in Eq. (1) of the main text are linked to data through the modeling framework as where ϕ(⋅) is the standard normal density function that comes from the rate of progress toward maturation. Again, the anomalies do not appear in this expression for trends because trends in the anomalies and year effects enter through γ.
The first four lines in Eq. (7) are, respectively, the effects of trends in spring minimum temperatures ΔT sp,j , summer mean temperature ΔT j , moisture deficit ΔD j , and size ΔG ij , where the latter comes from growth on inventory plots. The contribution of maturation to change in fecundity is the first term in the fourth line, α G ϕ(μ ij,t )/Φ(μ ij,t ). A map of this term in Fig. 7b shows the strong contribution to fecundity in the East due to the young (Fig. 7a) and/or small (Fig. 4b) trees there. The sum of these terms dominates the patterns in Fig. 3c.
All terms in Eq. (7) have units of mean change in proportionate fecundity, and these are mapped in figures of the main text. We focus on proportionate fecundity because it reflects the full effect of climate as opposed to total fecundity, which would often be dominated by one or a few trees of a single species. However, from proportionate fecundity we can obtain change in fecundity as ΔF ij,t = F ij,t × Δf ij . Stand-level effects on fecundity change at site j can be obtained from individual change as