Testing basic assumptions of species richness hypotheses using plant species distribution data

Dissertation zur Erlangung der naturwissenschaftlichen Doktorw¨ urde (Dr. sc. Nat.) vorgelegt der Mathematisch-naturwissenschaftlichen Fakult¨at der Universit¨at Z¨ urich von Katharina Steinmann von Waltenschwil AG Promotionskomitee Prof. Dr. H. Peter Linder (Vorsitz) Dr. Niklaus E. Zimmermann (Leitung der Dissertation) Dr. Felix Gugerli Z¨ urich, 2008

“Menon: Kann man suchen, was man nicht kennt?”

1

1

Einwand des Menon im Gespr¨ach mit Sokrates. In Platon: S¨amtliche Werke II, 1970. Hrsg. E. Grassi unter Mitarbeit von Walter Hess. Rowohlt.

Contents Summary

7

Zusammenfassung

9

1 Introduction 1.1 Why are there so many kinds of species? . . . . . . . . . . . . . . 1.2 Ecological niches . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Spatial scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Modelling the spatial distribution of species richness . . . . . . . . 1.5 Hypotheses and aims . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Climate gradients and species richness of functional groups 1.5.2 Habitat diversity and area effect on species richness . . . . 1.5.3 Climate gradients and historical effects on species richness

. . . . . . . .

1 1 2 4 5 6 6 7 7

. . . . . . . . . . . . . .

9 10 12 12 12 12 13 13 14 14 16 16 16 16 19

3 Niches and Noise - Disentangling habitat diversity and area effect on species diversity 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Material and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Study area . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23 24 26 26

2 Modelling species richness of functional groups 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Material and Methods . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Study site . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Response variables . . . . . . . . . . . . . . . . . . . . 2.2.3 Predictors . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.4 Regression Modelling and Model testing . . . . . . . . 2.2.5 Clustering Species into Functional Groups . . . . . . . 2.2.6 Habitat Affinity of the Functional Groups . . . . . . . 2.2.7 Niche Width and Model Quality . . . . . . . . . . . . . 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Model Performance . . . . . . . . . . . . . . . . . . . . 2.3.2 Spatial Prediction . . . . . . . . . . . . . . . . . . . . . 2.3.3 Habitat Affinity, Realized Niches and Model Evaluation 2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

1

3.3 3.4

3.2.2 Species richness data . . . . . . . . . . . . . . . . . . . . 3.2.3 Habitat classification . . . . . . . . . . . . . . . . . . . . 3.2.4 Constructing sample based rarefaction curves . . . . . . 3.2.5 Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Grain size, niche width and habitat heterogeneity effect . 3.4.2 Confounding effects . . . . . . . . . . . . . . . . . . . . . 3.4.3 Implications for the spatial modelling of species richness

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

26 27 27 28 29 32 32 33 34

4 Is the maximum extent of the Quaternary ice shield necessary to understand the current spatial patterns of alpine plant species richness? 4.1 Inroduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Material and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Study area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Model parameters . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Classification of the environment . . . . . . . . . . . . . . . . 4.2.4 Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 The effects of habitat type, area and refugial distance . . . . . 4.4.2 The response shape to the refugial distance and hypotheses . . 4.4.3 Migration, establishment and equilibrium state . . . . . . . . 4.4.4 Prediction uncertainty . . . . . . . . . . . . . . . . . . . . . . 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35 38 39 39 39 40 41 42 47 47 48 48 49 49

5 Synthesis 5.1 Implications for conservation management . . . . . . . . . . . . . . . 5.2 Implications for ecological theory . . . . . . . . . . . . . . . . . . . . 5.3 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51 52 53 54

References

55

Appendix I

75

Appendix II

85

Acknowledgements

87

2

List of Figures 1.1

Conceptional niches in a 3-dimensional hyper-space . . . . . . . . . .

4

2.1

Spatial distribution of the total species richness and the species richness of the four plant groups with differing longevity. Grey regions represent areas with climatic parameters extending the range of values used for calibration. . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.1

3.2

4.1

4.2 4.3

5.1

Three types of sample rarefaction curves are compared on three different grain scales, while the extent was held constant, including the area of Switzerland. . . . . . . . . . . . . . . . . . . . . . . . . . . . . The habitat effect calculated as the proportion of the difference between the species number sampled within dissimilar and similar habitats to the number of species found within dissimilar habitats. (a) Habitat effect calculated from the data directly. (b) Habitat effect calculated from curves fitted to an Arrhenius power law function. . . Observed (a) and predicted (b) species richness 1000 m.a.s.l. Green areas represent potential glacial refugia. While the red squares represent areas of high prediction precision (lower 5% of the prediction uncertainty), the yellow squares represent areas of low prediction precision (upper 5% of the prediction uncertainty). . . . . . . . . . . . . The predicted species number explains 85% of the observed species number (r2 = 0.85). . . . . . . . . . . . . . . . . . . . . . . . . . . . . Posterior probability of the distance coefficient. a) The probability for the hotspot hypothesis being true accounted for 15% when the habitats were defined by 3 classes per predictor only. b) The probability for the hotspot hypothesis being true accounted for 99% when habitats were defined by 5 classes per predictor. . . . . . . . . . . . . . . . . . Spatial distribution of the species richness of the 40 functional groups. Grey regions represent areas with climatic parameters extending the range of values used for calibration. . . . . . . . . . . . . . . . . . . .

30

31

44 45

46

76

3

4

List of Tables 1.1

Overview of general hypotheses of species richness . . . . . . . . . . .

2.1

Plant traits used for clustering the perennial herbs into functional groups. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Description of the habitat types used to test group affinities. . . . . . 15 Mean absolute error (MAE) and explained deviance (D2) of the models. 18

2.2 2.3 3.1

4.1

5.1

5.2

5.3

3

Comparison of the different types of fitted sample rarefaction curves. Sample rarefaction curves of similar habitat types and dissimilar habitat types respectively, were compared to randomly aggregated rarefaction curves (which were considered as null models). . . . . . . . . . .

29

Model output of selected models. Mean b1: mean of the posterior distribution of the distance effect. Mean b2: mean of the posterior distribution of the coefficient of the quadratic term of the distance. CI: credible interval of the estimated coefficients. MAE: mean absolute error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

Environmental predictors. After collinearity diagnostics only 25 (printed in italics) out of the 95 predictors remained for the model calibration. The indices [i] range from 1 to 12 and refer to the monthly mean predictor values. The months that remained after collinearity analysis and significantly improved the models, are indicated with their respective number in brackets (e.g. [4] for April). . . . . . . . . . . . 81 Plant traits (rosette type, shoot metamorphosis, lateral growth form, leaf anatomy) and subgroups clustered thereof are listed together with the cluster habitat affinity (where low p-values indicate that the null model is true). The abbreviations are defined in Table 2.1 of the main script. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Predictors and model quality for each plant group. Predictors labelled with an asterisk entered the model as linear and quadratic terms. MAE: mean absolute error. D2 : deviance explained. . . . . . . . . . . 83

5

6

Summary The earth’s species richness is the result of biological evolution over the last billions of years. Manifold processes interact together and influence the spatial distribution of species richness. Trying to answer the question of “why are there so many kinds of species”, researchers developed countless species richness hypotheses over the last two centuries. Different mechanisms influencing species richness patterns act on different spatial and temporal scales. As in natural systems space and diversity are correlated, it is difficult to disentangle the manifold factors, which are correlated either as a consequence of mechanistic relationships, or as a matter of stochasticity. In the present thesis we tested the following aspects of plant species diversity: 1. Climate gradients and species richness of functional groups (Chapter 2) 2. Habitat diversity and area effect on species richness (Chapter 3) 3. Climate gradients and historical effects on species richness (Chapter 4) In the first part (Chapter 2), it was hypothesised that accounting for the autecology of species improves the model quality of spatial predictions of species richness patterns. Therefore, two modelling approaches were compared: a direct versus a cumulative modelling approach, where the latter gives more weight to the ecology of functional species groups. In the direct modelling approach, species richness was predicted by a single model calibrated for all species. In the cumulative modelling approach, the species were partitioned into functional groups. Models were calibrated for each functional group separately. The estimated species richness of each group was cumulated to predict the total species richness. Climate and topographic gradients explained species richness by ca. 25%. However, depending on the functional group up to 67% of the variability in species richness could be explained by climate and topographic gradients. Even though both modeling approaches performed equally well on average, the models of the different functional groups highly varied in their quality and their spatial richness patterns. Part of this variability could be a result of the differing growth forms of plants belonging to different functional groups. The differing growth patterns could lead to a scale effect. In the second part (Chapter 3) the aim was to disentangle the habitat diversity and area effect on species richness. With increasing area, habitat diversity increases, and so does species richness. However, the species area curve is not only driven by habitat diversity, but also by immigration and extinction as well as other stochastic

7

