# Global patterns of tree density are contingent upon local determinants in the world's natural forests

Jaime Madrigal-González, Joaquín Calatayud, Juan Ballesteros-Cánovas, Adrián Escudero, Luis Cayuela, Laura Marqués, Marta Rueda, Paloma Ruiz-Benito, Asier Herrero, Cristina Aponte, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jaime Madrigal-González, Joaquín Calatayud, Juan Ballesteros-Cánovas, Adrián Escudero, Luis Cayuela, et al. Global patterns of tree density are contingent upon local determinants in the world's natural forests. Communications Biology, 2023, 6 (1), pp.1-6. 10.1038/s42003-023-04419-8. hal-03942333

HAL Id: hal-03942333

https://hal.science/hal-03942333

Submitted on 17 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Global patterns of tree density are contingent upon local determinants in the world's natural forests

Jaime Madrigal-González<sup>1,2✉</sup>, Joaquín Calatayud<sup>1</sup>,<sup>3</sup> Juan A. Ballesteros-Cánovas<sup>1,4</sup>, Adrián Escudero<sup>5</sup>, Luis Cayuela<sup>5</sup>, Laura Marqués<sup>6</sup>, Marta Rueda<sup>7</sup>, Paloma Ruiz-Benito<sup>8</sup>, Asier Herrero<sup>9</sup>, Cristina Aponte<sup>10,11</sup>, Rodrigo Sagardia<sup>12</sup>, Andrew J. Plumptre<sup>13</sup>, Sylvain Dupire<sup>14</sup>, Carlos I. Espinosa<sup>15</sup>, Olga V. Tutubalina<sup>16</sup>, Moe Myint<sup>1</sup>, Luciano Pataro<sup>17</sup>, Jerome López-Sáez<sup>1</sup>, Manuel J. Macía<sup>17,18</sup>, Meinrad Abegg<sup>19</sup>, Miguel A. Zavala<sup>8,20</sup>, Adolfo Quesada-Román<sup>1,21</sup>, Mauricio Vega-Araya<sup>22</sup>, Elena Golubeva<sup>23</sup>, Yuliya Timokhina<sup>23</sup>, Guillermo Bañares de Dios<sup>3</sup>, Íñigo Granzow-de la Cerda<sup>24</sup> & Markus Stoffel<sup>1,25,26</sup>

Previous attempts to quantify tree abundance at global scale have largely neglected the role of local competition in modulating the influence of climate and soils on tree density. Here, we evaluated whether mean tree size in the world's natural forests alters the effect of global productivity on tree density. In doing so, we gathered a vast set of forest inventories including >3000 sampling plots from 23 well-conserved areas worldwide to encompass (as much as possible) the main forest biomes on Earth. We evidence that latitudinal productivity patterns of tree density become evident as large trees become dominant. Global estimates of tree abundance should, therefore, consider dependencies of latitudinal sources of variability on local biotic influences to avoid underestimating the number of trees on Earth and to properly evaluate the functional and social consequences.

<footnote n="1">¹Climate Change Impacts and Risks in the Anthropocene (C-CIA), Institute for Environmental Sciences (ISE), University of Geneva, 66 Boulevard Carl Vogt, CH-1205 Geneva, Switzerland. ²EiFAB-iuFOR, Universidad de Valladolid, Campus Duques de Soria s/n, 42004 Soria, Spain. ³Departamento de Biología y Geología, Física y Química Inorgánica. ESCET, Universidad Rey Juan Carlos, C/Tulipán s/n, Móstoles, C.P. 28933 Madrid, Spain. ⁴Spanish Scientific Research Council, Museo Nacional de Ciencias Naturales, Madrid, Spain. ⁵Grupo de Investigación de Alto Rendimiento en Ecología de Comunidades, Departamento de Biología y Geología, Física y Química Inorgánica, ESCET, Universidad Rey Juan Carlos, C/Tulipán s/n, Móstoles, C.P. 28933 Madrid, Spain. ⁶Department of Environmental Systems Science, Swiss Federal Institute of Technology (ETH Zürich), Universitätstrasse 2, 8092 Zürich, Switzerland. ⁷Departamento de Biología Vegetal y Ecología, Universidad de Sevilla, C/Profesor García González s/n, 41012 Sevilla, Spain. ⁸Forest Ecology and Restoration Research Group, Departamento de Ciencias de la Vida, Universidad de Alcalá, ctra., Madrid-Barcelona, km 33.4, 28805 Alcalá de Henares, Spain. ⁹FisioClima CO2 Research Group, Department of Plant Biology and Ecology, Faculty of Pharmacy, University of the Basque Country, 01006 Vitoria-Gasteiz, Basque Country, Spain.</footnote>

