Competition contributes to both warm and cool range edges

Study area and species

We selected three sites across an elevation gradient in the western Swiss Alps (Bex, Canton de Vaud), situated at 890, 1400 and 1900 m above sea level (hereafter, the low, middle, and high sites; Supplementary Fig. 1). The three sites span a temperature gradient ranging from 2.5 to 9.6 °C (mean annual temperature from 1981 to 2015⁴⁴; Supplementary Table 1). With increasing elevation, soil moisture increased, and the growing season length was shortened by a longer snow-covered period, as measured from July 2019 to June 2020 (Supplementary Fig. 2). All sites were established on south-facing and shallow slopes in pasture and fenced to exclude livestock.

We included 14 herbaceous focal species that frequently occur in this region, half of which originated from low elevation (hereafter lowland species) and half from high elevation (highland species, Supplementary Table 2). Lowland species had upper range limits (defined as the 90th percentile of their elevation distribution) below 1500 m (with the exception of Plantago lanceolata, with a 90th percentile of 1657 m), while highland species had lower range limits (defined as the 10th percentile of their elevation distribution) above 1500 m, based on a dataset of 550 vegetation plots from the study area⁴⁵. These species consisted of 12 perennial and two biennial species, which are the dominant life histories in this region. Species were selected to include a range of functional types (7 forbs, 4 grasses, 3 legumes) and functional traits (based on plant height, specific leaf area and seed mass). Seeds were obtained from regional suppliers given the large quantities that were needed to establish the experiment (Supplementary Table 2).

Field competition experiment

We designed a field experiment to study the effects of elevation on population growth rates and competitive outcomes by growing focal plants either without competition or competing with a background monoculture of the same or another species (Supplementary Fig. 1). In spring 2017, we established 18 plots (1.6 × 1 m, 0.2 m deep) at each of the three field sites, lined with wire mesh to exclude rodents (except at the high site) and with weed-suppressing fabric on the sides to prevent roots growing in from outside. To control for soil effects, the beds were then filled with a silt loam soil that originated from a nutrient-poor meadow at 1000 m a.s.l. within the study area. Four plots were maintained as bare soil plots (non-competition plots). The other 14 plots received 9 g m⁻² of viable seeds of each species, which allowed the establishment of a monoculture of relatively high density (competition plots). We then periodically weeded the plots to maintain monocultures over the course of the experiment. All species except for two (Arnica montana and Daucus carota) successfully established monocultures, of which 11 species (including six lowland species and five highland species) were fully established by autumn 2017. We then resowed the other plots that failed to establish, which subsequently established either in spring 2018 (Poa trivialis and Poa alpina in the low site and Bromus erectus in the middle site) or autumn 2018 (Aster alpinus, P. trivialis and P. alpina in the middle site and Sesleria caerulea in the low and high sites). Species that failed to establish were included only as focal species for the calculation of invasion population growth rates (i.e. the density was low for A. montana and D. carota in all sites, Trifolium badium in the low site and S. caerulea in the middle site, probably due to high mortality rates caused by drought).

We first raised focal seedlings of each species in a greenhouse for six weeks on standard compost and then transplanted them into the field (Supplementary Fig. 1). To test for responses to elevation in the absence of competition, focal plants were transplanted into non-competition plots at 25 cm apart in autumn 2017 (n = 9 per species and site). To test for effects of competition, we transplanted focal individuals into established plots with 14 cm spacing (n = 9 per focal species, competitor and site). Focal plants that died within two weeks of transplanting were replaced (ca. 5%), assuming mortality was caused by transplant shock. Note that we transplanted focal plants into plots only when the background monocultures were fully established. In 2018 and 2019, we replaced dead focal individuals in spring and autumn (ca. 10% each time). The full design included 56 unique interspecific pairs in each site accounting for 61% of all 14 × 13 = 91 possible pairwise combinations. These pairs were selected to evenly sample differences in functional trait space based on a pilot analysis using plant height, specific leaf area and seed mass obtained from the LEDA dataset⁴⁶. Each focal species competed against four lowland and four highland species, yielding 14 lowland–lowland and highland–highland pairs and 28 lowland–highland pairs. Across all three sites, this design resulted in N = 3780 individuals in total ([56 interspecific pairs × 2 + 14 intraspecific pairs + 14 non-competition] × 9 individuals × 3 sites).