processes. In order to distinguish the relative importance of habitat heterogeneity and the area effect on total species richness, three different types of sample rarefaction curves were generated: (1) A randomly aggregated rarefaction curve, (2) a rarefaction curve where areas of similar habitat types were aggregated and (3) a rarefaction curve, where areas of dissimilar habitat types were aggregated. The classification of the habitat types was based on three environmental variables. In order to account for scale effects, these three types of curves were produced separately for three floristic surveys of different grain sizes. All of them had the spatial extent of Switzerland. The relative contribution of the habitat heterogeneity was then estimated from the difference of the dissimilar and similar rarefaction curve. The analyses showed that habitat heterogeneity contributes nearly 40% to species richness. The remaining 60% are likely related to effects caused by the area and by stochasticity. However, the habitat heterogeneity effect varied with grain size, contributing only 20% when a sample unit of 1.25ha was used instead of 10m2 or 28m2 , thus the detection of the habitat heterogeneity effect depends on the grain size. From the results, the prediction of biodiversity from habitat heterogeneity is expected to be most precise at medium grain sizes. While small grain sizes are affected by stochastic noise, large grain sizes lack in variability, as the heterogeneity is becoming uniform at larger scales of sampling units. In the third part (Chapter 4) the relation of climate gradients and historical effects on species richness was analysed. Using high mountain species data compiled over the Alpine system, the following hypotheses were tested: 1) Putative past glacial refugia still represent hotspots of plant species richness, presumably as speciation and diversification processes had a longer history in these areas during the Quaternary ice age. (2) During the Quaternary ice age, species differentiated in the isolated peripheral refugia. After the retreat of the ice shield, these species migrated and mixed in suture zones, which are characterized by high current species richness. (3) The alpine flora has already reached its equilibrium state on a large scale. Hence, the spatial pattern of species richness is no longer affected by historical distribution patterns. In order to test these hypotheses, the alpine environment was classified in areas of similar climatic and topographic conditions. The spatially replicated species richness observations within similar environments were used to estimate effects of the relative location to refugia on species richness. A Bayesian hierarchical mixed model of the type Poisson regression was used to test the species richness hypotheses related to the distance to refugia. Even though the distance had a weak effect on species richness, there is evidence for a negative correlation of species richness with distance to refugia, when accounted for similar environmental conditions. It might be that the hotspot hypothesis is still true but starts to be overridden by post glacial expansion of species adapted to a wide range of climatic condition.

8

Zusammenfassung Die heutige Artenvielfalt hat sich u ¨ber Milliarden von Jahren entwickelt. Aus einem komplexen Zusammenspiel von biotischen und abiotischen Faktoren entstehen r¨aumliche Muster von Artenvielfalt. Seit mehr als 2 Jahrhunderten versuchen Forscher immer wieder eine Erkl¨arung f¨ ur die Muster der Artenvielfalt zu finden. Zur Beantwortung der Frage weshalb es so viele Arten gibt wurden unz¨ahlige Hypothesen zur Erkl¨arung der Artenvielfalt entwickelt. unterschiedliche Erkl¨arungsmechanismen wirken auf ganz unterschiedlichen r¨aumlichen Skalen und stehen auch in Wechselbeziehungen untereinander. Gerade weil in nat¨ urlichen Systemen Vielfalt und Fl¨ache korreliert sind, ist es schwierig, die einzelnen Mechanismen isoliert zu betrachten. In dieser Arbeit wurden folgende Aspekte der Artenvielfalt von Pflanzen n¨aher betrachtet: 1. Klimagradienten und Artenvielfalt funktioneller Gruppen (Kapitel 2) 2. Habitatdiversit¨at und Fl¨acheneffekt auf die Artenvielfalt (Kapitel 3) 3. Klimagradienten und historische Einfl¨ usse auf die Artenvielfalt (Kapitel 4) Im ersten Teil (Kapitel 2) wurde angenommen, dass Modelle zur Vorhersage der Artenvielfalt verbessert werden k¨onnen, wenn man die Aut¨okologie der Pflanzenarten im Modellierprozess ber¨ ucksichtigt. um diese Hypothese zu testen, wurden zwei Modelliertypen verglichen. In einem direkten Modell wurde die Artenvielfalt mit einem Modell, welches f¨ ur alle Arten kalibriert wurde, vorhergesagt. In einem kumulativen Modell wurde die Artenvielfalt f¨ ur funktionelle Gruppen einzeln modelliert. Die gesch¨atzten Artenzahlen der einzelnen Modelle wurden anschliessend zusammengez¨ahlt um die Gesamtartenzahl vorherzusagen. Es zeigte sich, dass klimatischeund topographische Variablen im Durchschnitt ca. 25% der Artenvielfalt erkl¨aren k¨onnen. Je nach funktioneller Gruppe, konnten diese umweltvariablen jedoch beinahe 70% der Variabilit¨at der Artenzahl erkl¨aren. Die Modellqualit¨at der verschiedenen funktionellen Gruppen variierte betr¨achtlich. In der Vorhersage der Gesamtartenvielfalt erlangten die beiden Modelltypen jedoch dieselbe Pr¨azision. unterschiedliche Wachstumsformen von Pflanzen unterschiedlicher Gruppen k¨onnten einen Skaleneffekt in die Modelle eingebracht haben, der die unterschiedliche Modellqualit¨at erkl¨aren k¨onnte. Im zweiten Teil (Kapitel 3) war das Ziel, die korrelierten Effekte der Habitatvielfalt und der Fl¨achengr¨osse auf die Artenvielfalt voneinander zu trennen. Mit steigender Fl¨achengr¨osse steigt nicht nur die Artenvielfalt sondern auch die Habitatvielfalt,

9

die ihrerseits wiederum positiv mit der Artenvielfalt korreliert ist. Die Form von Artenarealkurven ist jedoch nicht nur von der Habitatvielfalt abh¨angig, sondern auch von Zufallsprozessen, wie auch Einwanderungs- und Sterberaten, welche ihrerseits von der Fl¨achengr¨osse abh¨angen. um die relative Einfl¨ usse der Habitatvielfalt und der Fl¨achengr¨osse auf die Artenvielfalt abzusch¨atzen, wurden 3 Typen von Artenarealkurven generiert: (1) Artenarealkurven, die sich aus zuf¨allig gew¨ahlten Fl¨achen zusammensetzten, (2) Artenarealkurven bestehend aus Fl¨achen von a¨hnlichem Habitattyp, und (3) Artenarealkurven, die sich aus Fl¨achen unterschiedlicher Habitattypen zusammensetzten. Die Klassifikation der Fl¨achen in Habitattype, basierte auf klimatischen und topographischen Variablen. um allf¨allige Skaleneffekte aufzudecken, wurden diese 3 Kurventypen separat mittels drei verschiedenen Datens¨atzen generiert. Die Vegetationsaufnahmen aller drei Aufnahmen u ¨berstreckten sich u ¨ber die ganze Schweiz, unterschieden sich jedoch in der Fl¨achengr¨osse, in der die Aufnahmen gemacht wurden. Der relative Beitrag der Habitatvielfalt zur Artenvielfalt wurde an Hand der Differenz der beiden Kurven mit a¨hnlichen und un¨ahnlichen Habitattypen gesch¨atzt. Den Berechnungen zu Folge kann die Habitatvielfalt bis zu 40% der Artenvielfalt erkl¨aren. Mehr als die H¨alfte der Variabilit¨at der Artenvielfalt wird jedoch sehr wahrscheinlich durch andere, Zufalls- und fl¨achenabh¨angige Prozesse verursacht. Je nach Gr¨osse der Aufnahmefl¨achen variiert jedoch der Effekt der Habitatvielfalt. Bei Aufnahmefl¨achen von 1.25 ha (im Vergleich zu Aufnahmen, die auf Fl¨achen von 10m2 bzw. 28m2 gemacht wurden), erkl¨art die Habitatvielfalt nur noch 20% der Artenvielfalt. Die Resultate lassen vermuten, dass Artenzahlen basierend auf Vegetationsaufnahmen auf Fl¨achen mittlerer Gr¨osse die pr¨azisesten Modelle liefern. Wenn Vegetationsaufnahmen auf zu kleinen Fl¨achen gemacht werden, werden die Modelle durch Zufallseffekte verschlechtert. Bei zu grossen Aufnahmefl¨achen wird die Heterogenit¨at uniform, weil eine Fl¨ache gewisser Gr¨osse potentiell alle regional vorhandene Arten enthalten kann. Im dritten Teil (Kapitel 4) wurden historische Effekte untersucht. Mittels Daten u ¨ber die r¨aumliche Verteilung von alpinen Pflanzen u ¨ber den ganzen Alpenbogen wurden folgende Hypothesen getestet: (1) Auf Fl¨achen, die w¨ahrend der quart¨aren Eiszeit eisfrei blieben, stand der Pflanzenwelt mehr Zeit zur Spezialisierung und Diversifizierung zur Verf¨ ugung. Deshalb wird vermutet, dass Regionen ehemaliger glazialen Refugien bis heute u ¨ber eine reichere Flora verf¨ ugen. (2) In den w¨ahrend der letzen Eiszeit isolierten Refugien entwickelten sich voneinander unabh¨angig neue Arten. Nach R¨ uckzug der Eisdecke wanderten diese isoliert gebildeten Arten ins Zentrum der Alpen, wo sie sich zu vermischen begannen. Auf Grund dieser Vorg¨ange sollten die Zentralalpen eine h¨ohere Artenvielfalt aufweisen als die Randalpen. (3) u ¨ber die Skala der Alpen gesehen, haben die Pflanzen einen Gleichgewichtszustand mit der klimatischen umgebung erreicht. Deshalb reichen klimatische und topographische Variablen aus, um die r¨aumlichen Muster der Artenvielfalt zu erkl¨aren. um diese Hypothesen zu u ¨berpr¨ ufen, wurde der Alpenraum in klimatisch und topographisch

10

