Evidence of co-limitation in global forest diversity gradients


See “Docs” for the author list and Supplemental Materials



The latitudinal diversity gradient (LDG) is one of the most recognized global patterns of species richness exhibited across a wide range of taxa1,2. Numerous hypotheses3-5 have been proposed in the last two centuries to explain LDG, but none have been sufficiently supported. This shortcoming is in part due to a lack of high-quality global species richness data that can be used to test existing hypotheses. Here, we developed a multi-variable, multi-scale analytical framework to organize a global forest inventory database with individual tree information and local biophysical characteristics from ~1.3 million sample plots. Using this framework, we quantified drivers of local tree species richness patterns across latitudes. Generally, annual mean temperature was a dominant predictor of tree species richness, which is most consistent with the metabolic theory of biodiversity (MTB)6,7. However, MTB underestimates LDG in the tropics where high species richness was moderated also by topographic, soil, and anthropogenic factors operating at local scales. We suggest a unique passage to MTB – the clause of co-limitation – which specifies the lack of a definite dominant factor at low latitudes such that local landscape variables operate synergistically with bioclimatic factors in shaping the global LDG pattern. Our analysis also produced a high-resolution (0.025°×0.025°) map of local tree species richness which is a significant step forward in mapping Earth’s biodiversity and underlying factors at scales that matter for conservation.

Identifying which mechanisms moderate global biodiversity patterns8,9 has perplexed the scientific community for more than two centuries3,4. The most noticeable pattern, LDG, is a trend of declining local species richness (alpha diversity) from low to high latitudes. This trend has been observed for many taxonomic groups and across land, freshwater, and marine environments1,2. More than 30 hypotheses have been proposed10 to explain LDG11, but few can be reconciled with existing observational data for predicting biodiversity decline towards the poles. To test these varied hypotheses, biodiversity data must be assembled that are global in scope with sufficient sample coverage across all ecoregions and biomes.

In addition to biodiversity data, testing these varied hypotheses also requires data on a wide spectrum of potential drivers that may moderate biodiversity at local scales11,12, such as climate, soil and land features, as well as anthropogenic factors. For instance, environment temperature (i.e., ambient temperature of the air, represented by annual mean temperature) is largely responsible for the generation and maintenance of biodiversity, through the effects of solar radiation on demographic rates (e.g., growth and mortality), ecological interactions (e.g., predation and competition) and evolutionary rates of change (e.g., speciation and extinction)13,14. Soil and topographic heterogeneity facilitate niche partitioning via inducing microclimatic variation, contributing to compositional variation15 and biodiversity maintenance16,17. Furthermore, humans have a long history of reshaping biodiversity through the selective use of natural resources and the modification of native species composition18. In addition, multiple subordinate factors jointly affecting biodiversity could potentially increase the diversity of niche opportunities, thereby resulting in species-rich assemblages.

Here, we quantified the relative contribution of a wide range of environmental factors across space on local tree species richness in forested areas around the world. To accomplish this, we standardized a global tree species richness (i.e., as alpha diversity) database (Figs. 1&2) and quantified the relative contribution of 47 explanatory variables including bioclimatic conditions (e.g., annual mean temperature), vegetation and survey attributes (e.g., sample plot size), topographic covariates (e.g., terrain roughness), soil covariates (e.g., bulk density), and anthropogenic spatial features (e.g., size of roadless areas) in an attempt to test whether local co-limitation exists when multiple subordinate drivers co-dominate (Fig. 3). We conducted a three-stage analysis (Fig. 1, see Methods in Supplementary Information for details) based on two independent ground-sourced forest inventory datasets (Phase-I and Phase-II, Fig. 2). The main dataset (Phase-I) consisted of 1,255,444 sample plots, while the validating dataset (Phase-II) consisted of 22,131 sample plots, most of which are located in unsampled and under-sampled regions of Phase-I dataset. Together, our sample data covered 424 of the 435 (97%) forested ecoregions worldwide (Fig. S1), with a total of ~55 million sample trees representing more than 32,000 species.

Our analyses confirmed, with a high level of accuracy, one general spatial trend in local tree species richness worldwide that has led us to three conclusions regarding the underlying mechanisms to tree species richness patterns. Firstly, we found that LDG for tree species richness was consistent with that of most other groups of organisms, with a decline from the tropics to the poles (Fig. 4). In the Northern Hemisphere, for instance, tree species richness dropped sharply from the equator (98 species·ha-1) to 10°N with an average rate of decline of 6 species·ha-1 per 1° increase in latitude, after which the decline diminished and stabilized at 4 species·ha-1 at 50°N. In the Southern Hemisphere, tree species richness declined from the equator to 25°S on average by 3 species·ha-1 per 1° increase in latitude, after which tree species richness fluctuated before another steep drop from 25 species·ha-1 (43°S) to 4 species·ha-1 (50°S). We were able to detect and map regional patterns and global peaks of tree species diversity, with a high spatial resolution (0.025°×0.025°). The Amazonian, Southeast Asian, and Melanesian rainforests are clearly the regions with the greatest local tree species richness worldwide, containing >200 tree species·ha-1 above the 5 cm diameter-at-breast-height (DBH) threshold, confirming previous findings19,20. Tropical African rainforests generally contain 50% fewer tree species per hectare than Amazonian rainforests. In the temperate forests of the Northern Hemisphere, the Changbai Mountains in Northeast Asia (up to ~28 species·ha-1) and the Central Appalachian forests in the Eastern United States (up to ~20 species·ha-1) display high local species richness. In the Southern Hemisphere, the sclerophyllous and Nothofagus-dominated forests in south-central Chile are among the most species-rich temperate communities (up to 50 species·ha-1). Boreal forest communities are consistently low in local tree species richness, with typically five or fewer tree species per hectare.  

