Higher-order interactions mitigate direct negative effects on population dynamics of herbaceous plants during succession

Plant succession is regulated by a combination of abiotic and biotic factors. However, previous studies of biotic drivers have focused overwhelmingly on direct pairwise species interactions, ignoring the likely prevalent higher-order interactions (HOIs) in natural systems. Climate also plays a significant role in determining successional dynamics with both direct effects and indirect effects via altered biotic interactions. Here we explored the relative effects of direct species interactions, HOIs, climate, and their interactions on population dynamics of herbaceous plants during 50 years of post-agricultural secondary succession and tested whether the inclusion of HOIs and climate data improved forecasts of population dynamics. Direct intraspecific interactions were competitive and prevalent across the 90 herbaceous plants examined, while direct interspecific interactions only affected populations of 29% species. HOIs, mainly arose from intraspecific HOIs of conspecifics, were mostly positive and thus largely mitigated the competitive effects of direct intraspecific interactions. Species with lower peak cover experienced stronger intraspecific competition and positive intraspecific HOIs of conspecifics. Direct interspecific interactions had neutral or facilitative effects on species with lower peak cover, and tended to have competitive effects on species with higher peak cover. Climate simultaneously influenced population dynamics both directly and indirectly via altered species interactions. Forecast performance was significantly improved with the inclusion of HOIs or climate for about half and one-third of species, respectively. Our study emphasizes the importance of HOIs, which largely mitigated direct competitive effects on population dynamics of herbaceous plants during succession. Teasing apart HOIs from direct species interactions substantially refined our understanding of successional dynamics of herbaceous plants and improved the accuracy of forecasting population dynamics during succession in a changing world.


