Higher-order interactions enhance the latitudinal tree diversity gradient

The systematic decline in species diversity from the equator to the poles, known as the latitudinal diversity gradient, is among the most widely observed biogeographic patterns^1,7. The latitudinal decline in tree species diversity is particularly prominent and well documented^2,8. CNDD, whereby conspecific neighbours exert more negative effects on the performance of a focal tree than heterospecific neighbours⁹, is a considered to be a primary mechanism that maintains local diversity¹⁰ and promotes the latitudinal tree diversity gradient^11,12. The lemma that CNDD regulates latitudinal diversity gradient posits that CNDD is stronger in tropical than in temperate forests⁴. However, empirical evidence for this proposition is mixed: some studies based on field data across forest plots have reported significant declines in CNDD with increasing latitude¹² or decreasing species richness¹¹, whereas others did not support this trend^3,13,14. These inconsistent findings have led to a longstanding debate on whether there is a latitudinal pattern of CNDD and whether it contributes to latitudinal tree diversity gradient^{4,15,16,17,18}.

Although negative density dependence has been the main focus of interest in previous studies of the latitudinal tree diversity gradient, positive density dependence (facilitation) generated by mutualists (for instance, mycorrhizal fungi)¹⁹ could also play an important part in maintaining local species diversity^20,21 and contribute to the latitudinal diversity gradient^22,23. Moreover, existing studies have mostly been built on the assumption that organisms interact in a pairwise fashion^5,24, neglecting the potential effects of HOIs that emerge when pairwise interactions between two neighbouring trees are modified by other neighbours (Fig. 1). Specifically, the effect of a neighbouring tree of species j on a focal tree of species i in the absence of other trees, denoted ({alpha }_{{ij},text{true}}) (Fig. 1a), can be competitive (({alpha }_{{ij},text{true}} < 0)) or facilitative (({alpha }_{{ij},text{true}} > 0)), and conspecific ((j=i)) or heterospecific ((jne i)). This pairwise effect can be either strengthened (Fig. 1b) or weakened (Fig. 1c) in the presence of another neighbouring tree of species k. The higher-order effect of k (initiator) on species i (receiver) through j (transmitter) is denoted ({beta }_{{ijk}}). As a consequence, the effect of j on i in the presence of k, denoted ({alpha }_{{ij},text{modified}}) (implicitly incorporating ({beta }_{{ijk}})), can be stronger (Fig. 1b) or weaker (Fig. 1c) than ({alpha }_{{ij},text{true}}), depending on whether ({alpha }_{{ij},text{true}}) and ({beta }_{{ijk}}) have the same or opposite signs. As an example, Eucalyptus urophylla (initiator) can inhibit the root growth of Cryptocarya concinna (transmitter) through allelopathic effects²⁵, thereby weakening the competitive effect of C. concinna on its neighbouring species (receiver).

There is accumulating evidence that HOIs regulate the survival and growth of trees^26,27 and thus play a substantial part in the structuring of community assemblages^28,29,30. Consequently, estimating CNDD (when ({alpha }_{{ii}} < {alpha }_{{ij}} < 0)) without explicitly considering HOIs may evoke divergent latitudinal changes in CNDD and obscure the true contribution of CNDD to the latitudinal tree diversity gradient. To better understand the contributions of pairwise interactions and HOIs, we examined how their cumulative effects changed with species abundance and latitude. Pairwise interactions and HOIs may act as stabilizing forces promoting species diversity when their cumulative effects are negatively associated with species abundance (that is, rare species are less negatively or more positively affected). Such stabilizing effects may further reinforce the latitudinal tree diversity gradient if this cumulative effect–abundance relationship weakens with increasing latitude. It has remained unclear, however, whether HOIs are common in forests and how they contribute to the latitudinal tree diversity gradient.