a¨hnliche Klassen eingeteilt. Mittels Artenzahlen aus Gebieten mit a¨hnlichen umweltbedingungen, aber unterschiedlicher geographischer Lage, wurde der Effekt der Distanz zu fr¨ uheren glazialen Refugien getestet. Es zeigte sich, dass die Distanz zu den Refugien einen schwachen, aber doch signifikanten Effekt auf die Artenvielfalt haben kann. Die Artenvielfalt aus Gebieten mit a¨hnlichen Umweltbedingungen scheint mit zunehmender Distanz abzunehmen. Regionen fr¨ uherer glazialer Refugien sind also bis heute durch eine hohe Artenvielfalt gekennzeichnet. Eury¨oke Pflanzenarten k¨onnten allerdings durch postglaziale Wanderungsprozesse das Muster abgeschw¨acht haben.

11

12

1 Introduction The geochemical and biological evolution over the last 4.7 billion years bore manifold life forms whose interaction together and with their environment result in a uniquely versatile place. Our planet harbours at least 1.4 Mio. animal species and about 400’000 plant species are described (Flindt, 2000). Why are there so many kinds of species? How are they distributed in space and time? What drives their distribution? How do they interact and coexist with each other? Will our ecosystems also be stable with less species? All these questions have fascinated biogeographers and ecologists for more than two centuries (Forster 1778, von Humboldt 1807, Hutchinson 1959). However, even for the simpler question of why are there so many kinds of plant species, no generally accepted unified theory exists. Understanding why some communities are richer in species than others, is not only fascinating for the postulation and testing of ecological theories, but has direct practical consequences in conservation planning and management. Once the underlying mechanisms causing species diversity are known, species richness can be modeled and predicted in space and time. In a further step, local and regional extinction risks of single species might be foreseen and avoided. In regard to the recently anthropogenic accelerated, extinction rate (exceeding the extinction in fossil records 50- to 500 fold) understanding the functioning of biodiversity has become more important than ever (IUCN, 2004).

1.1 Why are there so many kinds of species? When Hutchinson collected Corixidae (water boatmen) species in ponds near the last resting-place of Santa Rosalia in Palermo in Sicily, he observed that the bigger species C. punctata was breading earlier in the season than the smaller C. affinis (Hutchinson, 1959). This type of observation was made many times before by other ecologists. However his observation led him to ask, why the larger species should breed first, and then to the more general question as to why there should be two and not 20 or 200 species of the genus in the pond (Hutchinson 1957). Even though we will likely never find out, why there are so many kinds of species, not fewer of more, it is strongly assumed, that the limitation of available energy and resources prevent a small number of species to monopolize them (see Rohde 1992). Phenological differences between species that result in a temporal division of resources, differences in growth form, dispersal and life histories are all mechanisms which increase the amount of potentially suitable resources, necessary for life (Tokeshi, 1999). In other words, species divide up available resources. Each species can exploit some resources

1

well enough to monopolize them, but trade-offs prevent them to monopolize all resources. According to Gause’s principle no two species can permanently occupy the same niche (Gause, 1936). Therefore resource partitioning, niche differentiation and speciation are the key mechansims for the coexistence of multiple species (Wilson, 1990, Huston 1994). Hence they represent important aspects in many hypotheses posted during the last decades to explain species richness (for reviews see e.g. Brown 1981, Huston 1994, Ricklefs & Schluter 1993). A short overview of species richness hypotheses is given in Table 1.1. The list is by no means complete (see Palmer 1994, for a list of 120 named species richness hypotheses) and the proposed hypotheses are certainly not mutually exclusive. Moreover, most of the hypotheses can be ascribed to one common underlying ecological concept - the niche theory.

1.2 Ecological niches The concept of ecological niches was already formulated nearly a century ago. Grinnell was one of the first who described the niche as a place in an environment which is occupied by a specific species (Grinnell 1914, 1917). He defined the niche as the combination of necessary conditions for a species’ existence, including physiological tolerances, morphological limitations, feeding habitats, and interactions with other members of the community (see Chase and Leibold, 2003 for a detailed review). His niche definition arose from the logic of the competitive exclusion principle, an idea which was reflected later by Gause (1936), and further developed and expressed by Hardin (1960). Around the same time as Grinnell, Elton developed a similar niche theory. However, he accounted for the species’ functional role within the food chain and its impact on the environment (Elton 1927). While Grinnell’s niche concept focused on the effects of the environment on the species, Elton stressed the effect of a species on the environment. As the niche definition by Grinnell mainly considers the abiotic environment, it is closer related to what is known today as the “fundamental niche”. By including biotic interactions, Elton’s niche definition is closer to the “realised niche” concept (Ellenberg 1953, Hutchinson, 1957). In fact already in 1917 Tansley explicitely differentiated between conditions in which a species could theoretically exist (fundamental niche) and the actual conditions in which a species indeed does exist (realised niche; Tansley, 1917). Hutchinson quantified the niche concept and defined the term niche as the range of environmental and biotic conditions within which its population can persist without imigration (Hutchinson 1957). He thus defined the niche as a region of an n-dimensional “hyper-space” (Hutchinson 1944). An illustrative example for an 3-dimensional “hyper-space”, also called “hyper-volume” (Hutchinson 1957) is given in Fig. 1. Given this framework of niche theory, the energy hypothesis states simply than that with increasing energy, more potential niches are available and therefore species richness is supposed to increase. Similarly, the productivity richness hypothesis can

2

Table 1.1: Overview of general hypotheses of species richness

Hypothesis

Argument

Authors

Energy

Species richness increases with increasing available energy (i.e. the coexistence of species is limited by the energy supply).

Brown 1981, Wright 1983

Productivity

Species richness reaches its optimum at intermediate productivity. Species richness is limited by nutrients and light at the low and high end of a nutrient gradient respectively.

Tilman 1982

Climate

Species richness varies with temperature and or water availability

Currie & Paquin 1987

Physiological tolerance

Species richness in a particular area is limited by the number of species that can tolerate the local conditions

Woodward 1987, Woodward 1990, Root 1988

Gradual climate change

Species richness is higher in a system where climate is gradually changing, as the changing environment prevents the achievement of the equilibrium state.

Connell 1978, Huston 1979, Warner & Chesson 1985, Grubb 1988

Niche diversification

Species richness increases with the number of different available niches

Connell 1978, Aarssen 1983

Predator Pressure

More intense predation reduces competition and thus permits greater niche overlap. Hence by enhancing the coexistence of species predation promotes higher species richness.

Gillett 1962, Connell 1970, Shmida & Ellner 1984

Intermediate disturbance

Species richness reaches its optimum at intermediate disturbance. High disturbance excludes most of the species, while low disturbance gives way to the equilibrium state.

Grubb 1977, Connell 1978

Area

Larger areas harbor more species

MacArthur & Wilson 1967

Spatial mass effect

The local species richness increases with increasing chance of immigration from nearby habitats.

van Steenis 1972, Webb & Hopkins 1984

Rapoport’s rule

Species richness increases toward the equator as a consequence of decreasing mean geographical range sizes of species at low latitude.

Stevens 1989

Mid-domain effect

Species richness peaks in the center of a study region or domain simply as a consequence of placing data sets of species of varying range size (whether empirically or theoretically generted ) randomly within a bounded domain.

Colwell & Lees 2000

Speciation rate

Speciation rate varies with climate, due either to faster evolutionary rates or stronger biotic interactions increasing the opportunity for evolutionary diversification in some regions.

Rohde 1992, Allen et al. 2002

Historical factors

Glaciation effects, dispersal, speciation rates

Ricklefs and Schluter 1993

3

Figure 1.1: Conceptional niches in a 3-dimensional hyper-space be seen as the result of species occupying niches which are either favourable in the abiotic niche space (high nutrient availability), but less favourable in the biotic niche space (high competition) or vice versa. The hump shaped species richness response to productivity could be seen as an energy constrain at the stressful lower and upper end of the nutrient and competition gradients (niche axes) respectively (Brown 1996). At least half of the described hypotheses listed in Table 1.1 can be reduced to the niche theory. While some of the hypotheses give higher weight to abiotic niche axes, other focus more on biotic niche axes. However different hypotheses act over different spatial and temporal scales.

1.3 Spatial scale Biological diversity increases with the area sampled (Arrhenius 1921, Gleason 1922). Over a small scale, species richness increases with the total abundances of individuals. Therefore species richness on the local scale (including an area of ca. 10 m2 , see Weber et al. 2004) is mainly driven by the growth form and growing density of single species (Crawley & Harral, 2001) jointly with area. Over the scale of the landscape (including an area of ca. 106 m2 ), species richness depends on the habitat diversity (Williams 1964, Abele 1974, Deshaye & Morisset 1988, Kohn & Walsh 1994, Triantis et al. 2003). Extinction of local populations, colonisation rates and whole metapopulation dynamics act on the species richness pattern of the regional scale (containing an area of ca. 1010 m2 ; MacArthur & Wilson 1963, 1967). At the macro-

4

scale, when entire provinces are considered, the degree of evolutionary independence between species pools (addition of new floras or faunas) dictates the species richness (Preston 1960, Rosenzweig 1995). While climatic conditions affect the local species composition, it has less influence on the local (plot size) species richness. However on a meso- up to the macroscale, species richness has been shown to correlate with climatic variables (Currie & Paquin 1987, O’Brien 1993, Heikkinen & Birks 1996, Lobo et al. 2001). Thus the conspicuous latitudinal gradient in species richness is mainly driven by climate either directly (Francis & Currie 2003) or indirectly via other mechanisms as e.g. the Red Queen hypotheses (Anonymous correspondent 1973), which finally are also driven by climate.