<footnote n="10">¹⁰School of Ecosystem and Forest Sciences, The University of Melbourne, 500 Yarra Boulevard, Richmond, Victoria 3121, Australia. ¹¹Department of Environment and Agronomy, Centro Nacional Instituto de Investigación y Tecnología Agraria y Alimentaria, INIA-CSIC, Ctra. Coruña Km 7.5, 28040 Madrid, Spain. ¹²Instituto Forestal de Chile, Sucre 2397, Ñuñoa, Santiago de Chile, Chile. ¹³Key Biodiversity Area Secretariat, c/o BirdLife International, Cambridge, UK. ¹⁴Université Grenoble Alpes, Inrae, LESSEM, 38000 Grenoble, France. ¹⁵EcoSs_Lab, Departamento de Ciencias Biológicas, Universidad Técnica Particular de Loja, San Cayetano Alto, 110107 Loja, Ecuador. ¹⁶Scott Polar Research Institute, University of Cambridge, Lensfield Road, Cambridge CB2 1ER, UK. ¹⁷Departamento de Biología (Botánica), Facultad de Ciencias, Universidad Autónoma de Madrid, calle Darwin 2, Madrid, Spain. ¹⁸Centro de Investigación en Biodiversidad y Cambio Global (CIBC-UAM), Universidad Autónoma de Madrid, Calle Darwin 2, ES-28049 Madrid, Spain. ¹⁹Swiss Federal Institute for Forest, Snow and Landscape Research, WSL, Zürcherstrasse 111, 8903 Birmensdorf, Switzerland. ²⁰Instituto Franklin, Universidad de Alcalá, Calle Trinidad 1, Alcalá de Henares, 28801 Madrid, Spain. ²¹Escuela de Geografía, Facultad de Ciencias Sociales, Universidad de Costa Rica, Ciudad de la Investigación, Montes de Oca, 2060 San José, Costa Rica. ²²Instituto de Investigación y Servicios Forestales (INISEFOR), Universidad Nacional de Costa Rica, 86-3000 Heredia, Costa Rica. ²³Faculty of Geography, Lomonosov Moscow State University, Moscow, Russia. ²⁴Real Jardín Botánico - CSIC, Plaza Murillo 2, Madrid ES-28014, Spain. ²⁵Department of Earth Sciences, University of Geneva, 13 rue des Maraîchers, CH-1205 Geneva, Switzerland.</footnote>

<footnote n="26">²⁶Department F.-A. Forel for Environmental and Aquatic Sciences, University of Geneva, 66 Boulevard Carl Vogt, CH-1205 Geneva, Switzerland.</footnote>

<footnote>✉email: jaime.madrigal@uva.es</footnote>
The most recent assessment of global tree abundance suggests that more than one trillion trees inhabit planet Earth¹. These trees represent a massive stock of organic carbon and offer critical support for millions of species including animals, plants, fungi, lichens, and bacteria. Expressed as tree density (i.e., the number of trees per unit area), abundance represents a major structural component of natural forests linked to ecosystem functioning and energetics¹⁻³. Moreover, tree abundance is a major component of diversity and can have a direct contribution to population viability and species richness in natural forests under limiting climatic conditions⁴. Thus, improving our understanding of the drivers of tree density is imperative for Earth science, ecology and conservation of forest ecosystems.

Climate is among the most conspicuous determinants of tree density. In fact, correlational evidence supports a positive, combined influence of rising temperatures (i.e., more energy) and increased water availability on the number of assembled trees per unit area⁵. Whereas evidence of this influence is hitherto inconclusive in certain biogeographical regions⁶, soil fertility has been reported to be positively correlated with maximum stand density in other regions⁷. At fine scales, tree density is negatively correlated with tree size, following some density-size rules including the Yoda's law⁸. Specifically, this density-size rule implicitly depicts the critical role of competition in driving tree dynamics at the forest stand level through self-thinning constraints based on saturation of light demands by tree canopies over the course of secondary succession. It is well known that self-thinning dynamics in forests leads to a maximum carrying capacity for a given site⁹. Consequently, any biomass loss resulting from mortality of individual trees is eventually compensated by secondary growth in the remaining trees¹⁰. Importantly, the determinants of tree density have commonly been evaluated as independent drivers across contrasting spatial and temporal scales, such that interactive biotic influences have barely been addressed in the literature. Recent findings based on time-series of tree mortality data in a tropical forest highlighted the tight interplay between interannual climate variability and competitive interactions driving tree density¹¹. Yet, it remains poorly understood whether such controls of tree density can also be recognized in the well-known, large-scale productivity patterns of tree density through competition and self-thinning dynamics operating over a set of different pools of species and biogeographical contexts.

