Predicting hotspots for threatened plant species in boreal peatlands

Abstract

Understanding the spatial patterns of species distribution and predicting suitable habitats for threatened species are central themes in land use management and planning. In this study, we examined the geographic distribution of threatened mire plant species and identified their national hotspots, i.e. areas with high amounts of suitable habitats for threatened mire plant species. We also determined the main environmental correlates related to the distribution patterns of these species. The specific aims were to: (1) identify the environmental variables that control the distribution of threatened peatland species in a boreal aapa mire zone, Finland; and (2) to identify the richness patterns and hotspots of threatened species. Our results showed that the combination of individual species models offers a useful tool for identifying landscape-scale richness patterns for threatened plant species. The modeling performance was high across the modelled species, and the richness patterns generated by single models coincide with the expected richness pattern based on expert knowledge. The method is therefore a powerful tool for basic biodiversity applications. In cases where reliable models for species occurrences and hotspots can be produced, these models can play a significant role in land-use planning and help managers to meet different conservation challenges.

Keywords

Introduction

Peatlands are core ecosystems of biological diversity and are known for their wide range of ecosystem services (Ramsar Convention Secretariat 2013). As highly productive ecosystems, they are used increasingly to support economic development and human wellbeing. Drainage and resource exploitation of wetlands are the main reasons why they are among the most threatened ecosystems in the world. For example the area covered by peatlands (the most widespread wetland type) has reduced by 10–20% since 1800 (Joosten and Clarke 2002). In Finland, over half of peatlands have been drained for forestry (Finnish Forest Research Institute 2014), which has caused habitat degradation and increased the number of threatened peatland species. At present, there are 223 red-listed vascular plant and bryophyte species with peatlands as their primary habitats (4.5% of all red-listed species) and 420 red-listed species with peatlands as one of their habitats (Rassi et al. 2010). The ongoing bioeconomy development (Spatial Foresight, SWECO, ÖIR, t33, Nordregio, Berman Group, Infyde 2017) and high interest in Arctic countries for their mineral resources (Boyd et al. 2016) are increasing pressures on peatlands. These intense activities are expected to have strong, and mostly negative, impacts on peatland biodiversity.

Many of the most adverse effects resulting from peatland use can be avoided through careful planning. This requires an analysis of potential ecological values in an area before intensive and/or large-scale use is planned and carried out. Hotspots, or concentrations of threatened species, are important surrogates of biological diversity that have a significant role in conservation and management strategies (Gaston 1994). Although locations of hotspots should not be the guiding principles in land use planning, they can be used to avoid disturbing valuable sites with high numbers of rare species (Loiselle et al. 2003; Elith and Leathwick 2009). Predictive species-distribution modeling offers a cost-effective method of exploiting the limited empirical data for the evaluation of biodiversity. Statistics-based spatial models are valuable for generating biogeographical information that can be applied across a broad range of fields, including ecology, land use planning and climate change (e.g. Barbet-Massin et al. 2012; Bolliger et al. 2007; Thuiller et al. 2008). Habitat suitability models rely on the concept of niche conservatism (the tendency of species to retain ancestral ecological characteristics) and assume that environmental variables will play an important and consistent role in shaping species distributions (Wiens and Graham 2005). Predictive habitat suitability models of species’ geographical distributions and species richness are increasingly used as an alternative for incomplete or spatially biased survey data as a basis for conservation planning (Hirzel and Le Lay 2008; Elith and Leathwick 2009; Freeman et al. 2013; Lemes and Loyola 2013).

A traditional way to develop spatial projections of species richness is to directly measure numbers of species from surveyed sites and to relate this information to environmental variables derived from GIS data. The analysis produces models that yield predictions of species richness for unsampled sites. In this study, species are first modelled individually, and species richness is estimated by stacking individual habitat suitability models (see also Algar et al. 2009; Parviainen et al. 2009; Mateo et al. 2012, 2013). The top 5% of grid squares ranked by species richness can be classified as hotspots. Individual models are constructed by relating species occurrence data to environmental variables and projecting the modelled relationships onto geographical space (Elith et al. 2006). This method may provide some useful advantages, such as better control for poorly modelled species, and easier identification of the set of the most important explanatory variables and of the response shapes between species and their environment in certain subgroups of species.

The aim of this study was to provide a comprehensive picture of environmental prerequisites for a set of threatened peatland species, thereby helping to plan peatland use in a more ecologically sustainable way. The specific aims were: (1) to scrutinize how the occurrence of species is affected by different environmental correlates; and (2) to identify the richness patterns and hotspots of threatened species. The models were developed for the “aapa mire zone” in central and northern parts of Finland. It is worth noting that species richness per se was the target of our study. Predictive modelling was focused on the richness of species that characterize valuable environments that need specific attention in planning. The study thereby provides a new approach to focused biodiversity modeling.

Study area

The study was carried out in Finland, located between 63° and 68° latitudes in northern Europe (Fig. 1). Biogeographically the study area lies in the middle and northern boreal zones covering almost the entire aapa mire zone, where climate is more continental than in most other parts of northern Europe but with some humid, maritime effect (Ahti et al. 1968). The annual mean temperature declines from south (+ 5 °C) to north (− 2 °C) and the mean annual precipitation sum varies between 450 and 750 mm (Pirinen and Ruuhela 2012). Peatlands and pine and spruce-dominated forest are frequent, as well as numerous lakes and rivers characterizing the landscape of the study area.

The location of the study area in a boreal landscape in northern Finland, together with the distribution map with the observational points of the threatened plant species studied. Land cover classification is based on data about the drainage status of peatlands (SYKE)