1.4 Modelling the spatial distribution of species richness A wide range of statistical models are currently used to predict the spatial distribution of plant species (e.g. Hill 1991, Lenihan & Neilson 1993, Huntley et al. 1995, Franklin 1998, Guisan et al. 2002, Gelfand et al. 2005), communities (e.g. Brzeziecki et al. 1993, Zimmermann & Kienast 1999), functional types (e.g. Box 1996) or biodiversity (e.g. Heikkinen 1996, Wohlgemuth 1998). Among different techniques such as ordination and classification methods, neural networks and Bayesian statistics, generalised linear models (GLMs) and generalised additive models (GAMs) have become the most popular models in predicting species distributions (Elith et al., 2007). Models predicting the spatial distribution of single species directly rely on the niche theory and gradient analysis (Austin & Austin 1980, Austin, 2002). They assume that each species has its optimal growth along direct and indirect gradients. Bioclimatic predictors such as e.g. temperature, degree days, amount of precipitation, etc. are assumed to represent these gradients. With GIS layers of the bioclimatic parameters, the spatial distribution of the realised niches of species can then be predicted by using models calibrated with presence absence data of the single species. The ecological concept of modelling species richness is more complex. As described above, there are a vast number of species richness hypotheses. While some drivers of species richness can easily be represented by bioclimatic and topographic variables (such as e.g. degree days or solar radiation as a surrogate of energy), other hypotheses are more difficult to test. For example the habitat hypothesis implicitly relies on the niche theory. It assumes that with increasing variability of habitats, the number of niches increases. However, as we do not (and probably never will) know the niche of every single species it is difficult to predict the species richness by the number of different niches. Only approximations such as arbitrarily classifying the environment in hyper-cubes (as illustrated in Fig. 1.1) can be made. A further caveat in modelling species richness based on environmental conditions

5

relies in the fact that two different habitat types might be characterised by two completely different vegetation types, but harbour the same number of species. This may lead to imprecise models, as only the species composition but not the species number is directly related to the environmental conditions. Finally as the countless mechanisms driving species richness might interfere with each other and act on different scales (Pickett et al. 1994), more complex models, might be necessary. Due to progresses in computational statistics, Bayesian hierarchical models which accommodate multiple stochastic elements (Carlin & Louis 2000, Gelman et al. 2004, Clark & Gelfand 2006, McCarthy 2007) might provide a new framework within which species richness hypotheses can be tested and predictions can be done in a more realistic way (Clark et al. 2001, Clark 2005).

1.5 Hypotheses and aims In natural systems, space and diversity are correlated (Currie, 2007; Bahn & McGill 2007) . As trivial this issue seems, as difficult it actually is to disentangle the manifold factors, which are correlated either as a consequence of mechanistic relationships, or as a matter of stochasticity. The overall goal of the thesis therefore was to test basic assumptions of species richness hypotheses and their potential interference at different scales.

1.5.1 Climate gradients and species richness of functional groups In the second Chapter, the relation between climate, and the autecology of functional groups was tested. Species are known to be distributed along direct and indirect environmental gradients. With their specific morphological and physiological adaptations to their environment, they create their own niches, and hence plant types segregate along these axes. Each plant type has its optimum growth at a specific location along the axes. Theoretically, the species richness along a specific axis could remain constant. However the species composition along this axis is likely to change gradually. For example only frost tolerant species are found at the lower end of a temperature gradient, while at the upper end, rather dry adapted species will be found. A classical overall species richness model may therefore exclude the temperature as a predictor, even though the presence or absence of single species is highly sensitive to temperature. Therefore the following hypothesis was tested: • A model approach which takes the autecology of the species into account will result in improved species richness predictions. In order to test the hypothesis we compared two modelling approaches. (1) In the direct approach, species richness is predicted by one single model calibrated for all

6

species. (2) In the cumulative modeling approach, species are partitioned into functional groups. Models were calibrated for each functional group separately and the estimated species richness of each group was then cumulated to predict the total species richness.

1.5.2 Habitat diversity and area effect on species richness In the third Chapter, the relation of the ‘habitat diversity - species richness’ and ‘area - species richness’ hypotheses was tested. Therefore three types of sample rarefaction curves were generated. (1) A randomly aggregated rarefaction curve (serving as a null model); (2) a rarefaction curve where areas of similar habitat types are aggregated (similar curve); and (3) a rarefaction curve, where areas of dissimilar habitat types are aggregated (dissimilar). These three curve types were used to test the following hypothesis: • As a consequence of the habitat diversity hypothesis, the highest species accumulation rate is expected when areas of dissimilar habitat types are aggregated. • Assuming that rarefaction curves constructed with areas of similar habitat types catch the area effect only, their species accumulation rates and saturation levels are expected to be lower compared to randomly generated rarefaction curves. The different curve types were used to disentangle the correlated effects of area and habitat diversity (and their interaction) on species richness. The proportion of the difference of the saturation levels of the dissimilar curve and the similar curve to the total species richness (dissimilar curve) was assumed to indicate the contribution of habitat diversity to the total species richness. In order to test for potential scale effects, the effect of habitat heterogeneity was estimated for three different grain sizes.

1.5.3 Climate gradients and historical effects on species richness In the fourth Chapter, the relationship between the climate gradient hypothesis and historical effects was tested. Commonly it is assumed, that alpine species survived the Quaternary ice age either in peripheral refugia (the tabula rasa hypothesis) or on ice-free areas above the ice-shield (the nunatak hypothesis). In both cases, reimmigration was necessary for colonisation of the alpine area after the retreating ice shield. Whether this migration process is completed and species richness on a meso-scale has reached an equilibrium state with the current climate remains an open question. Therefore the following hypotheses were tested:

7

H0 : The alpine flora has already reached its equilibrium state on a large scale. Hence, the spatial pattern of species richness is no longer affected by historical distribution patterns. H1 : Putative past glacial refugia still represent hotspots of plant species richness, presumably as speciation and diversification processes had a longer history in these areas during the Quaternary ice age. H2 : During the Quaternary ice age, species differentiated in the isolated peripheral refugia. After the retreat of the ice shield, these species migrated and mixed in suture zones, which are characterized by higher current species richness. While the hypotheses of climate, ecological niches, habitat diversity and area were tested for the extent of Switzerland only, the interrelationship of climate and history was tested over the extent of the whole alpine mountain system.

8

2 Modelling species richness of functional groups 1

K. Steinmann, 2 H. P. Linder, 1 N.E. Zimmermann

1

Land Use Dynamics, Swiss Federal Research Institute WSL, Zuercherstrasse 111, CH-8903 Birmensdorf, Switzerland 2

Institute for Systematic Botany, Zollierstrasse 107, CH-8008 Zuerich, Switzerland

Abstract: Conservation biologists increasingly rely on empirical biodiversity distribution models for decision-making. Therefore, good model quality is required. While statistical techniques were optimized to improve model quality, less focus has been given to the question how the autecology of single species might affect the model quality. In the present study, two modelling approaches are compared: a direct versus a cumulative modelling approach, where the latter gives more weight to the ecology of functional species groups. In the direct modelling approach, species richness is predicted by one single model calibrated for all species. In the cumulative modelling approach, the species are partitioned into functional groups. Models were calibrated for each functional group separately. The estimated species richness of each group was cumulated to predict the total species richness. As we hypothesized that the model quality depends on the ecology of single species, we expected the cumulative modelling approach to predict species richness more accurately. In general the predictors explained plant species richness by ca. 25 %. However, depending on the functional group the deviance explained varied from 3 to 67%. While both modelling approaches performed equally well on average, the models of the different functional groups highly varied in their quality and their spatial richness pattern. This variability helps to improve our understanding on how functional groups respond to ecological gradients.

9

Keywords: clustering species, generalized linear models, meso-scale, mor- phological traits, niche width

2.1 Introduction Predictive modelling of species richness relies on niche theory and gradient analysis (Austin, 2002). The fundamental niche is the environmental range a species population can persistently occupy without immigration (Hutchinson, 1957). The realized niche of a species additionally includes biotic interaction and competitive exclusion (Hutchinson, 1957). The response of species along environmental gradients is often assumed to follow a Gaussian shape with identical width and height but individually distributed optima (Whittaker, 1956; Gauch and Whittaker, 1972; ter Braak, 1985; but see Abrams, 1995 and Austin, 2002). However, Austin and Smith (1989) pointed out that physiological processes and interactions between species may lead to skewed, bimodal or more complex response curves. As a species’ niche response is the result of multiple reactions to varying environmental gradients, it is unlikely to exclusively take a multivariate Gaussian shape (Minchin, 1989; Austin et al., 1994). Woodward and Kelly (1997) showed that already a simple temperature gradient results in highly skewed maximum photosynthetic rates in different biomes. It is now accepted that the physiological responses of species to environmental factors are skewed (Oksanen and Minchin, 2002). In addition, species with optima close to the end of a gradient seem to have narrower niche widths, compared to species with growth optima near the mid-point of the same gradient (Thuiller et al., 2004). Hence, according to its physiological and morphological traits, each species has its own niche, with its characteristic shape and optimum along a given environmental gradient. This expression of the realized niche is often calibrated in statistical models for predicting species spatial occurrence or richness patterns (Guisan and Zimmermann, 2000). While there is a wide range of studies concerning the optimization of different statistical techniques and comparing their quality (see Guisan et al., 2002; Thuiller et al., 2003; Rushton et al., 2004; Austin et al. 2006), comparably few studies attempt to optimize the selection of model predictors (but see Mac Nally; 2000, Austin et al. 2006). The question whether or to what degree the autecology of single species might affect the model quality, has been neglected even more. In the present study, we assume that the width, position and shape of realized niches of single species will influence the predictive power of species richness models. Conceptually, species richness can be modelled in two different ways: (1) by direct modelling of species richness from a single model (direct approach); (2) indirect prediction of species richness by cumulating the outputs from individual presence and absence modelling of each single species (complete cumulative approach). The first method is straightforward, but may show deficiencies as pooling of all species with