Here, we unveil the interdependency of the main determinants of tree density at local and global extents using forest inventory data covering a large network of well-conserved forests worldwide [Fig. 1]. We hypothesize that, beyond the direct contribution of primary productivity related to well-known local and latitudinal determinants on tree density, fine-scale competitive dynamics would significantly modulate tree density patterns. We anticipate that self-thinning would be less intense in the most productive tree assemblages (i.e., tropical forests) where competition, as a major driver of natural selection in evolutionary time, promotes niche segregation¹² which in turn results in a more efficient filling of ecological space above and below the ground¹³. On the contrary, we posit that self-thinning dynamics would reduce tree density more drastically due to niche overlap among tree individuals in the most stressed, less productive environments, where environmental constraints strongly sort species with similar or convergent functional adaptations¹⁴. An abiotic filtering in such limiting conditions would reduce the number of co-existing species, which would in turn trigger within and between species competitiveness due to a less efficient use of nutrients at the stand level¹⁵, ¹⁶. If true, a latitudinal pattern of tree density would be more evident for mature forests composed of large trees because maximum forest stand density is contingent on climate conditions¹⁷.

## Results and discussion

**Determinants of tree density at global to local scales.** We fitted a general additive mixed model (GAMM) to test for dependencies of the latitude-tree density relationship on plot-level mean diameter at breast height (DBH). Model selection using the Bayesian Information Criterion (BIC) supported the interaction between latitude and mean DBH (BIC<sub>interaction</sub> = 44,542.87; BIC<sub>no interaction</sub> = 44,575.27; BIC<sub>null</sub> = 44,748.20), indicating that, across the 23 studied regions, the highest tree densities are reached in forest plots dominated by small trees irrespective of latitude [Fig. 2a], and that the latitudinal gradient of tree density is more evident towards higher mean DBH values. Interestingly, plot size was not supported as a plausible determinant of tree density (BIC<sub>with plot size</sub> = 44,558; BIC<sub>without plot size</sub> = 44,542.87). Different combinations of irradiance through climatic net primary productivity (NPP) with certain values of soil cation exchange capacity (CEC) a surrogate of soil fertility, induce peaks of maximum density at different latitudes. Results of a GAMM in which NPP and CEC were evaluated as a function of latitude

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Geographic distribution of the 23 studied forests. This map (image and every element) was fully created by the authors using averaged coordinates of the studied forest regions and forest cover data (public domain) retrieved in the FAO's website³⁴.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Latitude and mean DBH as predictors of tree density. a Colored contour biplot for predictions of a generalized additive mixed model (GAMM) in which tree density (see the colored scale on the right) is expressed as a function of mean DBH (diameter at breast height) and latitude. Gray points are the data used in the GAMM. b predicted smoothing values of climatic Net Primary Productivity (NPP, blue lines) and Cation Exchange Capacity (CEC, red lines) as function of latitude. Solid lines represent the fitted mean values whereas dotted lines represent the 95% confidence intervals around the mean.</figcaption></figure>

separately support the idea that, while NPP can be critical driving tree density in tropical areas, CEC can have a prevalent role driving tree density in temperate and boreal forest biomes [Fig. 2b]. These findings are in line with results of previous studies reporting climatic controls on tree abundance at global¹ and regional¹⁸ scales. In general terms, the combination of moderate temperature and high precipitation seems to maximize the potential number of trees per unit area¹. Nonetheless, the interplay between soil fertility and climatic NPP seems to be critical to unveil this pattern. Results of a structural equation model (SEM), defined to evaluate determinants of tree density from latitudinal sources of variability to fine grain biotic factors, evidence opposing patterns of climatic NPP and CEC. Whereas NPP is negatively associated with latitude, the opposite is true for CEC [Fig. 3]. Interestingly, both climatic NPP and CEC had positive and rather similar influences on tree density. Thus, moderate to high limiting climatic conditions in cold temperate and boreal regions could compensate for soil fertility, likely promoting elevated tree densities at high latitudes as well.