The peatland habitats of the studied plant species are of three types. Mesotrophic fens are mainly open peatlands with deep peat deposits. The field layer vegetation is characterized by sedges and herbaceous plants, and the ground layer consists of sphagnum mosses or other bryophytes. Rich fens are open or sparsely wooded peatlands with high species diversity of vascular plant and mosses. They are typically found in areas where the bedrock and soil are calcium-rich. Approximately half of Finland’s threatened peatland species are primarily associated with rich fens (Rassi et al. 2010). Spruce swamp forests are wooded minerotrophic peatlands where the dominant tree species is usually Norway spruce (Picea abies), though deciduous trees may also grow abundantly in Spruce forests that are richer in nutrients. The presence of living and dead trees of different sizes and ages is an important structural feature for the species diversity of Spruce swamp forests. Finnish mires have been intensively drained in the last century, and more than half of the 10.0 million hectares of originally pristine mires have been drained to improve timber growth (Finnish Forest Research Institute 2014).

The study area was divided into grid cells of 25 ha (500 m × 500 m), and cells where peatlands covered less than 5% were excluded from the study. Thus, the study area consists of a total of 500,545 grid cells (125,136 km2).

Plant species data

We used the occurrence records of threatened mire plant species from the national database of red-listed species (Rassi et al. 2010) (Table 1). The field records produced by voluntary amateur and professional botanists are the most important data sources for this database, but information on species occurrences has also been gathered from the scientific literature and herbaria (Ryttäri et al. 2012; Rassi et al. 2010). Species data included detailed information on the geographical location of the occurrences (coordinates in the uniform grid system, Grid 27°E). In total, 48 species with ten or more records in the whole study area were used in the analyses (Fig. 1, Table 1). Only observations with accuracy better than 100 m and presence observations from 1990 or later, were selected for this study. As the databases of red-listed species do not include records for the absence of species, the assumption was made that the absence of a record from a sampled grid square corresponded to true absence of the species, because a quasi-exhaustive sampling could be assumed for most squares with presence records (Guisan and Zimmermann 2000).

Table 1

Number of presence records in the study area, red list category of the species, main habitats of the species and group for the studied plant species

We modelled the habitat requirements for all species by using the same environmental predictors for each species. Based on their different associations with the various environmental predictors, we grouped the species into five groups as follows: Mesotrophic fen species (n = 6), rich fen species (n = 10), calcareous species (n = 22), spruce swamp forest species (n = 3) and decaying wood species (n = 7). Rich fen species and calcareous species are partially overlapping. The reason for separating these two species groups was that rich fen species can also be found outside the calcareous areas whereas calcareous species are restricted only to calcareous habitats.

Environmental correlates

We selected a set of quantitative correlates that would reflect the main biophysical gradients with a recognized, physiological influence on plants. In total, 17 environmental variables were calculated for all of the studied grid squares of 25 ha and were then used to explain plant species distribution. Two correlates/variables indicated climate, one topography, one geology and 13 local habitat features (Table 2). Correlations among these variables were only moderate (Spearman correlation < 0.7) and thus, none of the variables was excluded a priori from the actual modelling.

Table 2

List of selected environmental variables used as explanatory variables in the modeling experiments

Environmental variable

Abbreviation

Unit

Mean [min–max]

Data source

Growing degree days (> 5 °C)

GDD

–

910 [517–1154]

FMI

Mean water balance

WAB

mm year−1

330 [137–442]

FMI

Mean topographical wetness index

TWI

–

14.11[6.85–25.22]

NLS, DEM

Proportion of undrained peatlands in grid square

UNDRAINED

%

16.9 [0–100]

SYKE

Proportion of drained peatlands in grid square

DRAINED

%

23.4 [0–100]

SYKE

Proportion of open peatlands in grid square

OPEN

%

6.6 [0–100]

SYKE

Proportion of herb-rich site type (peatlands) in grid square

KP1

%

0.5 [0–100]

Luke, MS-NFI

Proportion of Vaccinium myrtillus site type (peatlands) in grid square

KP2

%

2.1 [0–100]

Luke, MS-NFI

Proportion of Vaccinium vitis-idaea site type (peatlands) in grid square

Temperature and moisture requirements reflect the principal limitations on plant growth and survival (Skov and Svenning 2004). Thus, growing degree days (> 5 °C) (GDD) and water balance (mm) (WAB) were calculated for the years 1981–2010 from climate data with 1 km2 resolution (Finnish Meteorological Institute, Pirinen and Ruuhela 2012). Water balance was used because precipitation alone is not a good measure of the water available for plant growth. A simple water balance variable was calculated as the monthly difference between precipitation and potential evapotranspiration, as suggested by Skov and Svenning (2004). The potential evapotranspiration (PET) was calculated as:

Moisture conditions in peatlands can be related to many ecological processes across landscapes, e.g. species composition and distribution, peatland productivity (Iverson et al. 1997) and hydrology in terms of ombrotrophy and minerotrophy. Topographic wetness index (TWI) was used to describe local relative differences in moisture conditions (Gessler et al. 2000). High values represent lower catenary positions (wet) and small values upper catenary positions (dry). The moisture level of the study area was calculated by defining the wetness index (or compound topographic index) using the following formula (Burrough and McDonnell 1998):

where ɑ is the upslope contributing area per width orthogonal to the flow direction, and tanβ is the local slope in radians.

Moreover, data about the drainage stage of peatlands with a resolution of 25 m × 25 m (Finnish Environmental Institute 2009) were used to calculate the proportion of undrained (UNDRAINED) and drained (DRAINED) peatland area as percentage cover for each grid square.