Secondly, the LDG pattern of tree species richness was mostly attributable to bioclimatic conditions, but environmental temperature alone was not a good predictor of local tree species richness at low latitudes (Fig. 5). At higher latitudes, bioclimate is a dominant determinant of tree species richness, except in areas of Central Europe and East Asia that are also subject to strong topographic influences. Based on the partitioning of spatial variance (see Fig. 6 and §Variance Partitioning in Supplementary Methods) across all the forested continents, the observed LDG pattern was almost entirely attributable to the effect of environmental conditions on tree species richness. Among all the environmental attributes, annual mean temperature was the key predictor of tree species richness. This finding is consistent with the metabolic theory of biodiversity (MTB)6,7, which specifically postulates that environment temperature (represented by annual mean temperature) is largely responsible for the generation and maintenance of biodiversity14,21,22. According to MTB, the natural logarithm of species richness is linearly associated with 1000/T, where T is the environmental temperature in Kelvin (mean annual temperature + 273.15K), with a slope ranging from -7.5 to -9.0 K. Our global tree species richness gradient was consistent with MTB, with a slope of -8.0 K (p<0.001). However, despite a coefficient of determination of 0.82 (see §Metabolic Theory of Biodiversity in Supplementary Methods), MTB substantially underestimated LDG at low latitudes. In fact, near the equator where the actual LDG peaked (98 species·ha-1), tree species richness was almost twice as high as predicted by MTB (56 species·ha-1) (Fig. 4A). Our results suggest that within this low latitudinal range other factors are also important to the maintenance of biodiversity.

Our third conclusion, which marks a unique contribution to the understanding of LDG patterns, is that the under-estimation of local tree species richness by MTB at low latitudes is attributable, in part, to the lack of a definite dominant environmental factor, suggesting a co-limitation of multiple drivers at low latitudes (Fig. 5). In general, bioclimatic factors dominantly determined species richness in 82.6% of the forested areas, while co-limitation (i.e., absence of any dominating factor) occurred in 11.7% of forested areas globally. However, in the low-latitude range between 5°N and -15°S, the percentage area of co-limitation increased to 37.1%, more than three times the global average. Furthermore, forested areas under co-limitation contained on average 81.1 ± 0.1 species per hectare, much higher than the average local tree species richness of forested areas dominantly determined by topographic (43.9 ± 0.1), anthropogenic (35.6 ± 0.2), soil (33.9 ± 0.2), and bioclimatic (19.4 ± 0.02) factors (Fig. 5B). This suggests that the pattern of co-limitation is pervasive in species-rich tropical forests. In South America, for instance, transitional areas between Amazonia and savanna formations nearby are subject to co-limitation that is partly attributable to a dynamic equilibrium between closed forest and savanna23, edaphic conditions, and natural fire regimes24. In Africa, anthropogenic influences such as selective timber extraction and fuelwood collection, together with large-scale degradation25 affect local tree species richness (Figs. 5 & S7). In Central Africa, the evolution of anthropogenic influences from prehistoric to present times has imposed a substantial effect on species diversity26 and resulted in the development of a complex system of mixes with light-demanding and old-growth tree species.

In addition to an overall positive response of local tree species richness to the rise of annual mean temperature (see the partial dependence plot [PDP] of C1 in Fig.3 and Fig.S4), the importance of environmental temperature (2.7%) was topped by the total annual precipitation (C12, 7.6%) (Fig. 3). Our findings are consistent with previous discoveries of a joint role of water and temperature/energy – as a proxy for net primary productivity27 – on plant species richness, with water dominating particularly at warmer, lower latitudes22,28. Predicted tree species richness accelerated exponentially with temperature and rainfall, although independently, as shown in the cold-dry quadrant and the convex contours of the 2D PDP (Fig. S4), until each has reached its respective threshold (1500mm for total annual precipitation and 10°C for annual mean temperature). Beyond one of these thresholds, species richness is only limited by the predictor below its threshold (i.e., by annual mean temperature in the cold-wet quadrant, or by annual precipitation in the hot-dry quadrant). When both predictors have reached their thresholds, i.e., in the hot-wet quadrant, co-limitation predominates in most tropical forests. Net primary productivity in the tropics, thus, requires co-limitation of other factors besides only temperature and rainfall29. As carbon flux responses mirror the co-limitation clause for tree species diversity, the matching determinants for both diversity and productivity may explain both the similar latitudinal gradient in productivity and the positive diversity-productivity relationship30,31. Our findings also indicate that under climate change, intensified droughts coupled with increased annual mean temperature32 can potentially trigger declines of tree species richness, although possible increases in water-use efficiency from elevated CO2 and the dominance of highly contingent co-limiting factors may partially buffer this effect in the tropics33.

Here, we articulate for the first time, to our knowledge, evidence for a co-limitation clause to the MTB hypothesis. Resource co-limitation is a common concept in ecology (e.g., 34,35), often used to describe how the synergistic interactions of two or more factors limit ecological productivity36. Our use of the term co-limitation emphasizes the reduced significance of a globally predominant driver of species richness at low latitudes, recognizing that several local subordinate factors synergistically contribute to increased tree species richness in this latitudinal range. We thus argue that the explanatory power of MTB could be improved by considering multiple subordinate factors where single-factor dominance is lacking, especially in the tropics. We posit that the increasing intricacy of biodiversity maintenance mechanisms towards low latitudes constitutes the co-limitation clause to the MTB hypothesis. At high latitudes, bioclimatic conditions, particularly environmental temperature, are the major limiting factors, and thus the dominant drivers of tree species diversity. As the latitude declines especially between 5°N and 15°S, the influence of bioclimatic conditions dwindles, and the maintenance of tree species richness is moderated by many interacting drivers without a clear dominance (Fig. 5). This prevalence of co-limiting factors is thus not a mere coincidence as to why the observed LDG at low latitudes is almost double that predicted by MTB (Fig. S3). While each of the existing hypotheses underpinning LDG addresses a certain process10,12 (e.g., selection, drift, dispersal, or speciation), the co-limitation clause highlights synergistic interactions of local processes across the latitudinal gradient.