**Latitudinal patterns of tree density are contingent on plot-level mean tree size.** We also found that mean tree size affected tree density negatively, as stipulated by the Yoda's law, a pattern that has long been supported in forest science¹⁹. Importantly, the Bayesian information criterion supported the interactive effect between local (Mean DBH) and global (primary productivity) determinants in a linear mixed-effects model (LMM) (BIC<sub>interaction</sub> = 6980; BIC<sub>no interaction</sub> = 6987) [Fig. 4], yielding a model that explained a considerable proportion of forest density variability

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Global and local determinants of tree density. Results of the structural equation model (SEM) fitted to tree density values taking into consideration mean diameter at breast height (MDBH), Net Primary Productivity (NPPc), Cation Exchange Capacity (CEC), and plot size, an exogenous variable included as a statistical controlling factor. Both NPP and CEC depends on latitude and are considered global factors. Solid arrows represent significant effects while dotted arrows denote non-significant relationships. Red arrows represent negative effects whereas black arrows do so for positive effects. All estimates are included in the figure as standardized estimates.</figcaption></figure>

(conditional $R^2$ = 0.55). This interaction suggests that plots can only support a high density of large trees when productivity is high and there is a decline in density as productivity decreases. In the high-productivity plots, the number of trees per hectare are rather similar irrespective of climate/soil conditions (≈500–800 trees/ha on average). Hence, if forests mature and large trees become dominant, the self-thinning dynamics seemingly operate but with different intensity depending on where trees grow along the global productivity gradient. Specifically, self-thinning dynamics seem to be more intense under severely limiting conditions than in productive environments. Consequently, the latitudinal gradient of tree density was particularly conspicuous when considering plots dominated by large trees. This means that tree density was significantly higher towards highly productive forests at low latitudes, namely tropical and equatorial rainforests. At the other end of the latitudinal gradient, in the boreal biome, tree density yielded minimum values in plots dominated by large trees. In addition, tree density was also relatively low in seasonal dry tropical and subtropical forests, where water scarcity represents a major limiting factor to tree establishment and growth. Thus, self-thinning dynamics seem to be particularly strong in areas where climate and/or soil severity limit tree performance. In these regions, abiotic filtering can result in a functional clustering tied to a selection of functional traits conferring the ability to thrive with specific environmental constraints¹⁴. Such a high functional overlap can be expected to result in more intense and
<figure><img src="image_7.png" /><figcaption>**Fig. 4** The productivity-tree density relationship is contingent upon mean DBH. Graphical representation of tree density as a function of the interaction between mean tree size (Mean DBH) and productivity calculated as 0.232NPP + 0.237CEC (being NPP the climatic Net Primary Productivity, and CEC the Cation Exchange Capacity). Shadowed areas represent the 95% confidence intervals of model predictions. The right panels display the distributions of mean DBH and productivity, respectively. Results support the idea that the positive influence of productivity on tree density reinforces towards higher values of mean tree size.</figcaption></figure>

competitive interactions with the course of forest development and fine-scale (competitiveness) dynamics¹⁵. Likewise, temperate forests were reported to undergo important self-thinning events during forest stand development and succession according to density-size relationships²⁰. In tropical forests, self-thinning dynamics have been observed to affect tree density along fine-scale environmental gradients depending on dry season length and frequency of natural disturbances²¹. Nonetheless, tree density in tropical rainforests tends to be high even for big trees, which might be justified in part by high productivity²² but also by the huge functional trait dispersion, likely forced in evolutionary time by competition to induce niche segregation and natural selection. High niche segregation has been reported in these type-forests in relation to light and water niche axes²³.

**Implications and limitations.** Our study underlines the crucial role played by climate and soil fertility on global gradients of tree density in natural forests but also stipulates that global tree density is significantly modulated by competitiveness. Although a simplified pool of potential drivers of abundance was used, our findings point to the important role of fine-scale determinants, such as mean tree size, in modulating global patterns of tree density associated with climatic and soil productivity. Further research should focus on disentangling more specific and complex interactions between latitudinal sources of variability and local context variables so as to elucidate mechanisms through which biotic interactions (i.e., competition) can alter global tree density patterns in more detail. For instance, dealing separately with precipitation and soil nutrients might help to unveil the actual links between water availability and nutrients with competition and forest dynamics along the latitudinal gradient. Whereas water scarcity can play a major role in dry forests, nutrient availability becomes a key driver in tropical rainforests, where abundant precipitation is responsible for the leaching loss of nutrients. Similarly, potential signals of local disturbances and historical contingencies, but also methodological choices (e.g., the inclusion of trees with a DBH ≥

10 cm only) deserve future attention. For the moment, even if the data used in this study cover most of the size spectrum, and DBH cutoffs represent standardized methods in the national forest inventories (NFI), further efforts should be made to extend monitoring protocols to trees smaller than 10 cm DBH (as is the case in some NFI). Whereas these sources of uncertainty might potentially blur signals of tree density patterns, drivers used in this study accounted for a considerable portion of variability in tree density, and thus support the importance of the interplay between productivity and competition.