As many of the threatened plant species require calcareous substrate, the proportion of calcareous rock (CALC) as percentage extent for each grid square was calculated from digital maps of Quaternary deposit and pre-Quaternary rocks (Digital Base Map, NLS) using ArcGIS software (ESRI 1991). The presence of springs (SPRINGS) was used as many of the modelled species benefit from springs.

Information on the percentage cover of main peatland site types and site fertility in each 25-hectare grid square was derived from the Multi-source National Forest Inventory (MS-NFI) from 2011 (Natural Resources Institute Finland 2013) with a resolution of 20 m × 20 m. Site fertility classes were selected to match the habitat requirements of the studied plant species as accurately as possible: eutrophic peatlands and corresponding drained peatlands (namely herb-rich types, KP1), mesotrophic mires and fens and corresponding drained peatland forests, (Oxalis-myrtillus type, KP2), meso-oligotrophic natural and drained peatlands (Vaccinium myrtillus type, KP3), and Sphagnum fuscum-dominated (ombrotrophic) natural and drained peatlands (Cladina type, KP6). The proportions of open peatlands (OPEN MIRE) and spruce swamp forests (SPRUCE SWAMP FOREST) in grid squares were used to reflect general habitat patterns of peatland properties in each grid squares. Moreover, mean volumes (m3/ha) of four tree species—pine (PINE) (Pinus sylvestris), spruce (SPRUCE), birch (BIRCH) (Betula pendula and B. pubescens) and other broad-leaved trees (OTHER) were calculated from MS-NFI-data and employed in the modeling.

Habitat suitability modelling

The presence-only habitat suitability modelling method Maxent v3.3.3 k (Phillips et al. 2006) was used to predict species distributions across the aapa mire zone. The resulting habitat suitability model represents the relative probability of the species’ distribution over all grid squares in the defined geographic space, where a high probability value indicates that the location is predicted to have suitable environmental conditions for the species (Hirzel et al. 2002). Maxent has been utilized extensively to model species’ ranges using presence-only data, and it has been shown to perform well even with scarce and noisy presence data subsets collected by different researchers and methodologies (Elith et al. 2006; Franklin 2010). Maxent has also performed well in modelling other ecosystem services, such as the distribution of GHG-balances (Parkkari et al. 2017).

To be able to compare and combine or stack models for multiple species, the same environmental predictors and Maxent parameters were used for all species. Model calculations were made using the Maxent logistic output, rather than raw or cumulative output, in order to facilitate comparisons between species (Merow et al. 2013). Maximum iterations were set at an average of 5000, based on model performance across all target species. The remaining settings were left at the default setting. Moreover, response curves were created to show how the predicted relative probability of occurrence depends on the value of each environmental variable.

Model evaluation

The models and model predictions were evaluated using the area under the curve (AUC) of a receiver operating characteristic (ROC) plot based on a four-fold cross-validation (Fielding and Bell 1997), routinely calculated for each run with Maxent. Cross-validation was performed with subsets of the entire dataset, where each subset contained an equal number of randomly selected data points. Each subset was then dropped from the model, the model was recalculated, and predictions were made for the omitted data points. A combination of the predictions from the different subsets was then plotted against the observed data (Lehmann et al. 2002). Following Swets (1988), model accuracy was considered low if AUC was below 0.7, fair if it was between 0.7 and 0.8, good if between 0.8 and 0.9, and excellent if AUC was above 0.9.

To identify the relative importance of various environmental variables for species, we employed two outcomes of the Maxent model: percent contribution and permutation importance of each environmental variable. The percent contribution values are only heuristically defined: they depend on the particular path that the Maxent code uses to arrive at the optimal solution, and a different algorithm could give rise to the same solution via a different path, resulting in different percent contribution values. If there are highly correlated environmental variables, the percent contributions should be interpreted with caution. The permutation importance measure depends only on the final Maxent model, not the path used to obtain it. The contribution for each variable is determined by randomly permuting the values of that variable among the training points (both presence and background) and measuring the resulting decrease in training AUC. A large decrease indicates that the model depends heavily on that variable. Hence, permutation importance appears to be a better measure of a variable’s explanatory power—since it is path—(algorithm-) independent. Modelling performance was evaluated using the regularized training gain, which describes how much better the Maxent distribution fits the presence data compared to a uniform distribution (Phillips and Dudic 2008).

A jackknife test was also run to obtain alternate estimates of variable importance. Each variable was excluded in turn, and a model was created with the remaining variables. The model was then created using each environmental variable in isolation. In addition, a model was created using all variables. For the variables with the highest predictive values, response curves show how each of these environmental variables affects the Maxent predictions (Phillips and Dudík 2008). The curves illustrate how the logistic prediction changes as each environmental variable is varied, while keeping all other environmental variables at their average sample value. The curves thus represent the marginal effect of changing any single variable alone.

Hotspot maps

First, we produced projected distribution maps for individual species at a spatial resolution of 25 ha. The continuous Maxent output maps were reclassified into binary maps of suitable (1) and unsuitable (0), using the averaged species-specific logistic threshold value that “maximises training sensitivity plus specificity” (Liu et al. 2013). This threshold selection method has been shown to perform rather well with presence-only data (Liu et al. 2005, 2013), and is suitable for this study considering the goal is to predict where current suitable habitats are located. Choosing a relatively high threshold reduces the risk of choosing unsuitable sites by identifying only those areas with the highest suitability (Pearce and Ferrier 2000).