Demographic data

We followed each focal individual between 2017 and 2020 to monitor individual-based demographic performance (i.e. vital rates; Supplementary Fig. 4). Survival was monitored twice a year at the beginning and the end of the growing season. Towards the end of the growing season each year (August–September), we measured all individuals to record plant size, whether they flowered, and to estimate seed production on flowering individuals. To estimate focal plant size, we measured size-related morphological traits on all focal individuals at each census (i.e. the number and/or length of flowering stalks, leaves or ramets, depending on the species) and estimated dry aboveground biomass using regression models fitted using collected plant samples (mean R² = 0.871; Supplementary Data 1; Supplementary Methods). To estimate seed production, we counted the number and measured the size of fruits on reproductive individuals; we then estimated the number of seeds produced by each individual using regression models fitted using intact fruits of each species collected at the early fruiting stage on background plants (mean R² = 0.806; Supplementary Data 2; Supplementary Methods). We conducted a separate experiment to estimate the germination and recruitment of each species in each site (Supplementary Methods).

Population modelling

To estimate population growth rates (λ), we built integral projection models to incorporate multiple vital rates across the life cycle⁴⁷ (see Supplementary Table 3 for model structure and parameters). Separate IPMs were built to estimate intrinsic growth rates using plants growing in the absence of competition (in non-competition plots) and invasion growth rates using plants growing within the background monocultures (in competition plots), under the assumption that monocultures were at equilibrium (see Supplementary Fig. 5 for a test of this assumption) and that focal individuals did not interact with each other but only with the background species. We used plant size (i.e. estimated dry aboveground biomass, log scale) as a continuous state variable and fitted linear models to estimate vital rate parameters by combining multiple-year demographic data over three censuses (i.e. 2017–2018, 2018–2019, and 2019–2020; see Supplementary Methods for consideration of more complex models). We modelled the probability of survival, flowering, and seedling establishment using generalized linear models with a binomial error distribution, modelled growth and seed production using general linear models and described the offspring size distribution using Gaussian probability density functions. We modelled seed germination, seedling establishment and the seedling size distribution as size-independent functions, assuming they are unaffected by maternal size (Supplementary Fig. 4; Supplementary Table 3). For each vital rate of each species, we selected the best-fitted vital rate model by comparing all nested models of the full models using the Akaike information criterion corrected for small samples (AICc), which allowed us to avoid overfitting models and to borrow strength across competitor species and sites in cases where full models were outperformed by reduced models (Supplementary Methods; Supplementary Data 3 and 4).

We calculated population growth rates (λ) as the dominant eigenvalue of the IPMs, which represents the discrete per-capita growth rate (i.e. ({N}_{t+1}=lambda {N}_{t}))⁴⁷. To evaluate the uncertainty around λ, we performed parametric bootstraps for size-dependent vital rates (i.e. survival, growth, flowering, and fecundity). Specifically, we resampled the parameters of each vital rate model using multivariate normal distributions based on their means and covariance matrices⁴⁸. We then fitted all IPMs and estimated λ for each of the 500 bootstrap replicates (Supplementary Data 5).

Estimation of niche differences, relative fitness differences, and coexistence outcomes

We quantified niche and relative fitness differences and predicted coexistence outcomes following the method of Carroll et al.⁴⁹. This method is based on species’ sensitivity to competition defined as the proportional reduction of the population growth rate of a focal species i when invading a population of a competitor species j that is at its single-species equilibrium, and is mathematically equivalent to one minus the response ratio:

where λ_ij denotes the invasion growth rate of focal species i and λ_i is its intrinsic growth rate. The natural logarithm of discrete population growth rates λ estimated from IPMs are equivalent to per-capita growth rate in continuous population growth models⁵⁰, and this transformation makes sensitivities compatible with the coexistence analysis described below. Sensitivity is greater than 0 for antagonistic interactions, with higher values equating to stronger competition, while facilitative interactions lead to negative sensitivities.