In this study, we assembled census data from 32 large permanent forest plots spanning tropical to boreal forests (Fig. 2a and Supplementary Table 1) to address three key questions. (Q1) Are HOIs prevalent among trees across forest plots? (Q2) How do pairwise interactions (({alpha }_{{ij},{rm{true}}}) and ({alpha }_{{ij},{rm{modified}}})) and HOIs (({beta }_{{ijk}})) vary with latitude? (Q3) How do latitudinal changes in HOIs contribute to the latitudinal tree diversity gradient (Fig. 2b)? To answer these questions, we estimated both pairwise interactions and HOIs from demographic growth (for 1,543 tree species–plot combinations) and survival models (for 1,340 tree species–plot combinations). With these data, we built three types of growth and survival model: (1) null models with no biotic interactions, (2) pair-only models including only pairwise interactions and (3) HOI-inclusive models including both pairwise interactions and HOIs. The pairwise interactions estimated from pair-only models implicitly incorporate the effects of HOIs and thus represent ({alpha }_{{ij},mathrm{modified}}) (blue arrows in Fig. 1), whereas pairwise interactions estimated from HOI-inclusive models isolate HOIs and thus represent ({alpha }_{{ij},mathrm{true}}) (orange arrows in Fig. 1), the true interaction between two trees of species i and j in the absence of other neighbours. We then compared Akaike information criterion (AIC) values of the three types of model (Q1), tested whether the estimated pairwise interactions (({alpha }_{{ij},mathrm{true}}) and ({alpha }_{ij,mathrm{modified}})) and HOIs (({beta }_{{ijk}})) declined with latitude (Q2), and evaluated how the cumulative effects of pairwise interactions and HOIs on growth and survival changed with species abundance and latitude (Q3).

Our results demonstrate strong statistical support for the prevalence of HOIs across 32 large permanent forest plots, as evidenced by the AIC values of HOI-inclusive growth models being at least two units lower than those of the alternative models for 40% of the 1,543 species–plot combinations across all 32 forest plots (Fig. 3a). Moreover, HOI-inclusive growth models significantly improved predictions of tree growth rates compared with the alternative models (Fig. 3b). Although HOIs have been reported at local scales in forests^26,27 and other systems (for instance, bacteria–paramecium–protozoan⁷, amphibians³¹, microcrustaceans³², annual plants³³), our study shows that HOIs are ubiquitous in global forests. The widespread HOIs in forests highlight the importance and necessity of further exploration of their latitudinal patterns and ecological consequences.

Competitive ((alpha < 0,beta < 0)) and facilitative ((alpha > 0,beta > 0)) interactions were almost equally frequent across all the species–plot combinations (Extended Data Fig. 1). The average strength of true intraspecific pairwise interactions (({alpha }_{{ii},{rm{true}}})), both competitive and facilitative, declined rapidly with latitude (orange dots and lines in Fig. 4a). However, the average strength of true interspecific pairwise interactions (({alpha }_{{ih},{rm{true}}})), both competitive and facilitative, remains relatively constant across latitudes (orange dots and lines in Fig. 4b). These results are consistent with some previous work indicating that higher tree diversity is associated with stronger CNDD^11,12. HOI coefficients (except ({beta }_{{ihh}})), both competitive and facilitative, also declined with latitude (Fig. 4c–f).

Strong negative correlations were detected between ({alpha }_{{ii},{rm{true}}}) and ({beta }_{{iih}}) ((r=-0.96)) and between ({alpha }_{{ih},{rm{true}}}) and ({beta }_{{ihh}}) ((r=-0.89)) (Extended Data Fig. 2). This meant that both intraspecific (({alpha }_{{ii}})) and interspecific (({alpha }_{{ih}})) pairwise interactions were more strongly weakened by heterospecific neighbours (({beta }_{{iih}}:{alpha }_{{ii}}leftarrow h) and ({beta }_{{ihh}}:,{alpha }_{{ih}}leftarrow h)) than by conspecific neighbours. Consequently, the latitudinal signal of intraspecific pairwise interactions modified by HOIs (({alpha }_{{ii},{rm{modified}}})) became much weaker (blue dots and lines in Fig. 4a) or even disappeared (Supplementary Text 2.3), because ({beta }_{{iih}}) weakened ({alpha }_{{ii},{rm{true}}}) more strongly in low latitudes than in high latitudes (Fig. 4e). By contrast, interspecific pairwise interactions (({alpha }_{{ih},{rm{true}}})) were weakened uniformly by relatively constant ({beta }_{{ihh}}) across latitudes (Fig. 4f), resulting in relatively constant modified interspecific pairwise interactions (({alpha }_{{ih},{rm{modified}}})) (blue dots and lines in Fig. 4b). Our findings suggest that HOIs, by modulating pairwise interactions across latitudes, may reconcile the divergent observations of latitudinal changes in CNDD^3,11,12.