Next, to create richness maps, we combined the binary maps representing suitable habitats for individual species and used a simple summation of the predicted suitabilities using the Raster Calculator feature in ArcGIS v10.2 for each species group separately, and also for all 48 species. The reason for doing so was that this allowed us to investigate whether certain species groups are more intimately related to certain environmental predictors than other groups. Spruce swamp forest species were excluded from the species group analysis, as species richness and hotspot based on only three species is not particularly informative. However, they were included in the analyses of total species richness and summary hotspot based on individual species.

We then identified richness hotspots as the top 5% of grid squares ranked by species richness in each species groups (see Prendergast et al. 1993; Williams et al. 1996). Finally, summary hotspot map was produced by stacking the hotspot maps from all individual species.

Results

For all models, the AUC was excellent for the training data (mean AUC 0.924, ranging from 0.841 to 0.989) and good for the evaluation data (mean AUC 0.855, ranging from 0.607 to 0.959) (Table 3). The highest test AUC values were obtained with calcareous species (mean 0.88) and the lowest with decaying wood species (mean 0.80), perhaps because decaying wood environments exist extensively outside peatland habitats.

Table 3

Performance of the models applied to individual threatened plant species distribution, as assessed by the AUC values from the training (Training AUC) and evaluation data (Test AUC) and the stability of the models (Standard Deviation)

Species

Training AUC

Test AUC

AUC standard deviation

Rich fen

Carex heleonastes

0.92

0.90

0.02

Dactylorhiza incarnata subsp. cruenta

0.93

0.89

0.02

Dactylorhiza incarnata subsp. incarnata

0.84

0.84

0.01

Hamatocaulis vernicosus

0.92

0.91

0.01

Leiocolea bantriensis

0.94

0.82

0.09

Lophozia grandiretis

0.85

0.69

0.14

Meesia longiseta

0.97

0.94

0.03

Moerckia hibernica

0.94

0.87

0.04

Riccardia multifida

0.97

0.96

0.02

Sphagnum contortum

0.90

0.77

0.08

0.92

0.86

0.05

Mesotrophic fen

Carex laxa

0.88

0.81

0.04

Epilobium laestadii

0.94

0.88

0.05

Hamatocaulis lapponicus

0.98

0.89

0.08

Hammarbya paludosa

0.92

0.89

0.02

Lycopodiella inundata

0.86

0.69

0.12

Rhynchospora fusca

0.95

0.94

0.01

0.92

0.85

0.05

Calcareous

Amblyodon dealbatus

0.97

0.96

0.03

Botrychium virginianum

0.95

0.88

0.07

Bryum pseudotriquetrum var. neodamense

0.94

0.91

0.03

Calypso bulbosa

0.94

0.93

0.01

Campyliadelphus elodes

0.95

0.86

0.05

Carex appropinquata

0.92

0.87

0.03

Carex viridula var. bergrothii

0.93

0.88

0.04

Cypripedium calceolus

0.90

0.89

0.01

Dactylorhiza fuchsii

0.91

0.80

0.05

Dactylorhiza lapponica

0.97

0.95

0.02

Dactylorhiza traunsteineri

0.85

0.81

0.02

Dicranum acutifolium

0.95

0.91

0.02

Eriophorum brachyantherum

0.91

0.87

0.04

Malaxis monophyllos

0.99

0.88

0.08

Palustriella commutata

0.95

0.91

0.03

Palustriella decipiens

0.91

0.87

0.03

Palustriella falcata

0.92

0.89

0.02

Philonotis calcarea

0.95

0.90

0.06

Pseudocalliergon angustifolium

0.95

0.87

0.10

Pseudocalliergon lycopodioides

0.96

0.84

0.11

Saxifraga hirculus

0.90

0.89

0.01

Schoenus ferrugineus

0.94

0.86

0.07

0.93

0.88

0.04

Spruce swamp forest

Carex atherodes

0.93

0.84

0.06

Epipogium aphyllum

0.86

0.77

0.05

Poa remota

0.95

0.81

0.08

0.92

0.81

0.06

Decaying wood

Anastrophyllum hellerianum

0.89

0.83

0.03

Calypogeia suecica

0.96

0.81

0.08

Jungermannia leiantha

0.87

0.73

0.10

Lophozia ascendens

0.94

0.88

0.04

Lophozia ciliata

0.92

0.82

0.11

Lophozia longiflora

0.87

0.66

0.10

Riccardia palmata

0.90

0.85

0.06

0.91

0.80

0.07

All species on average

0.924

0.855

0.050

On average, the proportions of drained peatland area, open peatland and undrained peatland area had the highest permutation importance (15.2, 14.7, and 10.2%, respectively) across the modelled species groups (Table 4). When considering importance between species groups, proportions of undrained (mainly positive responses) and drained peatland area (negative responses) showed the greatest impact (18.0% and 18.3%) followed by the proportion of open peatland with positive association (8.8%) in rich fen species models (Tables 4, 5).

Table 4

Relative importance (%) of single environmental variables for predicting the habitat suitability of five species groups based on permutation importance (in Maxent)

Environmental variable

Rich fen

Mesotrophic fen

Calcareous

Spruce swamp forest

Decaying wood

Mean

CAL

3.17

0.38

8.87

2.15

0.52

3.02

GDD

7.38

13.51

8.36

2.42

2.40

6.81

KP1

3.29

4.98

7.26

4.27

7.03

5.37

KP2

0.77

3.56

2.76

1.64

5.74

2.89

KP3

2.93

4.74

1.83

4.79

2.10

3.28

