Drivers of large-scale spatial demographic variation in a perennial plant

To understand how the environment drives spatial variation in population dynamics, we need to assess the effects of a large number of potential drivers on the vital rates (survival, growth and reproduction), and explore these relationships over large geographical areas and long environmental gradients. In this study, we examined the effects of a broad variety of abiotic and biotic environmental factors, including intraspecific density, on the demography of the forest understory herb Actaea spicata between 2017 and 2019 at 40 sites across Sweden, including the northern range margin of its distribution. We assessed the effect of potential environmental drivers on vital rates using generalized linear mixed models (GLMMs), and then quantified the impact of each important driver on population growth rate (λ) using integral projection models (IPMs). Population dynamics of A. spicata were mostly driven by environmental factors affecting survival and growth, such as air humidity, soil depth and forest tree species composition, and thus those drivers jointly determined the realized niche of the species. Soil pH had a strong effect on the flowering probability, while the effect on population growth rate was relatively small. In addition to identifying specific drivers for A. spicata’s population dynamics, our study illustrates the impact that spatial variation in environmental conditions can have on λ. Assessing the effects of a broad range of potential drivers, as done in this study, is important not only to quantify the relative importance of different drivers for population dynamics but also to understand species distributions and abundance patterns.


44
Abiotic and biotic environmental conditions such as climate and competition vary across time 45 and space and drive the population dynamics of species through effects on vital rates such as 46 survival, growth and reproduction (Bruna andOli 2005, Doak andMorris 2010). These vital 47 rates vary both with regards to how sensitive they are to different environmental factors and to 48 how much they influence population dynamics (Silvertown et al. 1993, Pfister 1998, Nicolè et 49 al. 2011, Ehrlén and Morris 2015. Yet, underlying demographic mechanisms are not 50 accounted for in standard models of species distributions and abundances (Guisan and Thuiller 51 2005, Araújo andRahbek 2006, Elith andLeathwick 2009; but see Merow et al. 2017). 52 Moreover, species distribution models based on occurrence patterns are typically based on the 53 assumption that species are in equilibrium with their environment, and thus do not capture 54 ongoing or recent changes in environmental drivers. To achieve an in-depth understanding of 55 the population dynamics of a species, and to describe and predict distributions and abundances, 56 we need to assess how environmental variation is linked to demographic variation (Ehrlén and 57 Morris 2015). Using environmentally explicit demographic models allows us to disentangle 58 which vital rates are most affected by the environment, and how those vital rates contribute to 59 population growth rate (λ). Such a more mechanistic understanding of the processes underlying 60 a species' distribution and abundance allows for more accurate predictions of species responses 61 to environmental and climate changes (Doak andMorris 2010, Csergő et al. 2017). 62 Environmentally explicit demographic models also provide a means to identify the drivers of 63 the short-term dynamics and of the realized niches of species. 64 To avoid preconceived notions about which drivers are important, assessments of 65 environmental effects on population dynamics should include a broad range of potential 66 drivers, ideally including climatic and other abiotic factors, as well as biotic and anthropogenic 67 survival and growth would have a larger effect on population growth rate than those affecting 93 reproduction, because such dynamics are described for other long-lived herbs (Silvertown et 94 al. 1993). We used generalized linear mixed models (GLMM) to link the environmental drivers 95 to vital rates describing survival, individual growth, flowering probability and fruit number. 96 We then incorporated the GLMMs into an integral projection model (IPM,Easterling et al. 97 2000) to assess the effects of environmental drivers on the population growth rate of A. spicata.

100
Baneberry (Actaea spicata L., Ranunculaceae) is distributed over most of Europe and parts of 101 Asia and North America (Anderberg A. and Anderberg A.-L. 2017). In northern Europe, it 102 occurs in shady, well-drained and nutrient rich forests, often on limestone (Pellmyr 1984). The 103 plant is common throughout Sweden but does not occur in the far north (Mossberg and 104 Stenberg 2014). Actaea spicata has a morphology typical of many early summer flowering 105 forest herbs, with a greater height than spring ephemerals and an umbrella-like leaf display (cf. 106 Givnish 1987). Previous studies have found that individuals can produce several shoots, each 107 typically bearing up to four inflorescences (Eriksson 1995), with the first, top inflorescence, 108 typically bearing most flowers and fruits. Each black berry contains 8-16 seeds (Zeipel et  We collected a range of environmental variables including 15 abiotic factors related to climate, 151 topography and soil properties, and three biotic factors concerning plant community structure 152 and intraspecific population density. We collected weather data throughout the whole study 153 period, with air-temperature and humidity data loggers (EasyLog EL-USB-2, Lascar 154 electronics) placed at the center of each plot and covered with a plastic cup to protect them 155 from direct sunlight and rain. We used the logger data to calculate (i) growing degree days 156 (base temperature: 5 °C) for the spring period 15th May-15th June.06, (ii) the start of spring, 157 using the definition from the Swedish meteorological office (SMHI 2011) and (iii) vapor 158 pressure deficit (hereafter VPD), the difference between the actual vapor pressure in the air and 159 that of saturated air, where vapor pressure is a temperature corrected measure for moisture 160 (Anderson 1936. VPD was calculated for the summer period (15th June -15th 161 August) using the R package plantecophys . We measured slope inclination 162 and slope aspect and used them to calculate the extent to which the ground was slanted towards 163 the sun, as an additional indicator of the exposure to light and moisture (hereafter "slope"). In  as soil potassium has previously been shown to have important effects on A. spicata (Dahlgren 190 and . Final models thus included the following environmental variables in 191 the vital rate regression models: VPD, soil potassium concentration, soil pH, soil depth,

192
percentage of coniferous trees, intraspecific density, canopy gap fraction and slope.

193
We assessed effects of environmental variables on vital rates (probability of survival, growth, including site and plant ID as random effects. We also analyzed the two years separately, to 205 investigate the consistency of effects. In all models, we included plant size as a fixed effect. 206 We included all environmental variables and intraspecific population density as both linear and 207 quadratic fixed-effect terms. To ease interpretation and avoid effects of overfitting (Harrell 208 2001) we then created a reduced model, which only included those quadratic terms that where 209 statistically significant following strict cut-off points (for the GLMMs: P < 0.05; for the LMMs: 210 no P-values are provided by the lme4 package, we therefore used -1.96 < t < 1.96 which 211 corresponds to a P-value of < 0.05 in models with high degrees of freedom). We used this 212 reduced model to assess the statistical significance of remaining linear and quadratic terms, and 213 we evaluated the statistical significant of effects according to the same criteria as above. We 214 conducted all modelling, data exploration and visualization using R version 3.5.1 (R Core 215 Team 2018).

217
An IPM is a population model with similar properties as a matrix population model but with 218 continuous instead of discrete stage classes and a projection kernel instead of a projection 219 matrix (Easterling et al. 2000). We constructed the IPM in this study as described in Dahlgren and fn the number of fruits. fs is the number of seeds per fruit, which was assumed to be 9.61, 236 as in

245
Over the two observed transitions the average yearly survival of individuals was 91%. The Individual plant size was correlated with all four investigated vital rates; survival, growth, 251 probability of flowering, and number of fruits (Appendix S1: Fig. S4 and table S1). Soil depth, 252 soil pH, VPD, intraspecific density and percentage of coniferous trees were significantly 253 related to one or several vital rates (Fig. 2). Most environment-vital rate relationships were 254 linear, but density and percentage of coniferous trees had significant quadratic relationships 255 with fruit number (Fig. 2d). In our regression models, the percentage of coniferous trees affected survival positively (Fig. 2a) whereas fruit number was highest at average values of 257 percentage of coniferous trees (Fig. 2d). VPD negatively influenced growth and flowering 258 probability ( Fig. 2b and c). Soil depth had a positive relationship on growth (Fig. 2b). Soil pH 259 was negatively related to flowering probability (Fig. 2c) and the effect of intraspecific density 260 on number of fruits was u-shaped with higher fruit production at lowest and highest density 261 values (Fig. 2d). The effects of environmental variables were similar in both transitions 262 (Appendix S1: table S4). rates also varied between sites (Fig. 1) and ranged from 0.84 to 1.05 over the study period 277 (sd = 0.05, Appendix S1: table S2). We found no significant relationship between population 278 growth rate and latitude (P = 0.77).

279
Our IPMs showed that coniferous trees and VPD had the largest influence on population growth 280 rate, closely followed by the effect of soil depth (Fig. 3, Appendix S1: Fig. S5 and table S3).

281
While the effects of increasing proportion of coniferous trees and soil depth were positive, 282 increasing VPD had a negative association with λ. The influence of increasing pH was also 283 negative but smaller than the effect of VPD. The effect of density was negligible.

286
Our results demonstrate substantial effects of climatic, non-climatic abiotic, and biotic 287 environmental drivers on the population dynamics of the perennial herb Actaea spicata. 288 Specifically, we showed that population growth rate is strongly influenced by vapor pressure 289 deficit (VPD), soil depth, and tree composition. There was also a substantial difference in the  The negative effects of VPD on several vital rates and population performance is particularly 301 interesting in the context of ongoing climatic changes. If humidity in Sweden will increase at 302 the same rate as it has during the last 50 years (Wern 2013), and temperature increases as and predicting how environmental change will affect species. We conclude that 343 environmentally explicit demographic models, fitted to data over large geographical scales, 344 and including a broad variety of environmental drivers are a valuable tool for identifying 345 species' niches. We also assert that these analyses can be valuable for designing management 346 strategies for species of special interest. Finally, combining these demographic models with 347 classical species distribution models (e.g. Greiser et al. 2020) has the potential to lead to a 348 deeper understanding of the drivers of changes in species abundances and distributions.      Table S1: Significance of environmental variables in the vital rate regression models over all transitions. Estimates (with standard deviations in parentheses) are from models including all tested environmental factors, but with non-significant polynomial terms omitted. • = p < 0.1, * = p < 0.05, ** = p < 0.01, *** = p < 0.001, NS = quadratic term was not significant in the model where all non-linear terms were included and has therefore not been included in the reduced model.    Table S4: Significance of environmental variables in the vital rate regression models for the individual transitions. Estimates (with standard deviations in parentheses) are from models including all tested environmental factors, but with non-significant polynomial terms omitted. • = p < 0.1, * = p < 0.05, ** = p < 0.01, *** = p < 0.001, NS = quadratic term was not significant in the model where all non-linear terms were included and has therefore not been included in the reduced model.

Appendix S2
Section 1: Demographic data collection

Size data
The size of individuals of Actaea spicata was calculated from plant height and stem diameter following  and as given in Where height is defined by the distance from the ground to the horizontal plane formed by the largest leaves.
To measure this, we used a ruler starting with 0 and measure from ground to horizontal plane, ignoring the inflorescences (Fig. S1). Stem diameter was measured on the vegetative stem, 1 cm below the point, where the different leaves branch. The diameter of A. spicata is not round but oval. We therefore measured the point with the largest width as shown in Fig. S2.

Fecundity data
Actaea spicata has up to four inflorescences, sometimes even more. The inflorescences do not start growing at the same time and therefore were in various stages during the census (e.g. flowering or already having fruits). The order in which inflorescences appear is shown in Fig. S3. We recorded the number of flowers or fruits in each inflorescence separately. Figure S3: Flowering-order of Actaea spicata

Section 2: Environmental variable collection
In addition to the demographic data we collected a range of environmental variables regarding climate, soil nutrients, topography and plant community for each site. We collected weather data throughout the whole study period. Air-temperature and humidity loggers (EasyLog EL-USB-2 from Lascar electronics) were placed in the center of each plot, tied to a wooden pole and protected from rain with a plastic mug equipped with holes to provide air circulation (as described in Terando et al. 2017). We used the logger data to calculate growing degree days with a base temperature of 5 °C for the spring period 15.05-15.06 and first spring date by the definition from the Swedish meteorological office (SMHI 2011), which is seven days in a row with mean temperatures above 0 °C after the 15th of February. In addition, we calculated Vapor Pressure Deficiency, a temperature corrected measure for moisture in the air , during summer (15.06-15.08, hereafter VPD) using the R package plantecophys  since it is highlighted as more important for plants compared to relative humidity (Ashcroft and Gollan 2013).
In 2018 we took soil samples from five different spots within a site. The soil was taken from around 10 cm below ground (after scraping away litter from the surface). The spots were chosen so that they are spread out evenly over the site (e.g. like on a dice). If possible, the soil was taken from inside the patch. In case of a very small patch or a very dense population, we took samples from the closest possible location. Soil samples from one site were mixed and stored in a plastic-zip bag where we pressed out the air before storing them in a cooler. The samples were then analyzed in the laboratory at the University of Southern Denmark. In 2018 we collected soil samples and analyzed the concentrations of several important plant nutrients. Nitrate (NO3 -) and exchangeable ammonia (NH4 + ) was extracted from fresh soil by a 2M KCl solution (Keeney 1982). Both NO3and NH4 + were analyzed spectrophotometrically on a SKALAR SAN plus Analyzer. Plant available Phosphorus (P) was extracted from fresh soil by a 0.5M NaHCO3 solution at pH 8.5 (Olsen 1954) and phosphate (PO4 3-) was determined spectrophotometrically, according to Koroleff and Grasshof (1983). Plant available potassium (K) was extracted from air dried soil by a mixture of 0.10M ammonium lactate and 0.40M acetic acid at pH 3.75 for 90 minutes (Egnér et al. 1960) and analyzed on an ICP OES (Perkin Elmer Optima 2100 DV instrument). Soil pH was measured pH on a suspension of fresh soil and water. In addition, we measured soil depth using a simple metal stick which was put into the ground towards the center of earthdisregarding potential slopes of the side. As with the soil samples, we took at least five soil depth samples, spread out like on a dice, minimizing the distance from each plant to measurement-points.
The collected topographic information including exact latitude and longitude positions of the populations, taken with a Garmin GPS map 64s when standing in the middle of the plot as well as the slope inclination of the site measured using a digital mechanic's level in the android app Bubble level galaxy. Slope inclination was measured from the highest point of the site to the lowest in the way that one person stood at the highest point of the site, the other person on the lowest point. The upper person will look at a point on the lower person that matches its own eye height. Using the mechanic's level, we can write down the approximate slope of the site (Fig. S4 and S5). Aspect direction of the slope was determined using a compass. Slope inclination and aspect direction were used to determine the variable southwards slope, calculated as standardized aspect of the site multiplied with the standardized aspect direction as distance to South, were 0 indicated a slope towards South while 190 indicates a slope aspect towards North. This variable describes the extent to which the ground was slanted towards the south and thereby acts as an additional indicator for the exposure to light and moisture. In addition, we measured site area and calculated intraspecific plant density.
In 2018 we took canopy cover pictures using the back camera of a Sony Xperia L1 with an attached fisheyelens (180° Supreme Fisheye Lens, Model MFE4 by MPOW Inc). Pictures were taken pointing away from the earth center perpendicular to the ground, disregarding the potential slope of the site, from 60 cm above ground using a stick to ensure the same distance for every picture. Those images were then processed in ImageJ ) using the plugin Hemispherical 2.0  to calculate the gap fraction. If different parts of the site had substantially different canopy covers, several pictures from different spots within the site were taken and then averaged.
Lastly, we collected data on the plant community. We defined tree species and their abundance using a relascope, a forestry instrument to measure the basal area, for each tree visible from the center of the site to determine the percentage of coniferous vs. broadleaf trees. Looking through the relascope, a tree that is smaller than the relascope opening was counted, one that exactly fits the opening counted as 0.5, and trees bigger than the opening are counted as one (Fig. S6). Each stem was counted individually at the height of the relascope (≙ eye height). One exception for this procedure were hazel trees since they often grow in a cluster with stems being too small to be measured in the traditional way. Instead, we measured the stem diameter, as we expected it to be if all stems were tightly layered together. Trees that were smaller than 5 cm in diameter were not counted. Figure S4: Upper and Lower are standing on even ground. Figure S5: Upper stands at the highest point, Lower at the lowest. Figure S6: Looking through a relascope the left tree would not count; the right tree would be counted as 0.5; trees bigger are counted as 1.