Previous studies have primarily focused on competitive pairwise interactions (that is, CNDD) and their contribution to the latitudinal tree diversity gradient^3,4,12. By contrast, our results show that facilitative pairwise interactions and HOIs are also common and strong (Fig. 4 and Extended Data Fig. 1), highlighting the necessity of considering their overall contributions (including both competitive and facilitative effects) to the latitudinal tree diversity gradient (Fig. 2b). Here we assessed the relative changes in growth rate caused by cumulative effects (both competitive and facilitative effects) of pairwise interactions and HOIs separately and examined how they changed with species abundance and latitude. The relative changes in growth rate caused by cumulative effects of pairwise interactions and HOIs both declined with species abundance, shifting from beneficial for rare species to detrimental for common species (Fig. 5). This pattern implies that both the cumulative effects of pairwise interactions and HOIs promote the growth rates of rare species but suppress those of common species, thereby potentially favouring species diversity at local scales. Moreover, the stabilizing effect of pairwise interactions changed little across latitudes (the P value for the interaction between abundance and latitude for the relative changes caused by pairwise interactions was 0.736; Table 1), whereas the stabilizing effect of HOIs became weaker (marginally significant change) towards higher latitudes (the P value of the interaction between abundance and latitude for the relative change caused by HOIs was 0.093; Table 1). Taken together, these findings suggest that latitudinal changes in HOIs enhance the latitudinal tree diversity gradient, whereas latitudinal changes in pairwise interactions contribute little to this gradient.

Our study documents the prevalence of HOIs and their potential contribution to the latitudinal tree diversity gradient. For accurate interpretation of our results, we note some limitations of our analyses. First, the present study estimated pairwise interactions and HOIs separately for tree growth and survival (Extended Data Figs. 3–5 and Extended Data Table 1), thereby overlooking their covarying effects on population dynamics. Future studies should integrate multiple demographic processes (such as survival, growth and recruitment) to estimate overall species interactions influencing population growth^34,35. Second, pairwise interactions and HOIs estimated at the neighbourhood scale should be upscaled to the community scale before their contribution to species diversity can be rigorously assessed^36,37. To stimulate further investigation, we conducted a preliminary upscaling analysis and explored the potential contribution of HOIs to species coexistence and the latitudinal diversity gradient based on structural stability and found that forest plots with higher species richness were less susceptible to loss of species (Supplementary Text 3). Third, we assumed that pairwise interactions would be modified by HOIs in a linear form, consistent with the findings of a broad range of previous theoretical^28,30,38 and empirical studies^33,39. However, a recent study explored alternative functional forms and found that nonlinear formulation (using exponential and sigmoid functions) more accurately predicted population trends of annual plants than the linear form⁴⁰; this warrants further investigation in forests.

How the latitudinal diversity gradient is maintained is a longstanding question in ecology and biogeography. We provide empirical evidence that HOIs are widespread in global forests. HOIs decline with latitude and weaken both the intraspecific (CNDD) and interspecific pairwise interactions, and thus enhance the latitudinal tree diversity gradient. Our findings help to resolve the debate on the role of pairwise CNDD in maintaining latitudinal diversity gradient and contribute to a shift from a paradigm of classical community ecology theory based solely on pairwise interactions towards a framework that the incorporates facilitations and HOIs.

Study sites and census data

