Determinants of fine-scale plant diversity in dry calcareous grasslands within the Baltic Sea region

Determinants of fine-scale plant diversity in dry calcareous grasslands within the Baltic Sea region

Agriculture, Ecosystems and Environment 182 (2014) 59–68 Contents lists available at ScienceDirect Agriculture, Ecosystems and Environment journal h...

2MB Sizes 0 Downloads 53 Views

Agriculture, Ecosystems and Environment 182 (2014) 59–68

Contents lists available at ScienceDirect

Agriculture, Ecosystems and Environment journal homepage: www.elsevier.com/locate/agee

Determinants of fine-scale plant diversity in dry calcareous grasslands within the Baltic Sea region Triin Reitalu a,∗ , Aveliina Helm b , Meelis Pärtel b , Karin Bengtsson c , Pille Gerhold b,d , Ejvind Rosén e , Krista Takkis b , Sergey Znamenskiy f , Honor C. Prentice g a

Institute of Geology, Tallinn University of Technology, Ehitajate tee 5, 19086 Tallinn, Estonia Department of Botany, Institute of Ecology and Earth Sciences, University of Tartu, Lai 40, 51005 Tartu, Estonia c School of Culture, Energy and Environment, Gotland University, 62167 Visby, Sweden d Department of Experimental Plant Ecology, Institute for Water and Wetland Research, Radboud University Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands e Department of Ecology and Genetics, EBC, Uppsala University, Norbyvägen 18D, 75236 Uppsala, Sweden f Institute of Biology, Karelian Research Centre RAS, Pushkinskaya 11, 185910 Petrozavodsk, Karelia, Russian Federation g Department of Biology, Lund University, Ecology Building, 22362 Lund, Sweden b

a r t i c l e

i n f o

Article history: Received 4 July 2012 Received in revised form 12 October 2012 Accepted 12 November 2012 Available online 16 January 2013 Keywords: Alvar Biogeographic Habitat history Immigration histories Phylogenetic clustering Species density Standardized effect size

a b s t r a c t We used an extensive dataset (1220 vegetation plots of 1 m2 ) to study vegetation gradients and fine-scale plant diversity in dry calcareous grasslands (including alvar grasslands) in the Baltic Sea region. The study area covers the entire European distributional range of alvar habitats: Sweden (Öland, Gotland, Götaland), Estonia (Saaremaa, Hiiumaa, north Estonia, west Estonia), and western Russia (Izhora, Izborsk). Fine-scale plant diversity was characterized by species richness and standardized phylogenetic diversity (comparing the observed mean pairwise phylogenetic distance (MPD) with MPD values from random communities). Ordination techniques (DCA) were used to characterize the main vegetation gradient. Variables describing local environment, climate, the biogeographic composition of the plant communities, and geographic location were related to fine-scale species richness and phylogenetic diversity using variation partitioning techniques and linear mixed models. The main vegetation gradient in the dry calcareous grasslands in the Baltic Sea region had a strong geographic component, was associated with soil depth, species’ stress- and disturbance-tolerance and the age of the grassland habitat. Fine-scale phylogenetic diversity and species richness were negatively associated suggesting that these two diversity components are influenced by different sets of environmental and historical parameters. Fine-scale species richness was unimodally associated with the main vegetation gradient, and the highest levels of species richness were found under intermediate environmental (disturbance, light conditions and temperature) conditions where there was a mixture of species from different biogeographic regions. In contrast to species richness, fine-scale phylogenetic diversity was negatively associated with the main vegetation gradient. The highest phylogenetic diversity was found in the extremely thin-soiled alvar grasslands in Götaland and on the Baltic islands (especially on Öland) where the high phylogenetic diversity is likely to be a reflection of a long history of continuous openness that has allowed time for the “collection” of phylogenetically different species within these unique habitats. © 2012 Elsevier B.V. All rights reserved.

1. Introduction In Europe, the climax community for most of the boreo-nemoral region is forest. The majority of European grasslands have an anthropogenic origin and require at least a moderate level of management (grazing or mowing) for their long-term persistence. Naturally open habitats only persist where conditions are unfavourable for trees; for example, on floodplains, in saline coastal

∗ Corresponding author. Tel.: +372 530321; fax: +372 6203011. E-mail address: [email protected] (T. Reitalu). 0167-8809/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.agee.2012.11.005

areas and on thin soils overlying bedrock (Laasimer, 1965). Primary open habitats with sparse vegetation also characterize regions that have relatively recently emerged from the sea as a result of neotectonic land uplift (Zobel and Kont, 1992). Within the Baltic Sea region, dry calcareous grasslands occupy a range of habitats: from those where the vegetation is more-orless naturally open (for example, in sites on limestone outcrops), to those on deeper soils and calcareous moraines where the grassland vegetation has a clearly secondary origin. Calcareous grasslands that occur on thin soils (generally <20 cm) on Ordovician or Silurian limestone are referred to as “alvars” (Laasimer, 1965; Ekstam and Forshed, 2002). The word “alvar” originates from

60

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