10

their different realized niches may result in information loss (Minchin, 1989). For example two environmentally completely different sites may share the same number of species, yet the species composition at these two sites usually differs completely (Sarr et al., 2005). While in one habitat species may primarily be limited by temperature, in the other habitat species may be constrained by precipitation (Broennimann et al., 2006). An overall richness model may therefore include (or exclude) both parameters, hence the predicted species richness response may be fuzzy. In contrast, predicting the spatial distribution of each species separately might result in more precise cumulative predictions of species richness. However, depending on the targeted ecosystem and the spatial scale, the total number of species can be huge (> 1000 species), while presence data of single species may be scarce (and may thus hamper the calibration of accurate models). Therefore, overlaying simulated presence absence maps of each single species may accumulate to large prediction uncertainties for species richness, and thus lower the strength of the full cumulative model. A compromise between the two described approaches consists in dividing all species into functional groups of similar eco-physiology, predicting the species richness within each group separately and cumulating richness values group-wise to total species richness (partially cumulative approach).Thus both penalties - i.e information loss and weak model calibration with an inflated confidence interval - may be reduced. The challenge of the intermediate approach is to find appropriate criteria to define functional groups representing similar eco-physiological behaviour (Smith, 1997). Adaptations to specific environmental constraints require specific traits, which are reflected in morphological, physiological or life history characters of each single species (Parkhurst and Loucks, 1972). Therefore morphological as well as life history traits seem suitable to determine the niche of a specific plant species (Lavergne et al., 2004) and are often used to divide plant species into functional groups. The principal goal of our study was to evaluate whether the modelling of species richness of functional groups is more accurate than modelling the total species richness. The second goal was to evaluate whether the result differs depending on the level of functional differentiation. Therefore, functional groups were built on a coarse level (4 life form groups) and on a finer level (40 groups based on morphological traits). As more than 70% of the species used are perennial herbs, only the latter group was partitioned into smaller groups using a clustering method applied on morphological characters. In order to evaluate the ecological relevance of the clusters, the habitat affinity of the species of each cluster was tested. By this, we aimed at testing the following two hypotheses: (1) The model quality for species richness can be improved (over direct modelling of species richness) if species of similar autecology (functional groups) are modelled separately and aggregated a posteriori. (2) The model quality for the species richness prediction of a functional group depends on its average niche width. Higher model quality is expected for groups with narrow niches.

11

2.2 Material and Methods 2.2.1 Study site The study was carried out in Switzerland ( 6-10.5◦ N and 46◦ -47.5◦ E). The climate is generally humid temperate with rather mild winters and moderately warm summers. Annual mean temperature ranges from ca. 1.5 C◦ (in high alpine valleys) to 11.5 C◦ (in the southern part of Switzerland). The mean annual precipitation sums vary between 600 mm and 1700 mm.

2.2.2 Response variables The species richness data are derived from vegetation recordings at 784 sites in Switzerland. The data originate from the biodiversity monitoring program of the Swiss federal office of environment (Plattner et al., 2004). Each record contains all vascular plant species growing in a circle plot of 10 m2 (r = 1.78 m). The sampling sites are organised in a regular grid with mesh size of 6x4 km. The initial point of the grid was generated randomly. Urban areas and inappropriate growing conditions such as lakes and glaciers were excluded. Therefore only 421 out of the 784 sites were kept for further analyses.

2.2.3 Predictors For the predictive modelling of species richness, we used climatic (12), topographic (5) and edaphic (4) parameter maps in combination with a coarse habitat classification (8 classes) as environmental predictors (see table S1 in the electronic appendix). Most of the climatic parameters were available as monthly variables. The climatic parameters were generated according to Zimmermann and Kienast (1999) using a digital elevation model (DEM) with a spatial resolution of 25 m and data from meteorological measurements. Climate variables represent monthly Normals of the period 1961-1990. Degree days of the growing season are based on a threshold of 0◦ C. Site water balance was calculated as the monthly difference of precipitation and potential evapotranspiration over the water holding capacity of the soil (see Maggini et al., 2006 for details). The topographic wetness index expresses the lateral water flow. It was calculated according to Beven and Kirkby (1979). The topographic position is a measure of convexity of the terrain and was calculated from the DEM directly. As a surrogate of potential stand productivity the two edaphic layers coarse fragment content and nutrient availability were included. Finally, the presence of limestone bedrock was derived from the digitally geotechnical map of Switzerland (1:200’000, de Quervain et al., 1963-1967). While all continuous data were available at a 25 m resolution, the categorical data had a resolution of 100 m only.

12

2.2.4 Regression Modelling and Model testing To avoid model calibration bias, the predictive parameters were first checked for collinearity. After discarding collinear parameters according to the method described in Belsley (1991), 25 predictors out of the original 95 were kept for modelling species richness pattern. The plot species richness showed a negative binomial distribution, thus we used generalized linear models with a negative binomial model family. In two cases the group of modelled species was too small (0 or 1 species per plot), thus a binomial model family was used. We employed a stepwise forward variable selection to build models. The number of parameters entering the final model was set to minimize the Bayesian information criterion (BIC), a comparably strict selection algorithm (Schwarz, 1978). Continuous predictors were allowed to enter the model as linear and quadratic terms. Model calibration was calculated from the deviance explained (D2 ) and model performance was evaluated from the mean absolute error (MAE) and from the prediction bias (mean of the error) originating from a 10fold-cross validation. GIS layers of each predictor were used for spatial predictions of species richness throughout Switzerland. Areas outside the calibrated predictor space were excluded from the predictions (safe predictions, Guisan et al., 2006). Maps were generated for all calibrated models. The two modelling approaches (partially cumulative vs. direct prediction) were compared by testing the MAE of both approaches for significant differences (paired Wilcoxon test). All statistical analyses were performed in the statistical software environment R (R Development Core Team, 2006).

2.2.5 Clustering Species into Functional Groups For the different models, the total list of plant species was divided into four functional groups (annuals & biannuals, perennial herbs, shrubs and trees). The information for this classification of life form was taken from Aeschimann et al. (2004). Perennial herb species were further classified into 40 subgroups according to a set of morphological characters. The characters included rosette type, lateral growth, shoot metamorphosis and leaf anatomy with 3, 4, 7 and 6 classes, respectively (see Table 2.1 for definitions). We used Klotz et al. (2002) for defining the morphological characters and we compiled the characters for the species from Hess et al. (1976) and Klotz et al. (2002). We applied a cluster analysis by partitioning around medoids (Kaufmann and Rousseeuw, 1990) to divide the perennial plant group into subgroups. We then used silhouette coefficients (sc: defined as the maximal average silhouette width for the entire data set) to judge the appropriateness of the cluster structure. A sc > 0.5 indicates a reasonable structure in the data (Kaufmann and Rousseeuw, 1990). Hence the number of morphological species groups (k) built was set so that k was minimized while sc was not allowed to drop below a value of 0.5. This procedure resulted in 40 subgroups.

13

2.2.6 Habitat Affinity of the Functional Groups To test whether clustering based on morphological characters resulted in ecologically meaningful groups, we applied a correlation analysis for habitat affinity. We assigned the species to nine habitat types (Table 2.2) according to Aeschimann et al. (2004). The different proportions of species per habitat type (considering all species) served as null model. For each cluster group, the species proportions per habitat type were compared against the null model using Spearman’s rank correlation coefficient. Low correlations (with associated high p-values) indicate that the two species lists are not equally distributed across the habitats. Thus, it indicates that the respective cluster group contains species, more associated to a specific habitat and hence the group reflects ecological idiosyncrasies.

2.2.7 Niche Width and Model Quality The variance of the normalized predictor values among all observations per group was used as a measurement of niche width. The niche width of a single cluster group was then defined as the mean variance along the predictor gradients of the observed presence data points of the respective group. Only predictors that remained in the respective models were used to calculate the average species niche width per group. As the mean absolute errors increased with the maximal observed species richness, they were expressed relative to the maximal observed richness number per group. Finally the niche widths were correlated with the relative MAEs using Pearson’s product-moment correlation.

14

Table 2.1: Plant traits used for clustering the perennial herbs into functional groups.

Rosette type

Lateral growth

Shoot metamorphosis

Leaf anatomy

fr: full rosette cushion b: bulbill hel: helomorph hr: half rosette standard n: standard shoot hyd: hydromorph nr: no rosette turf p: pleiocorm hyg: hygromorph tussock rh: rhizome m: mesomorph rht: rhizome and tiller building scl: scleromorph rp: rhizome pleiocorm suc: succulent rt: root tuber

Table 2.2: Description of the habitat types used to test group affinities.

Code Habitat type 1 2 3 4 5 6 7 8 9

Water courses Nitrophilous vegetation Screes, gravel and rocky area Creeks Swamps Meadows, pasture, turfs, snow beds Dwarf-shrub heathland, tall forb meadow Scrubland sensu lato Forest