More research is needed to fully elucidate patterns of LDG driven by climatic and other influences, especially those outlined in competing hypotheses. First, our analyses lack explicit consideration of some evolutionary, ecological and historical factors. These include mid-domain stochastic effects37, the legacies of the poleward expansion of tree species after the Last Glacial Maximum38,39, and recent human colonization. Alternative hypotheses, such as niche conservatism or climatic history, are more difficult to test due to data limitations. In addition, long-term effects at geological and millennial time scales also play a role, but it is difficult to disentangle these effects due to collinearity40. A major source of uncertainty in our results (Fig. 4B) came from an uneven sample coverage between developed and developing countries (Fig. 2). To address this gap, we argue that there needs to be a shared responsibility among forestry agencies at various levels of government, scientists, indigenous communities, and other biodiversity monitoring groups to improve sample coverage of forest inventories in developing countries. Innovative biodiversity funding mechanisms, e.g., forest inventories funded by carbon initiatives such as REDD+, should be incorporated into a comprehensive global forest biodiversity database. Meanwhile, the severe shortage of experts and database management infrastructures, especially in developing countries, poses another major challenge to address this gap41. The education and training of new generations of forest scientists, taxonomists, and foresters can bring tangible benefits to biodiversity monitoring while improving local economies as well.

The co-limitation clause to MTB enables a refined description of the biogeographic distribution of biodiversity and mechanisms underlying LDG. Our analysis has resulted in the production of a high-resolution map of tree species richness across the global forest range, along with visuals of those factors responsible for the moderation of local tree species richness. Such tools are necessary for conservation management which requires assessments of factors responsible for biodiversity patterns at multiple scales that matter – from local, regional to global scales. Patterns of local tree species richness and associated drivers may provide insights into how and why the diversity of other forest flora, fauna, and microbes42,43 vary across space and time. Furthermore, the high-resolution map of local tree species richness presented here provides a benchmark for evaluating the impact of biodiversity loss on the productivity and functioning of forest ecosystems31,44. Finally, aligned with current international calls for spatially explicit monitoring of ecosystem attributes45, this study delivers detailed biogeographic information to support international endeavors46 focused on valuing natural capital and advancing global conservation.