Swedish, where it applies in a broader sense to “unproductive ground” (Pärtel et al., 1999a). In Europe, alvars are found within an area that stretches from southeastern Sweden, through Estonia and into north-western Russia (Laasimer, 1965; Ekstam and Forshed, 2002) and the vegetation composition and environmental conditions of European alvar habitats show considerable variation (Dengler and Löbel, 2006; Znamenskiy et al., 2006; Helm et al., 2007). Similar plant communities on thin calcareous bedrock have also been described in southeastern Canada (Belcher et al., 1992; Schaefer and Larson, 1997). Dry calcareous grasslands have a substantial nature conservation value because they contain a high diversity of specialized species belonging to a range of different taxonomic groups of organisms (plants, insects, lichens, birds) (Schaefer and Larson, 1997; Veen et al., 2009). European grasslands also represent an important historical and cultural heritage. The history of human influence on the alvar grasslands in the Baltic Sea region dates back to the Neolithic (7000–4000 years ago) (Königsson, 1968; Poska et al., 2004). Alvars have interested biologists since 1741 when Carl von Linné visited Öland and Gotland (Linnaeus, 1745). The largest area of alvar is found on the Baltic island of Öland where the Great Alvar, with its 255 km2 of grassland habitat, has become the main centre of alvar research in Europe. Research on the Great Alvar has focused on wide range of different topics – from habitat management and restoration (e.g. Rosén, 1982; Bakker et al., 2007) to vegetation patterns (e.g. Bengtsson et al., 1988; Rusch and Fernández-Palacios, 1995), determinants of species diversity (e.g. van der Maarel and Sykes, 1993; Löbel and Dengler, 2007) and genetic diversity (e.g. Prentice et al., 1995; Lönn and Prentice, 2002). Estonian alvar grasslands have also been the subject of a number of plant ecological studies focusing on community classification (Pärtel et al., 1999a), environmental conditions (Pärtel and Helm, 2007), and plant diversity patterns (e.g. Pärtel and Zobel, 1999; Zobel et al., 2000; Helm et al., 2006, 2009; Pärtel et al., 2007). In addition, a few studies have compared the vegetation in different alvar regions (Znamenskiy et al., 2006; Helm et al., 2007) or examined relationships between the vegetation composition of alvar grasslands and other calcareous grasslands in Europe (Dengler and Löbel, 2006). Even though alvars and associated dry calcareous grasslands have been the subject of a large number of plant ecological and phytosociological studies, there are no vegetation studies that have covered the whole distributional range of European alvar grasslands. Environmental conditions in the most extreme, thin-soiled alvar sites are characteristically harsh: the soils typically dry out completely during the summer but are also often waterlogged during the spring or autumn. Frost heaving may cause severe physical disturbance during the winter months (Rosén, 1982; Pärtel et al., 1999a). Extremely thin-soiled alvar areas, such as those that are characteristic of the southeastern part of the Great Alvar on Öland, are thought to have remained more-or-less open during most of the Holocene (Königsson, 1968). Some of the taxa that are presently common in the alvar areas are known to have been present in the open habitats around the margins of the retreating ice-sheet during the Late Glacial (e.g. Gypsophila, Helianthemum, Artemisia, Saxifraga) (Berglund, 1966; Königsson, 1968; Veski et al., 2012) and it is likely that the populations of these light-demanding alvar species have been able to persist throughout the Holocene in fragments of open habitat in the otherwise forested landscape (Königsson, 1968; Bengtsson et al., 1988). Many of the characteristic alvar species, especially on the Great Alvar on Öland, are on the limits of their geographic ranges or represent disjunct occurrences outside their typical distributional ranges (Bengtsson et al., 1988). These species belong to a range of phytogeographical groups, for example, Artemisia rupestris and Astragalus danicus have their main distributional areas in steppe regions in southern Siberia and southeastern

Europe, Poa alpina and Draba incana in northern European mountain areas and Sedum album and Gypsophila fastigiata in southern and central Europe (Hultén and Fries, 1986). Fine-scale plant species richness in semi-natural grasslands within local landscapes has been shown to depend on a range of local variables: present-day management and disturbance regime (Rosén, 1982; Schnoor and Olsson, 2010), management history (Pärtel et al., 1999b; Reitalu et al., 2010), nutrient status (Diekmann et al., 2004; Löbel et al., 2006), and environmental and microclimatic heterogeneity (Rosén, 1982; Pärtel and Helm, 2007). Extending the extent of a study area from a local (within an individual landscape) to a regional (extending over several thousands of square kilometres) scale is expected to increase the importance of large-scale factors, such as climatic variation and species’ immigration histories, as determinants of fine-scale species richness (Pärtel and Zobel, 1999; Harrison and Cornell, 2008). Both patterns of immigration from different geographic areas and species’ bioclimatic tolerance limits are expected to have a strong influence on the composition of a particular regional species pool – which, in its turn, will influence levels of local and finescale species richness within plant communities (cf. Pärtel et al., 1996). We can, therefore, predict that levels of fine-scale species richness are likely to be enhanced in regions where species have a wide range of geographic origins, and that the highest levels of fine-scale species richness will be found in communities that have diverse biogeographic composition. Biodiversity is most commonly characterized by measures of taxonomic diversity, such as species richness. Biodiversity can also be characterized using measures of phylogenetic diversity, based on the variety of different evolutionary lineages (Faith, 1992). Phylogenetic diversity has a high conservation value because it reflects both the evolutionary history of a community and its potential to persist in a changing environment (Sechrest et al., 2002; Forest et al., 2007). Local phylogenetic diversity depends on largescale and long-term processes such as speciation and immigration (Proches¸ et al., 2006). For example, Gerhold et al. (2008) showed that in northern Europe, phylogenetic diversity in local communities is associated with the size of the regional species pools that have been built up by species’ immigration since the last glaciation. The aim of the present study is to characterize the fine-scale taxonomic and phylogenetic diversity in alvar grasslands (and associated dry calcareous grasslands) over the entire European distribution of alvar habitats (Sweden, Estonia, northwestern Russia). We have assembled a large set (1220 vegetation plots) of fine-scale (1 m2 ) vegetation data and we ask the following questions: (1) How is the vegetation of dry calcareous grasslands related to large scale environmental and geographic gradients across the Baltic Sea region? (2) What is the relationship between fine-scale species richness and phylogenetic diversity in dry calcareous grasslands in the Baltic Sea region? (3) What are the main environmental and climatic factors associated with fine-scale grassland species richness and phylogenetic diversity? (4) How is the biogeographic composition of the grassland vegetation related to fine-scale species richness and phylogenetic diversity? 2. Materials and methods 2.1. Study area The study area covers a ca. 200 000 km2 region (13–30◦ E, 56–60◦ N) (Fig. 1) that can be divided into nine distinct grassland regions in Sweden (Öland, Gotland, Götaland), Estonia (Saaremaa,

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

61

Finland

Sweden Estonia

Latvia

Götaland

Gotland

Hiiumaa

N Estonia

Öland

Saaremaa

W Estonia

Izborsk

Russia

0

100 km

Izhora

Fig. 1. Locations of the studied alvar grasslands.

Hiiumaa, north (N) Estonia, west (W) Estonia) and western Russia (Izhora, Izborsk). While the Götaland, Öland, N Estonia and Izhora regions are situated on Ordovician limestone, the grasslands in W Estonia, Saaremaa and Gotland occur on Silurian limestone outcrops. The Izborsk region is underlain by Devonian limestone. Mean annual temperature ranges from 4.3 ◦ C in Izhora to 7.6 ◦ C on Öland and annual precipitation ranges from 509 mm on Öland to 675 mm in Izhora (Appendix A). Comparable information on the current management status was not available for all the grasslands from which data were collected. At the time of data collection, the majority of the sites were either currently under grazing management or had been grazed within the last decades. The definition of “alvar” is usually based on soil depth (<20 cm) and not on any biological criterion. However, individual alvar sites often contain grassland patches with different soil depths, including areas where the soil is deeper than 20 cm. In the present study, we did not want to rely only on a soil-depth-based definition of alvar, and our intention was to compile the dry grassland data from the whole distribution area of alvars including vegetation data from calcareous grasslands that are in close connection to the alvars sensu stricto. In the text, we therefore refer to “alvars and associated dry calcareous grasslands”. 2.2. Vegetation data Vegetation data were assembled from earlier published studies in which the authors of the present study have been involved. In regions where there were no available existing data, we carried out additional fieldwork. In each of the nine study regions, the vegetation was characterized in several individual grassland sites (3–36 sites per region) (Appendix A). In each grassland site, the presence of all herbaceous vascular plant species was recorded within 1 m × 1 m plots (3–15 plots per site). The vegetation data from Öland (collected in 1984) consists of 254 plots from 28 sites (Bengtsson et al., 1988). In Izhora 75 plots in 5 sites were described in 2005 (Znamenskiy et al., 2006). In Estonia, the fieldwork was carried out over longer period (1992–1999) and the dataset consists of 394 plots in 36 sites in Saaremaa, 45 plots in 3 sites in Hiiumaa, 180 plots in 12 sites in N Estonia, 135 plots in 9 sites in W Estonia (Pärtel et al., 1999a, 2001; Helm et al., 2007). On Gotland, the vegetation data were collected in 1999 from 45 plots in 9 grassland sites by