Our study was based on multiple censuses of 32 large permanent forest dynamic plots, with an average plot size of 24.5 ha (range: 9–50 ha). Data were sourced from the Forest Global Earth Observatory (http://www.forestgeo.si.edu) and Chinese Forest Biodiversity Monitoring Network (http://www.cfbiodiv.org) (Fig. 2a and Supplementary Table 1). Plots span tropical to boreal terrestrial biomes with latitude ranging from 1.92° S to 61.30° N. All plots were established and censused several times following a standardized protocol⁴¹. In each census, all free-standing woody stems with a diameter at breast height (DBH) greater than 1 cm were tagged (unique ID), mapped (coordinates), identified (species identity), measured (DBH) and recorded (alive, dead or recruit). The census was repeated every 5 years to monitor forest dynamics (for instance, survival, growth and recruitment). Most plots were only censused twice. For a few plots with three or more censuses, we selected two consecutive censuses between 1998 and 2022 for analysis (Supplementary Table 1). Overall, we compiled data for more than 3 million trees of 5,000 species across the 32 study plots.

Growth and survival models with HOIs

We estimated species interactions from demographic growth and survival models using the compiled dynamic forest census data^3,27. For trees with more than one stem (that is, with multiple branches at height less than 1.3 m), we considered the survival and growth of the main stem. The growth of a focal tree f of a species i (({{rm{Growth}}}_{{i}_{f}})) was modelled as a function of its potential growth rate in the absence of neighbours (({G}_{i})), size (({{rm{DBH}}}_{{i}_{f}})), and neighbourhood pairwise (({mathrm{Pair}}_{{i}_{f}})) and higher-order effects (({{rm{HOI}}}_{{i}_{f}}))^27,42:

Similarly, the survival probability of a focal tree f of a species i (({{rm{Survival}}}_{{i}_{f}})) is modelled as²⁷:

where ({lambda }_{i}) is the intrinsic survival probability in the absence of neighbours. The inverse of diameter (({{rm{DBH}}}_{{i}_{f}}^{-1})) is included to model rapid decline in mortality rate with increasing diameter, whereas the terms ({{rm{DBH}}}_{{i}_{f}}) and ({{rm{DBH}}}_{{i}_{f}}^{2}) model the U-shaped senescence effect⁴³.

The pairwise effects of all neighbours on the focal tree ({i}_{f}) (({mathrm{Pair}}_{{i}_{f}})) can be decomposed into conspecific and heterospecific effects:

where ({alpha }_{{ii}}) and ({alpha }_{{ih}}) are intraspecific and interspecific pairwise interaction coefficients, and ({n}_{i,{i}_{f}}) and ({n}_{h,{i}_{f}}) are conspecific and heterospecific pairwise neighbourhood crowding indices, respectively. Here we take a mean-field approximation⁴⁴ by replacing the species-specific interaction coefficients with a constant (average) heterospecific coefficient (({alpha }_{{ih}}))^36,37. The higher-order effects of all neighbours and all neighbours’ neighbours on focal tree ({i}_{f}) (({{rm{HOI}}}_{{i}_{f}})) can be decomposed into four components:

where ({beta }_{{iii}}), ({beta }_{{iih}}), ({beta }_{{ihi}}) and ({beta }_{{ihh}}) are HOI coefficients associated respectively with intraspecific pairwise interactions modified by conspecific neighbours (({beta }_{{iii}}:{alpha }_{{ii}}leftarrow i)) or heterospecific neighbours (({beta }_{{iih}}:{alpha }_{{ii}}leftarrow h)), and interspecific pairwise interactions modified by conspecific neighbours (({beta }_{{ihi}}:{alpha }_{{ih}}leftarrow i)) or heterospecific neighbours (({beta }_{{ihh}}:{alpha }_{{ih}}leftarrow h)); and ({n}_{{ii},{i}_{f}}), ({n}_{{ih},{i}_{f}}), ({n}_{{hi},{i}_{f}}) and ({n}_{{hh},{i}_{f}}) are the corresponding higher-order neighbourhood crowding indices. The conspecific and heterospecific pairwise neighbourhood crowding indices ({n}_{i,{i}_{f}}) and ({n}_{h,{i}_{f}}) add up the crowding contributions of all conspecific and heterospecific neighbours, respectively, around a focal tree within a given radius. The contribution of a neighbour is directly proportional to its size and inversely proportional to the distance to the focal individual (equations S3 and S4 in Supplementary Text 1). The contribution of a neighbour to the corresponding higher-order crowding indices ({n}_{{ii},{i}_{f}}) and ({n}_{{hi},{i}_{f}}) is further multiplied by its conspecific crowding index (equations S7 and S9 in Supplementary Text 1), and that to ({n}_{{ih},{i}_{f}}) and ({n}_{{hh},{i}_{f}}) is further multiplied by its heterospecific crowding index (equations S8 and S10 in Supplementary Text 1). The calculation of pairwise and higher-order neighbourhood crowding indices is summarized in Supplementary Text 1 according to our previously developed method²⁷.

Model fitting and evaluations

In each forest plot, we fitted demographic growth and survival models for each species with more than 100 trees (and also with at least 20 surviving and 20 dead trees for the survival model³) to ensure model performance and robustness. Overall, we fitted growth and survival models for 1,543 and 1,340 tree species–plot combinations, respectively. We log-transformed equation 1 to linearize the growth model and reduce model heteroscedasticity and residuals. We used logistic regression for tree survival (0 for death and 1 for alive). Overall, three classes of demographic models were constructed to determine the importance of tree size, pairwise interactions and HOIs: (1) null models without inclusion of biotic interactions, (2) pair-only models including pairwise interactions only, and (3) HOI-inclusive models including both pairwise interactions and HOIs. The pairwise interaction coefficients estimated from HOI-inclusive models were ({alpha }_{{rm{true}}}) (orange arrows in Fig. 1), because the effects of HOIs were isolated. By contrast, the pairwise coefficients estimated from pair-only models were ({alpha }_{{rm{modified}}}) (blue arrows in Fig. 1), because they implicitly incorporated effects of HOIs. Then, we compared the AIC values of the three models. The inclusion of HOIs was statistically supported when the AIC of the HOI-inclusive model was at least two units lower than those of the alternative models (Q1).

Latitudinal changes in biotic interactions

To explore the latitudinal trend in biotic interactions (Q2), we fitted exponential relationships between the estimated coefficients of pairwise interactions (({alpha }_{{ii}}) and ({alpha }_{{ih}})) and HOIs (({beta }_{{iii}}), ({beta }_{{iih}}), ({beta }_{{ihi}}) and ({beta }_{{ihh}})) and absolute latitude for competitive ((alpha < 0,beta < 0)) and facilitative ((alpha > 0,beta > 0)) interactions separately, given that they were nearly equally frequent (Extended Data Fig. 1). We further evaluated how HOIs (({beta }_{{ijk}})) modified the true pairwise interactions (({alpha }_{{ij},{rm{true}}})) by examining their correlations.

Contribution of biotic interactions to latitudinal diversity gradient

To further assess how the latitudinal change in pairwise interactions and HOIs might contribute to the latitudinal tree diversity gradient (Q3), we first calculated the relative change in growth rate and survival probability caused by the neighbourhood cumulative effects of true pairwise interactions (({{rm{e}}}^{{{rm{Pair}}}_{{i}_{f}}})) and HOIs (({{rm{e}}}^{{{rm{HOI}}}_{{i}_{f}}})) separately for each focal tree (equations 1–4)^27,33 and then averaged the relative change across trees for each species (({{rm{RC}}}_{{rm{Pair}}}={sum }_{f=1}^{{N}_{i}}{{rm{e}}}^{{{rm{Pair}}}_{{i}_{f}}}/{N}_{i}), ({rm{R}}{{rm{C}}}_{{rm{HOI}}}={sum }_{f=1}^{{N}_{i}}{{rm{e}}}^{{{rm{HOI}}}_{{i}_{f}}}/{N}_{i})). Finally, we explored how the relative changes (({mathrm{RC}}_{mathrm{Pair}}) and ({rm{R}}{{rm{C}}}_{{rm{HOI}}})) varied with species abundance per hectare and absolute latitude using linear regression. The relative changes and species abundance were log-transformed to improve normality.

Results from the survival models (Extended Data Figs. 3–5 and Extended Data Table 1) were qualitatively consistent with those from the growth models. Therefore, we only reported the growth model results in the main text. Moreover, as demonstrated in Supplementary Text 2, the growth model results were robust to variations in parameter settings (for instance, for different neighbourhood radii), inclusion of spatial autocorrelation (with or without quadrats as random effects), uncertainty in estimation of interaction coefficients, and the distinction between small (DBH < 10 cm) and large (DBH ≥ 10 cm) trees. All analyses were conducted in R v.4.4.1 (ref. ⁴⁵). The map in Fig. 2 was generated using R packages ggplot2 (v.4.0.0) and ggrepel (v.0.9.6).

Reporting summary

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