When history matters: The overlooked role of priority effects in grassland overyielding

This 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. © 2019 The Authors. Functional Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society 1Ecosystem Functioning and Services, Institute of Ecology, Leuphana University, Lüneburg, Germany 2Plant Sciences, Institute for Bio and Geosciences, IBG‐2, Forschungszentrum Jülich GmbH, Jülich, Germany


| INTRODUC TI ON
Long-term biodiversity-ecosystem functioning (BEF) experiments have shown that communities with a greater plant species or functional group richness are often more productive above-ground (Hector et al., 1999;Marquard et al., 2009;Tilman et al., 1997) and below-ground (Oram et al., 2018;Ravenek et al., 2014). Several mechanisms such as multitrophic interactions, resource partitioning and abiotic facilitation have been proposed to explain these positive biodiversity-productivity relationships, but their relative contributions to grassland overyielding remain unclear Eisenhauer, 2012;Weisser et al., 2017).
Over the years, the use of statistical methods developed to partition the net effect of biodiversity on ecosystem functioning into two (Loreau & Hector, 2001) or three (Fox, 2005) additive components has allowed researchers to quantify the contribution of niche differences and/or interspecific interactions (complementarity effect) as well as dominance of highly productive species (dominance/selection effect) to the increased functioning of diverse plant communities. Although these additive partitioning methods do not allow a direct identification of the biological processes driving grassland overyielding , they largely contributed to a better understanding of the mechanisms behind the patterns observed in BEF experiments (Cadotte, 2017;Cardinale et al., 2007;Fox, 2005;Loreau & Hector, 2001Marquard et al., 2009;Oram et al., 2018;Roscher et al., 2005).
Plant species and functional group richness, however, are not the only drivers of ecosystem functioning in natural habitats. Both the order and timing of species arrival during community assembly can also have long-lasting impacts on community structure and functioning Körner, Stöcklin, Reuther-Thiébaud, & Pelaez-Riedl, 2008;Švamberková, Doležal, & Lepš, 2019;Weidlich et al., 2017Weidlich et al., , 2018Wilsey, Barber, & Martin, 2015), as well as on the shape of the relationship between biodiversity and productivity (Fukami & Morin, 2003). This phenomenon is referred to as a priority effect and is a biotic component of historical contingency (Fukami, 2015;Grainger, Letten, Gilbert, & Fukami, 2019;Ke & Letten, 2018).
Priority effects occur when early arrival of species affects the establishment, growth or reproduction of species arriving later (Eriksson & Eriksson, 1998) and can thus lead to alternative states in vegetation (Fukami, 2010(Fukami, , 2015Fukami & Nakajima, 2011).
Despite the importance of priority effects for community assembly, we lack an understanding of their importance in influencing the direction and magnitude of the relationship between biodiversity and ecosystem functioning. At a given level of plant species and functional richness, however, it is probable that different sequences of plant species arrival could affect grassland overyielding via its effects on complementarity and dominance effects. For instance, an early arrival of N 2 -fixing species (legumes) in the community could favour the establishment of late-arriving species such as grasses and non-N 2 -fixing forbs via nitrogen facilitation mechanisms (Temperton, Mwangi, Scherer-Lorenzen, Schmid, & Buchmann, 2007), thus leading to the creation of positive priority effects. This would be in line with the greater net biodiversity and complementarity effect values usually observed in grassland communities containing legumes (Loreau & Hector, 2001;Marquard et al., 2009). Because larger complementarity effect values are expected when species facilitate one another (Fox, 2005;Loreau & Hector, 2001), positive priority effects would then be associated with greater complementarity effect values. An early arrival of species performing well in monoculture plots (e.g. grasses), however, could lead to negative priority effects and larger dominance effect values because early-arriving species might dominate mixtures at the expense of species arriving later during assembly. Therefore, we hypothesize that different sequences of plant species arrival during community assembly would lead to the creation of priority effects affecting the magnitude of net biodiversity effects as well as the relative contributions of complementarity and dominance effects to grassland overyielding.
To test this hypothesis, we used species-specific plant biomass data collected in 2013 and 2014 in a subset of the plots of the Jülich Priority Effect experiment located in Germany (Weidlich et al., 2017(Weidlich et al., , 2018. In this field experiment, the order of arrival of three plant functional groups (PFG: legumes, grasses and non-N 2 -fixing forbs) was manipulated to investigate how priority effects affect plant community structure and ecosystem functioning in temperate grasslands. Each PFG arrived either 6 weeks earlier (legumes, grasses or forbs sown first) or at the same time (synchronous) as the other PFGs. For each experimental plot, the net biodiversity effect was quantified as in Loreau and Hector (2001) and was partitioned using the method of Fox (2005) into three additive components: trait-independent complementarity effect, trait-dependent complementarity effect and dominance effect (Table 1). These data were analysed using a two-step approach. First, we investigated whether PFG order of arrival affected overyielding as well as its drivers in assembling grassland communities. Second, we investigated whether the magnitude of net biodiversity, complementarity and dominance effects was dependent on the strength and direction of priority effects. This study provides strong evidence that manipulating plant order of arrival during community assembly can lead to the creation of priority effects of various strengths and directions affecting the magnitude of net biodiversity effects in grassland communities.

| The Jülich Priority Effect experiment
We used above-ground plant biomass data collected at the spe- Detailed meteorological data measured in Jülich from 2012 to 2014 are provided in Figure S1. A detailed description of the experiment can be found in Weidlich et al. (2017). Briefly, this experiment was set up in 2012 using a full factorial and randomized design to study how PFG order of arrival and sown species richness affect ecosystem functioning and plant community structure in temperate grasslands (Weidlich et al., 2017(Weidlich et al., , 2018. Priority effects in community assembly | 3 Functional Ecology DELORY Et aL. TA B L E 1 Interpretation of the terms of the two-way additive partitioning method of Loreau and Hector (2001) and the tripartite additive partitioning method of Fox (2005)

Negative values
Net biodiversity effect (NBE) NBE The observed yield in mixture is greater than the weighted average of the monoculture yields (overyielding) The observed yield in mixture is equal to the weighted average of the monoculture yields (no effect) The observed yield in mixture is lower than the weighted average of the monoculture yields (underyielding)

Complementarity effect
Trait-independent complementarity  were created by manipulating the order of arrival of three PFGs: N 2 -fixing forbs (legumes), non-N 2 -fixing forbs (forbs) and grasses. In synchronous communities, all plant species were sown at the same time during the first sowing event. In forbs-first (F-first), grassesfirst (G-first) and legumes-first (L-first) communities, however, all the species from one PFG were sown 6 weeks before the others. Even though this experiment was initially set up using two sown species richness levels (9 or 21 species), we only used data collected from plots in which nine species (three species per functional group) were sown, as monoculture plots were only available for the species present in the 9 species mixtures (see Table S1). In addition, because the plots were set up on two different areas characterized by two different soil types, we only used plant biomass data collected from the plots located on the same area as the monocultures to ensure comparability. In total, data collected from 16 mixture plots (n = 4 for each PFG order of arrival treatment; surface: 16 m 2 /plot) and 18 monoculture plots (n = 2 for each species; surface: 4 m 2 /plot) were used for this study. Both monoculture and mixture plots were established at the same time when the experiment was set up in 2012.  Table 1 summarizes information from the papers of Fox (2005) and Loreau and Hector (2001) to calculate and help interpreting the terms of their additive partitioning methods. All calculations were performed using the apm function of the bef r package developed for the purpose of this study. This r package is available on GitHub (https ://github.com/Benja minDe lory/bef).

| Quantification of priority effects
In our field experiment, we created priority effects by sowing a group of N − p (N minus p) species 6 weeks after a group of p earlyarriving species (N is the total number of species sown in the plots).
The cost of arriving late during plant community assembly (P, priority effect) for the N − p late-arriving species was calculated using Oi and Y Sync Oi are the observed yields of species i when it arrived later or at the same time as the early-arriving species, respectively. This priority effect index has the same mathematical properties as the additive neighbour-effect intensity index developed by Díaz-Sierra, Verwijmeren, Rietkerk, Dios, and Baudena (2017): it is standardized, symmetric (additive symmetry) and bounded between −1 (competitive exclusion of late-arriving species) and +2 (obligate facilitation of late-arriving species). The direction and strength of the priority effect are given by the sign and absolute value of P, respectively (Figure 1). As we had four replicates for the Synchronous treatment, we calculated four values of P for each F-first, G-first and L-first plot. In Figure 3 and Figure S4, we reported the mean value of P (n = 4) calculated for each plot with a priority effect treatment as well as its 95% confidence interval computed by bootstrapping using the percentile method (10,000 iterations).

| Statistical analyses
We used two-way ANOVA models to determine whether the PFG order of arrival in assembling communities, the sampling year (1 or 2 years post-seeding) or their interaction affected overyielding and its three additive components (DE, TICE and TDCE). ANOVA assumptions were systematically checked by looking for any pattern in a plot showing the values fitted by the linear model against model residuals. Detailed ANOVA tables are provided in Table S2.
Pairwise comparisons using Tukey contrasts were performed on estimated marginal means computed with the emmeans (Lenth, 2018) r package. p-values and 95% confidence intervals were adjusted for multiple comparisons using Tukey adjustment. Sidak's adjustment method was used when Tukey's method was not appropriate. We tested whether the estimated marginal means were significantly different from zero by examining whether their 95% confidence interval contained zero or not. If the confidence interval did not contain zero, we considered that the estimated marginal mean significantly deviated from zero. (1) When a significant linear relationship was found between two variables (p < .05), standardized major axis regression models were fitted using the smatr r package (Warton, Duursma, Falster, & Taskinen, 2012). Correlation coefficients and model parameters with their 95% confidence interval are provided in Table S3.
Results for TICE and TDCE are provided as supplementary material ( Figure S4).
Statistical analyses were performed in R 3.5.2 (R Core Team, 2018) with an alpha value of .05.

| Plant order of arrival during community assembly affects complementarity and dominance effects
One forbs-first (F-first) and grasses-first (G-first) communities were simultaneously driven by positive dominance and trait-independent complementarity effects (Figure 2b,c), overyielding in plots where legumes were sown first (L-first) was mainly caused by a dominance effect ( Figure 2b). Contrary to our expectations, the trait-independent complementarity effect was the lowest (in fact, not significantly different from zero) when legumes were the first to arrive in the community (Figure 2c). We did not find any significant difference in trait-dependent complementarity effect between PFG order of arrival treatments. As shown in Figure 2d, this effect was either close to (G-first) or not significantly different from zero (Sync, F-first and L-first). The strong dominance effect observed in communities where legumes were sown first was due to the fact that one legume species (Trifolium pratense) with higher-than-average monoculture yield dominated the mixtures at the expense of all the other species sown in the plots, except for Medicago sativa. In these plots, the yield achieved by T. pratense was on average 23% lower than its yield in monoculture (see Figure S2). In plots where grasses and forbs were sown first, all the species that arrived first performed better or as good as what would be expected under the null hypothesis (i.e. for each species i, its yield in mixture Y Oi equals its yield in monoculture M i divided by sown species richness), despite the fact that T. pratense F I G U R E 1 Framework for the quantification of priority effects. In our experiment, we created priority effects by manipulating plant order of arrival. To do so, one plant functional group was sown 6 weeks before the two others. A scenario without any order of arrival manipulation (synchronous) was also included in the experimental design. Y Late is the total yield of late-arriving species in the treatment with order of arrival manipulation. Y Sync is the total yield of the late-category species in the synchronous treatment. For each order of arrival scenario, both Y Late and Y Sync are calculated using the same species pool. Y Sync was fixed at 1 unit. The priority effect index (P) shares the same mathematical properties as the additive intensity index proposed by Díaz-Sierra et al. (2017): it is standardized, symmetric around zero, and is bounded between −1 and +2 was also dominating the plots (see Figure S2). In plots where all PFG were sown at the same time, however, at least one species of each PFG performed better or as good as what would be expected based on monoculture yields (see Figure S2).
In 2014 Figure   S3). These two species performed remarkably well in legumesfirst plots, with D. glomerata even having a yield in mixture not significantly different from that obtained in monoculture (see Figure S3). Except for plots where legumes were sown first, the biomass of T. pratense measured in 2014 was lower than what would be expected by the null model. The species with the greatest monoculture yield (Festuca pratensis) had a low productivity in all experimental plots, except in those where grasses were sown first (see Figure S3), thus explaining the positive dominance effect values measured in grasses-first plots (Figure 2b). In plots where they were the first to arrive, all forb species performed at least as well as what would be predicted by the null model (see Figure S3).

| Moving from negative to positive priority effects increases grassland overyielding via increased complementarity effects
Overyielding in our grassland experiment was positively correlated to the priority effect index (Figure 3a,b, Table S3). Interestingly, this increase in net biodiversity effect observed when moving from negative to positive priority effects was solely due to an increase in complementarity effects (Figure 3c,d). Both trait-independent and trait-dependent complementarity effects were positively correlated to P (see Figure S4 and Table S3), but no relationship was found between dominance effect and P (Figure 3e,f). The same pattern was F I G U R E 2 PFG order of arrival alters overyielding drivers in the Jülich Priority Effect experiment. The tripartite method of Fox (2005) was used for the partitioning. For each sampling year and each PFG order of arrival treatment, the panels show the net biodiversity effect (a) and its three additive components: dominance effect (b), trait-independent complementarity effect (c) and trait-dependent complementarity effect (d). Values are estimated marginal means ± 95% confidence intervals (n = 4). Individual data points are displayed as grey dots on the left side of each group. For each sampling year, PFG order of arrival treatments that do not share a common letter are significantly different from each other (p < .05). Mean values that are significantly different from zero are shown with a filled dot (p < .05). Means values are otherwise shown with an empty dot (p > .05). Detailed ANOVA tables are available in Table S2 found in both sampling years despite the fact that the priority effect index values calculated for each experimental plot were very different in 2013 and 2014, particularly for plots where legumes or grasses were sown first ( Figure 3).
In legumes-first plots, strong negative priority effects (P close to −1) were measured in 2013 but, one year later, 75% of these plots were characterized by strong positive priority effects (P close to +1) and had the greatest overyielding values. These results strongly suggest that sowing legumes first can lead to the creation of positive priority effects.
In plots where grasses were sown first, however, 25% of the plots were characterized by positive priority effects (P close to 0.5) in 2013, with only 50% of the plots having negative priority effect values (P close to −0.5) on the same year. One year later (2014), all grasses-first plots were characterized by strong negative priority effects (P close to −1).

| D ISCUSS I ON
Linking BEF research with the field of community assembly is one of the next important steps in ecology, since natural communities experience assembly and are not weeded as BEF experiments are (Bannar-Martin et al., 2018). In order to do this, we believe that assembly processes that are important for the structure and functioning of plant communities, such as historical contingency in plant species order of arrival, have to be considered alongside plant species and functional group richness. This study enabled an important step in this direction within a priority effect experiment that includes natural assembly as well as monocultures. We show that plant order of arrival can affect overyielding drivers, namely complementarity and dominance effects, in the first years of assembly of a temperate grassland. We also provide evidence that the magnitude of complementarity and net biodiversity effects is dependent on the strength and direction of priority effects. More specifically, we showed that the greatest overyielding values were achieved in plots characterized by positive priority effects, and that the main reason for this was increased complementarity effects.
Overall, the net effect of biodiversity on above-ground productivity in the Jülich Priority Effect experiment markedly decreased from 2013 to 2014. This drop in overyielding was paralleled by a F I G U R E 3 Moving from negative to positive priority effects increases complementarity and overyielding in a temperate grassland. The relationship between the priority effect index and overyielding (a, b), complementarity (c, d) or dominance (e, f) is shown separately for each sampling year (left panels, 2013; right panels, 2014). When two variables were significantly correlated (p < .05), the regression line (solid line) is shown. Values of P are shown as mean ± 95% confidence interval (n = 4). The confidence intervals were computed by bootstrapping using the percentile method (10,000 iterations). Pearson's product-moment correlation coefficients (r) and regression parameters (slope and intercept) can be found in Table  S3. The symbol used for each individual observation refers to the PFG order of arrival treatment (see legend in panel a) decrease in complementarity (both trait-independent and trait-dependent) and dominance effects. Even though grassland overyielding usually tends to strengthen over time because of increased species complementarity (Cardinale et al., 2007), year of the experiment in plots where grasses were sown first. In these plots, grasses competitively excluded forbs and legumes (the total above-ground biomass production of forbs and legumes was on average 94% lower in grasses-first plots than in synchronous plots) (see Figure S3). The reasons behind this strong grass dominance are still unclear, but knowledge gained in other ecosystems, such as Mediterranean grasslands, can help to identify possible environmental drivers favouring grass-dominated transient states. In a field experiment testing the importance of year and site effects on the structure of plant communities in Californian grasslands, Stuble, Fick, and Young (2017) also found good 'grass' years leading to different vegetation states and identified both mean annual temperature and total number of rainy days as likely drivers of plant community dissimilarities. In a different study, Clary (2008) found that the relative abundance of annual and perennial grass species in Mediterranean grasslands was mainly determined by rainfall seasonality, with low summer precipitation levels favouring annual grass dominance. Although we cannot confirm it, the overall warmer and drier conditions during the third growing season of our experiment might have driven the convergence of plant communities towards a grass-dominated state. Our study makes abundantly clear that more research is now needed to improve the predictability of transient community dynamics with regard to weather conditions during plant establishment and plant order of arrival during assembly (Fukami & Nakajima, 2011;Temperton, Baasch, Gillhaussen, & Kirmer, 2016).
Contrary to our expectations, sowing leguminous species before the other PFGs led to a strong dominance of T. pratense and  Wright, Wardle, Callaway, & Gaxiola, 2017). Following the framework proposed by Barry et al. (2018), we present three non-mutually exclusive mechanisms that could explain increased complementarity between early and late-arriving species: below-ground resource partitioning, biotic feedbacks (both negative and positive) and abiotic facilitation (physical stress buffering).
Below-ground resource partitioning in space and/or time has been one of the most prevalent hypotheses to explain the positive biodiversity-productivity relationships found in grassland ecosystems Loreau & Hector, 2001). However, results from experiments that manipulated plant species richness without manipulating plant order of arrival often did not support this hypothesis (Jesch et al., 2018;Mommer et al., 2010;Oram et al., 2018;Ravenek et al., 2014), thus suggesting that mechanisms other than resource partitioning drive above-ground and below-ground grassland overyielding. If different species arrive at different times during plant community assembly, however, one can expect belowground niche partitioning to occur as a consequence of soil resource pre-emption and/or niche modification (sensu Fukami, 2015) by early-arriving species. In the Jülich Priority Effect experiment, we found that the root length density in the topsoil layer depended on the sequence of arrival of PFGs. In 2014, the standing root length density was indeed lower in plots where legumes were sown first in comparison with synchronous and grasses-first plots (Weidlich et al., 2018). Whether this pattern was due to changes in vertical root distribution, total root productivity or both is still unknown, but it is a strong indication that plant order of arrival during assembly can have important consequences for ecosystem functioning, particularly below-ground. Future research using imaging techniques to non-destructively follow root development over time (e.g. minirhizotrons) (Rewald & Ephrath, 2013) as well as molecular or spectral techniques to disentangle the relative contribution of individual plant species to biomass production in different soil layers (Meinen & Rauber, 2015;Mommer, Wagemaker, De Kroon, & Ouborg, 2008) hold much potential to investigate how the sequence of arrival of different species or functional groups affect vertical root distribution and below-ground productivity in temperate grasslands.
Because species arriving first during plant community assembly can alter the biotic and abiotic soil conditions that will be experienced by species arriving later (Baxendale, Orwin, Poly, Pommier, & Bardgett, 2014;Bezemer et al., 2006;Hu et al., 2018), historical contingency effects can arise as a consequence of plant-soil feedbacks, thus affecting plant community assembly (Kardol, Cornips, Kempen, Bakx-Schotman, & Putten, 2007;van der Putten et al., 2013). When the relative abundance of a species is high, as is typically the case in monoculture plots, the accumulation of species-specific pathogens such as bacteria, fungi or nematodes in the rhizosphere can result in negative biotic feedbacks leading to negative frequency dependence (Guerrero-Ramírez, Reich, Wagg, Ciobanu, & Eisenhauer, 2019;Hendriks et al., 2013;Mommer et al., 2018). According to modern coexistence theory (Chesson, 2000;Fukami, Mordecai, & Ostling, 2016), such negative feedbacks act as a stabilizing mechanism allowing species coexistence. They are also thought of as one of the primary mechanisms (alongside N facilitation) behind the increased ecosystem functioning (overyielding) observed in speciesrich grassland communities compared to monocultures Mommer et al., 2018). Although we did not verify this in our experiment, the build-up of species-specific pathogens in monocul- Positive biotic feedbacks can result from the accumulation of symbiotic mutualists in the rhizosphere such as N 2 -fixing rhizobia and mycorrhizal fungi (Eisenhauer, 2012;Semchenko et al., 2018;van der Putten et al., 2013;Wright et al., 2017). Because these mutualists are able to increase the amount of resources that can be taken up by plants ( Barry et al., 2018;Wright et al., 2017), their accumulation in the rhizosphere of early-arriving species could lead to increased establishment of species arriving later during plant community assembly. Although increased soil N availability can be expected as a consequence of legume presence in plant communities (Temperton et al., 2007), we do not think that it played an important role in our experiment for at least two reasons: (a) there was no effect of time and PFG order of arrival on the soil N content (Weidlich et al., 2017) and (b) we did not find evidence for N transfer using N content and δ 15 N natural abundance data (see Figure S5). Positive biotic feedbacks favouring the establishment of late-arriving species could also arise via the accumulation of plant growth-promoting rhizobacteria in the rhizosphere of early-arriving species. Although the link between nonresource mutualists and increased ecosystem functioning has been far less studied in comparison with rhizobia and mycorrhizal fungi , studies have shown that plant growth-promoting rhizobacteria positively impact on plant performance by inhibiting soil-borne pathogens, particularly in species-rich plant communities (Eisenhauer, 2012). In addition, some rhizobacterial strains are known to modulate root development and root system architecture as well as promoting root and shoot growth (Delaplace et al., 2015;Verbon & Liberman, 2016). Despite evidence showing that plant species richness drives the structure and activity of the root-associated microbiota, notably via root biomass and root exudate-dependent mechanistic pathways Lange et al., 2015;Steinauer, Chatzinotas, & Eisenhauer, 2016), the role of historical contingency related to plant order of arrival is currently unclear despite its obvious relevance for the cyclical feedback loops between plants, microbes and soils. In order to gain a better understanding behind the mechanisms driving plant-plant interactions in naturally assembling communities (such as complementarity), we argue that a better understanding of how the sequence of arrival of different plant species or functional groups affects root exudation patterns as well as the soil and plant-associated microbiota is much needed.
Next to resource partitioning and biotic feedbacks, abiotic facilitation via physical stress buffering (or microclimate amelioration) is a mechanism that could have also contributed to increase species complementarity and grassland overyielding in plots characterized by a positive priority effect . When some species arrive earlier than others, they modify the local environment by providing shade, thus reducing temperature and evapotranspiration as well as increasing air relative humidity and soil water content (Bruno, Stachowicz, & Bertness, 2003;Wright et al., 2017). These Because our experimental design did not include monocultures for 12 of the species used in the plots sown with 21 species (high diversity plots in Weidlich et al., 2017), we were not able to measure net biodiversity effects for these plots. Therefore, our results are based on overyielding values measured at one species richness level only (with 9 sown plant species). Whether the findings presented in this study hold true across a species richness gradient still needs to be investigated. Because our results can have strong implications for restoration settings where the sequence of arrival of different plant species or functional groups can be manipulated to create priority effects that alter community trajectories (Temperton et al., 2016;Wilsey et al., 2015;Young, Stuble, Balachowski, & Werner, 2017), we believe that future research combining both assembly and biodiversity approaches are needed. BEF findings can be applied in the 'real world' either by comparing ecosystems undergoing natural assembly or by including an element of intervention such as is commonly done in ecological restoration where species are added to a system (but usually at the same time) (Jochum et al., 2019;Manning et al., 2019). Further research combining our PFG order of arrival and BEF approach seems very promising, although challenging to design. This is due to the high number of possible treatment combinations associated with the creation of two orthogonal gradients (plant order of arrival × species richness). In addition, these experiments should be designed in such a way that species-specific responses can be separated from functional group responses (Weisser et al., 2017), which was not the case in our study. Nevertheless, we believe that such experiments are now needed to improve our understanding of the functioning of grassland ecosystems and increase the predictive power of community ecology in conservation and ecological restoration.

ACK N OWLED G EM ENTS
We thank Marlene Mueller, Edelgard Schoelgens, Agnes Höltkemeier, and all the students and friends who helped to collect data from the Jülich Priority Effect experiment. We are grateful to Axel Knaps (Jülich Forschungszentrum) for providing the meteorological data presented in Figure S1. We also thank Andreas

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

AUTH O R S ' CO NTR I B UTI O N S
E.W.A., P.v.G. and V.M. performed the Jülich Priority Effect experiment and collected the data; B.M. analysed the data, wrote the bef r package and led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.