E. Rosén, R. Kalamees, N. Ingerpuu and M. Otsus but have not been included in any publication. We carried out additional fieldwork in Götaland (13 plots, 3 sites (Kinnekulle, Österplana and Högstena)) in 2006 and in Izborsk (79 plots, 6 sites) in 2007–2008. Altogether, 1220 vegetation plots were included in the present study. The vegetation data were checked for synonyms and species with possible identification problems were pooled leaving 292 species/taxa for the analyses (Appendix B). Vegetation in dry calcareous grasslands has been shown to respond relatively slowly (on a time scale ranging from half a century to several centuries) to environmental changes (Helm et al., 2006; Reitalu et al., 2010). Therefore, we can assume that the plant data that are included in the present study are comparable – despite having been collected over two decades. The majority of the vegetation plots were described either early or late summer when the moisture conditions were suitable for the growth of annual species. 2.3. Climate Climate data for each of the studied grassland sites were extracted from the WorldClim global climate database (“bioclim” database with 2.5 arc minutes resolution, www.worldclim.com, Hijmans et al., 2005). The bioclim database includes 19 bioclimatic variables that are calculated from monthly temperature and rainfall values (Hijmans et al., 2005). Most of the bioclimatic variables are highly intercorrelated in our study region (Appendix C) and we used only a subset of the variables excluding variables with high pairwise correlations (Pearson’s r >0.85). The seven bioclimatic variables included in the analyses were: minimum temperature of the coldest month, maximum temperature of the warmest month, mean temperature of the wettest quarter (3-month period), mean temperature of the driest quarter, precipitation of the driest month, precipitation of the warmest quarter and isothermality (mean of monthly temperature range divided by annual temperature range). 2.4. Biogeographic composition To characterize the biogeographic composition of each of the vegetation plots, we used the Atlas of North European Vascular Plants (Hultén and Fries, 1986) to visually classify each of the species into a biogeographic distribution category. When most of

62

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

a species’ distributional area lies to the west of our study region, the species was classified as “western”; “eastern”, “southern” and “northern” were defined in a similar way (Appendix B). Additional distribution classes were “central” (species with their distributions centred on our study area) and “local” (species that are endemic to our study area) (Appendix B). In case of the taxa that were pooled because of identification problems, the pooled taxon was only used in the characterization of biogeographic composition when all the species that it included had similar distributions. Each vegetation plot was characterized by proportions of all six biogeographic distributional categories. For example, the proportion of western species was calculated as the number of species in the western distributional category divided by the total number of species in the plot. For statistical analyses, the proportions of different distribution categories were log-ratio transformed in relation to the proportion of species with central distribution. For example, in the case of the western distributional category, the log-ratio was calculated as: log10 (the proportion of species with western distribution/the proportion of species with central distribution). In addition to the proportional data, we calculated biogeographic richness and biogeographic evenness. Biogeographic richness was defined as the number of different biogeographic distributional categories within each plot. Biogeographic evenness was calculated as Simpson’s measure of evenness: E = 1/sum(p2i )]/S where pi is the proportion of ith biogeographic category and S is the biogeographic richness. 2.5. Local environment We used Ellenberg indicator values for “light” (L), “temperature” (T), “moisture” (F), “reaction” (R), “soil fertility” (N) and “continentality” (K) (Ellenberg et al., 1991) to characterize the environmental preferences of the species within the study plots. The local environment within each plot was characterized by the mean of the indicator values for all the species that were present in the plot. 2.6. Spatial variation To describe spatial variation, we used a set of nine geographic variables constructed from the basic geographic coordinates of the grassland sites – latitude (x) and longitude (y) and their quadratic (x2 , y2 ) and cubic (x3 , y3 ) terms: x, y, x2 , y2 , x*y, x2 *y, x*y2 , x3 , y3 (Borcard et al., 1992). 2.7. Grassland diversity measures We used two measures to characterize the fine-scale diversity in the grassland communities. Species richness was defined as the total number of vascular plant species/taxa within the 1 m × 1 m vegetation plots. Phylogenetic diversity was calculated using the index “mean pairwise distance” (MPD, Pavoine and Bonsall, 2011), which was computed as the mean phylogenetic distance (measured as branch length) between all pairwise combinations of species on the phylogenetic tree. The phylogenetic tree was constructed for the 250 species/taxa in the present study (the regional species pool) (excluding 42 species/taxa which occurred in only one vegetation plot or had no available phylogenetic information) on the basis of the phylogeny for higher plants of Central Europe (BiolFlor; Kühn et al., 2004). We used a standardized measure of phylogenetic diversity – the standardized effect size (SES) of MPD that has been shown to be intrinsically independent of species richness (Pavoine and Bonsall, 2011). The observed MPD values (MPDobs ) were compared with the MPD values from 999 random communities (drawing the same number of species as that in the plot randomly from the set of 250 species). SESMPD was calculated as (MPDobs − MPDrand.mean )/MPDrand.sd , where MPDrand.mean