Modulation of global tree density by mean tree sizes will have major implications on the way carbon is stored in forests through the formation of different size structures under specific climatic and/or soil productivity gradients. Thus, our findings also call for such dependencies to be considered more thoroughly in global assessments of tree density and the ensuing functional consequences related to forest functions and the provisioning of services to human societies.

**Materials and methods**

**Forest inventory data and selection criteria.** We used information on tree density recorded in more than 3000 well-conserved forest plots distributed across protected natural areas in 23 regions worldwide. Specifically, data were obtained from North (United States), Central (Costa Rica), and South (Ecuador, Brazil, Bolivia, Peru, and Chile) America, Africa (Uganda), Oceania (Australia), Asia (eastern Russia, Bhutan, and Myanmar), and Europe (Sweden, Switzerland, France, and Spain) (see Table 1 for details). All forest plots used here were retrieved from national forest inventories or obtained from research projects. Sampling plots can be either circular or rectangular (depending on the region considered) with identical sizes within each forest region. We only considered forest regions for which a minimum of 30 forest plots were available to ensure robustness in the statistical analysis. To exclude anthropogenic impacts, we only considered natural forests within protected areas lacking evidence of recent disturbance. Further information on the database can be found in Madrigal-González et al.⁴.

**Biotic and abiotic data.** Every tree exceeding 10 cm in DBH from each sampling plot was identified at the species level or classified as a morphospecies (in the case of Bhutan, Myanmar, Ecuador, and Peru). For each region, tree density was estimated as the number of standing trees per hectare and mean tree size as plot-level
<table><caption>Table 1 Descriptive geographic and sampling information for each of the 23 forest regions.</caption><thead><tr><th>Country</th><th>Region</th><th>Longitude (°)</th><th>Latitude (°)</th><th>Elevation (m.a.s.l.)</th><th>No. Plots</th><th>Plot Size (m²)</th></tr></thead><tbody><tr><td>Russia</td><td>Kola</td><td>39.19</td><td>67.58</td><td>199.045</td><td>28</td><td>400</td></tr><tr><td>US</td><td>Sequoia National Park</td><td>-118.41</td><td>36.20</td><td>2274.69</td><td>132</td><td>168.33</td></tr><tr><td>US</td><td>Great Canyon National Park</td><td>-112.12</td><td>36.26</td><td>2250.94</td><td>68</td><td>168.33</td></tr><tr><td>Sweden</td><td>Northern Sweden</td><td>22.31</td><td>66.17</td><td>141.40</td><td>101</td><td>1256</td></tr><tr><td>Spain</td><td>Sierra Nevada National Park</td><td>-3.72</td><td>38.10</td><td>1733.92</td><td>56</td><td>1963.49</td></tr><tr><td>Switzerland</td><td>Alps</td><td>8.89</td><td>46.39</td><td>1496.74</td><td>234</td><td>500.34</td></tr><tr><td>Bhutan</td><td>Toepisa</td><td>89.78</td><td>27.52</td><td>2371.11</td><td>160</td><td>500.34</td></tr><tr><td>US</td><td>Alaska</td><td>-133.16</td><td>56.30</td><td>874.33</td><td>491</td><td>168.33</td></tr><tr><td>US</td><td>New York</td><td>-74.47</td><td>43.66</td><td>1815.22</td><td>197</td><td>168.33</td></tr><tr><td>Brazil</td><td>Bahia</td><td>-40.90</td><td>-9.99</td><td>884.37</td><td>106</td><td>400</td></tr><tr><td>Ecuador</td><td>Western Ecuador</td><td>-80.19</td><td>-4.22</td><td>580.62</td><td>48</td><td>400</td></tr><tr><td>Australia</td><td>Victoria</td><td>146.80</td><td>-37.13</td><td>468.11</td><td>35</td><td>900</td></tr><tr><td>France</td><td>Mercantour National Park</td><td>7.12</td><td>44.13</td><td>1806</td><td>61</td><td>706</td></tr><tr><td>Chile</td><td>Northern Patagonia</td><td>-72.63</td><td>-41.59</td><td>361.22</td><td>109</td><td>500</td></tr><tr><td>France</td><td>Cévennes National Park</td><td>3.65</td><td>44.26</td><td>1123.82</td><td>98</td><td>706</td></tr><tr><td>Spain</td><td>Fuentes Carrionas Natural Park</td><td>-4.54</td><td>42.93</td><td>1316</td><td>117</td><td>1963.49</td></tr><tr><td>US</td><td>Klamath Forest</td><td>-123.87</td><td>41.78</td><td>787.43</td><td>74</td><td>168.33</td></tr><tr><td>Ecuador</td><td>Podocarpus National Park</td><td>-79.01</td><td>-4.12</td><td>1959</td><td>30</td><td>1000</td></tr><tr><td>Peru</td><td>Río Abiseo National Park</td><td>-77.38</td><td>-7.64</td><td>1948.80</td><td>30</td><td>1000</td></tr><tr><td>Myanmar</td><td>Wetphuyay</td><td>96.54</td><td>20.72</td><td>1095.66</td><td>62</td><td>500.34</td></tr><tr><td>Uganda</td><td>Forest Reserve</td><td>31.54</td><td>1.76</td><td>1042.00</td><td>622</td><td>100</td></tr><tr><td>Bolivia</td><td>Madidi National Park</td><td>-67.94</td><td>-14.21</td><td>656.25</td><td>44</td><td>1000</td></tr><tr><td>Costa Rica</td><td>Costa Rica</td><td>-84.15</td><td>9.87</td><td>933.97</td><td>44</td><td>1000</td></tr></tbody></table>