For a pair of species, modern coexistence theory predicts that niche differences (ND) promote coexistence by reducing the intensity of interspecific competition experienced by both species. Therefore, a pair of species with a large niche difference should display small mean sensitivities to competition from each other. Consequently, niche differences can be calculated as one minus the geometric mean of the two sensitivities (i.e. niche overlap). In contrast, relative fitness differences (RFD) quantify the degree of asymmetry in species’ competitive abilities. Therefore, a pair of species with a large fitness difference should display large differences in their sensitivities to competition from each other, as quantified as the geometric standard deviation of sensitivities⁴⁹:

There are three possible outcomes of competition between a given pair of species: stable coexistence, a priority effect, and competitive exclusion. These can be quantified based on either invasion criteria or the relative magnitude of niche differences versus relative fitness differences^15,51. Stable coexistence is only possible when both species are able to invade each other’s equilibrium populations; this condition is met when ND > 0 and ({{{{{rm{RFD}}}}}} , < , frac{1}{1-{{{{{rm{ND}}}}}}})⁴⁹, which is equivalent to (frac{1}{{{{{{rm{RFD}}}}}}(1-{{{{{rm{ND}}}}}})} > 1), with greater values indicating more stable coexistence and providing a metric for the strength of coexistence (i.e., coexistence metric²⁶). When neither species can invade when rare, then priority effects occur, meaning that whichever species is initially established within a community has an advantage and excludes the other. This could happen when a species pair has a small niche difference and a small relative fitness difference, that is ND < 0 and ({{{{{rm{RFD}}}}}} , < , frac{1}{1-{{{{{rm{ND}}}}}}})⁵¹. Otherwise, if only one species can invade, then it is predicted to competitively exclude the inferior competitor; this occurs when relative fitness differences are big enough to override the stabilizing effect of niche differences, that is when ({{{{{rm{RFD}}}}}} , > , frac{1}{1-{{{{{rm{ND}}}}}}}). We quantified competitive outcomes and coexistence metrics for each of the 500 bootstrap replicates of the dataset (Supplementary Data 6).

Note that we excluded facilitative interactions that were present in 13% of all pairs because the equations for niche differences and relative fitness differences are not compatible with negative values of sensitivity (Eq. 2 and 3); we did not exclude facilitative interactions for other analyses. We quantified the coexistence determinants of species pairs in cases where either one or both of the species were predicted to be unable to persist in the absence of neighbours (i.e. ln(λ_intrinsic) < 0; 8% of all pairs). This was because the quantification of niche and fitness differences were based on the competitive impacts on population growth (i.e. sensitivity to competition) and independent of the actual estimates of population growth rates.

Statistical analysis

We tested the effect of elevation on species’ intrinsic and invasion growth rates, and whether these effects differed between lowland and highland species, by fitting mixed-effects models with elevation (continuous), species origin (i.e. lowland or highland) and their interaction as fixed effects and focal species identity as a random effect. To test whether coexistence was strengthened at high versus low elevation, we fitted a mixed-effects model of the coexistence metric as explained by elevation and competitor identity (i.e. the three types of pairs: lowland–lowland, lowland–highland, and highland–highland pairs), with species pair as a random effect. We fitted similar models for niche differences and relative fitness differences (absolute values) to gain insights into whether changes in coexistence were driven predominantly by variation in niche or relative fitness differences. Additionally, we tested the effect of elevation on relative fitness differences for each type of pair separately and used the actual values of lowland–highland pairs since we were interested not only in the magnitude of relative fitness differences but also in the direction of effects (i.e. whether lowland or highland species were dominant competitors). The response variables in these models were log-transformed to meet model assumptions. Note that we therefore had to use niche overlap (1-ND) for the log-transformation instead of niche differences, because niche differences included negative values (Eq. 2). Firstly, we performed these analyses based on the mean of the 500 bootstrap replicates for a given focal species/background competitor combination and used likelihood ratio tests to determine the significance of fixed effects. To account for the uncertainty around estimated λ, we also performed individual tests for each bootstrap replicate of the dataset and determined the significance of fixed terms based on whether the 95% confidence interval of estimated coefficients, or contrasts of interest, included zero. All population modelling and statistical analysis were conducted in R version 4.0.3⁵².

Reporting summary

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

Source: Ecology - nature.com