and MPDrand.sd are the mean and the standard deviation of MPD in the randomized communities. Negative SESMPD values indicate phylogenetic clustering and positive values indicate phylogenetic over-dispersion. SESMPD values lower than −2 and higher than 2 indicate significant (p < 0.05) clustering and over-dispersion, respectively (Gotelli and Graves, 1996). We used the R statistical programming environment (R Development Core Team, 2009) and the package “picante” for calculation of the phylogenetic diversities. 2.8. Statistical analyses We used an ordination approach to identify and interpret the main environmental and geographic gradients in the calcareous grasslands of the study region. We first carried out a detrended correspondence analysis (DCA) on the full dataset of 1220 vegetation plots. DCA was used because simple correspondence analysis showed a clear arch effect (Hill and Gauch, 1980). To visualize the geographic gradients, we plotted the first two ordination axes and displayed the vegetation plots from different geographic regions using different symbols. The ecological distributions of the individual species were also summarized by displaying their occurrences on the first two axis of the plot-based ordination (cf. Bengtsson et al., 1988). Linear models were used to describe the associations between species richness, phylogenetic diversity (SESMPD ) and the plot scores on the first DCA axis (DCA1). Significance testing included tests for quadratic terms, to account for possible non-linear relationships. Linear mixed effects (LME) models were used to explore the relationships between grassland diversity measures (species richness and phylogenetic diversity (SESMPD )) and explanatory variables describing climate, biogeographic composition and local environment, with grassland region (nine categories) and grassland site (111 categories) included as random variables. For both response variables (species richness and phylogenetic diversity), separate LME models were constructed for each of three following groups of fixed explanatory variables: (1) climate (minimum temperature of the coldest month, maximum temperature of the warmest month, mean temperature of the wettest quarter, mean temperature of the driest quarter, precipitation of the driest month, precipitation of the warmest quarter, and isothermality); (2) biogeographic composition (log-ratio transformed proportions of western, southern, northern, eastern and local species, biogeographic richness, and biogeographic evenness); (3) local environment (“light”, “temperature”, “soil fertility”, “continentality”, “moisture”, and “reaction”). The quadratic terms of the fixed variables were included in the LME models, to account for the possible non-linear relationships. All the explanatory variables were standardized to zero mean and unit variance prior to the analysis. In each LME model, we followed the variable selection procedure for mixed effects models proposed by Zuur et al. (2009). We started the model selection procedure with full model – including all the fixed variables and their quadratic terms, together with both the random variables. First, we selected the best random structure (area and/or site) by choosing the model that gave the lowest AIC value. We then used a stepwise backward selection procedure for the fixed variables where, at each step, the variable (or quadratic term) with the highest pvalue for the t-statistic was left out until only significant (p < 0.01) terms remained in the model. Variation partitioning (Borcard et al., 1992) was used to quantify the shared and unique effects of the three groups of explanatory

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

63

Öland

Gotland

Saaremaa

Hiiumaa

W Estonia

N Estonia

Izborsk

Izhora

0 −2

−1

0

1

2

−2

−1

DCA2

1

2

−2

−1

0

1

2

Götaland

−3

−2

−1

0

1

2 −3

−2

−1

0

1

2 −3

−2

−1

0

1

2

DCA1 Fig. 2. Plot of the first two DCA axes with alvar locations superimposed.

variables – climate, biogeographic composition and local environment – on species richness and phylogenetic diversity. Spatial variation (described by nine spatial variables derived from latitude and longitude) was included, as a fourth group of variables, in the variation partitioning analysis. In each group, only the significant variables derived from the model selection procedure were included in the variation partitioning. All statistical analyses were carried out in the R environment (R Development Core Team, 2009) using the packages “nlme” for LME models and “vegan” for variation partitioning. 3. Results The species richness within the 1 m × 1 m plots varied between 5 and 49 (mean = 21.1) and phylogenetic diversity ranged between −4.2 and 1.8 (mean = −0.7). The highest average species richness (28.4 ± 7.6) was found in W Estonia and the lowest on Gotland (15.4 ± 5.9) (Appendix A). The highest average phylogenetic diversity was found on Öland (−0.24 ± 0.58) and the lowest in W Estonia (−1.24 ± 0.84) (Appendix A). The phylogenetic diversity was skewed towards phylogenetic clustering, with 8% of the studied plots (95 plots) showing significant clustering (SESMPD < −2) and none of the plots showing significant over-dispersion. The Eigenvalues for the first two DCA axes were 0.47 and 0.27, respectively. The distribution of the different grassland regions on the plot of the first two DCA axes shows a strong geographic

clustering (Fig. 2), with the vegetation plots from Öland, Gotland, and Hiiumaa situated in the left part of the ordination plot, and the vegetation plots from N Estonia, W Estonia, Izborsk, and Izhora grouped on the right side of the ordination plot. Both latitude and longitude showed a highly significant (p < 0.001) correlation with the first DCA axis (Pearson’s c 0.69 and 0.83, respectively). Among the local environmental variables (mean Ellenberg indicator values), “light” had the highest correlation with DCA1 (r = −0.68, p < 0.001). Species scores on the first DCA axis are given in Appendix B and the ecological distributions on the site ordination, for a chosen set of species, are presented in Fig. 3. Species richness and phylogenetic diversity were highly significantly negatively associated (Fig. 4). While species richness showed a unimodal association with DCA1, with the highest values in the middle of the axis, phylogenetic diversity was negatively associated with DCA1. The results of LME models show that species richness is negatively associated with most of the individual biogeographic categories but positively associated with biogeographic richness and evenness (Table 1) indicating that species richness is highest in grasslands that contain species with a mixture of different biogeographic origins. Phylogenetic diversity, on the other hand, was positively associated with the proportion of western species and unimodally associated with the proportion of local species. Among the climate variables, isothermality was negatively associated with species richness: grasslands in areas with relatively even annual temperatures have lower species richness, and