Values for longitude, latitude, elevation, number of plots and plot size are shown next to names of regions and countries.

mean DBH. We interpret the negative relationship between mean DBH and tree density as a static realization of self-thinning dynamics over time in line with expectations of Yoda's law⁸. Geographical coordinates were obtained for each plot using global positioning systems (GPS). The climatological net primary productivity (NPP) index²⁴ was retrieved as a nonlinear combination of temperature and precipitation following the equations of the Miami model²⁵ (see(1) for details) where NPP increases with rising temperature and precipitation up to a saturation of 3000 g dry matter m⁻² year⁻¹. Soil cation exchange capacity (CEC, cmolc/kg) was used here as a surrogate of soil fertility and retrieved from the SoilGrids system, which provides global gridded soil information at 250 m resolution accounting for climate, land cover and topography variability²⁶. We calculated an all-inclusive productivity variable combining NPP and CEC proportionally to their estimated effects on tree density in the structural equation model (0.232, 0.237 respectively; unstandardized coefficients; see statistical analyses below).

**Statistics and reproducibility.** To explore the potential contribution of global and local drivers to tree density, we first fitted a generalized additive mixed model (GAMM²⁷) where tree density was the response variable and mean tree size and absolute latitude were included as potential covariates. GAMM is a flexible method that allows modeling complex non-linear relationships of a response variable and one or several covariates. To account for site-to-site variability, a nested random factor of plots within tracks within regions was included in the model. Tracks are spatial clustering of sampling plots within regions. For instance, in the US, sampling plots are spatially clustered in groups of four plots around each fixed point of the grid considered across the country.

We then applied a structural equation model (SEM²⁸) to test the relative contribution of each environmental determinant on tree density. Following the theoretical framework defined in Fig. 3, latitude was included as an exogenous variable controlling variability in NPP and CEC, which in turn are considered to affect tree density at such a broad scale. Mean tree size and plot size were included in the SEM as exogenous local variables that could affect tree density. We applied mixed-effect models with random intercepts following the previous criterion (i.e., tracks within regions) for each causal relationship inside our SEM. For the final model selection, we used Fisher's C information criterion.

To evaluate how local mean tree size modulates the global productivity effect on tree density, we used linear mixed-effects models with the following structure:

$$Y = Xa + Zb + \varepsilon \quad (1)$$

where Y represents tree density across forest regions, a is the vector of parameters in the fixed effects term, that includes plot size and the interaction between the all-inclusive productivity and mean tree size, b is the vector of parameters of the random effects (i.e., nested structure of tracks within a forest region), X and Z are regression matrices of fixed and random effects, respectively, and ε is the within-group error component. For all the above-mentioned analyses, tree density was log-transformed to meet the assumption of linear models of normally distributed errors and homoscedasticity. All continuous predictor variables were standardized to

improve the interpretability of effect sizes and interactions. We built contrasting models differing in their fixed effects using a backward selection procedure starting with a full model that includes all hypothesized variables and the interaction (productivity × mean_DBH). In all cases, model comparison was conducted with the BIC. We used BIC as it performs better than the Akaike Information Criterion when dealing with large datasets²⁹. First, we eliminated the interaction between productivity and mean tree size and compare this model with the initial full model. Second, we kept the interaction between productivity and mean tree size and removed plot size. Lastly, we compared this model with a null model. All the analyses were performed with the R statistical software version 4.1.1³⁰. Generalized additive mixed models were computed using the *gamm* function from the *mgcv* R package³¹. SEM model were analyzed using the *psem* function from the *piecewiseSEM* R package³². Linear mixed-effects model were conducted using the *lme* function from the *nlme* R package³³.

