Biomass-density relationships of plant communities deviate from the self-thinning rule due to age structure and abiotic stress

1 –––––––––––––––––––––––––––––––––––––––– © 2020 Nordic Society Oikos. Published by John Wiley & Sons Ltd is is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Sa Xiao Editor-in-Chief: Dries Bonte Accepted 20 May 2020 00: 1–11, 2020 doi: 10.1111/oik.07073 A pertinent debate in plant ecology centers around the generality of the self-thinning rule. However, studies focused on highly simplified settings such as even-aged monospeci c populations or optimal conditions. is neglects the fact that most natural communities, to which the classical self-thinning slope is often applied, are age-structured, composed of multiple species and exposed to various types of abiotic stress. With the help of an individual-based model, we relax these simplified assumptions and systematically test for changes in the biomass–density relationships of unevenaged, functionally diverse plant communities across a complete stress gradient, using excessive to insu cient soil water as a case study. We show that frequent recruitment, which resulted in an uneven-aged community, and stress intensity caused predictable changes in the entire biomass–density trajectory. Increasing stress resulted in steeper (more negative) slopes and increased the intercept in the classical self-thinning section irrespective of excessive or insu cient soil water as a stress type. Recruitment steepened the slope, too and enabled a novel section in the biomass–density trajectory. is novel section represented a quasi-steady state of the density-dependent dynamics of new generations which occurred locally within patches of recruitment. At the community level, the slope of the biomass–density relationship at quasi-steady state had a signi cantly atter slope of −1.1 under optimal soil water conditions. Functional diversity showed little impact on density-dependent mortality. Namely, it resulted in an earlier onset of mortality but not in changes in the values of the slope and intercept. We conclude that the classical −3/2 slope is not useful to describe the biomass–density relationship in natural and semi-natural plant communities. The magnitude and direction of variation in the slope are related to the age-structure and abiotic stress intensity in the plant community.