64

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

Antennaria dioica

Dactylis glomerata

Sedum album

Helianthemum nummularium

Centaurea scabiosa

Saxifraga tridactylites

Phleum phleoides

Primula veris

Artemisia rupestris

Potentilla tabernaemontani

Leucanthemum vulgare

Festuca ovina

Linum catharticum

Poa compressa

0 −2

−1

0

1

2

−2

−1

0

1

2

−2

−1

DCA2

1

2

−2

−1

0

1

2

−2

−1

0

1

2

Allium schoenoprasum

−3

−2

−1

0

1

2 −3

−2

−1

0

1

2 −3

−2

−1

0

1

2

DCA1 Fig. 3. Examples of DCA ordination plots visualizing the ecological amplitudes of species. Plots in which the species of interest are present are indicated by black dots, plots without the species are indicated by grey dots. The species were chosen to show examples of different distributions along the first DCA axis were chosen to represent the whole distributional range along the first DCA axis.

2 1 0 −1 −2

Phylogenetic diversity

65

−4

−3

50 30

Richness

10

20

30 20 10

Richness

c

40

b

40

a

50

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

−4

−3

−2

−1

0

1

2

−2

−1

Phylogenetic diversity

0

1

2

−2

−1

DCA1

0

1

2

DCA1

Fig. 4. Associations between fine-scale diversity measures and first DCA axis in alvar grasslands. (a) Species richness and phylogenetic diversity, (b) species richness and first DCA axis, and (c) phylogenetic diversity and first DCA axis. An R2 value is given for each association, all three associations were highly significant (p < 0.0001).

grasslands in areas with larger annual temperature fluctuations have higher species richness (Table 1). The maximum temperature of the warmest month was unimodally associated with species richness and negatively associated with phylogenetic diversity. In addition, phylogenetic diversity was positively associated with the mean temperature of the wettest quarter. Among the local environmental variables (mean Ellenberg indicator values), “light”, “temperature”, “continentality” and “moisture” were unimodally associated with species richness. The only environmental variable that showed a positive association with richness was “reaction” (Table 1). Phylogenetic diversity was only associated with “temperature” (U-shaped association). In the variation partitioning analysis, the four sets of explanatory variables (biogeographic composition, environment, climate and spatial variation) explained 53.7% of the total variation in species richness and 33.2% of the total variation in phylogenetic diversity. The largest fractions of the variability in species richness were accounted for by local environmental variables (9.2%), biogeographic composition (8.3%) and the joint effect of local environmental variables and biogeographic composition (9.5%) (Fig. 5). In the case of phylogenetic diversity, the largest proportion of variation was accounted for by the joint effect of biogeographic composition and spatial variables (8.2%) followed by the joint effect of climate and spatial variables (6.9%) (Fig. 5). The largest individual contributions to the total variation were made by local

environmental variables (for species richness) and by biogeographic composition (for phylogenetic diversity).

4. Discussion The present study is the first in which vegetation gradients, fine-scale species richness, and phylogenetic diversity, have been analyzed for alvar grasslands (and associated dry calcareous grasslands) from the entire European distributional range of alvar habitats. Earlier, more local, studies of alvar grasslands have shown that the main gradient in vegetation composition is related to soil depth (Bengtsson et al., 1988; Pärtel et al., 1999a; Löbel and Dengler, 2007). However, while the comparison of the ecological distributions of different species on the ordination plot (Fig. 3) in the present study reveals an association between DCA1 and soil depth, other factors are also associated with the dominant axis of variation in grassland community composition. The first DCA axis can also be interpreted in terms of a disturbance-competition gradient, with species known for their tolerance of stress and disturbance (e.g. S. album, Saxifraga tridactylites) (cf. Grime, 2001) having lower scores, and species with high competitive ability (e.g. Dactylis glomerata, Leucanthemum vulgare) (cf. Grime, 2001) having high scores (Fig. 3). DCA1 also has a strong geographic component. The extremely

Table 1 Summary of linear mixed effects (LME) models for species richness and phylogenetic diversity. Site and region are included as categorical variables in the random part of the models. Three different models with explanatory variables characterizing (1) biogeographic composition, (2) climate, and (3) environment are presented for both species richness and phylogenetic diversity. The fixed variables within each model were selected with the help of backward selection with 0.01 significance level. Quadratic terms are marked with ∩ for hump-back associations and with ∪ for U-shaped associations. Model Biogeographic composition No of random groups: 9 regions, 111 sites

Climate No of random groups: 111 sites

Environment (mean Ellenberg values) No of random groups: 9 regions 111 sites

Explanatory variable

Species richness

Proportion of spp. with E distribution Proportion of spp. with N distribution Proportion of spp. with local distribution Proportion of spp. with S distribution Proportion of spp. with W distribution Biogeographic richness Biogeographic evenness

–*** –***, ∩*** –***, ∪*** +***, ∪*** –**, ∪*** +*** +***

Isothermality Maximum temperature of warmest month Mean temperature of wettest quarter

–***, ∪** + n.s., ∩***

“Light” “Temperature” “Continentality” “Moisture” “Reaction”

–*, ∩** –*, ∩*** –n.s., ∩*** –***, ∩*** +***

Phylogenetic diversity

+ n.s. ∪*** +***, ∪*** –*

–*** +***

+ n.s., ∪***

66

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

a

Residuals 46.3 %

Biogeographic composition 8.3 %

0.4%

3.4 %

5.3 %

2.2 %

0.9 %

1.2 %

0.3 % Environment 3.0 %

7.0 %

Climate 1.4 %

8.2 %

3.0 %

Space 4.1 %

0.8 % 0%

-0.6 %

1.4 %

0.2 %

-2.3 %

0.4 % Space 4.1 %

Residuals 66.8 %

Biogeographic composition 6.0 %

-0.8 %

9.5 %

Environment 9.2 %

b

6.9 %

0.7 % -0.5 %

Climate 3.2 %

Fig. 5. Variation partitioning of (a) fine-scale species richness and (b) phylogenetic diversity in terms of biogeographic composition, environment, climate and space.