**Reporting summary.** Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

## Data availability

The forest data used can be accessed in https://figshare.com/articles/dataset/Forest_csr/13072211. Source data for Fig. 2 can be found in Supplementary Data 1. Source data for Fig. 4 can be found in Supplementary Data 2.

Received: 23 May 2022; Accepted: 4 January 2023;

Published online: 13 January 2023

## References

1. Crowther, T. W. et al. Mapping tree density at a global scale. *Nature* 525, 201–205 (2015).

2. Asner, G. P. et al. A universal airborne LiDAR approach for tropical forest carbon mapping. *Oecologia* 168, 1147–1160 (2012).

3. Walker, A. P. et al. Predicting long-term carbon sequestration in response to CO₂ enrichment: How and why do current ecosystem models differ? *Glob. Biogeochem. Cy.* 29, 476–495 (2015).

4. Madrigal-González, J. et al. Climate reverses directionality in the richness-abundance relationship across the World's main forest biomes. *Nat. Commun.* 11, 1–7 (2020).

5. Stephenson, N. L. Climatic control of vegetation distribution: the role of the water balance. *Am. Nat.* 135, 649–670 (1990).
6. Weiskittel, A. R., Maguire, D. A., Monserud, R. A. Development of a hybrid model for intensively managed Douglas-fir in the Pacific Northwest. *In Forest Growth and Timber Quality*, 49 (USDA, Portland, 2009).

7. Paoli, G. D., Curran, L. M. & Slik, J. W. F. Soil nutrients affect spatial patterns of aboveground biomass and emergent tree density in southwestern Borneo. *Oecologia* **155**, 287–299 (2008).

8. Yoda, K., Kira, T., Ogawa, H. & Hozami, K. Self-thinning in overcrowded pure stands under cultivated and natural conditions. *J. Biol. Osaka City Univ.* **14**, 107–129 (1963).

9. Westoby, M. The self-thinning rule. *Adv. Ecol. Res.* **14**, 167–225 (1984).

10. Weiner, J. & Freckleton, R. P. Constant final yield. *Annu. Rev. Ecol. Evol. Sci.* **41**, 173–192 (2010).

11. Pillet, M. et al. Disentangling competitive vs. climatic drivers of tropical forest mortality. *J. Ecol.* **106**, 1165–1179 (2018).

12. Schluter, D. Experimental evidence that competition promotes divergence in adaptive radiation. *Science* **266**, 798–801 (1994).

13. Pacala, S.W. & Levin, S.A. Biologically generated spatial pattern and the coexistence of competing species. *Spatial Ecology: The Role of Space in Population Dynamics and Interspecific Interactions*. (Princeton University Press, Princeton, NJ, p. 204-232, 1997).

14. Asefa, M., Cao, M., Zhang, G., Ci, X. & Li, J. Yang Environmental filtering structures tree functional traits combination and lineages across space in tropical tree assemblages. *Sci. Rep.* **7**, 1–10 (2017).

15. Pretzsch, H. & Biber, P. Tree species mixing can increase maximum stand density. *Can. J. For. Res.* **46**, 1179–1193 (2016).

16. Escudero, A. et al. Every bit helps: The functional role of individuals in assembling any plant community, from the richest to monospecific ones. *J. Veg. Sci.* **32**, e13059 (2021).

17. Jump, A. S. et al. Structural overshoot of tree growth with climate variability and the global spectrum of drought-induced forest dieback. *Glob. Chang. Biol.* **23**, 3742–3757 (2017).

18. M. Takyu, Y. Kubota, S.I. Aiba, T. Seino, T. Nishimura. Pattern of changes in species diversity, structure and dynamics of forest ecosystems along latitudinal gradients in East Asia. In *Forest Ecosystems and Environments* (Springer, Tokyo, 2005), pp. 49–58.

19. Rivoire, M. & Le, G. A. Moguedec, generalized self-thinning relationship for multi-species and mixed-size forests. *Ann. Sci.* **69**, 207–219 (2012).

20. Salas-Eljatib, C. & Weiskittel, A. R. Evaluation of modelling strategies for assessing self-thinning behaviour and carrying capacity. *Ecol. Evol.* **8**, 10768–10779 (2018).

21. Schietti, J. et al. Forest structure along a 600 km transect of natural disturbances and seasonality gradients in central-southern Amazonia. *J. Ecol.* **104**, 1335–1346 (2016).

22. Vanclay, J. K. & Sands, P. J. Calibrating the self-thinning frontier. *For. Ecol. Manag.* **259**, 81–85 (2009).