Introduction
Plant succession has been a central theme in ecology for more than a century, providing a fertile area for testing and developing multiple fundamental ecological theories (Clements 1916, Connell and Slatyer 4 These authors contributed equally to this work. 5 Author to whom any correspondence should be addressed. 6 Mailing address: Xingang Xi Road 135, Sun Yat-sen University, Guangzhou, 510275, China; Tel: (+86)2084111541. 1977, Tilman 1985, Koffel et al 2018, Png et al 2019, Poorter et al 2019. Refining our knowledge of plant succession is essential for understanding, managing and predicting anthropogenic impacts on natural systems (Clements 1935, Meiners et al 2015a, Chang and Turner 2019.
The trajectory of plant succession can be influenced by both abiotic and biotic factors. Climate change is one of the most important abiotic factors affecting contemporary biodiversity and ecological processes including plant succession.
Several studies have found associations between climate and secondary succession rates, leading to potential alterations of successional plant dynamics under climate change. For instance, increasing temperature accelerated the rate of secondary succession (Anderson et al 2006, Prach et al 2007, Fridley and Wright 2018, while the effects of precipitation varied between systems (Otto et al 2006, Martínez-Ramos et al 2018. Moreover, the individual responses of species to climate change vary considerably (Classen et al 2010, Clark et al 2011, Martínez-Ramos et al 2018, which may generate contingent impacts of climate on population dynamics during succession. In addition to the direct effects of climate on population dynamics during plant succession, there may also be indirect impacts of climate via altered biotic interactions (Klanderud 2005, Adler et al 2009, Ockendon et al 2014, Chu et al 2016. This potential suggests that it is necessary to simultaneously test the effects of biotic and abiotic factors on long-term species dynamics to predict plant community succession and development (Mantgem et al 2003, Clark et al 2011, Zhang et al 2015.
Intra-and interspecific interactions are often presumed to be responsible for species transitions during succession. Intraspecific competition can speed up species replacement during succession by reducing the performance of the plant or its offspring (van de Voorde et al 2011). Interspecific interactions may have either positive, negative or neutral effects on species replacement (Connell and Slatyer 1977, Hils and Vankat 1982, Armesto and Pickett 1986, Chang and Hillerislambers 2016. However, previous studies have focused overwhelmingly on direct pairwise interactions between species, with little attention on indirect species interactions. As indirect interactions have been empirically demonstrated to be prominent in natural communities (Miller 1994, Padilla et al 2013, Yu and D'Odorico 2015, Levine et al 2017, Mayfield and Stouffer 2017, Bartomeus and Godoy 2018, Gallien et al 2018, the knowledge gap of indirect interactions during succession must be addressed. Indirect interactions may arise in ecological systems through interaction chains or higher-order interactions (HOIs; Wootton 1993, Billick and Case 1994, Levine et al 2017. HOIs have been historically used to describe the phenomenon that the per capita effect of one species on another depended on intermediary species Peacor 2003, Levine et al 2017). Such HOIs have been well explored in food-webs, where HOIs usually occur between trophic levels (Benedetti-Cecchi 2000, Werner and Peacor 2003, Trussell et al 2006, Polis and Winemiller 2013, Greeney et al 2015. Recently, theoretical and empirical evidence of HOIs has been reported within one trophic level in plant communities (Padilla et al 2013, Bairey et al 2016, Grilli et al 2017, Mayfield and Stouffer 2017. HOIs can emerge when one of the two competing species has a plastic morphological or behavioural response to a third species, analogous to the trait-mediated indirect interactions described in the trophic literature Peacor 2003, Levine et al 2017). For example, plantain (Plantago lanceolata) may suppress root growth of red fescue (Festuca rubra, an efficient competitor for soil nutrients), and this plastic response will weaken F. rubra's per capita effect on other competitors (Levine et al 2017). Moreover, HOIs were extended to a wider phenomenon of non-linear density dependence (Bairey et al 2016, Mayfield andStouffer 2017, Letten andStouffer 2019), which were termed as soft HOIs (Kleinhesselink et al 2019). In this generalized definition, HOIs might be driven by interactions between individuals belonging to just one or two species as they were driven by interactions between individuals belonging to three separate species. In our study, we adopted the non-linear density dependence definition (soft HOIs) and for the first time tested their prevalence and importance in population dynamics of herbaceous plants during succession. Moreover, little attention has been paid to variation among species in the strength of direct and higher-order species interactions and the possible consequences of that variation in determining species abundances during succession. In this study, we aim to explore how the strength of species interaction varies across species and test whether that variation is related to differences in species' peak abundance, which would refine our understanding of how communities are structured during succession.
Due to the logistic difficulty of obtaining longterm data from a single system, most studies have concerned only the initial stages of succession (e.g. Stoll et al 1994, Fridley andWright 2018), or have relied on chronosequence-based inferences which may differ markedly from longitudinal dynamics within a site (Feldpausch et al 2007, Johnson andMiyanishi 2008). Therefore, studies of long-term successional responses to climate and species interactions are undoubtedly needed and urgent (Curzon et al 2017). In the present study, we used over 50 years of succession data from the Buell-Small Succession Study (BSS; Pickett et al 2001, Cadenasso et al 2009, to (1) quantify the relative effects of biotic factors (direct species interactions and HOIs), climate and their interactions on population dynamics of herbaceous plants, (2) test whether variation in the strength of species interactions among species is related to differences in species' peak cover in succession, and (3) test whether the inclusion of HOIs or climate improves forecasts of population dynamics of herbaceous plants during succession.

Data collection
Our study is based on the data collected from the abandoned agricultural land at the Hutcheson Memorial Forest Center (HMFC) in the Piedmont of New Jersey, USA (40 • 30 ′ N, 74 • 34 ′ W), as a part of the long-term Buell-Small Succession Study (figure S1, available online at: https://stacks.iop.org/ERL/15/074023/mmedia); Meiners et al 2015b, Duffin et al 2019. The study site comprises ten fields and 480 permanently-marked 0.5 m × 2.0 m plots (Meiners et al 2015b). Our analyses included 90 herbaceous plants with sufficient samples (the number of transitions with non-zero cover in year t and/or t + 2 was more than 200), consisting of 17 annuals, 15 biennials and 58 perennials (table S1, taxonomy follows Gleason and Cronquist 1991). Woody plants were not included in the focal species list, because the 0.5 m × 2.0 m plots are not large enough to quantify their demographic parameters, though they were included in the models of herbaceous species as competitors. The 90 taxa constituted 32% of the total richness and 98% of the total cover of herbaceous species that occurred during this period.
Climate data for the Buell-Small Succession Study were derived from the Climatic Research Unit (CRU) TS 4.02 gridded dataset (Harris et al 2014). Specifically, we retrieved monthly climate data from 1958 to 2013 and calculated ten climatic variables: annual mean cloud cover (cld), annual mean diurnal temperature range (dtr), annual frost day frequency (frs), annual mean potential evapotranspiration (pet), annual precipitation (pre), annual mean temperature (tmn), annual mean minimum temperature (tmp), annual mean maximum temperature (tmx), annual mean vapour pressure (vap), and annual wet day frequency (wet).

Statistical analyses 2.2.1. Fitting statistical models of population growth.
Population growth of each species was modelled as a function of direct species interactions, higher-order interactions and climate variables as well as the interactions between the climates and biotic interactions. The model takes the general form: (1) The response variable G i,q,g,t is the average annual population growth in terms of cover change of focal species i in quadrat q of field g at year t and calculated as the difference between the natural log cover of focal species i in quadrat q of field g at year t + 2 and year t divided by the time interval 2. If a species was not present at time t or time t + 2, the cover is assigned a small value of 10 −6 for log-transformation. In the population growth model (equation 1), γ is a time-dependent intercept, φ is the random effect of group defined as a set of plots located within one field to account for spatial location.
Direct intraspecific and interspecific interactions are quantified by conspecific (D i,q,g,t ) and heterospecific density dependence indices (D j,q,g,t ), with the strength of these effects given by α ii and α ij , respectively. D i,q,g,t is quantified as the cover of focal species i in the focal quadrat q at year t, and D j,q,g,t is quantified as the total cover of all heterospecifics, including woody and non-woody species, in the focal quadrat q at year t. It should be noted that the total quadrat cover (sum of D i,q,g,t and D j,q,g,t ) could be greater (more overlaps than gaps) or smaller (less overlaps than gaps) than 100%.
HOIs are defined as the quadratic densitydependent effects on the focal species (Bairey et al 2016, Mayfield andStouffer 2017, Letten andStouffer 2019). HOIs include intraspecific higher-order interaction indices (H i,ii , H i,jj ) and the interspecific higherorder interaction indices (H i,ij ), with the strengths of these effects given by β-terms. Specifically, H i,ii , the intraspecific HOIs of conspecifics, capture nonadditive effects of conspecific cover on population growth of focal species i: H i,jj , the intraspecific HOIs of heterospecifics, capture non-additive effects of heterospecific cover on population growth of focal species i: H i,ij , the interspecific HOIs, capture non-additive effects of both conspecific and heterospecific cover on population growth of focal species i: q,g,t D j,q,g,t. (4) Annual mean temperature (tmn) and annual precipitation (pre) are two key factors influencing oldfield succession (Fridley and Wright 2018, Martínez-Ramos et al 2018) and were included in the population growth model. Two other climate variables, annual mean diurnal temperature range (dtr) and annual mean vapour pressure (vap), that were not strongly correlated to annual mean temperature (tmn) and annual precipitation (pre) (r < 0.7; figure  S2), were also included in the model. Annual mean temperature (tmn) and annual precipitation (pre) increased significantly during the focused period of community succession, while diurnal temperature range (dtr) and vapour pressure (vap) did not change significantly (figure S3). To match the 2 year interval population growth, we calculated the daily mean temperature, average diurnal temperature range, total precipitation and mean vapour pressure for each 2 year sampling interval (from August of year t to July of year t + 2 where t to t + 2 is the transition of interest). In the population growth model (equation 1), η are the effects of the climate variables C, ω and ζ are interactions between the climate variables and direct interactions, λ, ψ and Γ are the effects of interactions of climate variables with HOIs (H i,ii , H i,jj , and H i,ij ), respectively. All climate variables were centered to zero mean and scaled to unit variance before analyses. The full model (equation 1) for each species was fitted using the 'lmer' in the R package 'lme4 ′ (Bates et al 2015). Year and group variables were treated as random factors allowing intercepts to vary among years and spatial locations. To avoid problem of overfitting, we compared the full model with a series of candidate models including only subsets of terms in full model using the Akaike's Information Criterion (AIC). The best-fit model was selected as the most complex model that reduced AIC by more than two units from the next simplest model, because the penalty for one additional parameter is +2 AIC units (Burnham and Anderson 2003). It should be noted that the best model selected for different species may contain different subsets of terms of the full model (table S2). For the best model of each species, the proportion of variance explained by fixed and random terms were estimated using marginal and conditional R 2 (table S2).
Following the same procedures, we also explored the effects of direct intraspecific and interspecific interactions, higher-order interactions, climate variables and their interactions with direct interactions and HOIs on population growth of herbaceous plants when grouping heterospecifics by life form (woody vs. non-woody). Detailed descriptions of the population growth model when grouping heterospecifics by life form are available in supplementary information (SI).

Quantifying model effects.
We used modified partial-residual plots to assess the magnitude and direction of the impact of any given set of explanatory variables on population growth (Mayfield and Stouffer 2017). To do this, we multiplied the observed values of explanatory variables by their corresponding regression coefficients in the best model. Detailed information is available in supplementary information.

Relating the strength of biotic interactions to species' peak cover.
We assessed the relationship between the peak cover achieved in succession for each species and the strength of direct species interactions and HOIs across species using linear regression models. For each species, we first arranged the cover data by field age (number of years since abandonment) and then calculated the mean cover of the species in all ten fields at each field age. This approach generated composite trajectories of species abundance during succession. Species' peak cover was the maximum percent cover of each species across all possible ages (Caplan et al 2019). Species' peak cover ranged from 0.09%-29.55%. Values of species' peak cover were log-transformed before regression. The strength of direct species interactions and HOIs was given by their corresponding regression coefficients in the best-supported model. We only included species for which the corresponding effects were significant.

Model validation.
For the species whose best population growth model includes HOIs or climate, we tested whether the inclusion of HOIs or climate significantly improved forecast of population growth. To test the ability of forecasting population growth of model, we performed leave-one-year-out cross-validation (Tredennick et al 2017) and calculated the forecast accuracy (correlation between forecasts and observations) and forecast error (mean absolute error between forecasts and observations). A model with greater forecast accuracy (or smaller forecast error) has better forecast ability. Therefore, we can test whether the model including HOIs or climate had greater forecast accuracy (or smaller forecast error) than the model without HOIs or climate, to determine if the inclusion of HOIs or climate improved forecast of population growth. Statistical tests for comparing forecast accuracy and error between models were conducted using the algorithms from (Ye et al 2015). Specifically, we compared forecast accuracy (ρ) between models with and without HOIs/climate using one-sided t test with standard errors calculated using the HC4 estimator from (Cribari-Neto 2004) and adjusted degrees of freedom following (Wilcox 2009). Comparisons of forecast error between models with and without HOIs/climate use a one-sided paired t test for the difference, treating each forecast as an independent sample. We assumed that inclusion of HOIs or climate did not improve forecast of population growth for species whose best supported model did not include HOIs or climate. All analyses were conducted in R (R Core Team 2018).

Effects of biotic interactions and climate on population growth
Direct intraspecific interactions had negative (competitive) effects on population growth for 87% of the species (figure 1(a)). In contrast, direct interspecific interactions had significant effects on population growth for only 29% of the species, and these effects were negative for the majority (except for 9%) of those species (figure 1(a)). The intraspecific HOIs of conspecifics (β i,ii ) were positively correlated to population growth for 93% of the species (figure 1(b)), while the intraspecific HOIs of heterospecifics (β i,jj ) and interspecific HOIs (β i,ij ) did not have significant and non-significant (ns), respectively. Within each population growth model, α ii represents direct intraspecific interactions; α ij represents direct interspecific interactions; β i,ii and β i,jj represent intraspecific HOIs, β i,ij represents interspecific HOIs; η1, η2, η3, and η4 represent the effects of mean temperature, precipitation, diurnal temperature range, vapour pressure, respectively; ω1, ω2, ω3, and ω4 represent the effects of interactions between direct intraspecific interactions and the four climate variables; ζ1, ζ2, ζ3, and ζ4 represent the effects of interactions between direct interspecific interactions and the four climate variables; λ1, λ2, λ3, and λ4 represent the effects of interactions between intraspecific HOIs of conspecifics and the four climate variables; ψ1, ψ2, ψ3, and ψ4 represent the effects of interactions between intraspecific HOIs of heterospecifics and the four climate variables; Γ 1, Γ 2, Γ 3, and Γ 4 represent the effects of interactions between interspecific HOIs and the four climate variables.
effects on population growth of most species (figure 1(b)). Annual mean temperature had significantly negative direct effects for 43% of the species, while the direct effects of other climate variables were relatively weak ( figure 1(c)). Interactions between the four climate variables and direct intraspecific interactions or intraspecific HOIs of conspecifics had significant positive or negative effects for 34% to 53% of the species, while interactions of climate with direct interspecific interactions, intraspecific HOIs of heterospecifics, or interspecific HOIs only affected a few species (figures 1(d) and (e)). The results based on the full model were qualitatively similar to those based on the best model ( figure S4). In addition, these results were qualitatively similar when grouping heterospecifics by life form ( figure S5). Cumulative effects of direct biotic interactions were predicted to generate substantial reductions in population growth rates for 77% of the species (tables 1 and S3). Cumulative effects of HOIs were positive for 58% of the species and greatly mitigated direct competitive effects on population dynamics, but the net biotic effects usually remained negative (tables 1 and S3). Cumulative effects of climate and their interactions with biotic interactions were usually mixed (tables 1 and S3). Table 1. Summary of the observed impacts of biotic effects, climatic effects and their interactions on the population growth rates of all 90 species. Each column shows percentages of species for which intervals between 2.5% and 97.5% quantiles of the observed impacts were positive, negative or overlapped with zero (Mixed). Mixed impacts indicate that the observed cumulative effects were either positive or negative across observations during succession. Total effects were the sum of biotic effects, climatic effects and their interactions. Net biotic effects were deconstructed into direct effects (Direct) and higher-order effects (HOI). Full details for each species were provided in

Relationships between the strength of biotic interactions and species' peak cover
Direct intraspecific competition (α ii ) had more negative and intraspecific HOIs of conspecifics (β i,ii ) had more positive effects on species with lower peak cover (figures 2(a) and (c); p < 0.001). Direct interspecific interaction (α ij ) was negatively related to species' peak cover, with facilitative effects on species with lower peak cover and competitive effects on species with higher peak cover (figure 2(b); p = 0.007). No significant relationships were detected between species' peak cover and the strength of intraspecific HOIs of heterospecifics (β i,jj ; figure 2(d)) and interspecific HOIs (β i,ij ; figure 2(e)).

Importance of HOIs and climate to population growth forecasts
The model including HOIs had greater forecast accuracy for 44% species and smaller forecast error for 47% species (figure 3(a)). Inclusion of climate variables was not as important as HOIs for forecast of population growth, but still improved forecast accuracy for 30% species and reduced forecast error for 32% species ( figure 3(b)).

Direct competitive effects on population growth
We demonstrated that the cumulative effects of direct species interactions were generally competitive (table 1). Direct intraspecific competition, i.e. con-specific negative density dependence (CNDD), was prevalent across herbaceous plants ( figure  1(a)). Moreover, the negative relationship between intraspecific competition and species' peak cover indicates stronger CNDD for rare species, which is consistent with the findings of other studies (Comita et al 2010, Mangan et al 2010, Johnson et al 2012, Xu et al 2015. CNDD commonly results from intense competition for shared resources and density-dependent, host-specific natural enemies (also known as the Janzen-Connell hypothesis; For each comparison, p-values are from one-sided t tests designed to assess whether the best model had higher accuracy or lower error than the model without HOIs or climate. Orange bars represent percentages of species for which the best models have higher accuracy or lower error than the best model without HOIs or climate. Blue bars represent percentages of species for which the best models have lower accuracy or higher error than the best model without HOIs or climate. Gray bars represent percentages of species for which no significant difference was detected between the best model and the best model without HOIs or climate. Detailed results for each species are provided in figures S6-7. Janzen 1970, Connell 1971. A study that included some of the species analyzed here, found that species' relative abundance was negatively related to the rate and degree that species accumulated host-specific soil pathogens (Klironomos 2002), which might be the underlying mechanism of CNDD detected in the present analyses. Direct interspecific interactions had a weaker effect on population dynamics relative to intraspecific competition ( figure 1(a)). In our study, direct interspecific interactions had primarily negative or neutral effects on population dynamics during succession ( figure 1(a)). Based on priority effects, (Connell and Slatyer 1977) distinguished three alternative mechanisms of succession: facilitation, tolerance and inhibition. Our results agree with removal experiments during old-field succession that favored the tolerance and inhibition models over the facilitation model (Hils andVankat 1982, Armesto andPickett 1986). Nonetheless, by grouping interspecific data as the total cover of all heterospecifics, we could miss strong positive or negative effects that might otherwise have appeared if the effect of each heterospecific was analyzed. Given computational limitations, parameterizing the effect of each heterospecific would have been intractable. Therefore, we also grouped heterospecifics by life form (woody vs. non-woody). As expected, significant interspecific interactions were detected for about half of the species (figure S5). Positive direct interspecific effects on population growth were detected for eight species with lower peak abundances (figures 1(a) and 2(b)). This finding, together with the strong intraspecific competition across species with lower peak cover, supports the hypothesis that rare species experience particularly strong negative frequency dependence, or stronger negative density dependence than common ones (

Positive higher-order effects of species interactions on population growth
We found that HOIs played a substantial role in driving population dynamics of herbaceous plants during old-field succession. To our knowledge, our study for the first time empirically tested the effects of HOIs on population growth during long-term plant succession. The intraspecific HOIs of conspecifics doative effects of HOIs were mostly positive, largely mitigating the direct competitive effects on population dynamics of herbaceous minated the contribution of HOIs to variance in population growth rates ( figure 1(b)). Of particular interest, the cumulplants (tables 1 and S3). Analytical studies have found that exploitative competition for resources with nonlogistic growth would generate non-linear density dependence (Abrams 1983(Abrams , 1987. As (O'Dwyer 2018) put it, 'to obtain something as familiar as logistic growth … we must assume logistic growth "all the way down".' However, this seems particularly unlikely for abiotic resources consumed and competed for by plants. A recent theoretical work on the mechanistic basis for HOIs found that for a single species utilizing a single resource under constant resource supply, the functional relationship between per capita growth rates and consumer densities is concave-up, leading to positive quadrative HOIs (Letten and Stouffer 2019), consistent with the results presented here. Light is the main limiting resource for herbaceous species within the BSS except for the first few years after abandonment (Meiners et al 2015b). About 90% of the focal species persisted for more than 30 years (table S1). Because the total incident light from outside the system remains constant every year, we infer that competition for light is a possible explanation for the positive HOIs observed. Moreover, we found that the positive intraspecific HOIs of conspecifics (β i,ii ) also had stronger effects on species with lower peak cover (figure 2(c)). The strength of intraspecific HOIs of conspecifics (β i,ii ) was negatively correlated with the strength of direct intraspecific interactions (α ii ; r = − 0.87, p < 0.001). This indicated that species that experienced stronger direct intraspecific competition generally tended to be affected by stronger positive effects of intraspecific HOIs of conspecifics. One possible explanation is that the strong positive intraspecific HOIs of conspecifics might result from strong intraspecific competition. Consider two individuals of a focal species. If the performance of these two individuals was negatively affected by other individuals of the focal species, the competitive direct interactions between these two individuals might be weakened. Similarly, the direct interaction between each pair of individuals of the focal species might be weakened by other individuals of the focal species, implying that intraspecific HOIs of conspecifics are positive. Further empirical studies are needed to justify this speculation.
Comparisons of the best model and the best model without HOIs showed that the best model outperformed (higher accuracy/lower error) the one without HOIs for almost half of the species (figure 3(a)), suggesting that the inclusion of HOIs significantly improved forecast performance for these species. Including HOIs decreased the forecast performance for only a small fraction of species (figures 3(a) and S6). Thus, accounting for HOIs substantially increased our capabilities to understand and forecast the population dynamics of herbaceous plants during succession. Previous experiments found that the three models proposed by (Connell and Slatyer 1977) would operate simultaneously during succession (Hils and Vankat 1982, Armesto and Pickett 1986, Pulsford et al 2016, Ulrich et al 2016. However, these studies have focused overwhelmingly on the direct pairwise interactions between successional species. Our results clearly demonstrated that positive HOIs, in addition to the direct pairwise interactions, are of critical importance to the course of succession. We argue that the three succession models are in fact not mutually exclusive, and likely to jointly function during succession (Meiners et al 2015b), even within the dynamics of individual species.

Climate influences population growth both directly and indirectly
Our results demonstrate that climate simultaneously influenced population dynamics during succession both directly, and indirectly via altered direct intraspecific interactions and the intraspecific HOIs of conspecifics. Temperature had direct negative effects on 43% of the species, implying that the population growth of these species was inhibited by increasing temperature (figure 1(c)). This result was consistent with biomass growth of herbaceous plants being inhibited in warmer sites (Fridley and Wright 2018). Tree seedling growth and survival are often reduced by the competitive effects of herbaceous plants during secondary succession (Wright andFridley 2010, Fridley andWright 2018). Therefore continued climate warming may promote tree establishment and accelerate the rate of transition from an herbaceous to woody-dominated ecosystem, as suggested by others (Anderson et al 2006, Prach et al 2007, Fridley and Wright 2018. Besides direct effects, climatic variables also exerted strong indirect impacts on population dynamics across succession for 34% to 53% of the species (figures 1(d) and (e)). Specifically, climate either strengthened or weakened direct intraspecific interactions and the intraspecific HOIs of conspecifics, and consequently modified the biotic interactions that affected species dynamics across succession. Indirect effects of climate have been a great source of uncertainty in attempts to predict the responses of communities to climate change Similarly, model comparisons showed that the best model outperformed (higher accuracy/lower error) the one without climate for about one-third of the species (figure 3(b)). This finding suggested that the inclusion of climate significantly improved forecast performance for these species. However, including climate variables decreased the forecast accuracy for 19% of the species and increased forecast error for 39% of the species, despite the significant climate effects detected ( figure 3(b)), possibly indicating idiosyncratic responses to climate change.
In conclusion, we teased apart the effects of HOIs from other variables on population dynamics of herbaceous plants during old-field succession, and found that the positive intraspecific HOIs of conspecifics were prevalent and important, largely counteracting the direct competitive effects of intraspecific interactions. Our results encourage additional studies to explore the underlying mechanisms driving HOIs during succession, such as trait plasticity, resource competition or mutualisms (Peacor and Werner 2001, Werner and Peacor 2003, Letten and Stouffer 2019. In addition, we found that climate simultaneously influenced population dynamics during succession both directly and indirectly. Forecast performance was significantly improved with the inclusion of HOIs or climate, which indicates that we should explicitly take indirect effects of biotic and abiotic factors into account to better understand and more accurately predict the processes of population dynamics.

Acknowledgments
We thank Serguei Saavedra for the constructive suggestion on the model validation. This research was funded by the National Natural Science Foundation of China (31925027, 31622014 and 31570426 to CC, and 31901106 to YL), the China Postdoctoral Science Foundation (2018M643295 to YL) and the Fundamental Research Funds for the Central Universities (17lgzd24 to CC). The BSS was supported by the National Science Foundation of the USA (most recently DEB-0424605 to SM). DBS acknowledges the support of a Marsden Grant (16-UOC-008), administered by the Royal Society of New Zealand Te Apārangi.

Data availability statement
The data that support the findings of this study are openly available at https://doi.org/10.5061/dryad.  Abrams P A 1983