15

2.3 Results 2.3.1 Model Performance Clustering the 789 perennial species from four morphological traits resulted in 40 distinct groups (table S2 in the electronic appendix). All 25 predictors (which remained after testing for collinearity) were significant in at least one of the 45 models (see table S3 in the electronic appendix). Direct prediction of the total species richness resulted in a D2 of 0.28 (2.3), whereas predicting functional group richness at the plot level revealed D2 values between 0.06 (bi-annual plants) and 0.68 (trees), with an average of 0.31. The model which predicted the total species richness directly showed a weak negative bias (-0.05, i.e. a marginal underestimation) and an MAE of 7.9 species. The accumulation of the modelled species richness from four functional groups resulted in a small positive bias (0.08) and a similar MAE of 7.97 species. The direct richness model of the perennial group showed a positive bias of 0.08 and an MAE of 7.02, whereas the cumulated predictions of perennial species richness from the 40 functional groups reduced the bias to 0.02 and the MAE to 6.69. All 40 perennial sub-group models were unbiased and had an MAE < 1.2 species. The comparison of the model residuals for total plot species richness revealed no significant difference between the predictions of the two modelling approaches (p = 0.32 for absolute residual values). Similarly, the comparison of the perennial model residuals did not reveal significant differences (p = 0.08 for absolute residuals). It indicates that the richness models from functional groups did not significantly increase prediction accuracies.

2.3.2 Spatial Prediction The sum of the predicted species richness of the four different life forms resulted in very similar spatial pattern compared to the directly modelled species richness (Fig. 2.1 a, b). Likewise we found a high similarity between the predicted perennial plant richness and the sum of all 40 richness models of perennial herbs clustered from functional traits (Fig. 2.1 c, d). Since roughly 75% of the recorded species were perennial herbs, it is not surprising that the richness pattern of this group shares a high similarity with the pattern of the total species richness. However, the spatial patterns of the different functional groups varied considerably (Fig. 2.1 e-h and A1 in the appendix 1).

2.3.3 Habitat Affinity, Realized Niches and Model Evaluation Ten out of the 40 perennial groups consisted of species clearly belonging to a specific habitat type (see Table S2 in the electronic appendix). However, groups that significantly reflect specific habitats did not necessarily result in better models (a t-test

16

comparing the explained deviances of the two groups was insignificant; p = 0.43). The realized niches, defined by the mean variance of the normalized predictors entering the model for a specific functional group, varied between 0.3 (Crassulaceae species, group c29) and 1.8 (species of the type Potentilla reptans, group c3). The model fit could not be explained by the average niche width per perennial group (r = 0.08, p = 0.63). Only a weak correlation was found between the number of plots containing species belonging to a specific functional group and the explained deviance of the respective model (r = 0.26, p = 0.1).

17

Table 2.3: Mean absolute error (MAE) and explained deviance (D2) of the models.

Group All species Trees Shrubs Perennial herbs (Bi-)annuals

MAE 7.9 0.75 1.44 7.02 0.73

Sum of 4 groups Sum of 40 groups

7.97 6.69

D2 0.28 0.68 0.17 0.34 0.06

Summary statistics of the 40 functional groups Mean Variance

18

0.36 0.06

0.23 0.01

2.4 Discussion Both the selected predictors and the model quality varied considerably between functional groups. In line with other studies (Sarr et al., 2005; Broennimann et al. 2006) the model for tree species richness performed best. This may be explained by the dominance of this life form, which is highly competitive after establishment. Thus we likely find trees more steadily compared to other life forms once the habitat is suitable. In contrast, short lived-species showed rather poor model fit, which was also found by Broennimann et al. (2006). Since (bi-)annual plants require open space to regenerate, they are found less often under dense vegetation, which reduces their detection probability and consequently their model quality. Such plants belong to the ”rural” plant type (according to Collins et al., 1993), a plant type that has been identified as being difficult to predict (Edwards et al., 2005; Zimmermann et al., 2007). Within the perennial species groups, wetland species were best predicted. Even though their habitat is restricted, they show a wide geographic distribution and high local abundance. Thus, they belong to predictable habitat specialists (Rabinowitz et al., 1986). Yet we did not find a general agreement between niche width and model fit. We believe that this is the main reason why the two different model approaches tested here did not yield differing model accuracies for the total species richness, and only marginal differences for the perennial species richness models. The influence of both the ecological information loss (mixture of niches) and the cumulated uncertainty (many group models accumulated) might cancel each other out when modelling species richness at a plot level where randomness is high, a conclusion that was also found by Guisan and Theurillat (2000) when comparing individual species distribution vs. species richness models. We did not test a complete cumulative model because 94% of our species occurred in fewer than 50 plots, which is below a threshold that allows sound model calibration. Effects of geographic range size and varying population densities of single species may have strongly influenced the different models (Venier et al., 1999; Stockwell and Peterson, 2002; McPherson et al., 2004). By grouping the perennial herbs according to their lateral growth and shoot metamorphosis, we differentiated between clonal species and species reproducing sexually. As shown by Kolb et al. (2006) the local abundance of a species is related to such life histories. The lateral spread of clonal species results in a high ramet density within a given habitat patch. Therefore groups containing clonal species are expected to result in higher model accuracy, compared to groups with similar habitat width and lower local abundance. Hence, independently of the according niche width, we generated models of different quality. In summary, we found no clear improvements for species richness modelling when splitting species into functional groups. The direct modelling of plot species richness performed equally well. However, in line with Bruun et al. (2006) our analysis showed a clearly differing response of species groups to environmental gradients. In particular, the spatial patterns of the 40 classes of the perennial herbs differed highly

19

Figure 2.1: Spatial distribution of the total species richness and the species richness 20 of the four plant groups with differing longevity. Grey regions represent areas with climatic parameters extending the range of values used for calibration.

from each other. Even though our hypothesis of improving predictive species richness models by considering the autecology of the species was not substantially supported by the data, it is not necessarily falsified either for the following reasons: (1) The sampling unit of 10 m2 is small and does not reflect the local flora of a site well enough. Thus, random effects and non equilibrium processes (Huston, 1979; Huston, 1994) introduce noise and impede the detection of ecological patterns, which depend on the spatial resolution (Willis and Whittaker, 2002; Rahbek, 2005; Luoto et al., 2007). To circumvent this species package problem, Crawley and Harral (2001) suggested that areas >100 m2 should be used to characterize alpha diversity. (2) In the data set used, more than 50% of all sampling sites were covered by forests. On the other hand the presence of dwarf-shrubs, boulder and gravel sites, mires, peat land or riverside vegetation had low frequencies (< 1% each). Therefore we had incomplete coverage of ecological gradients and some species were not well represented in the data set. This resulted in less precise models for groups containing several such species (Venier et al., 1999; Edwards et al., 2005). A stratified random sampling design may overcome this type of deficit (Till´e, 2001). (3) The selected plant traits may be irrelevant for partitioning of species into ecological meaningful groups. Yet, since some groups attained high model qualities, they seem to show a clear ecologically preference. Also, as one fourth of the functional groups showed a significant habitat affiliation, the selected plant traits seem to carry ecological information. Although the individual errors of the class models averaged to comparable model quality at the overall species richness level, the individual models allow us to make better informed predictions of subgroups at the landscape scale. This clearly improves our understanding of how functional groups are distributed spatially, and how these patterns and their associated species richness may alter in the future under scenarios of global changes. We see this as the most important reason for distinguishing functional groups for such future projections.

Acknowledgements This research was funded by the 6th Framework Programme of the European Union (Contract Number GOCE-CT-2003-505376). We thank the Swiss Federal Agency SAEFL and Hintermann and Weber Consultancy for access to the data from the biodiversity monitoring programme.

21

22

3 Niches and Noise - Disentangling habitat diversity and area effect on species diversity 1

K. Steinmann, 2 S. Eggenberg, 3 T. Wohlgemuth, 4 H.P. Linder, 1 N.E. Zimmermann

1

Land Use Dynamics, Swiss Federal Research Institute WSL, Zuercherstrasse 111, CH-8903 Birmensdorf, Switzerland 2

Atelier fr Naturschutz und Umweltfragen UNA AG, Mhleplatz 3, CH-3011 Berne

1

Forest Dynamics, Swiss Federal Research Institute WSL, Zuercherstrasse 111, CH8903 Birmensdorf, Switzerland 4

Institute for Systematic Botany, Zollierstrasse 107, CH-8008 Zuerich, Switzerland

Abstract: The species-area curve of Arrhenius is generated by habitat heterogeneity, immigration and extinction, and various other stochastic processes. Even though the use of environmental variables is widespread to predict the spatial distribution of species richness, it remains difficult to distinguish the relative importance of habitat heterogeneity and the area effect on total species richness. In our study we used different types of species area curves to disentangle the habitat heterogeneity effect and the area effect on species richness. We generated three types of sample rarefaction curves. (1) A randomly aggregated rarefaction curve, (2) a rarefaction curve where areas of similar habitat types were aggregated and (3) a rarefaction curve, where areas of dissimilar habitat types were aggregated, based on three floristic surveys from Switzerland, using different grain sizes. The classification of the habitat types was based on three environmental variables. From the difference of the dissimilar and similar rarefaction curve, we estimated that habitat heterogeneity contributed nearly 40% to the species richness. The remaining 60% are likely related to effects caused by the area and by stochasticity. However, the habitat heterogeneity effect varied with grain size, contributing only 20% when a sample unit of 1.25 ha was