23. Sapijanskas, J., Paquette, A., Potvin, C., Kunert, N. & Loreau, M. Tropical tree diversity enhances light capture through crown plasticity and spatial and temporal niche differences. *Ecology* **95**, 2479–2492 (2014).

24. Lieth, H. Modeling the primary productivity of the world. In H. Lieth & R. H. Whittaker, eds. Primary Productivity of the Biosphere (Springer-Verlag, New York, New York, USA, 1975), pp. 237–264.

25. Grieser, J., Gommes, R., Cofield, S., Bernardi, M. World Maps of Climatological net Primary Production of Biomass, NPP. Food and Agriculture Organization of the United Nations. (GEONETWORK. FAO, Rome, Italy, 2006).

26. Hengl, T. et al. SoilGrids250m: global gridded soil information based on machine learning. *PLoS One* **12**, e0169748 (2017).

27. Wood, S. N. Stable and efficient multiple smoothing parameter estimation for generalized additive models. *J. Am. Stat. Assoc.* **99**, 673–686 (2004).

28. J.B. Grace. Structural Equation Modeling and Natural Systems. (Cambridge University Press, Cambridge, 2006).

29. Aho, K., Derryberry, D. & Peterson, T. Model selection for ecologists: the worldviews of AIC and BIC. *Ecology* **95**, 631–636 (2014).

30. R Core Team (2021). R: A language and environment for statistical computing. (R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/).

31. Wood, S., Scheipl, F. & Wood, M. S. Package 'gamm4'. *Am. Stat.* **45**, 339 (2017).

32. Lefcheck, J. S. piecewiseSEM: Piecewise structural equation modelling in r for ecology, evolution, and systematics. *Methods Ecol. Evol.* **7**, 573–579 (2016).

33. Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D. R. Core Team nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-148 (2020).

34. Latham J, Cumani R, Rosati I, Bloise M. FAO Global Land Cover (GLC-SHARE) Database Beta-Release 1.0, Land and Water Division. 2014. http://www.fao.org/uploads/media/glc-share-doc.pdf.

## Acknowledgements

JMG was funded by the Research Commission of the Swiss National Park, Swiss Academy of Sciences SCNAT (project FLUCTREE- 02/2020). J.C. thanks to Spanish Ministerio de Ciencia e Innovación (Project UNIPER, PID2020-114851GA-100). J.B.C. was funded by the Talento fellowships, Comunidad de Madrid (project INOVA-RISK, 2020-T1/AMB-19913). AE thanks the European Social Fund managed by the Regional Government of Madrid (Remedial TE-CM: S2018/EMT-4338). L.C. and M.J.M. were partly funded through projects CGL2013-45634-P, CGL2016-75414-P, and PID2019-105064GB-100. L.M. was funded by the Swiss National Science Foundation (PCEFP2_181115) and the Margarita Salas Postdoctoral Fellowship from Universidad de Alcalá (MIU, Spain). A.H. was supported by the Basque Country Government funding support to FisioClimaCO₂ (IT1682-22) research group. GBD was funded by the Spanish Ministry of Education (MINEDU; FPU14/05303) and the Regional Government of Madrid (Remedial TE-CM: S2018/EMT-4338).

## Author contributions

Conceptualization: J.M.G., J.C., M.S., A.E. Data acquisition: J.M.G., J.C., J.B.C., A.E., L.C., L.M., M.R., P.R.B., A.H., C.A., R.S., A.P., S.D., C.E., O.T., M.M., L.P., J.L.S., M.J.M., M.A., M.A.Z., A.Q.R., M.V.A., E.G., Y.T., G.B.D., I.G.C., M.S. Methodology: J.M.G., J.C., A.E., L.M. Investigation: J.M.G., J.C., M.S. Visualization: J.M.G., M.S. Funding acquisition: J.M.G., J.C., M.S., A.E. Project administration: J.M.G., J.C., M.S., A.E. Supervision: J.M.G., M.S. Writing – original draft: J.M.G., J.C., J.B.C., A.E., L.M. Writing – review & editing: J.M.G., J.C., J.B.C., A.E., L.C., L.M., M.R., P.R.B., A.H., C.A., R.S., A.P., S.D., C.E., O.T., M.M., L.P., J.L.S., M.J.M., M.A., M.A.Z., A.Q.R., M.V.A., E.G., Y.T., G.B.D., I.G.C., MS.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s42003-023-04419-8.

**Correspondence** and requests for materials should be addressed to Jaime Madrigal-González.

**Peer review information** *Communications Biology* thanks Eoin O’Gorman and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Luke R. Grinham. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023