Introduction
A key principle in plant population ecology is the selfthinning rule which relates changes in mean plant biomass to plant density through time (Westoby 1984, Weiner andFreckleton 2010). The rule predicts that this relationship forms a straight thinning line with a slope of −3/2 when the logarithm of mean surviving biomass w is plotted against the logarithm of density N (Yoda et al. 1963, Westoby 1984, Weiner and Freckleton 2010. Usually, the rule is expressed as an allometric power law w = kN −ɑ , where k and ɑ are the self-thinning intercept and slope, respectively. The rule was originally designed and validated for the simplest case, i.e. crowded even-aged monospecific populations undergoing density-dependent mortality, and a diversity of alternative values for k and ɑ may emerge if any of these assumptions is violated (Weller 1989, Lonsdale 1990, Pretzsch 2006, Strigul et al. 2008). Yet, ɑ is frequently used for management decisions, productivity assessments and carbon budgeting in natural and semi-natural habitats, i.e. uneven-aged multispecies communities (Cao 1994, Pretzsch and Biber 2005, Stankova and Shibuya 2007, Luyssaert et al. 2008, Newton 2012. There is broad consensus that the general process of self-thinning holds in these communities, too (Kohyama 1992, Schulze et al. 2005, Vospernik and Sterba 2015. However, a conceptual unification of the aspects concerning the specifics of the self-thinning rule, especially the values of k and ɑ, in natural plant communities is still missing. For example, the thinning slope may be steeper (Schulze et al. 2005), equal to the classical −3/2 (Luyssaert et al. 2008) or vary unpredictably (Kohyama 1992) in uneven-aged multispecies communities.
Another fundamental aspect differentiating natural and semi-natural communities from even-aged monospecific populations is the ubiquity of abiotic stress (Grime 1977). Abiotic stress can impede key performances such as survival and growth (Grime 1977), and consequently influence selfthinning. For example, the intercept and slope of the thinning line may vary under resource deficit (Morris 2003). The existing theoretical and empirical studies on variations due to abiotic stress are based on even-aged monospecific populations (Morris 2003, Chu et al. 2010, Lin et al. 2014. Interestingly, even in this simplified setting, the intensity and direction of variations in the intercept and slope appear essentially unpredictable. Specifically, six different patterns of responses of the intercept k and slope ɑ to resource deficit have been observed (reviewed by Morris 2003): 1) both k and ɑ stay constant; 2) both k and ɑ decrease; 3) k decreases but ɑ stays constant; 4) k increases but ɑ stays constant; 5) k increases but ɑ decreases; or 6) k decreases but ɑ increases. The first pattern is referred to as the 'common line' pattern, whereas the other patterns are 'different line' patterns (Morris 2003). However, resource deficit is not necessarily the only abiotic stress type in nature, making it likely that predicting k and ɑ may be even more complex than the existing studies suggest. For example, plant physiological responses along below-and aboveground resource gradients usually follow a unimodal relationship, i.e. both a deficiency as well as a surplus of resources hamper plant performance (Austin 1990). Soil water is an excellent example for such a unimodal relationship since both insufficient and excessive water availability reduce plant performance (Feddes et al. 1978, Moeslund et al. 2013, Silvertown et al. 2015. Namely, water is a key resource and can simultaneously act as a disturbance (flooding) (Grime 1977, Silvertown et al. 2015. Furthermore, water is a measure of other soil conditions such as oxygen concentration and nutrient availability (Yang andJong 1971, Sauter 2013). Thus, studies along a complete water stress gradient could yield important insights into responses of the intercept and slope to abiotic stress and the generality of the responses irrespective of the precise mechanism which impedes plant performance. A previous experiment indicated no significant difference in intercept and slope between moderate water deficit and well-watered conditions (Liu et al. 2006). However, there is no study which systematically examines changes in the intercept and slope along a complete stress gradient.
Another shortcoming of using the allometric coefficient ɑ for management decisions or carbon budgeting in natural habitats is that ɑ is only valid for the straight thinning line, the self-thinning section (Yoda et al. 1963, Westoby 1984. A focus on the self-thinning section misses important processes along the entire time course of the biomass-density trajectory (hereafter: 'bdt') (Begon et al. 1991, Berger andHildenbrandt 2003). Classical bdts additionally include a fast growing section prior to self-thinning and a post-thinning section (Begon et al. 1991, Berger andHildenbrandt 2003). Even moderate water deficit can delay and retard the self-thinning process indicating that abiotic stress alters the duration of all sections of the bdt (Liu et al. 2006).
Our study systematically relaxes a number of assumptions of previous studies with the aim to derive predictable patterns for the entire biomass-density trajectory. Specifically, we study the duration and characteristics, including the values of k and ɑ, of different sections of the entire biomassdensity trajectory within an uneven-aged plant community along a complete soil water stress gradient using an individual-based model.

Material and methods
The individual-based model PLANTHeR PLANTHeR (PLAnt fuNctional Traits Hydrological Regimes) is a generic spatially explicit individual-based model designed to evaluate the relative importance of soil water for the functional and structural composition of plant communities (Herberich et al. 2017). Despite its generic design, PLANTHeR is capable of reproducing well excepted theoretical and experimental patterns such as the theory of competitive exclusion (Gause 1934) and hydrological niche segregation (Silvertown et al. 2015, Herberich et al. 2017).
Here we only present the ODD (Overview, Design concepts, Details, Grimm et al. 2006Grimm et al. , 2010 protocol elements 'Overview', and 'Initialization'. The 'Design concepts' and 'Details' can be found in Supplementary material Appendix 1, the parameter values in Supplementary material Appendix 2. PLANTHeR comprises three entities: grid cells, and two types of individuals (adult plants and viable seeds). Grid cells have only a single state variable -soil water potential Ψ (mm). Adult plants are characterized by the state variable biomass (mg). Both types of individuals are classified by their position, age and plant functional type (PFT) which is based on the parameterization of six functional traits. These are general, composite functional traits which can represent several alternative traits depending on the specific ecosystem. Two traits are related to competitive ability (seedling competitive ability, maximum adult growth rate), two traits represent spatiotemporal colonization abilities (seed dormancy, seed dispersal distance), one trait affects longevity (life form) and one trait affects the response to abiotic stress (water stress tolerance). Each functional trait is represented by two opposing strategies: strong(C)/weak(c) seedling competitive ability, high(G)/low(g) maximum adult growth rate, long-(S)/shortterm(s) seed dormancy, long(D)/short(d) seed dispersal distance, perennial(P)/annual(a) life form and high(T)/low(t) water stress tolerance (Seifan et al. 2012(Seifan et al. , 2013. PLANTHeR assumes no a priori defined trait tradeoffs so that full factorial combinations of the two opposing strategies of the six traits result in 64 PFTs (Supplementary material Appendix 2 Table  A3). Except for the functional traits life form and water stress tolerance (Supplementary material Appendix 1 Sect. A2.1, A2.7), the functional trait strategy indicated by the capital letter stands for 80% of the maximum trait value whereas the strategy indicated by the lowercase letter stands for 20%, respectively (Supplementary material Appendix 2 Table A1). The exact trait values are drawn for each individual separately from truncated normal distributions (range: 0-100% of trait values) with the PFT strategy as mean and 10% of the maximum trait value as standard deviation (SD). Seed dormancy, seedling competitive ability and maximum adult growth rate are probabilistic traits with a maximum of 100%. Maximum seed dispersal distance is 20 m which corresponds to the maximum dispersal distance of most plant species dispersed by abiotic dispersal vectors (Thomson et al. 2011).
The cell size of 5 × 5 cm corresponds to the size of an average adult plant typical for many temperate herbaceous plant communities (Schippers et al. 2001). The model landscape is a grid of 800 × 800 cells with periodic boundary conditions. One time step equals one year in which individuals can potentially finish a complete life cycle (Fig. 1).

Soil water scenarios
Soil water scenarios were chosen based on the reduction functions f T (Ψ) and f t (Ψ) to represent gradients of water stress intensities for high and low water stress tolerant PFTs (Supplementary material Appendix 1 Sect. A2.1, Supplementary material Appendix 2 Table A1). Here, only cases with negligible seasonal and inter-annual soil water potential variability are considered. The growing season is thus assumed to be statistically homogenous and modelled in 19 soil water potential scenarios Ψ (mm) ( Table 1). Each of the scenarios was repeated ten times.
At soil water potentials between −300 and −2000 mm, the values of both reduction functions f T (Ψ) and f t (Ψ) equal one, i.e. no PFT experiences stress. Therefore, Ψ = −1900 mm was chosen to represent a 'base model'.

Initialization
Similar to previous self-thinning simulation studies (Chu et al. 2010, Lin et al. 2014, the model is initialized with a random spatial distribution of individuals whereby each grid cell has a 50% chance of being occupied by an adult plant. The specific individual identity is then determined randomly such that each of the 64 PFTs has the same probability to occupy a cell. This resulted in nearly even initial distributions of the 64 PFTs and functional trait values. We initialized ten random communities, i.e. one for each replicate, which were then duplicated 19 times (resulting in 19 identical communities) so that the outcomes of soil water scenarios could be compared. Irrespective of the PFT, the initial size of the individuals' zone of influence (ZOI) is one cell. ZOI (cm 2 ) is a circular area from which individuals could potentially acquire water (Supplementary material Appendix 1 Sect. A2.1). ZOI is allometrically related to the individuals' aboveground dry biomass [mg], B, as ZOI ~ μB (2/3) with the proportionality constant μ = 1 (Weiner et al. 2001) resulting in 125 mg initial biomass.

Composition of plant communities and simulation time
A detailed analysis of the functional composition of the simulated plant communities can be found in Herberich et al. (2017). Here, we only outline those aspects which may affect the biomass-density relationship.
Communities were initialized with all 64 PFTs at high evenness that is the highest possible diversity. Diversity decreased over time, i.e. in the course of the biomass-density trajectory, so that a maximum of two PFTs coexisted at the end of the simulation time (Supplementary material Appendix 3 Fig. A1). This simulated diversity dynamic resembles the dynamics in species richness in natural plant communities following a large-scale high intensity disturbance in otherwise spatially and temporally relatively homogenous habitats (Mackey and Currie 2001). The decrease in species richness (PFT diversity) with increasing time since last disturbance (initialization) results from competitive exclusion (Gause 1934, Huston 1979. Adding temporal and/or spatial heterogeneity to the simulations may promote PFT coexistence and thus, maintain diversity. In line with field observations of hydrological niche segregation (Silvertown et al. 2015), successful PFTs and functional trait strategies segregated into specific ranges of water stress intensities. Notably, none of the PFTs was successful under all scenarios, thus no functional trait combination was consistently beneficial. Instead, the strength of functional trait correlations and tradeoffs was a function of water stress intensities. The rate of diversity loss was generally comparable across the scenarios with PFT richness and Shannon diversity index stabilizing latest after 1500 years (Supplementary material Appendix 3 Fig. A1). We therefore chose 2000 years as simulation time.

Statistical analysis
In this study, we investigate the biomass-density relationship through time ('time as the independent variable'). This is not to be confused with interspecific scaling studies which correlate biomass and density at different sites ('no independent variable', Enquist et al. 1998, West et al. 2009, Weiner and Freckleton 2010. Each bdt section corresponds to a specific phase in the development of size structure such as the emergence of size hierarchies in the fast growing section (Begon et al. 1991). An objective selection of data points for each bdt section may thus be based on measures for size structure such as the skewness (third-order moment) of a size distribution (Gates et al. 1983, Berger andHildenbrandt 2003). The sequential calculation of the skewness over time provides markers to distinct the bdt sections (Berger et al. 2002, Berger andHildenbrandt 2003). We divided the bdt a posteriori into different sections by means of the skewness of the biomass distribution. Thereby, bdt sections during which communities underwent self-thinning are characterized by a decay of positive skewness (sensu Berger et al. 2002). For each self-thinning section identified with that approach, regression lines were fitted using a reduced major axis analysis (RMA) (Warton et al. 2006). We tested with a one-sample t-test whether the slopes differed from −3/2 and whether the slopes and intercepts of each soil water scenario differed from those of the base model. Pearson correlations were used to investigate the relationship between water stress intensity and the duration of the bdt sections or the values of k or ɑ. All statistical analyses were performed in R ver. 3.1.2 (<www.r-project.org>). Table 1. 19 simulated soil water potential scenarios Ψ (mm) representing different stress intensities for high water stress tolerant (f T (Ψ)) and low water stress tolerant (f t (Ψ)) PFTs. The targeted stress intensity of the respective PFT strategy is highlighted in bold lettering.

Results
The biomass-density trajectories (bdts) could be subdivided into different sections indicated by the skewness of the biomass distribution for all 19 soil water scenarios (Fig. 2, Supplementary material Appendix 4). When the soil water potential Ψ was within the range where f T (Ψ) and f t (Ψ) differed, the duration and characteristics of bdt sections including the values of k and ɑ were dominated by the strategy for water stress tolerance with the higher f(Ψ), from now on called dominant reduction function f d (Ψ) (Supplementary material Appendix 4 Fig.  A3). That is because independent of similar initial densities among the initial 64 PFTs, the best competitor for soil water dominated the model landscape (Herberich et al. 2017).
When Ψ was within the range where f d (Ψ) = 1.0, i.e. the dominant PFT strategy did not experience any stress, the duration of bdt sections as well as the values of k and ɑ were not significantly different from the base model (Fig. 2, 3, Supplementary material Appendix 4 Fig. A8-A11). In these scenarios, the first maximum and both x-intercepts of the skewness subdivided the bdts into four distinct sections. Section I was characterized by an increase of positive skewness and ended with the first maximum skewness. In this section, bdts showed a non-linear relationship between log w and log N (Fig. 2b, Supplementary material Appendix 4 Fig. A8b-A11b). Section II exhibited a constant decay of positive skewness between the first maximum and the first x-intercept. In this section, log w increase and log N decrease were linearly related, corresponding to the classical thinning line. The fitted slopes 6 were marginally but significantly flatter than the classical −3/2 slope (Fig. 3). Section III showed a negative skewness between the two x-intercepts. Here, the bdts deviated from the linear thinning line. Section IV started at the second x-intercept with a positively increasing skewness and was characterized by periodic changes in skewness. Over time, this periodicity was less pronounced and less synchronized among simulation replicates. An increase in density and a decrease in mean biomass was accompanied by an increase in skewness and vice versa. Bdts in this section had a linear relationship between log w and log N with values of the slopes that were significantly flatter than −3/2 (Fig. 3). The transition between section III and IV corresponded approximately to the point in time of maximum yield (base model 235 ± 25 years) (Supplementary material Appendix 5).
With increasing stress intensity for both excessive and insufficient water availability; i.e. decreasing f d (Ψ), the duration of all bdt sections as well as the values of k and ɑ changed in a consistent manner (Supplementary material Appendix 4 Fig. A3). Namely, there were significant negative correlations between f d (Ψ) and the duration of section I (Pearson's r = −0.53, p = 0.020, n = 19), the duration of section II (Pearson's r = −0.95, p < 0.001, n = 15) and the duration of section IV (Pearson's r = −0.51, p = 0.026, n = 19). In contrast, the duration of section III decreased with decreasing f d (Ψ) Figure 3. Changes in mean slope ± SE and mean intercept ± SE of the biomass-density trajectory sections II (classical self-thinning section) and IV (density-dependent dynamics of new generations) with increasing water stress for both excessive and insufficient water availability, i.e. decreasing f d (Ψ), for both high (f T (Ψ)) and low (f dt (Ψ)) water stress tolerant PFTs. p-values indicate differences to the values of the base model, i.e. the scenario in which no PFT experienced any stress: × p ≥ 0.01. • p < 0.01 (one-sample t-test).
(Pearson's r = 0.98, p < 0.001, n = 7) so that if f d (Ψ) < 0.9 this section disappeared entirely (Supplementary material Appendix 4). In section II, with decreasing f d (Ψ) the values of ɑ significantly decreased (Pearson's r = 0.76, p = 0.001, n = 15) whereas the values of k significantly increased (Pearson's r = −0.76, p = 0.001, n = 15) (Fig. 3). The values of k and ɑ of section II were not significantly different from those of the base model if f d (Ψ) > 0.8 (common line pattern) but k and/or ɑ were significantly different from the base model if f d (Ψ) ≤ 0.8 (different line patterns) expect for Ψ = −41.0 E + 3 (Fig. 3). In section IV, there was no significant correlation between f d (Ψ) and the values of either ɑ (Pearson's r = 0.00, p = 0.987, n = 19) or k (Pearson's r = 0.41, p = 0.080, n = 19) (Fig. 3). The values of k and/or ɑ of section IV were significantly different from the base model if f d (Ψ) < 1.0 (different line patterns) except for f d (Ψ) = 0.4 (Fig. 3).

Discussion
Our overall results indicate that the allometric coefficient ɑ of the self-thinning rule is not generally applicable to represent the biomass-density relationship of natural plant communities. Instead the entire biomass-density trajectory (bdt), including the intercept and slope of the classical thinning line, changed as a function of abiotic stress intensity and the inclusion of recruitment. Thus, there is 'no slope that fits them all', i.e. a fundamentally valid thinning slope cannot be specified. While previous studies have discussed the generality of the classical thinning slope at length (Weller 1989, Lonsdale 1990, Osawa and Allen 1993, West et al. 1997, Enquist et al. 1998, White et al. 2007), our study aimed at systematically addressing the consequences of relaxing a number of assumptions unrepresentative of natural plant communities.
In the following, we first briefly highlight the similarities between our model results and that of previous studies on variations in the thinning line due to abiotic stress and then discuss the deviations and their reasons.
First, the growth reduction in stressed communities delayed the onset of self-thinning (longer duration of section I) as well as retarded the process itself (longer duration of section II) (Liu et al. 2006, Zhang et al. 2017. Second, under minor stress, the values of k and ɑ were not significantly different from those under optimal conditions, i.e. in these scenarios the common line pattern emerged (Kays and Harper 1974, Lonsdale and Watkinson 1982, Liu et al. 2006. Third, different line patterns appeared when stress by both excessive and insufficient water availability was serious (Morris 2003).
Interestingly, this shift from the common line pattern under minor stress to different line patterns under serious stress was not restricted to resource deficit as in previous studies (Morris 2003), but we could show that it also holds for other types of abiotic stress such as excessive water availability. Thereby, the intensity and direction of variation in k and ɑ as well as in the duration of bdt sections under serious stress could be predicted by the stress intensity irrespective of the stress type, a result that has not been observed before. Namely, the studies that addressed abiotic stress gradients, either focused on unidirectional stress gradients, e.g. resource deficit (Morris 2003, Liu et al. 2006, Lin et al. 2014 or nonresource abiotic stress (Chu et al. 2010, Zhang et al. 2017). Our finding resulted from the design of soil water scenarios by means of defined reduction in plant performance. This design integrates differences in physiological mechanisms underlying avoidance or tolerance strategies as well as differences in response rates towards changes in Ψ (Feddes et al. 1978, Feddes andRaats 2004). Thus, while the variations in k and ɑ were mirrored at optimal conditions if plotted against stress intensity (f(Ψ)), they were not mirrored if plotted against Ψ. The values of k and ɑ at both ends of the gradient resulted from the extreme growth reduction under serious abiotic stress which lead to virtually density independent community dynamics. The density independence is further indicated by a decrease of the density range (log 10 (N min /N max )) with increasing abiotic stress, i.e. less change in density over time. Previous studies on even-aged monospecific plant populations stated that populations with a small density range lack sufficiently long growth periods and can thus be considered in an early stage of self-thinning (Westoby 1984, Lonsdale 1990, Zhang et al. 2005. Our findings suggest that growth reduction due to abiotic stress can cause a similarly small density range in plant communities with corresponding effects for k and ɑ even over long growth periods.
Serious stress by both excessive and insufficient water availability steepened the slope of the classical self-thinning section while the direction of the variations in the slope in section IV changed depending on the stress intensity. A set of interspecific scaling studies across a moisture gradient showed that with decreasing precipitation the slope of the biomass-density relationship flattens (Deng et al. 2006, Dai et al. 2009, Bai et al. 2010, Zhu et al. 2015. However, lack of information, for example on actual drought stress intensity or positive plant interaction, precludes us from further comparisons. Here, we studied patterns of responses to abiotic stress in the entire bdt of uneven-aged plant communities using excessive and insufficient soil water as an example. We therefore assumed other abiotic factors such as nutrients or light to be optimal. In natural plant communities, however, deviations from optimality for plant growth are the norm (Grime 1977, Körner 2003. We suggest that the above described patterns of responses of the bdt to soil water stress might generally hold for deviations from optimality regardless of the stress type(s) causing the deviation. Two main findings strengthen this suggestion. First, several of the emergent variations in the classical thinning line due to soil water stress corroborated experimental evidence on variations due to nutrient (Morris 2003) and light deficiency (Lonsdale and Watkinson 1982). Second, the variations in the entire bdt sections including the values of the thinning intercept k and slope ɑ were a function of the stress intensity and irrespective of the type of soil water stress.
In addition to the prominent effects of stress, the inclusion of recruitment and functional diversity, i.e. key characteristics of natural plant communities, caused several deviations from traditional bdts. Most interestingly, the bdts included four sections instead of the traditional three without and under minor water stress (Lonsdale 1990, Begon et al. 1991, Berger and Hildenbrandt 2003. Section I and II delineated the beginning and the end of the self-thinning process of the initial cohort which was indicated by a characteristic rise and decay of positive skewness (Berger et al. 2002). In section I, the bdts showed an immediate onset of mortality which contradicts traditional bdts where the first section indicates fast growth without or little mortality (Lonsdale 1990, Begon et al. 1991, Berger and Hildenbrandt 2003. This inconsistency was directly related to the functional diversity included in our model and independent of density. Specifically, early mortality was mainly observed in annuals and perennials with a PFT-specific low adult growth rate. Section II corresponded to the classical self-thinning section with slopes marginally less negative than the expected −3/2. This is the outcome of an averaging process between the opposing effects of size-symmetric competition and recruitment on the thinning slope. Specifically, size-symmetric competition flattened (Stoll et al. 2002, Chu et al. 2010, Lin et al. 2014 whereas recruitment steepened the values of the slope (Schulze et al. 2005) (Supplementary material Appendix 6). In section III, the bdts deviated from the linear thinning line because individuals of the initial cohort approached their maximum size (Berger et al. 2002). Therefore, mortality due to senescence increased and subsequently, more gaps formed into which new individuals recruited. This new recruitment straightened the bdt compared to even-aged studies in which the bdts instead bent parallel to the x-axis (Kohyama 1992, Zeide 1995, Berger et al. 2002, Monserud et al. 2005 (Supplementary material Appendix 6). The novel section (IV) detected in our study was a result of allowing recruitment, too (Supplementary material Appendix 6). This section represented the density-dependent dynamics of new generations as a sequence of alternating phases of size diversification (increase in density with decreasing in mean biomass) and homogenization (density-dependent mortality). Over time, these phases desynchronized between simulation runs and decreased in amplitude and period despite similar longevities and maximum sizes of the PFTs. This was caused by a shift from the initial even-aged cohort to an uneven-aged community with large individuals and aggregated patches of recruitment. Consequently, density-dependent mortality occurred locally within patches of recruitment resulting in a spatially desynchronized self-thinning process (Moeur 1997). Section IV has not been observed in field studies which relate changes in biomass to density through time. However, the observed changes in the total biomass (community yield) through time (Supplementary material Appendix 6) in section IV resemble the temporal dynamics of the biomass development in uneven-aged multi-species forest landscapes in the so-called transition and quasi-equilibrium phase (Bormann andLikens 1979, Shugart Jr. 2014). Specifically, after a constant initial increase, the community yield reached a maximum at a point in time approximately at the transition between section III and IV. The yield then dropped and settled to a quasiequilibrium. Theory predicts that even-aged monospecific populations thinning along slopes of −3/2 eventually reach a maximum yield which then stays constant over time. At maximum yield, the thinning slope shifts to −1 (Lonsdale andWatkinson 1982, Westoby 1984). Surprisingly, our results indicate that in uneven-aged communities a corresponding shift in the thinning slope to −1.1 occurs at maximum yield even though then their yield is not constant over time but at quasi-equilibrium.
Future studies should examine whether sections I-IV or the shift in slope may be observed in natural plant communities. This could be achieved by analyzing the biomass-density relationship in existing long-term data sets of forests for which the above described yield dynamics has been described (Bormann andLikens 1979, Shugart et al. 2010). Another option would be to add additional information on the age distribution, yield and site conditions to interspecific scaling studies so that biomass-density relationships from communities on similar sites but of different ages could be pieced together ('space for time' substitution).
The respective influence of recruitment and functional diversity on the above-described deviations from traditional bdts without and under minor water stress is evidenced by the simulations without recruitment (Supplementary material Appendix 6) and the temporal dynamics of the skewness of age and biomass distributions (Supplementary material Appendix 7), PFT richness and Shannon diversity index. Specifically, the absence of influence of functional diversity on the bdt sections II-IV is reflected by the independence of the skewness of the biomass distribution from both PFT richness and Shannon diversity index over time. For example, despite pronounced differences in the skewness in section IV, PFT richness was similar across soil water scenarios (Supplementary material Appendix 4). On the other hand, the temporal dynamics of the skewness of the age and biomass distribution correspond in section II-IV but not in section I (Supplementary material Appendix 7). This supports that functional diversity was underlying the deviations from traditional bdts in section I and recruitment in section II-IV.
The fact that we relaxed the even-aged assumption underpinning the self-thinning rule by allowing frequent recruitment steepened the slope of the bdt relatively in both section II and IV. Thus, for a given mean surviving biomass, the potential density would be higher in uneven-aged compared with even-aged communities. This is in line with forestry studies on Reineke's stand density rule (Sterba and Monserud 1993, Gül et al. 2005, Monserud et al. 2005). Reineke's rule is a special case of the self-thinning rule which predicts a slope of −1.605 when the logarithm of mean stem diameter is plotted against the logarithm of tree density (Reineke 1933, Pretzsch 2002). In contrast to the self-thinning rule, Reineke's rule was derived by evaluating forest inventory data across different sites (indirect time series) whereby density was the dependent and not the predictor variable (Pretzsch 2002). A flattened Reineke's maximum density slope in uneven-aged forest stands (Sterba and Monserud 1993, Gül et al. 2005, Monserud et al. 2005 corresponds to the here detected steeper self-thinning slope in unevenaged plant communities. This implies that assuming an even age for management decisions, productivity assessments or carbon budgeting would underestimate density for a given mean stem diameter (Sterba and Monserud 1993) or mean surviving biomass.
Functional diversity neither influenced the classical selfthinning section nor the density-dependent dynamics of new generations. This is in contrast to forestry studies which reported a flattened Reineke's maximum density slope for multi-species forest stands (Sterba and Monserud 1993). Our findings can be explained with the fact that all individuals were competing for the same limiting resources and interand intra-PFT competition were of similar strength (Loreau and de Mazancourt 2008). Thus, density-dependence operated at an aggregated community level as if all PFTs were a single population (Loreau andde Mazancourt 2008, West et al. 2009).

Conclusions
The thinning slope ɑ was historically developed for even-aged monospecific populations (Yoda et al. 1963). Since then, it has been suggested that the self-thinning rule also holds for natural and semi-natural plant communities (Schulze et al. 2005, Vospernik andSterba 2015). The inclusion of a few important characteristics of plant communities into our model, specifically recruitment, functional diversity and abiotic stress, enabled us to derive general principles that broaden the applicability of the self-thinning rule to a community level.
Our results support a wider range of biomass-density relationships than predicted by the self-thinning rule. An interesting and new finding was that the changes in the biomass-density trajectory (bdt) sections including the values of the thinning intercept k and slope ɑ were a function of the stress intensity and irrespective of the stress type. Specifically, increasing stress increased k and decreased ɑ in the classical self-thinning section. Furthermore, the age structure resulting from frequent recruitment decreased ɑ and enabled a novel fourth bdt section. Vice versa, the impact of functional diversity on the bdt was minor. Namely, early mortality was density independent due to the coexistence of different PFTs with various growth rates and longevities. In summary, each relaxation of assumptions, led to a consistent and predictable change in different sections of the bdt. Consequently, to predict the biomass-density relationships in natural communities, the entire bdt should be considered rather than excluding data points or entire sections of the bdt (Liu et al. 2006, Chu et al. 2010, Lin et al. 2014.
Our study highlights that rather than trying to pinpoint fundamentally valid values to the thinning slope, studies should focus on mechanisms and patterns behind densitydependent dynamics. The fact that the values of the thinning intercept and slope, as well as other measures within the entire bdt, changed in a consistent manner across a gradient of water stress and with age structure does neither generalize nor dissipate the rule but instead increases its applicability to conditions beyond even-aged monospecific populations in optimal environments.