23

used instead of 10m2 or 28m2 , thus the detection of the habitat heterogeneity effect depends on the grain size. In conclusion, we expect the prediction of biodiversity from habitat heterogeneity to be most difficult at small grain sizes, while at medium to large grain sizes the predictive power of environmental gradient and heterogeneity variables is likely to increase. Yet, habitat heterogeneity is only a powerful predictor at intermediate grain sizes, since heterogeneity is becoming uniform at larger scales of sampling units. Keywords: Grain size, habitat heterogeneity, spatial scale, species area response, species richness, vascular plants

3.1 Introduction Modelling the spatial distribution of species richness is a powerful tool for testing hypotheses in ecology (Currie, 1991; Curry et al., 2004; Moreno-Rueda & Pizarro, 2007; Whittaker et al., 2007) and has gained an increased importance in conservation planning (Margules et al., 1994; Pearson & Carroll, 1998; Cabeza et al., 2004; Moilanen, 2005; Pawar et al., 2007). Most predictions for species richness are based on climate and topographic variables. These variables are more easily obtained than species lists from the relevant areas. However, the predictive power of climate and topographic variables depends on the spatial scale (as reviewed by Willis & Whittaker, 2002; Ricklefs, 2004 and Rahbek, 2005) and the type of organism (Boone & Krohn, 2000; Whittaker et al., 2007; Moreno-Rueda & Pizarro, 2007). For example, climate and topography explain up to 95% of species richness of amphibians, but they explain barely 60% of bird’s species richness (Boone & Krohn, 2000). These values are reduced to 72% and 15%, respectively, when a larger grain scale is used (Whittaker et al., 2007). For plant species richness, values from less than 15% up to 85% are reported (Currie & Paquin, 1987; Heikkinen & Birks, 1996; Birks, 1996 and Lobo et al., 2001), indicating that species richness is not driven by climate alone. In fact, it seems likely that species richness is determined by multiple processes operating at different spatial scales (Shmida & Wilson, 1985; Huston, 1994; Turner, 2005). The number of species found in small areas (micro-scale) depends on the size and growth form of the species (Crawley & Harral, 2001). On a meso-scale, species richness is thought to be affected by equilibrium processes such as immigration and local extinction (Mac Arthur & Wilson 1963, 1967), as well as habitat diversity (Williams, 1964; Abele, 1974; Deshaye & Morisset, 1988; Kohn & Walsh, 1994; Triantis et al., 2003). When entire provinces, e.g. at the macro-scale, are considered, the degree of evolutionary independence between species pools (addition of new floras or faunas) dictates the species richness (Preston, 1960; Rosenzweig, 1995).

24

The interplay of these processes leads to the well known pattern of the species area relationship, which was first described by Arrhenius (1921) and Gleason (1922). On a small scale, species richness increases rapidly with area as a result of increasing sampling effort. Once the area is big enough to contain all potentially detectable species the species richness increases further with area as the chance of inclusion of different habitat types carrying different species compositions increases. In addition, with increasing area of suitable habitats, the extinction probability of single populations decreases and the chance of immigration of new species increases, thus leading to a further increase of the species area curve. On a large scale, species richness increases further with area, as different species pools, which evolved independently, are aggregated (Ricklefs et al., 2004a). It is generally accepted that the species area curve follows a power law function (Arrhenius, 1921; Preston, 1960, 1962) rather than an exponential function (as proposed by Gleason, 1922, but see Tjrve, 2003, for a review of possible functions of species area curves). There is also strong evidence that the slope of log log curves depends on the habitat where species richness was observed (Crawley, 2001). However, when the area includes more than one habitat, we do not know what proportion of the shape of species area curves is explained by processes like habitat diversity, dispersal ability, and the spatial distribution of sink- and source populations. According to the unified neutral theory (Hubbell, 2001), the species area relationship is independent of environmental heterogeneity and solely the consequence of ecological drift. Hubbell shows that the shape of the species area curve depends only on the fundamental biodiversity number and the immigration rate, i. e. the shape of the species area curve is a function of the area occupied by the metacommunity, the mean density of individuals per unit area, and the speciation- and immigration rate only. One could argue that the mean density of individuals itself is a function of climate and its variability. However Hubbell’s model is not the only one resulting in a power-law species area relation without including habitat heterogeneity (see Leitner & Rosenzweig, 1997; Harte et al., 1999; Mc Gill & Collins, 2003). On the other hand empirical studies reported a positive correlation between habitat heterogeneity and species richness (Kohn & Walsh, 1994; Rosenzweig, 1995; Ricklefs & Lovette, 1999). The main issue in the debate about the effects of area and habitat heterogeneity on species richness is that area and habitat diversity are strongly correlated (Simberloff, 1976; Kohn & Walsh, 1994; Ricklefs & Lovette, 1999). Consequently species richness can readily be predicted by either of the two parameters. Area effect per se cannot be easily disentangled from the habitat diversity effect and so the interrelationship cannot be quantified. In addition, the area effect is often confounded with the sampling effort (Connor & McCoy, 1979; Cam et al., 2002). The detection probability of new species increases with the size of area, and therefore independently of the habitat heterogeneity and equilibrium processes, a higher number of different species is found within larger areas. Here, we assume that species richness is affected by habitat diversity. The aim of

25

our study was to disentangle the contributions of habitat diversity and area effect on vascular plant species richness at different grain sizes. Therefore we constructed three types of sample rarefaction curves (Gotelli & Colwell, 2001): (1) a randomly aggregated rarefaction curve (random curve); (2) a rarefaction curve where areas of similar habitat types are aggregated (similar curve); and (3) a rarefaction curve, where areas of dissimilar habitat types are aggregated (dissimilar curve). If habitat diversity is indeed a driver of species richness, then we expect that the shape of the three constructed curves differ significantly. Based on the classical niche theory (Hutchinson, 1957), we expect the highest accumulation rate for aggregation of dissimilar habitat types (Palmer, 2002). Low saturated species richness is expected when similar habitats are aggregated, whereas the random aggregation of habitat types (considered as a null model) should result in an intermediate species area curve. The proportion of the difference of the saturation levels of the dissimilar curve and the similar curve to the total species richness (dissimilar curve) is assumed to indicate the contribution of habitat diversity to the total species richness. To detect scale dependent patterns, we estimated the effect of habitat heterogeneity on three different grain sizes.

3.2 Material and Methods 3.2.1 Study area The study area encompassed Switzerland. The climate is generally humid temperate with rather mild winters and moderately warm summers. Annual mean temperature ranges from ca. 1.5C◦ (in high alpine valleys) to 11.5C◦ (in the southern part of Switzerland). The mean annual precipitation sums vary between 600mm and 2900mm. The topographic relief ranges from ca. 190m to 4634 m.a.s.l.

3.2.2 Species richness data For our analysis we used three different data sets encompassing different grain sizes. Two data sets are from the biodiversity monitoring program of the Swiss federal office of environment (Plattner, Birrer & Weber, 2004). In this monitoring program, vascular plant composition was recorded in circular plots with an area of 10m2 (plot data) and on a transect basis with length of 2.5km and width of 5m (transect data). In both data sets the sampling design is a regular grid with mesh size 6x4km (n = 422) and 12(24)x8km (n = 270) respectively (see Fig. 3.1). The initial points of the grids were generated randomly. The third data set contains grasslands only and originates from the dry meadows and pastures survey of the Swiss Federal Office of Environment (dry meadow data, Eggenberg et al., 2001). The plots (n = 13911) were recorded based on phytosociological methods. They include vegetation types of four alliances: Mesobromion, Xerobromion, Stipo-Poion and Cirsio-Brachypodion

26

(according to Delarze et al., 1999). The number of plant species was recorded within a circular area of ca. 28m2 (radius = 3m). The spatial distribution of the plots is given in Fig. 3.1.

3.2.3 Habitat classification Temperature, water and light as well as nutrients are the most important factors for a plant to grow. Therefore we used the following three environmental parameters as surrogates for the habitat classification: Mean annual temperature, mean moisture index of May, and the slope of the terrain. The slope was derived from a digital elevation model (DEM) with a spatial resolution of 25m. The same DEM was used to generate the climatic layers from meteorological measurements (for details see Zimmermann & Kienast, 1999). The mean annual temperature was derived from monthly Normals of the period 1961-1990. The moisture index was derived from the difference of the monthly precipitation sum and the potential evapotranspiration. The environmental parameters for the dry meadow and plot data were directly obtained by the climatic layers with a resolution of 25m. For the analysis of the transect data, where each transect covers an area of 12500m2 , we took the focal mean of all 25x25m2 pixels within windows of 100x100 meters and resampled the predictor grids to a resolution of 100m. Intersecting these layers at the center-coordinates from each transect field resulted in the environmental parameters for the transect data. Each of the three environmental parameters was classified into five classes of equal interval length, i.e. the range of a parameter (the difference of the maximal and minimal value) was equally split in five classes. This way, the 3-dimensional environmental parameter space spanned by the three parameters was divided into 125 (53) hyper-cubes. Each hyper-cube was considered as one habitat type. Based on the combination of a site’s temperature, moisture and slope, it was assigned to a specific environmental habitat type. By this procedure we produced a habitat map where each pixel was assigned to one of the possible 125 environmental habitat types. The spatial analyses were made with Arc Info Version 9.2 (Environmental Systems Research Institute, Inc., 380 New York Street, Redlands, CA 92373-8100, USA).