thin-soiled alvar areas occur mainly within the Swedish alvar regions in Götaland and on Öland and Gotland. These types of extreme alvar communities are also found in Estonia (on Saaremaa and Hiiumaa) whereas the calcareous grasslands on the Estonian mainland (W and N Estonia) and in Russia occur on slightly deeper soils (but usually still <20 cm). Human management has played an important role in the development of dry calcareous grassland vegetation in the entire Baltic Sea area, and the grasslands on the deeper soils require more continuous management than those on the extremely thin-soiled alvars – which can remain open for longer periods of time without grazing management (Rosén, 1982). It seems likely that some fragments of thin-soiled alvar habitat that are unsuitable for shrub and tree growth may have remained more-or-less open since the Late Glacial period (Königsson, 1968). The main vegetation gradient (DCA1) in the present study may, therefore, also be related to the overall age of the habitat, because anthropogenic management only became widespread in the area around 5000–4000 years ago (Berglund, 2003; Poska et al., 2004). The fact that the fine-scale species richness in dry calcareous grasslands in the Baltic Sea region was unimodally related to the main vegetation gradient (DCA1) (Fig. 4) is consistent with the intermediate disturbance hypothesis (Connell, 1978) and the unimodal association between species richness and the Ellenberg variable “light” is also consistent with the same hypothesis. Both the open, highly disturbed, thin-soiled alvar sites (where plants are unshaded), and the poorly grazed grasslands on deeper soils (where competition for light becomes more important) contain fewer species than the intermediate sites in the middle of the vegetation gradient (cf. Löbel and Dengler, 2007). In addition, soil reaction and soil moisture – both important determinants of species richness and vegetation composition (Bengtsson et al., 1988; Tyler, 1996; Löbel and Dengler, 2007) – are also among the variables that show significant associations with species richness in the present study. A positive association between soil reaction and fine-scale species richness is often observed in temperate and boreal regions (Pärtel, 2002). The first DCA axis is also closely related to the west–east geographic gradient (Fig. 2), and the unimodal association between species richness and the first DCA axis shows that the highest levels of fine-scale (within 1 m2 ) species richness are found in sites that are geographically central within the study area, where the western Estonian sites have an average of 28 species/m2 (Appendix A). We predicted that local species richness would also show a relationship with the biogeographic composition of the grassland vegetation and that the highest fine-scale species richness would be found in grassland sites containing a mixture of species belonging to different biogeographic categories. The vegetation analyses confirmed that, whereas most of the individual biogeographic

categories were negatively associated with richness, high finescale species richness was associated with a mixture of species from many different biogeographic regions. This association is consistent with the unimodal association between species richness and climatic variables, and suggests that species from different biogeographic groups (with overlapping bioclimatic requirements) can coexist under intermediate climatic conditions. In addition to the species’ bioclimatic tolerances, the biogeographic categories also reflect species’ immigration histories. Population genetic studies have shown that plant populations in the southern Baltic Sea region, which represents a confluence of different post-glacial immigration routes, show high levels of genetic diversity (Nordström and Hedrén, 2008; Prentice et al., 2011). Analogously, areas where the immigration routes for species from different glacial refugia merge are also likely to be characterized by species-diverse local plant communities. While variation in fine-scale species richness was more-or-less equally well explained by local environmental variables and biogeographic composition (Fig. 5a), phylogenetic diversity was best explained by the biogeographic composition of the communities (Fig. 5b). Proches¸ et al. (2006) related phylogenetic diversity to the evolutionary age of plant communities in South Africa, where young vegetation types (such as grasslands) are phylogenetically clustered whereas the phylogenetic structure in evolutionarily older vegetation types does not deviate from random expectations. All grasslands in the Baltic Sea region are “young”, given the relatively short time since deglaciation, and phylogenetic diversity in northern Europe is likely to be mainly associated with immigration histories after the last Ice Age (Gerhold et al., 2008). In the present study, fine-scale phylogenetic diversity is highest (and corresponds to a random phylogenetic structure) in the highly disturbed sites on the shallowest alvar soils, where the continuity of open habitats is likely to have been the longest (Königsson, 1968; Bengtsson et al., 1988). The long history of openness of the thin-soiled alvars is the most likely factor underlying their high phylogenetic diversity, with the long habitat-continuity having allowed time for the “collection” of phylogenetically different species within these extreme habitats. The occurrence of species with highly disjunct distributions belonging to a range of phytogeographic categories on the Great Alvar on Öland supports this interpretation (cf. Pigott and Walters, 1954). Dry calcareous grassland sites on deeper soils where habitat continuity is dependent on management by humans, are phylogenetically less diverse and also show phylogenetic clustering. The anthropogenic disturbance (cropping and trampling) associated with grazing management, differs from that in the naturally open, thin-soiled alvar vegetation, and has a relatively recent

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