KP6

3.88

2.93

2.67

3.50

5.65

3.73

BIRCH

4.81

5.32

4.76

3.34

3.91

4.42

SPRUCE

3.28

9.71

4.84

6.69

7.03

6.31

SPRINGS

2.70

1.17

1.86

2.15

1.08

1.79

PINE

5.59

3.85

8.89

5.76

4.82

5.78

OTHER

5.09

0.94

3.74

2.18

2.01

2.79

SPRUCE SWAMP FOREST

1.75

4.28

4.81

1.29

0.60

2.55

OPEN MIRE

8.82

12.81

8.15

36.40

7.06

14.65

UNDRAINED

18.04

13.49

8.20

4.31

7.19

10.24

DRAINED

18.33

11.99

13.59

12.98

18.96

15.17

TWI

2.23

4.20

6.86

3.96

11.45

5.74

WAB

7.95

2.15

2.56

2.18

12.47

5.46

Total

100

100

100

100

100

100

Three most important variables for each species group and on average are in bold. Abbreviations are explained in Table 2

Table 5

Summary of the shape of the response curves between modelled plant species (n = 48) and each environmental variable

In calcareous species models, the proportion of drained peatland area was the most influential variable with positive response (13.6%), followed by proportion of calcareous rock (8.9%, mainly positive responses), and volume of pine (8.9%, mainly negative responses). Increasing proportions of undrained peatland (13.5%), growing degree days (13.5%), and proportion of open peatlands (12.8%) contributed most and positively to mesotrophic fen species models.

Spruce swamp forest species were mainly negatively associated with proportion of open peatlands (36.4%), drained peatland area (13.0%), and volume of spruce (6.7%). In decaying wood species, the proportion of drained peatland area (19.0%, mainly negative associations), water balance (12.5%), and topographical wetness index (11.5%) showed the greatest impact.

Looking to the jackknife evaluations across all modelled species, we found that the proportion of undrained peatland, topographical wetness index, and volume of spruce were the three most effective predictors when used individually (Table 6). In addition, proportions of calcareous rock, open peatlands, and drained peatland area decreased the gain most when they were omitted, and therefore contained information that was not present in any other variable (Table 7). However, species groups differed from each other according to their habitat preferences. In mesotrophic fen species models, the proportion of undrained peatland area had the highest gain when used in isolation (Table 6) and the largest decrease in gain when omitted (Table 7). These environmental variables contain information that is useful on its own and not present in other variables. Likewise, the proportion of calcareous rock provided the most useful and unique information on the distribution of calcareous species. In spruce swamp forest species models, the greatest change occurred when volume of spruce was used in isolation.

Table 6

Results of jackknife evaluations of relative importance of environmental variables when used in isolation, with only the corresponding feature

Rich fen

Mesotrophic fen

Calcareous

Spruce swamp forest

Decaying wood

Mean

CAL

0.14

0.01

0.48

0.17

0.00

0.16

GDD

0.14

0.18

0.23

0.02

0.01

0.12

KP1

0.19

0.08

0.12

0.01

0.03

0.09

KP2

0.12

0.04

0.11

0.02

0.09

0.08

KP3

0.29

0.12

0.11

0.04

0.04

0.12

KP6

0.09

0.22

0.07

0.01

0.04

0.09

BIRCH

0.16

0.29

0.09

0.15

0.01

0.14

SPRUCE

0.13

0.44

0.06

0.21

0.09

0.19

SPRINGS

0.16

0.02

0.13

0.08

0.02

0.08

PINE

0.32

0.30

0.21

0.01

0.01

0.17

DECIDUOUS

0.08

0.15

0.05

0.04

0.01

0.06

SPRUCE SWAMP FOREST

0.03

0.07

0.07

0.13

0.02

0.06

OPEN MIRE

0.28

0.27

0.09

0.20

0.06

0.18

UNDRAINED

0.59

0.64

0.20

0.02

0.01

0.29

DRAINED

0.23

0.11

0.20

0.07

0.19

0.16

TWI

0.22

0.33

0.20

0.14

0.36

0.25

WAB

0.16

0.09

0.05

0.02

0.27

0.12

The most important variable for each species group and on average is in bold. Abbreviations are explained in Table 2

Table 7

Results of jackknife evaluations of relative importance of environmental variables when used in isolation, without the corresponding feature

Rich fen

Mesotrophic fen

Calcareous

Spruce swamp forest

Decaying wood

Mean

CAL

1.30

1.45

1.19

0.94

0.93

1.16

GDD

1.37

1.37

1.47

1.08

0.93

1.24

KP1

1.37

1.40

1.47

1.05

0.93

1.24

KP2

1.40

1.43

1.52

1.08

0.91

1.27

KP3

1.38

1.42

1.52

1.06

0.93

1.26

KP6

1.39

1.39

1.51

1.05

0.93

1.25

BIRCH

1.40

1.44

1.52

1.07

0.92

1.27

SPRUCE

1.39

1.43

1.50

1.02

0.88

1.24

SPRINGS

1.30

1.44

1.47

1.04

0.93

1.24

PINE

1.38

1.42

1.49

1.02

0.91

1.24

DECIDUOUS

1.40

1.46

1.52

1.05

0.92

1.27

SPRUCE SWAMP FOREST

1.40

1.45

1.51

1.08

0.94

1.28

OPEN MIRE

1.37

1.38

1.50

0.95

0.89

1.22

UNDRAINED

1.36

1.35

1.51

1.06

0.89

1.24

DRAINED

1.34

1.41

1.48

1.00

0.87

1.22