3.2.4 Constructing sample based rarefaction curves Three types of sample based rarefaction curves were constructed: a randomly aggregated rarefaction curve which served as a null model (random curve), a rarefaction curve of similar habitat types (similar curve) and a rarefaction curve of dissimilar habitat types (dissimilar curve). For each dataset the sites used for building the random curve were randomly selected from the whole pool of records. To construct the second curve of dissimilar habitats, only one record per habitat type was randomly sampled, i.e. the aggregation level of the curve was limited by the number of habitat types. As not all of the 125 defined hyper-cubes (habitat types) were realized in

27

the recorded sites, the number of aggregated areas was lower than 125. The third type of curve was composed of single curves representing sample rarefaction curves of similar habitats. Therefore the pool for random selection of records contained sites with similar habitats only. Groups of 3x3 neighbouring hyper-cubes were defined as similar habitat types. In total we defined 8 bigger cubes containing part of the small hyper-cubes. For each of the eight bigger cubes, we built a sample rarefaction curve analogous to the random curves described above. The mean of the eight curves represented the species area relationship of similar habitats. All three types of curves were iterated 100 times. Calculating the mean of the 100 iterations resulted in a smoothed species area curve and allowed to estimate a confidence interval.

3.2.5 Statistics The different types of curves were compared by using an F-test. For that purpose we assumed that the curves follow an Arrhenius power law function. Therefore all types of constructed curves were fitted to the Arrhenius function, using the non linear minimization procedure provided by the statistical software environment R. The fitted curves of similar and dissimilar habitat types respectively were tested against the relevant random curve. The contribution of habitat diversity to species richness was calculated as the ratio of the difference of the species richness between the rarefaction curves of similar and dissimilar habitats to the species richness indicated by the latter curve. All the algorithms were written in the statistical software environment R (version 2.4.1; R development Core Team, 2006), which was also used for all statistical analyses performed.

28

3.3 Results A consistent pattern of sample rarefaction curves was found with all three data sets (see Fig. 3.2). As expected from the niche theory, the aggregation of dissimilar habitat types resulted in higher species accumulation rates and saturation levels, compared to the randomly aggregated sample rarefaction curves, while the aggregation of similar habitat types resulted in lower species accumulation rates and saturation levels. In all three datasets both types of curves, i.e. aggregated similar as well as dissimilar habitat types, significantly differed in their species accumulations compared to the random sample rarefaction curves which were considered as null models (see Tab. 3.1). The contribution of habitat diversity to total species Table 3.1: Comparison of the different types of fitted sample rarefaction curves. Sample rarefaction curves of similar habitat types and dissimilar habitat types respectively, were compared to randomly aggregated rarefaction curves (which were considered as null models).

Data set

Type of curve tested against the null model

F-value

Degrees of freedom

p-value

Plot data

Similar Dissimilar Similar Dissimilar Similar Dissimilar

462.62 31.76 193.19 53.38 148.70 973.02

2, 2, 2, 2, 2, 2,

< < < < < 0 C◦ ) Precipitation days per growing season pday Summer - Winter precipitation Site water balance Water holding capacity Continentality index Number of frost days

Abbreviation

Unit

prcp[i] etpt[i];[10] index mind[i];[4];[7] srad[i];[12] tave[ i];[11] ddeg[i] # d prcp swb bucket gams sfroy

1/10 mm (1/10 mm) / day 1/10 mm kJoule/m2 /day 1/10 ◦ C 1/10 ◦ C days 1/10 mm (1/10 mm) / year mm / m3 index # day/100

Topographic parameters Elevation Aspect Slope Topographic position Topographic wetness index

elev aspval slope topos twi

m a.s.l. 0 to 100 degrees - inf to + inf index

Edaphic parameters Soil depth Course fragment content Nutrient content Presence of limestone

depth cfc nutri ca

cm % mval / cm2 binomial

Vegetation coverage Jura meadow Unproductive land Deciduous forest Coniferous forest Mixed deciduous forest Mixed coniferous forest Unspecified forest Remaining coverage

jura unprod decf conf mdecf mcf forest remain

binomial binomial binomial binomial binomial binomial binomial binomial

81

Table 5.2: Plant traits (rosette type, shoot metamorphosis, lateral growth form, leaf anatomy) and subgroups clustered thereof are listed together with the cluster habitat affinity (where low p-values indicate that the null model is true). The abbreviations are defined in Table 2.1 of the main script.

82

Group

Rosette type

Shoot metamorphosis

Lateral growth

Leaf anatomy

Habitat Type

p-value

c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20 c21 c22 c23 c24 c25 c26 c27 c28 c29 c30 c31 c32 c33 c34 c35 c36 c37 c38 c39 c40

hr nr fr hr hr hr hr hr hr hr hr fr nr fr nr nr hr hr hr hr hr nr hr hr hr nr nr nr nr hr nr nr fr hr hr hr fr hr nr fr

t t t t n rt rh p rp rh p rh t rh rh rh rht n rh rp rh t rh t t p t p t n t n rh t t rh n n rh b

standard standard standard standard standard standard standard standard standard standard standard standard standard standard standard standard standard standard standard standard tussock standard standard standard standard standard turf standard standard cushion standard standard standard tussock tussock tussock tussock tussock standard standard

he hy me me me me me sc me hy me me sc hy me sc me sc sc sc me me he hy sc sc me me ls sc he me sc sc me sc sc me hy ls

4 9 6 6 6 6 6 6 6 2; 6 6 6 3; 9 6; 2 6 6 6 6 6 6 6 2 6 6 6 3 3; 1 2 6 6 6 6 6 6 2; 6

0.44 0.00 0.016 0.00 0.004 0.00 0.001 0.007 0.003 0.00 0.007 0.001 0.002 0.005 0.009 0.00 0.00 0.00 0.005 0.22 0.00 0.005 0.003 0.008 0.036 0.356 0.114 0.052 0.222 0.709 0.029 0.001 0.074 0.00 0.002 0.021 0.287 0.003 0.001 0.186

6

7; 9 7

6

4; 9

Table 5.3: Predictors and model quality for each plant group. Predictors labelled with an asterisk entered the model as linear and quadratic terms. MAE: mean absolute error. D2 : deviance explained.

Group

Predictors

All species Trees Shrubs (Bi)annuals Perennials c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20 c21 c22 c23 c24 c25 c26 c27 c28 c29 c30 c31 c32 c33 c34 c35 c36 c37 c38 c39 c40

diffprcp, tave11*, ca, jura, remain, unprod mind7, diffprcp*, etpt10*, sfroy*, nutri, tave11, topos*, decf, mdecf gams*, cfc, sfroy, nutri*, tave11*, remain, unprod diffprcp*, jura, remain sfroy, tave11, topos, ca, jura, remain, unprod cfc*, tave11, jura, mcf gams*, mind7, diffprcp*, mind4*, aspval, nutri*, tave11*, mcf tave11*, jura, remain, unprod diffprcp, tave11*, jura, remain, unprod, forest forest tave11*, ca mind7, diffprcp, cfc, nutri, tave11, jura gams*, srad12, remain slp, diffprcp, jura, remain, unprod srad12, tave11*, twi25s*, remain mind7, diffprcp, sfroy*, tave11, ca, jura, remain, unprod mind7, etpt10*, sfroy, jura, unprod mind4, mcf, remain gams, mind7, bucket, aspval, tave11* etpt10, bucket, sfroy, nutri, tave11*, topos, conf diffprcp*, tave11*, jura etpt10*, srad12*, mind4*, sfroy*, topos slp*, mind7 tave11*, jura, remain, unprod gams*, ca, jura, remain, unprod aspval*, tave11*, jura, unprod mind7*, cfc*, tave11*, topos* mind7, srad12, tave11*, twi25s, jura, remain nutri*, twi25s* slp, mind7, depth*, forest tave11, ca, remain tave11, unprod gams, aspval*, sfroy, jura, remain, unprod nutri sfroy* diffprcp*, depth, bucket*, sfroy, tave11*, twi25s*, topos*, remain nutri* gams*, depth*, mind4*, tave11*, topos*, ca, conf cfc*, aspval* mind7, diffprcp, jura, remain, unprod tave11, jura, mcf, unprod gams*, aspval, jura cfc*, mind4*, jura srad12, mind4*, tave11* slp*, diffprcp, mind4, sfroy, jura, unprod

MAE 7.9 0.75 1.44 0.73 7.02 0.3 0.5 0.3 0.9 0.2 0.1 1 0.2 0.4 0.4 0.8 0.7 0.4 0.4 0.7 0.5 0.2 0 0.4 0.2 0.7 0.7 0.1 0.1 0.3 0.2 0.2 0.4 0.2 0.1 0 0.1 0.6 0.1 0.5 0.3 0.5 0.4 0.2 0.1

D2 0.28 0.68 0.17 0.06 0.34 0.19 0.40 0.33 0.29 0.03 0.22 0.25 0.20 0.19 0.27 0.37 0.31 0.13 0.24 0.17 0.11 0.29 0.12 0.42 0.23 0.21 0.12 0.32 0.10 0.17 0.15 0.09 0.25 0.11 0.41 0.67 0.09 0.18 0.19 0.21 0.34 0.17 0.24 0.18 0.35

83

84

Appendix II WinBugs code used for the models in Chapter 4 model { # priors b0Mean ~ dnorm(0,0.01) b1Mean ~ dnorm(0,0.01) b0SD ~ dunif(0,100) b0Prec