history in the Baltic Sea region. Species in managed grasslands share a range of traits that enable them to withstand the grazing disturbance (Kahmen et al., 2002). If these traits are phylogenetically conserved, human management and the consequent filtering of management-tolerant species into the grasslands would lead to phylogenetic clustering (Webb et al., 2002; Emerson and Gillespie, 2008). In addition, this kind of management-related filtering may interact with competition for light in more productive (or temporarily abandoned) managed grasslands: if competitively stronger species (which displace the other species) share similar phylogenetically conserved traits, competition for light should promote phylogenetic clustering (Grime, 2006; Mayfield and Levine, 2010). The study by Carboni et al. (2014) shows that phylogenetic clustering in grassland vegetation is related to the productivity gradient. The negative association between fine-scale phylogenetic diversity and species richness reflects the fact that these two diversity components are influenced by different sets of environmental and historical parameters and the two components of diversity often showed opposing relationships with the environmental and biogeographic variables in the present study. Role of the funding source The funding sources had no part in planning or performing this study. Acknowledgements We thank Rein Kalamees, Nele Ingerpuu and Merit Otsus for sharing the vegetation data from Gotland. The field-trips that initiated the compilation of the manuscript were financed by grants (to HCP) from the Swedish Institute Visby Programme (01391/2005) and The Royal Swedish Academy of Agriculture and Forestry. The project was supported by Estonian Science Foundation (Mobilitas Programme MJD4 to TR; ETF8323 to MP; ETF8613 to PG; ETF9223 to AH) and by the European Union through the European Regional Development Fund (Centre of Excellence FIBIR). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.agee.2012.11.005. References Bakker, J.P., Rosén, E., Steg, K., 2007. Re-establishment of alvar plant species in abandoned arable fields on Öland. Acta Phytogeogr. Suec. 88, 73–82. Belcher, J.W., Keddy, P.A., Catling, P.M., 1992. Alvar vegetation in Canada: a multivariate description at two scales. Can. J. Bot. 70, 1279–1291. Bengtsson, K., Prentice, H.C., Rosén, E., Moberg, R., Sjögren, E., 1988. The dry alvar grassland of Öland: ecological amplitudes of plant species in relation to vegetation composition. Acta Phytogeogr. Suec. 76, 21–46. Berglund, B.E., 1966. Late Quaternary vegetation in eastern Blekinge. Opera Bot. 12, 1–180. Berglund, B.E., 2003. Human impact and climate changes – synchronous events and a causal link? Quatern. Int. 105, 7–12. Borcard, D., Legendre, P., Drapeau, P., 1992. Partialling out the spatial component of ecological variation. Ecology 73, 1045–1055. ˇ Doleˇzal, J., Hornik, J., Lepˇs, J., Reitalu, T., Carboni, M., de Bello, F., Janeˇcek, S., Klimeˇsova, J., 2014. Trait divergence and convergence along a productivity gradient in wet meadows: implications for biodiversity management. Agric. Ecosyst. Environ. 182, 96–105. Connell, J.H., 1978. Diversity in tropical rain forests and coral reefs. Science 199, 1302–1310. Dengler, J., Löbel, S., 2006. The basiphilous dry grasslands of shallow, skeletal soils (Alysso-Sedetalia) on the island of Öland (Sweden), in the context of North and Central Europe. Phytocoenologia 36, 343–391. Diekmann, M., Dupre, C., van der Maarel, E., 2004. Fine-scale species associations in alvar limestone grasslands. Nord. J. Bot. 23, 115–128. Ekstam, U., Forshed, M., 2002. Svenska Alvarmarker: Historia och Ekologi (Swedish Alvars: History and Ecology). Naturvårdverkets Förlag, Stockholm.

67

Ellenberg, H., Weber, H.E., Düll, R., Wirth, V., Werner, W., Paulissen, D., 1991. Zeigerwerte von Pflanzen in Mitteleuropa. Scripta Geobot. 18, 1–248. Emerson, B.C., Gillespie, R.G., 2008. Phylogenetic analysis of community assembly and structure over space and time. Trends Ecol. Evol. 23, 619–630. Faith, D.P., 1992. Conservation evaluation and phylogenetic diversity. Biol. Conserv. 61, 1–10. Forest, F., Grenyer, R., Rouget, M., Davies, T.J., Cowling, R.M., Faith, D.P., Balmford, A., Manning, J.C., Proches¸, S., van der Bank, M., Reeves, G., Hedderson, T.A.J., Savolainen, V., 2007. Preserving the evolutionary potential of floras in biodiversity hotspots. Nature 445, 757–760. Gerhold, P., Pärtel, M., Liira, J., Zobel, K., Prinzing, A., 2008. Phylogenetic structure of local communities predicts the size of the regional species pool. J. Ecol. 96, 709–712. Gotelli, N.J., Graves, G.R., 1996. Null Models in Ecology. Smithsonian Institution Press, Washington and London. Grime, J.P., 2001. Plant Strategies, Vegetation Processes, and Ecosystem Properties. John Wiley and Sons, Chichester. Grime, J.P., 2006. Trait convergence and trait divergence in herbaceous plant communities: mechanisms and consequences. J. Veg. Sci. 17, 255–260. Harrison, S., Cornell, H., 2008. Toward a better understanding of the regional causes of local community richness. Ecol. Lett. 11, 969–979. Helm, A., Hanski, I., Pärtel, M., 2006. Slow response of plant species richness to habitat loss and fragmentation. Ecol. Lett. 9, 72–77. Helm, A., Oja, T., Saar, L., Takkis, K., Talve, T., Pärtel, M., 2009. Human influence lowers plant genetic diversity in communities with extinction debt. J. Ecol. 97, 1329–1336. Helm, A., Urbas, P., Pärtel, M., 2007. Plant diversity and species characteristics of alvar grasslands in Estonia and Sweden. Acta Phytogeogr. Suec. 88, 33–42. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. Hill, M.O., Gauch, H.G., 1980. Detrended correspondence analysis: an improved ordination technique. Vegetatio 87, 47–58. Hultén, E., Fries, M., 1986. Atlas of North European Vascular Plants: North of the Tropic of Cancer. Koeltz Scientific Books, Königstein. Kahmen, S., Poschlod, P., Schreiber, K.-F., 2002. Conservation management of calcareous grasslands. Changes in plant species composition and response of functional traits during 25 years. Biol. Conserv. 104, 319–328. Königsson, L.-K., 1968. The Holocene history of the Great Alvar of Öland. Acta Phytogeogr. Suec. 55, 1–172. Kühn, I., Durka, W., Klotz, S., 2004. BiolFlor: a new plant-trait database as a tool for plant invasion ecology. Divers. Distrib. 10, 363–365. Laasimer, L., 1965. Eesti NSV taimkate (Vegetation of the Estonian S.S.R.). Valgus, Tallinn. Linnaeus, C., 1745. Carl Linnaei – Öländska och Gothländska Resa – Förrättad Åhr 1741. Stockholm and Uppsala. Löbel, S., Dengler, J., 2007. Dry grassland communities on southern Öland: phytosociology, ecology, and diversity. Acta Phytogeogr. Suec. 88, 13–31. Löbel, S., Dengler, J., Hobohm, C., 2006. Species richness of vascular plants, bryophytes and lichens in dry grasslands: the effects of environment, landscape structure and competition. Folia Geobot. 41, 377–393. Lönn, M., Prentice, H.C., 2002. Gene diversity and demographic turnover in central and peripheral populations of the perennial herb Gypsophila fastigiata. Oikos 99, 489–498. Mayfield, M.M., Levine, J.M., 2010. Opposing effects of competitive exclusion on the phylogenetic structure of communities. Ecol. Lett. 13, 1085–1093. Nordström, S., Hedrén, M., 2008. Genetic differentiation and postglacial migration of the Dactylorhiza majalis ssp. traunsteineri/lapponica complex into Fennoscandia. Plant Syst. Evol. 276, 73–87. Pavoine, S., Bonsall, M.B., 2011. Measuring biodiversity to explain community assembly: a unified approach. Biol. Rev. 86, 792–812. Pigott, C.D., Walters, S.M., 1954. On the interpretation of the discontinuous distributions shown by certain British species of open habitats. J. Ecol. 42, 95–116. Poska, A., Saarse, L., Veski, S., 2004. Reflections of pre- and early-agrarian human impact in the pollen diagrams of Estonia. Palaeogeogr. Palaeoclimatol. Palaeoecol. 209, 37–50. Prentice, H.C., Lönn, M., Lefkovitch, L.P., Runyeon, H., 1995. Associations between allele frequencies in Festuca ovina and habitat variation in the alvar grasslands on the Baltic island of Öland. J. Ecol. 83, 391–402. Prentice, H.C., Andersson, S., Månsby, E., 2011. Mosaic variation in allozyme and plastid DNA markers in the European ranges of Silene vulgaris and its partially sympatric relative S. uniflora (Caryophyllaceae). Bot. J. Linn. Soc. 166, 127–148. Proches¸, S., Wilson, J.R.U., Cowling, R.M., 2006. How much evolutionary history in a 10 × 10 m plot? Proc. Biol. Sci. 273, 1143–1148. Pärtel, M., 2002. Local plant diversity patterns and evolutionary history at the regional scale. Ecology 83, 2361–2366. Pärtel, M., Helm, A., 2007. Invasion of woody species into temperate grasslands: relationship with abiotic and biotic soil resource heterogeneity. J. Veg. Sci. 18, 63–70. Pärtel, M., Zobel, M., 1999. Small-scale plant species richness in calcareous grasslands determined by the species pool, community age and shoot density. Ecography 22, 153–159. Pärtel, M., Helm, A., Reitalu, T., Liira, J., Zobel, M., 2007. Grassland diversity related to the Late Iron Age human population density. J. Ecol. 95, 574–582. Pärtel, M., Kalamees, R., Zobel, M., Rosén, E., 1999a. Alvar grasslands in Estonia: variation in species composition and community structure. J. Veg. Sci. 10, 561–570.