TWI

1.40

1.43

1.47

1.06

0.86

1.24

WAB

1.33

1.43

1.51

1.07

0.80

1.23

The most important variable, the exclusion of which decreased the gain most, is in bold. Abbreviations are explained in Table 2

Predictions of threatened plant species richness, based on the summation of single-species predictions, and hotpots, are shown in Figs. 2 and 3. Suitable habitats for mesotrophic fen and rich fen species were predicted in the western part of the study area, whereas rich fen species had high suitabilities also in northern parts. Eastern and southeastern parts of the study area had high suitability for presence of decaying wood species. Suitable habitats for calcareous species were, not surprisingly, mainly concentrated in the areas with much calcareous rock.

Spatial predictions of threatened plant species hotspots: by a species groups, and b total species richness (n = 48 species)

Discussion

Regional-scale biodiversity patterns, i.e. spatial resolutions ranging from 0.5 to 2 km, are an important component of the diversity that occurs in a landscape or region, and also represent the scale at which land use decisions are often made. Sustainable planning of peatland use requires an analysis of potential ecological resources and the effects of their utilization on an area. Without knowledge of the biodiversity values of peatlands, their unsustainable utilization continues to degrade their biodiversity and threaten the remaining valuable habitats and species. Predictive habitat suitability modelling may considerably increase the efficiency of biodiversity mapping schemes and incorporate un-surveyed regions into decision-making. Using models to predict potential species distributions is also likely to become increasingly important as environmental change and other dynamic processes are incorporated into land use planning efforts (Rondinini et al. 2006; Underwood et al. 2010). The aim of this study was not to reflect the full reality, but to construct and evaluate simple and ecologically significant habitat suitability models that approximate this reality and constitute useful tools for land use planning. By targeting on the richness of species of valuable environments also provide a new approach to focused biodiversity modelling.

Modelling performance: uncertainty issues

Although the predictive performance of the models was rather high across the species, it is important to be aware that the ecological significance of the observed relationships between environmental variables and species occurrences are not always obvious. In this study, species richness is simply predicted by stacking presence–absence predictions of all species. This method hence relies on our ability to model the distributions of individual species, a field that has greatly matured over the last two decades (see Guisan and Thuiller 2005; Elith and Leathwick 2009; Franklin 2010). Thus, the main factors are those controlling individual species distributions, and often purely abiotic variables are used.

One of the main caveats in using individual species models to generate species richness patterns is that it tends to overestimate actual species richness (i.e., commission errors; Algar et al. 2009; Trotta-Moreu and Lobo 2010; Guisan and Rahbek 2010). The individual species models created in this study do not include all environmental, ecological (particularly competition), and historical factors that affect species distributions. By solely using environmental variables at rather coarse resolution to construct predictions of a species’ suitable habitat, models fail to incorporate biological, geographical or historical influences on species distributions (e.g., Guisan and Thuiller 2005; Heikkinen et al. 2006). This, in turn, can lead to an overestimation of species suitable habitats, as only areas of suitable habitats, not their realized distributions, are projected. This kind of overestimation may increase when combining multiple individual models to create species richness maps, as was the case in this study. In the light of these facts, we propose that the habitat suitability models influence on potential over-prediction is the result of the inherent nature of habitat suitability models.

Moreover, response curves are just simplifications of reality, and their shape may be strongly dependent on the setting of the study and the variable selection criteria used. First, a local model is fit to a particular region of the geographical space, but the model can differ in different regions of the sample space. Certain abiotic factors, such as topography and land cover, may be important locally, but they generally can be applied only within a limited geographical extent (Thuiller et al. 2003). Thus, conclusions about the response curve of species may only be made within the context of the study area. Second, it is possible that species may respond to a combination of a different set of variables in different parts of its distributional range, as the shape of responses in a multivariate model may depend on the nature of the correlations between the indirect variable and the causal gradients (Franklin 1995; Guisan and Zimmermann 2000). The use of different geographical extents and spatial resolutions could provide contradictory answers to the same ecological question.

It should also be kept in mind that an area of suitable habitat is not occupied by a species if the species is unable to disperse there (Pulliam 2000; Newbold 2010). Dispersal limitation (Kadmon and Shmida 1990), source-sink dynamics (Pulliam and Danielson 1991) and metapopulation dynamics (Hanski 2005) will result in spatial patterns in species distributions that are at least partly independent of the environment. These spatial patterns are referred to as “endogenous spatial autocorrelation” (Legendre 1993). Several papers have discussed the importance of measuring spatial autocorrelation when evaluating the importance of different factors to explain species distributions (e.g. Dormann et al. 2007; Hawkins et al. 2007). However, as the 25 ha grid cells in the model setup were distributed rather sparsely across the whole study area we assumed that the effect of spatial autocorrelation was rather small. Moreover, Parviainen et al. (2008) carried out autocorrelation assessment in a similar environment with a similar grid-based approach at the same 25-ha resolution. They found that inclusion of the effect of spatial auto-correlation as autocovariate term reflecting the species occurrences in the surroundings of the focal grid cell, had only a minor effect on the importance of the environmental variables and the shapes of predictor-response curves.

The observed distribution of threatened species is also affected by historical facts that may have restricted current distribution patterns of species (Guisan and Thuiller 2005; Svenning et al. 2008). This may, at least partly, explain the true absence of species in many areas where the environmental conditions are apparently suitable. Moreover, populations of threatened plant species may be extremely small and thus prone to local extinctions arising from stochastic processes in areas with appropriate environmental conditions. In light of this fact, some plant species may actually have had populations in earlier periods of significantly better availability of suitable habitat, and current land use may have no relevance for their potential occurrence (Lindborg and Eriksson 2004; Helm et al. 2006; Wisz et al. 2007). Thus, the individual habitat suitability modelling approach is limited because without adding a dispersal filter, it may incorrectly predict species in areas that appear environmentally suitable but that are outside their colonizable or historical range.