We thank the following agencies, initiatives, teams, and individuals for data collection and other technical support: the Global Forest Biodiversity Initiative (GFBI) for establishing the data standards and collaborative framework; United States Department of Agriculture, Forest Service, Forest Inventory and Analysis (FIA) Program; University of Alaska Fairbanks; The SODEFOR, Ivory Coast; University Félix Houphouët-Boigny (UFHB, Ivory Coast); the Queensland Herbarium and past Queensland Government Forestry and Natural Resource Management Departments and staff for data collection for over seven decades; the National Forestry Commission of Mexico (CONAFOR). We thank Marc Baker (Carbon Tanzania), together with a team of field assistants (Valentine and Lawrence); all persons who made the Third Spanish Forest Inventory possible, especially the main coordinator, J. A. Villanueva (IFN3); The French National Forest Inventory (NFI campaigns [raw data 2005 and following annual surveys, were downloaded by GFBI at, site accessed on January 1st, 2015); the Italian Forest Inventory (NFI campaigns Raw data 2005 and following surveys were downloaded by GFBI at, site accessed on April 27th, 2019); Swiss National Forest Inventory, Swiss Federal Institute for Forest, Snow and Landscape Research WSL and Federal Office for the Environment FOEN, Switzerland; The Danish National Forestry, Department of Geosciences and Natural Resource Management, UCPH; Coordination for the Improvement of Higher Education Personnel of Brazil (CAPES, grant no. 88881.064976/2014-01); Rafael Ávila and Sharon van Tuylen, Instituto Nacional de Bosques (INAB), Guatemala for facilitating Guatemalan data; The National Focal Center for Forest condition monitoring of Serbia (NFC), Institute of Forestry, Belgrade, Serbia; The Thünen Institute of Forest Ecosystems (Germany) for providing National Forest Inventory data; the Food and Agriculture Organization of the United Nations (FAO) and the United Nations High Commissioner for Refugees (UNHCR) undertaking the SAFE (Safe Access to Fuel and Energy) and CBIT-Forest projects; the Amazon Forest Inventory Network (RAINFOR), the African Tropical Rainforest Observation Network (AfriTRON), and the initiative for their contributions from Amazonian and African forests. The Natural Forest plot data collected between January 2009 and March 2014 by the LUCAS programme for the New Zealand Ministry for the Environment, are provided by the New Zealand National Vegetation Survey Databank The International Boreal Forest Research Association (IBFRA). The Forestry Corporation of New South Wales, Australia. The National Forest Directory of the Ministry of Environment and Sustainable Development of the Argentine Republic (MAyDS) for the plot data of the Second National Forest Inventory (INBN2). The Sabah Biodiversity Council and the staff from Sabah Forest Research Centre. All TEAM data are provided by the Tropical Ecology Assessment and Monitoring (TEAM) Network, a collaboration between Conservation International, the Missouri Botanical Garden, the Smithsonian Institution, and the Wildlife Conservation Society, and partially funded by these institutions, the Gordon and Betty Moore Foundation, and other donors, with thanks to all current and previous TEAM site manager and other collaborators that helped collecting data; The people of the Redidoti, Pierrekondre and Cassipora village who were instrumental in assisting with the collection of data and sharing local knowledge of their forest; and the dedicated members of the field crew of Kabo 2012 census. This research was supported in part through computational resources provided by Information Technology at Purdue, West Lafayette, Indiana.



This work is supported in part by the NASA Grant #12000401 “Multi-sensor biodiversity framework developed from bioacoustic and space based sensor platforms;” the USDA National Institute of Food and Agriculture McIntire Stennis projects 1017711 and 1016676; São Paulo Research Foundation, #2014/14503-7; CBIT-Forest project (GCP/GLO/882/CBT); Proyecto FONACIT No. 1998003436 and UNELLEZ No. 23198105; EU, Sumforest – REFORM, Risk Resilient Forest Management, FKZ: 2816ERA02S; German Research Foundation (DRG), KROOF Tree and stand-level growth reactions on drought in mixed versus pure forests of Norway spruce and European beech, PR 292/12-1; Bavarian State Ministry for Food, Agriculture and Forestry, W07 long term yield experiments, 7831-26625-2017 and Project No E33; The Deutsche Forschungsgemeinschaft (DFG) Priority Program 1374 Biodiversity Exploratories; The International Tropical Timber Organization, ITTO-Project PD 53/00 Rev.3 (F); The State Forest Management Centre, Estonia, and the Environmental Investment Centre, Estonia; Natural Sciences and Engineering Research Council of Canada Discover Grant Project (RGPIN-2014-04181 and STPGP428641); European Structural Funds by FEDER 2014-2020 GY0006894; ICNF-Instituto de Conservação da Natureza e Florestas. Inventário Florestal Nacional (IFN); FCT – Portuguese Foundation for Science and Technology, under the project UID/AGR/04033/2020; Vietnam National Foundation for Science and Technology Development (NAFOSTED-106-NN.06-2016.10); Australian Centre for International Agricultural Research (Grant no: ASEM/2010/50); German Research Foundation (DFG, FOR 1246); The project LIFE+ ForBioSensing PL “Comprehensive monitoring of stand dynamics in Białowieża Forest supported with remote sensing techniques” co-funded by Life+ (contract number LIFE13 ENV/PL/000048) and the National Fund for Environmental Protection and Water Management in Poland (contract number 485/2014/WN10/OP-NM-LF/D); National Natural Scientific Foundation of China (31660055 and 31660074); The Polish State Forests National Forest Holding (2016); the National Science Center (Poland) grant 2011/02/A/NZ9/00108; The Dutch Ministry of Agriculture, Nature Management and Food Quality for funding the Dutch National Forest Inventory; The Danish Ministry of Environment for funding the Danish National Forest Inventory; The Grant 11-TE11-0100 from the U.S. National Space and Aeronautics Administration; the Tropical Ecology, Assessment, and Monitoring (TEAM) / Conservation International project for funding the data collection, and the Instituto Nacional de Pesquisas da Amazônia (INPA); The Ministère des Forêts, de la Faune et des Parcs du Québec (Canada); the European Union H2020 VERIFY project (contract 776810): the H2020 European Union project Resonate (contract  101000574); The Exploratory plots of FunDivEUROPE received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement 265171; DBT, Govt. of India through the project ‘Mapping and quantitative assessment of geographic distribution and population status of plant resources of Eastern Himalayan region’ (sanction order No. BT/PR7928/NDB/52/9/2006 dated 29th September 2006); DBT, Govt. of India through the project “Mapping and quantitative assessment of plant resources and its distribution in Madhya Pradesh, Central India” (Sanction Order no.: BT/PR12899/NDB/39/506/2015); The financial support from Natural Sciences and Engineering Research Council of Canada to S. Dayanandan; Czech Science Foundation Standard Grant (19-14620S) and European Research Council advanced grant (669609); RFBR #16-05-00496; The project implementation Demonstration object on the transformation of declining spruce forests into ecologically more stable multifunctional ecosystems, ITMS 26220220026, supported by the Research & Development Operational Program funded by the ERDF; The Swedish NFI, Department of Forest Resource Management, Swedish University of Agricultural Sciences SLU; The National Research Foundation (NRF) of South Africa (89967 and 109244) and the South African Research Chair Initiative; University Research Committee of the University of the South Pacific, and New Colombo Plan Funding through the Department of Foreign Affairs and Trade of the Australian government; The TEAM project in Uganda supported by the Moore foundation and Buffett Foundation through Conservation International (CI) and Wildlife Conservation Society (WCS); COBIMFO project funded by the Belgian Science Policy Office (Belspo), contract no. SD/AR/01A; The German Federal Ministry of Education and Research (BMBF) under Grant FKZ 01LL0908AD for the project “Land Use and Climate Change Interactions in the Vu Gia Thu Bon River Basin, Central Vietnam” (LUCCI); Programme Tropenbos Côte d’Ivoire : projet 04/97-1111a du “Complément d’Inventaire de la Flore dans le Parc National de Taï”; The Danish Council for Independent Research | Natural Sciences (TREECHANGE, grant 6108-00078B) and VILLUM FONDEN (grant 16549) to JCS; The Natural Environment Research Council of the UK (NERC) project NE/T011084/1 awarded to JAG and NE/ S011811/1; ERC Advanced Grant 291585 (“T-FORCES”) and a Royal Society-Wolfson Research Merit Award to OLP; RAINFOR plots supported by the Gordon and Betty Moore Foundation and the U.K. Natural Environment Research Council (NERC), notably NERC Consortium Grants ‘AMAZONICA’ (NE/F005806/1), ‘TROBIT’ (NE/D005590/1), and ‘BIO-RED’ (NE/N012542/1); Fundação de Amparo à Pesquisa e Inovação de Santa Catarina, FAPESC (2016TR2524), Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq [312075/2013-8); “Investissement d’Avenir” grant managed by Agence Nationale de la Recherche (CEBA, ref. ANR- 10-LABX-25-01); CIFOR’s Global Comparative Study on REDD+ funded by the Norwegian Agency for Development Cooperation (Norad), the Australian Department of Foreign Affairs and Trade (DFAT), the European Union (EU), the International Climate Initiative (IKI) of the German Federal Ministry for the Environment, Nature Conservation, Building and Nuclear Safety (BMUB), and the CGIAR Research Program on Forests, Trees and Agroforestry (CRP-FTA), and donors to the CGIAR Fund; The Nature and Biodiversity Conservation Union (NABU) under the project entitled “Biodiversity under Climate Change: Community Based Conservation, Management and Development Concepts for the Wild Coffee Forests”, funded by the German Federal Ministry for the Environment, Nature Conservation and Nuclear Safety (BMU) through the International Climate Initiative (IKI); The Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq); The institutional project “EXTEMIT – K”, no. CZ.02.1.01/0.0/0.0/15_003/0000433 financed by OP RDE; EC DG VIII grant BZ-5041 (ECOSYN), NWO-WOTRO (W84-204), and GTZ; AfriTRON network plots funded by the local communities and NERC, ERC, European Union, Royal Society and Leverhume Trust; BOLFOR (Proyecto de Manejo Forestal Sostenible- Bolivia); The Global Environment Research Fund F-071 and D-1006, and JSPS KAKENHI Grant Numbers JP17K15289; The National Institute of Biology (now Research Center for Biology), LIPI (Indonesian Institute of Sciences), Indonesia IFBN project (contract 4000114425/15/NL/FF/gp) funded by ESA; NSF grant DBI-1565046; NSF DEB-0424767, NSF DEB-0639393, NSF DEB-1147429, NASA Terrestrial Ecology Program,Swiss National Science Foundation (SNSF No. 130720, 147092); Projects D/9170/07, D/018222/08, D/023225/09 and D/032548/10 funded by the Spanish Agency for International Development Cooperation [Agencia Española de Cooperación Internacional para el Desarrollo (AECID)] and Fundación Biodiversidad, in cooperation with the Universidad Mayor de San Simón (UMSS), the FOMABO (Manejo Forestal en las Tierras Tropicales de Bolivia) project and CIMAL (Compañía Industrial Maderera Ltda.); The Agency for Economic and Environmental Development (DDEE) of the north province of New Caledonia (the projects Ecofor & Cogefor, 2011-2016); Russian Science Foundation (16-17-10284 “The accumulation of carbon in forest soils and forest succession status”); Norwegian Ministry of Food and Agriculture; A grant from the Royal Society and the Natural Environment Research Council (UK) to S.L.L.; The Spanish Agency for International Development Cooperation [Agencia Española de Cooperación Internacional para el Desarrollo (AECID)] and Fundación Biodiversidad, in cooperation with the governments of Syria and Lebanon; COBIMFO Project, Federal Science Policy, Belgium; Consejo Nacional de Ciencia y Tecnología, Mexico; Comisión Nacional Forestal, Mexico; BEF-China project (FOR 891) funded by the German Research Foundation (DFG); WWF Russell Train Fellowship to P.M.U. (Grant ST54); Wildlife Conservation Society DRC Program under CARPE Funding; Seoul National University Big Data Institute through the Data Science Research Project 2016, R&D Program for Forest Science Technology (Project No. 2013069C10-1719-AA03 & S111215L020110) funded by Korea Forest Service (Korea Forestry Promotion Institute); The European Union’s Horizon 2020 research and innovation program within the framework of the Multifunctionality Marie Skłodowska-Curie Individual Fellowship (IF-EF) under grant agreement 655815;  Tropenbos International-Suriname; The Institute for World Forestry, University of Hamburg; REMBIOFOR Project “Remote sensing based assessment of woody biomass and carbon storage in forests” funded by The National Centre for Research and Development, Warsaw, Poland, under the BIOSTRATEG program (agreement no. BIOSTRATEG1/267755/4/NCBR/2015); Project “Environmental and genetic factors affecting productivity of forest ecosystems on forest and post-industrial habitats” (2011-2015; no. OR/2717/3/11 ) funded by the General Directorate of State Forests, Warsaw, Poland; Project “Carbon balance of the major forest-forming tree species in Poland” (2007-2011; no. 1/07) funded by the General Directorate of State Forests, Warsaw, Poland; and the research professorship for “Ecosystem-based sustainable development” funded by Eberswalde University for Sustainable Development. Pontificia Universidad Católica del Ecuador supported fieldwork census in Yasuni National Park; National Forest Programme of the National Institute of Agricultural Research (INTA); the São Paulo Research Foundation (FAPESP; grants #2014/14503-7, 2017/05662-2, and 03/12595-7); CNPq/FAPESP-PELD (grants No. 403710/2012-0 and 2012/51509-8) MAUA group supported by FAPEAM-PRONEX (grant No. 1600/2006); CNPq/FAPEAM-PELD (grant No. 403792/2012-6); ATTO project (grant No. MCTI-FINEP 1759/10 and BMBF 01LB1001A); Czech Science Foundation Standard Grant (17-07378S and 17-19376S); the long-term research development project RVO 67985939 of Institute of Botany of the Czech Academy of Sciences; Slovak Research and Development Agency under the project APVV-20-0168; The Strategic Science Investment Fund of the New Zealand Ministry for Business, Innovation and Employment; FAPESP grants 2014/14503-7 and 2017/05662-2; and CNPq Universal (grant No. 479599/2008-4); research grants LTAUSA19137 (program INTER-EXCELLENCE, subprogram INTER-ACTION) provided by Czech Ministry of Education, Youth and Sports, 20-05840Y of the Czech Science Foundation and the Czech Academy of Sciences (grant nr. RVO 67985939). The Digital Environment for Enabling Data-driven Science (DEEDS) Project is funded by National Science Foundation CIF21 DIBBs: EI: #1724728; The Manchester Metropolitan University; Russian Science Foundation supported forest measurements in European Russia (project #19-77-30015) and forest plots established in Central Siberia (project #21-46-07002); the São Paulo Research Foundation (FAPESP grants 2003/12595-7 and 2012/51872-5); H.L. was supported by National Natural Science Foundation of China (31800374) and Shandong Provincial Natural Science Foundation (ZR2019BC083). Czech Science Foundation (21-26883S) and MSMT INTER-EXCELLENCE project (LTAUSA18007); Science and Engineering Research Board (SERB), Department of Science and Technology, New Delhi under National Post-Doctoral Fellowship Scheme (Ref. No.: PDF/2015/000447), project: Assessing the carbon sequestration potential of different forest types in Central India in response to climate change at Dr. Harisingh Gour Vishwavidyalaya (A Central University), Sagar, Madhya Pradesh, India; the São Paulo Research Foundation (FAPESP grants 2003/12595-7 and 2012/51872-5). T.J. was supported by a UK NERC Independent Research Fellowship (grant code: NE/S01537X/1). S. d-M benefited from a Serra-Húnter Fellowship provided by the Government of Catalonia (Spain). R.V. was supported by Russian Science Foundation project Project # 19-77-300-12. CNPq, grant 442640/2018-8, CNPq/Prevfogo-Ibama Nº 33/2018; Czech Science Foundation (20-16499S, 18-10781S, and 21-24186M) and Charles University (UNCE204069); The Tropical Plant Exploration Group 70 1-ha plots in Continental Cameroon Mountains are supported by Rufford Small Grant Foundation, UK and 4-ha in Sierra Leone are supported by the Global Challenge Research Fund through Manchester Metropolitan University, UK. The PilotMAB project (Frame Agreement Belgian Development Cooperation and Royal Museum for Central Africa); HERBAXYLAREDD project (Belspo – Belgian Science Policy). Creation of a Forest Information System Covering the Areas of the Sudety and the Western Beskidy Mountains within the Scope of Forest Condition Monitoring and Assessment” financed by the Polish State Forests-National Forest Holding. Plot sampling in northern Kenya by ACS was supported by Marie Curie Actions Intra-European Fellowships (IEF), number 328075, and by the Percy Sladen Fund (the Linnaean Society of London). Plot sampling in Cameroon by ACS was supported by Marie Skłodowska-Curie Actions Global Fellowships, number 74356. The Key Project of National Key Research and Development Plan, China (2017YFC0504005). Plot sampling in the DRC by ACS was supported by the National Geographic Explorer Grant, NGS-53344R-18.

Author contributions:

Conceptualization: JL, CH

Methodology: JL, CH, JGPG, NP

Data coordination: JL, MZ, SdM, TC, GJN, PBR, FS, KvG, JGPG, NP

Data collection: All

Writing, revision, & editing: All

Competing interests

The authors declare no competing interests.

Supplementary Information


Figs. S1 to S7

Tables S1 to S2

References (6195)


1          Hillebrand, H. On the generality of the latitudinal diversity gradient. The American Naturalist 163, 192-211 (2004).

2          Crame, J. A. Taxonomic diversity gradients through geological time. Diversity and Distributions 7, 175-189 (2001).

3          Pickering, C. The Geographical Distribution of Animals and Plants.  (Little, Brown & Co., 1854).

4          Humboldt, A. v. & Bonpland, A. Essai sur la géographie des plantes accompagné d’un tableau physique des régions équinoxiales. Paris: Levrault, Schoell & Co (1805).

5          Gough, L. & Field, R. Latitudinal diversity gradients. eLS (2007).

6          Allen, A. P., Brown, J. H. & Gillooly, J. F. Global biodiversity, biochemical kinetics, and the energetic-equivalence rule. Science 297, 1545-1548, doi:10.1126/science.1072380 (2002).

7          Stegen, J. C., Enquist, B. J. & Ferriere, R. Advancing the metabolic theory of biodiversity. Ecology Letters 12, 1001-1015 (2009).

8          Sutherland, W. J. et al. Identification of 100 fundamental ecological questions. Journal of Ecology 101, 58-67 (2013).

9          Gaston, K. J. Global patterns in biodiversity. Nature 405, 220-227 (2000).

10        Willig, M. R., Kaufman, D. M. & Stevens, R. D. Latitudinal gradients of biodiversity: pattern, process, scale, and synthesis. Annual Review of Ecology, Evolution, and Systematics 34, 273-309 (2003).

11        Pontarp, M. et al. The latitudinal diversity gradient: novel understanding through mechanistic eco-evolutionary models. Trends in Ecology & Evolution 34, 211-223 (2019).

12        Vellend, M. The theory of ecological communities (MPB-57). Vol. 75 (Princeton University Press, 2016).

13        Brown, J. H. Why are there so many species in the tropics? Journal of Biogeography 41, 8-22 (2014).

14        Currie, D. J. et al. Predictions and tests of climate-based hypotheses of broad-scale variation in taxonomic richness. Ecology Letters 7, 1121-1134 (2004).

15        Baldeck, C. A. et al. Soil resources and topography shape local tree community structure in tropical forests. Proceedings of the Royal Society B: Biological Sciences 280, 20122532 (2013).

16        Qian, H. & Ricklefs, R. E. Large-scale processes and the Asian bias in species diversity of temperate plants. Nature 407, 180-182 (2000).

17        Stein, A., Gerstner, K. & Kreft, H. Environmental heterogeneity as a universal driver of species richness across taxa, biomes and spatial scales. Ecology Letters 17, 866-880 (2014).

18        Sponsel, L. E. Human impact on biodiversity, overview. Encyclopedia of Biodiversity, 137-152 (2013).

19        Sullivan, M. J. et al. Diversity and carbon storage across the tropical forest biome. Scientific reports 7, 1-12 (2017).

20        Cazzolla Gatti, R. et al. The number of tree species on Earth. Proceedings of the National Academy of Sciences 119, e2115329119, doi:10.1073/pnas.2115329119 (2022).

21        Wright, D. H. Species-energy theory: an extension of species-area theory. Oikos 41, 496-506 (1983).

22        Hawkins, B. A. et al. Energy, water, and broad-scale geographic patterns of species richness. Ecology 84, 3105-3117, doi: (2003).

23        Staal, A., Dekker, S. C., Xu, C. & van Nes, E. H. Bistability, spatial interaction, and the distribution of tropical forests and savannas. Ecosystems 19, 1080-1091 (2016).

24        de L. Dantas, V., Batalha, M. A. & Pausas, J. G. Fire drives functional thresholds on the savanna–forest transition. Ecology 94, 2454-2463 (2013).

25        Bodart, C. et al. Continental estimates of forest cover and forest cover changes in the dry ecosystems of Africa between 1990 and 2000. Journal of Biogeography 40, 1036-1047 (2013).

26        Hubau, W. et al. The persistence of carbon in the African forest understory. Nature Plants 5, 133-140, doi:10.1038/s41477-018-0316-5 (2019).

27        Šímová, I. & Storch, D. The enigma of terrestrial primary productivity: measurements, models, scales and the diversity–productivity relationship. Ecography 40, 239-252 (2017).

28        Kreft, H. & Jetz, W. Global patterns and determinants of vascular plant diversity. Proceedings of the National Academy of Sciences 104, 5925-5930, doi:10.1073/pnas.0608361104 (2007).

29        Clark, D. A. et al. Net primary production in tropical forests: an evaluation and synthesis of existing field data. Ecological applications 11, 371-384 (2001).

30        Pärtel, M., Laanisto, L. & Zobel, M. Contrasting plant productivity–diversity relationships across latitude: the role of evolutionary history. Ecology 88, 1091-1097 (2007).

31        Liang, J. et al. Positive biodiversity-productivity relationship predominant in global forests. Science 354, aaf8957, doi:10.1126/science.aaf8957 (2016).

32        Hammond, W. M. et al. Global field observations of tree die-off reveal hotter-drought fingerprint for Earth’s forests. Nature Communications 13, 1761, doi:10.1038/s41467-022-29289-2 (2022).

33        Martens, C. et al. Large uncertainties in future biome changes in Africa call for flexible climate adaptation strategies. Global Change Biology 27, 340-358 (2021).

34        Chapin, F. S., Bloom, A. J., Field, C. B. & Waring, R. H. Plant Responses to Multiple Environmental Factors. BioScience 37, 49-57, doi:10.2307/1310177 (1987).

35        Tilman, D. Resource Competition and Community Structure.(MPB-17), Volume 17.  (Princeton university press, 1982).

36        Harpole, W. S. et al. Nutrient colimitation of primary producer communities. Ecology letters 14, 852-862 (2011).

37        Colwell, R. K., Rahbek, C. & Gotelli, N. J. The mid-domain effect and species richness patterns: what have we learned so far? The American Naturalist 163, E1-E23 (2004).

38        Feng, G. et al. Species and phylogenetic endemism in angiosperm trees across the Northern Hemisphere are jointly shaped by modern climate and glacial–interglacial climate change. Global Ecology and Biogeography 28, 1393-1402 (2019).

39        Fraser, C. I., Nikula, R., Ruzzante, D. E. & Waters, J. M. Poleward bound: biological impacts of Southern Hemisphere glaciation. Trends in Ecology & Evolution 27, 462-471 (2012).

40        Algar, A. C., Kerr, J. T. & Currie, D. J. A test of metabolic theory as the mechanism underlying broadscale speciesrichness gradients. Global Ecology and Biogeography 16, 170-178 (2007).

41        Wilson, E. O. A Global Biodiversity Map. Science 289, 2279-2279, doi:10.1126/science.289.5488.2279 (2000).

42        Castagneyrol, B. & Jactel, H. Unraveling plant-animal diversity relationships: a meta-regression analysis. Ecology 93, 2115-2124 (2012).

43        Steidinger, B. S. et al. Climatic controls of decomposition drive the global biogeography of forest-tree symbioses. Nature 569, 404-408, doi:10.1038/s41586-019-1128-0 (2019).

44        Bongers, F. J. et al. Functional diversity effects on productivity increase with age in a forest biodiversity experiment. Nature Ecology & Evolution 5, 1594-1603, doi:10.1038/s41559-021-01564-3 (2021).

45        UN Statistical Commission. System of Environmental-Economic Accounting—Ecosystem Accounting: Final Draft. 350 (United Nations, 2021).

46        The World Bank. Green Bond Impact Report. (Capital Markets Department, the World Bank Treasury, Washington, DC, 2019).


Fig. 1 | A conceptual diagram of the three-stage process employed in the study. (Stage 1) Two independent global forest biodiversity individual-based (GFBi) datasets (Phase-I and Phase-II, see Fig. 2 for details) were standardized into a global tree-level dataframe, and aggregated into a global species abundance matrix. Based on plot locations, we merged the abundance matrix with 47 explanatory variables (Fig. 3) into a standardized plot-level dataframe. (Stage 2) We compared three candidate models (RF: random forests, XGB: XGBoost, OLS: ordinary least squares) trained from the Phase-I plot-level dataframe, using random and spatial cross-validation based on Phase-I data, and post-sample validation based on Phase-II data (Fig. S2). The final model was then selected and re-calibrated with both Phase-I and Phase-II data. (Stage 3) Using the final model, we standardized and mapped local tree species richness per hectare across the global forest range. Based on this globally continuous map, we quantified the associated latitudinal diversity gradient (LDG, Fig. 4A), and tested for the metabolic theory of biodiversity (MTB, Fig. S2). We further developed the global map of co-limitation (Fig. 5A) based on model sensitivity analysis, and quantified the contribution of key factors to local species richness patterns using variance partitioning (Fig. 6). Dotted boxes represent processes or models, and dashed ones represent data or results. See Methods for details.


Fig. 2 | Distribution of sample plots in the global forest biodiversity individual-based database (GFBi). GFBi data consists of in situ measurements of approximately 1,255,444 forest sample plots (Phase-I sample plots, blue). Together with an independent and complementary dataset of 22,131 forest sample plots (Phase-II sample plots, yellow), our ground-based sample data cover 97% of forested ecoregions worldwide (Fig. S1). GFBi data were collected from 96 national, regional, and local forest inventories, spanning 110 countries and territories (Table S1). For the assessment of tree species richness in local assemblages and its latitudinal gradient, each plot is underpinned by exhaustive community-wide tree survey records including species and diameter-at-breast-height or above buttress if there is a deformity.

Fig. 3 | A total of 47 explanatory variables in five categories (bioclimatic, vegetation and survey, topographic, anthropogenic, and soil) were used in random forests models to predict local tree species richness and quantify LDG. According to standardized variable importance values (horizontal bar plots to the left), bioclimatic variables contributed the most to LDG, followed by vegetation and survey, topographic, anthropogenic, and soil variables. The correlogram to the right illustrates correlations between any two variables by the color (color ramp represents the correlation coefficient) and size of a disk. The partial dependence plots to the left (next to the variable names, see Fig. S4 for details) show the effect of each predictor variable on the species richness, while all the other predictors remained constant at their sample mean. See Table S1 for a detailed description of the explanatory variables.


Fig. 4 | Estimated tree species richness per hectare in forested areas worldwide. (A) Tree species richness per hectare were first derived for the ca. 1.3 million GFBi plots across the world, and then imputed to the global forest extent. Curves (Top left) represent the empirical latitudinal diversity gradient (LDG, black) of tree species diversity in comparison with LDG (red) predicted by the metabolic theory of biodiversity (MTB) based on local mean annual temperatures (see Fig. S3). (B) Width of the 95% confidence interval (C.I.) for the estimated tree species richness per hectare. All map layers are displayed at a 0.025°×0.025° resolution with an equirectangular projection (Plate-Carrée) for better illustration of the latitudinal gradients.

Fig. 5 | Dominant drivers of tree species richness in forested areas worldwide. (A) Driver dominance was derived for each pixel from four driver categories (i.e., bioclimatic, topographic, anthropogenic, and soil), as well as co-limitation which represents a lack of clear dominance among the four foregoing categories. The pixel-level drivers were then aggregated by 0.5° latitudinal bins to show the percentage prevalence of dominant drivers by latitude (Top left). (B) The violin charts show the kernel probability density of tree species richness per hectare for different drivers, with inside boxes indicating the medium and interquartile range. The numbers on top of the violin charts indicate the percentage of forested pixels globally that corresponds to each driver category. Red line represents the mean and 95% confidence interval of tree species richness per hectare (81.1 ± 0.1) for all the 0.025°×0.025° pixels of co-limitation. The vertical axis is on a logarithmic-10 scale for better illustration.


Fig. 6 | Patterns and variance of local tree species richness per hectare by continent.  The collage of maps shows the zoomed-in view of the distribution of predicted local tree species richness per hectare (Fig. 4A) by continent. Circular Venn diagrams (with the legend in the center) show, for each continent, the spatial variance in observed tree species richness partitioned as follows: [a] (mean=14.3%) represents the fraction of variance uniquely explained by environmental factors (i.e. bioclimatic, topographic, anthropogenic, and soil variables), after latitudinal effects had been accounted for. [b] (mean=68.2%) stands for the fraction of variance jointly explained by environmental factors and latitudinal effects. [c] (mean=0.3%) represents the fraction of variance explained by latitudinal effects after removing environmental effects. [d] (mean=17.2%) represents the fraction of unexplained variance in tree species richness. The fractions were based on contrasting the amount of local richness variations in sample data from ~1.3 million plots explained by the R2 statistics from the continental-scale random forest models with the full set of factors versus those with targeted factors removed