68

T. Reitalu et al. / Agriculture, Ecosystems and Environment 182 (2014) 59–68

Pärtel, M., Moora, M., Zobel, M., 2001. Variation in species richness within and between calcareous (alvar) grassland stands: the role of core and satellite species. Plant Ecol. 157, 203–211. Pärtel, M., Mändla, R., Zobel, M., 1999b. Landscape history of a calcareous (alvar) grassland in Hanila, western Estonia, during the last three hundred years. Landsc. Ecol. 14, 187–196. Pärtel, M., Zobel, M., Zobel, K., van der Maarel, E., 1996. The species pool and its relation to species richness: evidence from Estonian plant communities. Oikos 75, 111–117. R Development Core Team, 2009. R: A Language and Environment for Statistical Computing. R Foundation of Statistical Computing, Vienna, Austria. Reitalu, T., Johansson, L.J., Sykes, M.T., Hall, K., Prentice, H.C., 2010. History matters: village distances, grazing and grassland species diversity. J. Appl. Ecol. 47, 1216–1224. Rosén, E., 1982. Vegetation development and sheep grazing in limestone grasslands of south Öland, Sweden. Acta Phytogeogr. Suec. 72, 1–107. Rusch, G.M., Fernández-Palacios, J.M., 1995. The influence of spatial heterogeneity on regeneration by seed in a limestone grassland. J. Veg. Sci. 6, 417–426. Schaefer, C.A., Larson, D.W., 1997. Vegetation, environmental characteristics and ideas on the maintenance of alvars on the Bruce Peninsula, Canada. J. Veg. Sci. 8, 797–810. Schnoor, T.K., Olsson, P.A., 2010. Effects of soil disturbance on plant diversity of calcareous grasslands. Agric. Ecosyst. Environ. 139, 714–719.

Sechrest, W., Brooks, T.M., da Fonseca, G.A.B., Konstant, W.R., Mittermeier, R.A., Purvis, A., Rylands, A.B., Gittleman, J.L., 2002. Hotspots and the conservation of evolutionary history. Proc. Natl. Acad. Sci. U. S. A. 99, 2067–2071. Tyler, G., 1996. Soil chemistry and plant distributions in rock habitats of southern Sweden. Nord. J. Bot. 16, 609–635. van der Maarel, E., Sykes, M.T., 1993. Small-scale plant species turnover in a limestone grassland: the carousel model and some comments on the niche concept. J. Veg. Sci. 4, 179–188. Veen, P., Jefferson, R., de Smidt, J., van der Straaten, J. (Eds.), 2009. Grasslands in Europe of High Nature Value. KNNV Publishing, Zeist. Veski, S., Amon, L., Heinsalu, A., Reitalu, T., Saarse, L., Stivrins, N., Vassiljev, J., 2012. Lateglacial vegetation dynamics in the eastern Baltic region between 14,500 and 11,400 cal yr BP: a complete record since the Bølling (GI-1e) to the Holocene. Quaternary Sci. Rev. 40, 39–53. Webb, C.O., Ackerly, D.D., Mcpeek, M.A., Donoghue, M.J., 2002. Phylogenies and community ecology. Annu. Rev. Ecol. Syst. 33, 475–505. Znamenskiy, S., Helm, A., Pärtel, M., 2006. Threatened alvar grasslands in NW Russia and their relationship to alvars in Estonia. Biodivers. Conserv. 15, 1797–1809. Zobel, M., Kont, A., 1992. Formation and succession of alvar communities in the Baltic land uplift area. Nord. J. Bot. 12, 249–256. Zobel, M., Otsus, M., Liira, M., Moora, M., Möls, T., 2000. Is small-scale species richness limited by seed availability or microsite availability? Ecology 81, 3274–3282. Zuur, A.F., Ieno, E.N., Walker, N.J., Saveliev, A.A., Smith, G.M., 2009. Mixed Effects Models and Extensions in Ecology with R. Springer, New York.