Threshold selection is one of the many possible biases in habitat suitability modelling. As suggested by Trotta-Moreu and Lobo (2010) and Mateo et al. (2012), the selection of an appropriate suitability threshold can reduce over-prediction in species models. However, selection is not straightforward and the results can vary, sometimes dramatically, depending on the threshold chosen (Milanovicha et al. 2012; Liu et al. 2013). We chose to use the rather conservative maximum training sensitivity plus specificity threshold, as we found it to be a promising selection method for presence-only data. An interesting methodological line of future research would thus be to study the reliability of different thresholding approaches in modelling, as it may help to reduce over-predictions, at least in some cases. An additional problem in the selection of reliable and stable threshold values is the lack of real absences, as in the present study. When the modelling algorithm has no information on true absences, even small differences in the selected threshold value can have a substantial effect on the model outputs (see Jiménez-Valverde and Lobo 2007).

Finally, the lack of floristic data from remote areas may possibly lead to a partial bias in the modelling analyses. The accuracy of the model increases with increased amounts and accuracy of presence and absence data, and may be updated to include new information to further refine distribution predictions (Elith and Leathwick 2009), but we assume that the main drivers of habitat suitability (which were predominantly ecologically plausible) will remain.

Model transferability is one important feature in habitat suitability models and thus, developing models that are able to provide reliable predictions of species distributions in new areas or other times is a major challenge. As the results of this study revealed, plant distributions are often critically affected by certain local factors, such as soil or habitat properties or the occurrence of favourable microsites and microclimates (see also Parviainen et al. 2008; Elith and Leathwick 2009). Our models were generated for the boreal aapa mire landscape, where a relatively high proportion of the peatland landscape has remained in a semi-natural state despite intensive draining. On the other hand, the effect of draining is not limited to the drained peatlands only, but may extend over larger areas within the catchment (Holden et al. 2006). The models used in this study took into account only the local drainage effects within the grid cells. From the ecological viewpoint, our models may not be directly applicable to regions of highly fragmented, intensively used or cultural landscape typical of, for example, western and southern Europe. However, from the technical viewpoint, our approach is applicable over different ecosystems and habitats.

Importance of environmental variables to plant species

The distribution of plant species is limited by the availability of suitable habitats. For rare plants, especially those with limited geographic ranges, narrow habitat specificity can further limit distribution. While climate is an important driver of plant species distribution at the continental scale, soil properties and biotic interactions determine habitat availability at smaller scales (Pearson et al. 2004). Furthermore, competition may also affect species occurrence patterns and persistence capability (Virtanen et al. 2010). Variations in peatland vegetation are the result of many environmental factors at the landscape and local scales, including the origins of the water that feeds the peatland, acidity levels (pH), the availability of main nutrients (nitrogen and phosphorus), the water table level, and the depth of the peat. Seasonal variations in moisture levels have also been found to be related to the composition of vegetation communities (Laitinen 2008).

In this study, species responded differently to the analyzed habitat gradients. A mixture of unimodal and linear responses was typical of the gradients. As a general observation, with most of the species the importance of variables reflecting local-scale variation in the habitat and land cover was superior to climate variables operating at higher scales. Responses to growing degree days varied according to the geographic distribution of the species; for example, Eriophorum brachyantherum is a northern species, for which suitable habitats occur at low numbers of growing degree days.

Our model confirms the importance of particular environmental variables that influence the presence and quality of peatland habitat for the selected threatened plant species. In mesotrophic fen species models, a high amount of undrained open peatlands was the most powerful characteristic forcing distribution patterns of threatened species. Wetness, high variation in site fertility types, and microtopography are distinct characteristics for undrained peatlands. The wettest peatlands have not been used for intensive forestry and agriculture, and they therefore continue to be suitable habitats for peatland species. At drained sites, key hydrological characteristics have changed, which has led to the degradation of peatland vegetation (Similä et al. 2014). Based on our findings, species growing on wet surfaces are most susceptible to the effects of changes in the water table. These include Carex heleonastes, Dactylorhiza incarnata ssp. incarnata, Hamatocaulis vernicosus, Meesia longiseta, Carex laxa, Hamatocaulis lapponicus, Dactylorhiza lapponica, Dactylorhiza traunsteineri, and Saxifraga hirculus. Epilobium laestadii is a demanding northern species that is most closely associated with nutrient-rich fens with a sparse and low field layer, and often occurs around springs or in seepage areas. Rhynchospora fusca is a sedge species characteristic of open flark fens.

Amblyodon dealbatus, Cypripedium calceolus, and Dactylorhiza traunsteineri are calcium-demanding species of nutrient-rich, calcareous fens. As expected, the importance of calcareous bedrock was explicit to species demanding calcium, e.g. Amblyodon dealbatus, Malaxis monophyllos, and Pseudocalliergon lycopodioides. However, suitable habitats for calcareous species were also found in areas where there was not much calcareous rock. This kind of environments may contain for example small outcroppings, springs associated with distant calcareous deposits, or superficial calcareous deposits such as shell-rich sands. For example, Amblyodon dealbatus, Malaxis monophyllos, and Pseudocalliergon lycopodioides grow mainly in herb-rich fens, calcareous springs and wet calcareous rocky outcrops (Ulvinen 2001). The results of this study also revealed that numerous threatened species, such as Philonotis calcarea, Palustriella decipiens, and P. falcata, also occur most commonly in springs which have a neutral pH. Sphagnum contortum is a rich fen species that exists in a rather restricted region within the study area, and most of the variation can be explained by habitat, particularly by the increasing proportion of undrained mires. However, S. contortum does not occur throughout its potential habitat space because it performs best in nutrient- and calcareous-rich habitats and the presence of springs. Consequently, successful models require also other explanatory factors than undrained peatland alone.

The models of decaying wood species performed more weakly than those of other species groups. This may, at least partly, be explained by the fact that the models did not contain dead wood as an explanatory variable. Many decaying wood species are dependent on old-growth forest habitats with high amounts of decaying wood (Rassi et al. 2010). We can therefore expect that old spruce swamp forests are also important habitats for these species, although this could not be directly seen in the models. Nevertheless, there was a high contribution of spruce and pine in the models. Interestingly, topographical wetness index, with mainly negative association, was important in determining suitable habitats for decaying wood species. This kind of topographical variable can serve as proxies for other environmental variables, such as soil properties and plant-available water, which may drive plant distributions (Lassueur et al. 2006). The species associated with decaying wood require a continuing presence of deadwood at various stages of decay, as well as evenly moist microclimates and shady growth sites (Laaka-Lindberg et al. 2009). Despite favoring moisture, these species do not tolerate wet conditions, which prevail at low elevated sites where water accumulates.

Richness patterns and hotpots of species groups

The predicted species richness maps and the location of the most species-rich hotspots indicated important differences between species groups. High potential species richness occurred for rich fen species in northern parts of the study area, for mesotrophic fen species in southwestern parts, and for decaying wood species broadly throughout the central-southeastern part of the study area. These differences reflect the differing habitat requirements among the species groups. However, part of the differences may also arise from unbalanced species numbers within each group: mesotrophic fen (6), decaying wood (7), rich fen (19), and calcareous species (23).

Reliable identification of hotspot areas with a high number of potentially suitable habitats for threatened species has a central role in land use and conservation planning. The grouped species approach in the landscape makes the identification of potentially high-quality habitats for rare species more reliable and the argument for sympathetic management of these habitats more compelling. The species richness analysis indicates that management efforts, such as restoration, would provide the most benefit in northern parts for all threatened mire species as a whole. One reason may be that the amount of drainage decreases generally towards north (Finnish Forest Research Institute 2014), and the degradation of peatland habitats has not proceeded as intensively as in the heavily drained habitats further in the south.

It must be noted that the hotspots of threatened species are not the only habitats important for biodiversity. Complementary approaches are needed, whereby typicalness indicates that abundant habitats and species at the centre of their natural range are also important (Latimer 2009). These typical areas also need to be maintained and actions taken to mitigate, slow down or prevent their degradation. Our approach can also be used to present peatlands’ biodiversity “non-hotspots”. They may be either areas which are important because they hold characteristic assemblages of peatland species, or drained peatlands where forestry practices have degraded all but these ecologically valuable patches. Such areas are also important within the design of land management planning, before large scale peatland re-use is planned and carried out.

Management of specific areas on behalf of one species group may not equally benefit all species within the overall species assemblage, since each species has its own habitat requirements. In contrast, summed richness maps can be readily divided into different sub-categories, enabling land use planners to scrutinize the predictions for species with, for example, different endangerment status or species with different characteristics such as vascular plants and bryophytes. Species groups can be adaptively used to address the needs of both individual species and groups of species by first setting management targets for a group, then testing the benefits of that management for individual species, and thereafter adjusting management direction to best benefit all species of concern (Wisdom et al. 2001; Wiens et al. 2008).

Conclusions

Our results demonstrate that the combination of individual species models offers a useful tool for identifying landscape-scale richness patterns for threatened mire plant species. In conclusion, habitat suitability models can help in determining which aspects of the environment of a given species have a critical impact on its distribution, and thus advance our understanding of the ecological requirements of species, while also providing valuable information concerning where species are likely to be found in insufficiently surveyed landscapes. The modelling performance was high across the modelled species, and the richness patterns generated by single models coincide with the expected richness pattern based on expert knowledge. These generated richness patterns therefore offer a powerful tool for basic biodiversity applications (e.g., land use planning and conservation). Predictive habitat suitability models and the summed richness maps can provide a valuable means of delimiting potentially valuable geographic areas and focus survey and management efforts onto valuable geographic areas and focus such efforts towards ensuring the preservation of biological diversity in aapa mire landscapes. Thus, when examining larger landscape sites suitable for different land use, the models created in this study offer considerable scope for use as “first filters” for identifying potential locations of hotspots of threatened species in boreal peatland landscapes at the regional scale. It is, however, important to emphasize that areas should not be valued simply on the basis of model predictions of threatened species. Typical species and habitats are also important for the biodiversity. It is also vital that modelled valuations should be subject to ground-truth assessment before planning decision-making takes place.

Notes

Acknowledgements

Open access funding provided by Natural Resources Institute Finland (LUKE). A study of this nature would not have been possible without the hundreds of volunteers who contributed their data to the red-listed plant species database. We are thankful for two anonymous referees whose comments and suggestions greatly improved the manuscript. The study is part of the EU LIFE + project LIFEPeatLandUse (LIFE12 ENV/FI/000150). We also thank the Maj and Tor Nessling Foundation for providing a personal scholarship to M. Parkkari.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Copyright information

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.