Assessing the conservation status of European orchids in Germany

by Melanie H

Summary

The IUCN Red List inventory gives information about the conservation status of species worldwide and acts as an important tool for conservation management. For updating this database it is necessary to provide current assessment data of as many species as possible. With the R package ConR it is possible to calculate certain assessment parameters for multiple species for preliminary IUCN conservation assessments based on georeferenced distribution data. Based on free distribution data I investigate the conservation status of eight European orchid species in Germany. These species occur only on nutrient-poor grasslands and are threatened by land use change and climate change. Anacamptis pyramidalis, Himantoglossum hircinum, Ophrys apifera and Ophrys insectifera show positive trends with increasing occurrences and higher extent of occurrence as well as area of occupancy. This is probably mostly due to more occurrences in protected areas and involved conservation measures. However, the species Dactylorhiza viridis, Limodorum abortivum and Spiranthes spiralis show fragile or negative developments indicating a need for prioritized conservation efforts. Overall, the easy application of IUCN evaluation parameters with distribution data allows the assessment of occurrence trends over time and gives an estimate for prioritizing conservation efforts for different species.

Introduction

The IUCN Red List of Threatened Species provides an inventory of the conservation status of species worldwide representing an important tool for assessing biological diversity, extinction risks and conservation needs (IUCN 2019). Dauby et al. (2017) developed the ConR R package which allows preliminary conservation assessment of many species simultaneously with distribution data of species. This tool simplifies further Red List assessments and conservation efforts and gives an example for using geospatial objects in conservation management. In this project I use freely accessible distribution data and the ConR package to evaluate the conservation status of several plant species. I focus on the conservation status of some European orchids in Germany which occur only on nutrient-poor semi-natural grasslands. These habitats are threatened by land use change such as intensification, fragmentation and abandonment (Cousins & Eriksson 2008, Pykälä et al. 2005, Strijker 2005) and are also jeopardized by climate change (Araújo et al. 2011). Thus, many orchid species especially in Central Europe declined over the past decades and are listed in many national Red Lists (Kull et al. 2016). Conservation efforts rely on updated knowledge about current extinction risks as a whole or in a certain area. Consequently, the protection of most threatened species indicated by an increased extinction risk needs to receive priority. On the other hand, a decrease in extinction risk enables to allocate funds and efforts efficiently to the most critically endangered species. The aim of this project is therefore to assess the conservation status of the orchid species Himantoglossum hircinum, Orchis ustulata, Anacamptis pyramidalis, Limodorum abortivum, Dactylorhiza viridis, Ophrys insectifera, Ophrys apifera and Spiranthes spiralis in Germany.

Data and Methods

The distribution data of the chosen taxa used for this analysis comes from the GBIF database providing open access biodiversity data funded by the world’s governments (GBIF 2019). After registering on the website it is possible to download observation data of selected species and regions which is available as a .zip-file after a few minutes. The dataset I used is available here:

http://api.gbif.org/v1/occurrence/download/request/0042983-181108115102211.zip (GBIF.org (23rd January 2019) GBIF Occurrence Download https://doi.org/10.15468/dl.jfuurp)

After unzipping I read the data in R using the R package rgbif (Chamberlain et al. 2018). To analyze the extinction risk of the chosen orchid species, I used the IUCN criteria and the ConR package in R for a preliminary IUCN assessment and evaluation of endangerment based on occurrence data (Dauby et al. 2017). The IUCN extinction risk of species is assessed in nine categories: NE (Not Evaluated), DD (Data Deficient), LC (Least Concern), NT (Near Threatened), VU (Vulnerable), EN (Endangered), CR (Critically Endangered), EW (Extinct in the Wild), EX (Extinct) (IUCN 2019a). VU, EN and CR are categories described as threatened. Quantitative criteria (A-E) are used for classifying a taxon into a threatened category (IUCN Species Survival Commission 2012). These criteria include population size (i.e. total number of mature individuals), subpopulations (geographically distinct groups in the population), reduction (decline in number), extreme fluctuations, extent of occurence and others (IUCN Species Survival Commission 2012). Criterion B incoporates geographical components and it is particularly widely used for evaluation because it allows the assessment by using georeferenced distribution data of a taxon (Schatz 2002). This is important because more detailed data (e.g. abundance data) is often not available. Criterion B has the two subcriteria B1 (extent of occurrence (EOO)) and B2 (area of occupancy (AOO)) and three additional conditions (a, b and c) considering potential decline (Dauby et al. 2017). ConR can calculate the subcriteria B1 and B2 and give an estimate of the number of locations (condition a) (Dauby et al. 2017). Moreover, subpopulations can be calculated as required for the submission of assessments to the IUCN Red List (Dauby et al. 2017). Measures of continuing decline are not directly implemented, but I relate to that by comparing development of these aspects over the last 20-30 years. First I will look at numbers of occurrences in 5-year intervals and over the whole time period where data exists and then I will asses the endangerment of the chosen orchids by comparing developments between 1998, 2008 and 2018 because 10 years sufficiently reflect the life cycle if many orchids. For displaying occurrence in the world I used the R package rnaturalearth for mapping background (South 2017). I used R Studio with R Markdown for presenting code as well as visualizations and text.

Results

# Install and load necessary packages
library(alphahull)
library(sf)
library(sp)
library(dplyr)
library(rnaturalearth)
library(ggplot2)
library(broom)
library(ggiraph)
library(tidyr)
library(raster)
library(tidyverse)
library(ConR)
library(rgbif)
library(gganimate)
library(gifski)
library(png)
library(maptools)
library(readxl)
library(magick)

Data overview

After downloading and unzipping the data, I read them in R and looked at its contents.

orch_gbif <- occ_download_get(key = "0042983-181108115102211", overwrite = TRUE) %>% 
    occ_download_import(orch_gbif_download, na.strings = c("", NA))

orch_gbif
orch_gbif %>% count(genus, sort = TRUE) # gives number of occurrences of individuals of each genus
## Warning: package 'bindrcpp' was built under R version 3.5.2
## # A tibble: 7 x 2
##   genus              n
##   <chr>          <int>
## 1 Ophrys          9211
## 2 Neotinea        4246
## 3 Dactylorhiza    3854
## 4 Himantoglossum  2068
## 5 Spiranthes      2007
## 6 Anacamptis      1922
## 7 Limodorum         68
orch_gbif %>% count(species, sort = TRUE) # gives number of occurrences (68 to 5487) of individuals of each species (8)
## # A tibble: 8 x 2
##   species                     n
##   <chr>                   <int>
## 1 Ophrys insectifera       5487
## 2 Neotinea ustulata        4246
## 3 Dactylorhiza viridis     3854
## 4 Ophrys apifera           3724
## 5 Himantoglossum hircinum  2068
## 6 Spiranthes spiralis      2007
## 7 Anacamptis pyramidalis   1922
## 8 Limodorum abortivum        68
# or another way
# orch_gbif %>% filter(taxonRank == "SPECIES") %>% count(species) %>% arrange(desc(n))
orch_gbif %>% filter(species == "Spiranthes spiralis") # see the number of rows = occurences of S. spiralis
## # A tibble: 2,007 x 45
##    gbifID datasetKey occurrenceID kingdom phylum class order family genus
##     <int> <chr>      <chr>        <chr>   <chr>  <chr> <chr> <chr>  <chr>
##  1 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
##  2 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
##  3 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
##  4 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
##  5 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
##  6 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
##  7 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
##  8 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
##  9 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
## 10 1.95e9 64dabd3c-~ http://id.s~ Plantae Trach~ Lili~ Aspa~ Orchi~ Spir~
## # ... with 1,997 more rows, and 36 more variables: species <chr>,
## #   infraspecificEpithet <chr>, taxonRank <chr>, scientificName <chr>,
## #   countryCode <chr>, locality <chr>, publishingOrgKey <chr>,
## #   decimalLatitude <dbl>, decimalLongitude <dbl>,
## #   coordinateUncertaintyInMeters <dbl>, coordinatePrecision <dbl>,
## #   elevation <dbl>, elevationAccuracy <dbl>, depth <lgl>,
## #   depthAccuracy <lgl>, eventDate <chr>, day <int>, month <int>,
## #   year <int>, taxonKey <int>, speciesKey <int>, basisOfRecord <chr>,
## #   institutionCode <chr>, collectionCode <chr>, catalogNumber <chr>,
## #   recordNumber <chr>, identifiedBy <chr>, dateIdentified <chr>,
## #   license <chr>, rightsHolder <chr>, recordedBy <chr>, typeStatus <chr>,
## #   establishmentMeans <lgl>, lastInterpreted <chr>, mediaType <chr>,
## #   issue <chr>

Plot an overview

orch_gbif %>% count(species, sort = TRUE)  %>% filter(n > 50) %>% 
    ggplot(aes(x = reorder(species, n), y = n)) + geom_bar(stat = "identity", 
    show.legend = FALSE) + labs(x = "Species", y = "Number of Occurrence Records (observations)") + 
    coord_flip()+
  theme_minimal()

This graph shows the total number of observations of each selected orchid species in Germany between 1700 and 2019. Limodorum abortivum has the fewest recorded occurrences and Ophrys insectifera the most.

orch_2015 <- orch_gbif %>% filter(year == "2018")  %>% count(species, sort = TRUE) 
orch_2010 <- orch_gbif %>% filter(year == "2010")  %>% count(species, sort = TRUE)
orch_2005 <- orch_gbif %>% filter(year == "2005")  %>% count(species, sort = TRUE)
orch_2000 <- orch_gbif %>% filter(year == "2000")  %>% count(species, sort = TRUE)

# now some steps to create a dataframe with counts of occurrences and years of interest
# create year columns
orch_2000$year <- 2000
orch_2005$year <- 2005
orch_2010$year <- 2010
orch_2015$year <- 2015
# save as vectors
species <- c(orch_2000$species, orch_2005$species, orch_2010$species, orch_2015$species)
n <- c(orch_2000$n, orch_2005$n, orch_2010$n, orch_2015$n)
years <- c(orch_2000$year, orch_2005$year, orch_2010$year, orch_2015$year)
# put vectors in new dataframe
orch_combined <- data.frame(species, years, n)
str(orch_combined)
## 'data.frame':    30 obs. of  3 variables:
##  $ species: Factor w/ 8 levels "Anacamptis pyramidalis",..: 2 7 5 8 6 3 1 2 5 7 ...
##  $ years  : num  2000 2000 2000 2000 2000 ...
##  $ n      : int  33 31 30 20 18 6 5 71 34 27 ...
orch_combined$years <- as.integer(orch_combined$years)

# Plot
ggplot(orch_combined, aes(species, n, fill=n)) +
  geom_bar(colour="black", stat="identity")+
  coord_flip()+
  guides(fill=FALSE)+
  labs(title = 'Species occurrences between 2000 and 2015 \nYear: {closest_state}',
       subtitle=" ", x = 'Species', y ='Number of occurrences') +
  gganimate::transition_states(years,
                    transition_length = 2,
                    state_length = 1)+
  ease_aes('cubic-in-out') # Slow start and end for a smoother look

The comparison of occurrence records between the years 2000, 2005, 2010 and 2015 gives indication for changes in observation which might correspond to the frequency of the species in the wild. According to this, occurrences of Dactylorhiza viridis increased between 2000 and 2005 but declined again afterwards. Ophrys apifera, Anacamptis pyramidalis and Himantoglossum hircinum increased between 2005 and 2010, the latter increasing further in 2015. Limodorum abortivum was not observed in 2000 and 2005 but has few occurrences in 2010 and 2015. However, one has to bear in mind that especially orchids experience large population fluctuations between years depending e.g. on weather conditions (Pfeifer et al. 2006). Moreover, the occurrence data depends on recorded observations which must not necessarily reflect actual population developments but can also be due to different observation and registration efforts. A case study in a german nature reserve documented an actual increase of population size in Himantoglossum hircinum between 1975 and 2000 so that a further increase seems realistic. Also another Anacamptis species showed positive population trends between 1977 and 2010 (Hornemann et al. 2012), supporting the presented occurrence data. Number of Spiranthes spiralis individuals in a population in the Netherlands showed a decreasing trend until 2005 with an estimated time until extinction of 40 years, suggesting further endangerment (Jaquemyn et al. 2007). Especially the decreasing trend of Dactylorhiza viridis is concerning because as a short-living species it is particularly prone to extinction (Jaquemyn et al. 2005). Limodorum abortivum has been found to lose suitable niches in different climate change scenarios suggesting negative population developments (Kolanowska et al. 2017). Thus, the additional occurrences in the more recent years might be due to registered monitoring recordings but not necissarily due to establishment of this species.

orch_agg <- orch_gbif %>% filter(year %in% 1930:2018) %>% group_by(year ) %>% count(species, sort = TRUE)  %>% drop_na()

colvec <- c("#6161b1","#489b41", "#8fff48", "#1344d4", "#d49d13", "#d62d0a", "#17d6e4", "black") # choose hex colors for plotting

ggplot(orch_agg, aes(x = year, y = n, colour = species))+
  geom_line(size=1)+
  scale_color_manual(values=colvec)+
  labs(x = "Year", y = "Number of occurrences", title = "Occurrences of orchid species between 1900 and 2018") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        panel.background = element_rect(fill = "white", colour = "grey50"))

The patterns of observations over the last 120 years show an inreasing temporal resolution of data which supports the assumption that in earlier days less ocurrences were recorded or old records are not contained in the available data. However, there might have been a thorough (herbarium) investigation in the early 1900s presenting many occurrences of Ophrys insectifera and Neotinea ustulata. More frequent species like those show vast fluctuations over the years with no particular long-term trend. Himantoglossum hircinum, Anacamptis pyramidalis and Ophrys apifera show an overall increasing trend since the early 2000s with potential declines lately. Spiranthes spiralis and Limodorum abortivum have always been rare in Germany. Both are also distributed in southern Europe reducing their overall extinction risk, but german populations are important to preserve their genetic diversity and area of distribution.

IUCN evaluation with ConR

The following code refers to the ConR vignette (Dauby 2018).

Computing EOO (Extent of Ocurrence), subcriterion B1

To compare the EOO between 1998, 2008 and 2018, I prepared temporal subsets of the data and then calcualted the specific criteria.

# create dataset with all species and required structure
# ConR requires a data structure with 3 columns: 
#     ddlat: Latitude in decimal degrees
#     ddlon: Longitude in decimal degrees
#     Name of taxa
# field names do not matter but order of columns is mandatory.
orchids <- orch_gbif %>% dplyr::select(decimalLatitude, decimalLongitude, species, year)  # select the three necessary + one mandatory column
names(orchids)[4] <- "coly" # rename mandatory column

# create subsets 
orchids1998 <- orchids %>% dplyr::filter(coly == "1998")
orchids2008 <- orchids %>% dplyr::filter(coly == "2008")
orchids2018 <- orchids %>% dplyr::filter(coly == "2018")
# compute EOO for all species

EOO.computing(orchids1998, file.name = "orchids1998_EOO") # gives a csv file with EOO values for all species in 1998
##                           EOO
## Anacamptis pyramidalis   2330
## Dactylorhiza viridis     9034
## Himantoglossum hircinum   753
## Neotinea ustulata       62152
## Ophrys apifera          19999
## Ophrys insectifera      27984
## Spiranthes spiralis       573
EOO.computing(orchids2008, file.name = "orchids2008_EOO") # gives a csv file with EOO values for all species in 2008
## Warning in .EOO.comp(x, Name_Sp = ifelse(ncol(XY) > 2,
## as.character(unique(x$tax)), : EOO for Limodorum abortivum is not computed
## because there is only 1 unique occurrence
##                            EOO
## Anacamptis pyramidalis   62044
## Dactylorhiza viridis     27343
## Himantoglossum hircinum  43757
## Limodorum abortivum         NA
## Neotinea ustulata       107387
## Ophrys apifera          100120
## Ophrys insectifera      129818
## Spiranthes spiralis      59054
EOO.computing(orchids2018, file.name = "orchids2018_EOO") # gives a csv file with EOO values for all species in 2018
##                            EOO
## Anacamptis pyramidalis  134808
## Dactylorhiza viridis     33423
## Himantoglossum hircinum 123624
## Limodorum abortivum         10
## Neotinea ustulata       124610
## Ophrys apifera          141712
## Ophrys insectifera      149622
## Spiranthes spiralis      10509
# for single species if preferred
#EOO.computing(spiralis, file.name = "spiranthes_EOO") # gives the EOO of all occurrences of Spiranthes as output and as csv fil

Looking only at the EOO of the species in 2018, Limodorum abortivum would be listed as critically endangered (< 100 sqkm) and Spiranthes spiralis as vulnerable (< 20000 sqkm). 1998 Spiranthes spiralis, Himantoglossum hircinum and Anacamptis pyramidalis were regarded as endangered (< 5000 sqkm) and Dactylorhiza viridis as vulnerable according only to EOO.

Plotting the EOO of individual taxa as a shapefile in the map:

EOO.1998 <- EOO.computing(orchids1998, export_shp = T, alpha=5)
EOO.2008 <- EOO.computing(orchids2008, export_shp = T, alpha=5)
## Warning in .EOO.comp(x, Name_Sp = ifelse(ncol(XY) > 2,
## as.character(unique(x$tax)), : EOO for Limodorum abortivum is not computed
## because there is only 1 unique occurrence
EOO.2018 <- EOO.computing(orchids2018, export_shp = T, alpha=5)

par(mfrow=c(1,3))
sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Anacamptis pyramidalis\n(Convex Hull) 1998")
sp::plot(EOO.1998[[1]][[2]], col="darkred", add=TRUE)  # first species (in alphabetical order)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Anacamptis pyramidalis\n(Convex Hull) 2008")
sp::plot(EOO.2008[[1]][[2]], col="darkred", add=TRUE)  # first species (in alphabetical order)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Anacamptis pyramidalis\n(Convex Hull) 2018")
sp::plot(EOO.2018[[1]][[2]], col="darkred", add=TRUE)  # first species (in alphabetical order)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Dactylorhiza viridis\n(Convex Hull) 1998")
sp::plot(EOO.1998[[2]][[2]], col="darkred", add=TRUE)  # second species

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Dactylorhiza viridis\n(Convex Hull) 2008")
sp::plot(EOO.2008[[2]][[2]], col="darkred", add=TRUE)  # second species

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Dactylorhiza viridis\n(Convex Hull) 2018")
sp::plot(EOO.2018[[2]][[2]], col="darkred", add=TRUE)  # second species

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Himantoglossum hircinum\n(Convex Hull) 1998")
sp::plot(EOO.1998[[3]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Himantoglossum hircinum\n(Convex Hull) 2008")
sp::plot(EOO.2008[[3]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Himantoglossum hircinum\n(Convex Hull) 2018")
sp::plot(EOO.2018[[3]][[2]], col="darkred", add=TRUE)

par(mfrow=c(1,2))
sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Limodorum abortivum\n(Convex Hull) 1998")
sp::plot(EOO.1998[[4]][[2]], col="darkred", add=TRUE)

# sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
#         main="EOO Limodorum abortivum\n(Convex Hull) 2008")  # not enough data for 2008
# sp::plot(EOO.2008[[4]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
        main="EOO Limodorum abortivum\n(Convex Hull) 2018")  # not enough data for 2008
sp::plot(EOO.2018[[4]][[2]], col="darkred", add=TRUE)

par(mfrow=c(1,3))
sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Neotinea ustulata\n(Convex Hull) 1998")
sp::plot(EOO.1998[[5]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Neotinea ustulata\n(Convex Hull) 2008")
sp::plot(EOO.2008[[5]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Neotinea ustulata\n(Convex Hull) 2018")
sp::plot(EOO.2018[[5]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Ophrys apifera\n(Convex Hull) 1998")
sp::plot(EOO.1998[[6]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Ophrys apifera\n(Convex Hull) 2008")
sp::plot(EOO.2008[[6]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Ophrys apifera\n(Convex Hull) 2018")
sp::plot(EOO.2018[[6]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Ophrys insectifera\n(Convex Hull) 1998")
sp::plot(EOO.1998[[7]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Ophrys insectifera\n(Convex Hull) 2008")
sp::plot(EOO.2008[[7]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Ophrys insectifera\n(Convex Hull) 2018")
sp::plot(EOO.2018[[7]][[2]], col="darkred", add=TRUE)

par(mfrow=c(1,2))
# sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
#          main="EOO Spiranthes spiralis\n(Convex Hull) 1998") # not enough data for 1998
# sp::plot(EOO.1998[[8]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Spiranthes spiralis\n(Convex Hull) 2008")
sp::plot(EOO.2008[[8]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Spiranthes spiralis\n(Convex Hull) 2018")
sp::plot(EOO.2018[[8]][[2]], col="darkred", add=TRUE)

Beside the method convex hull there is also the generalization of an alpha hull. This method is less prone to problems with outliers especially if two hulls over a period of time shall be compared (IUCN Standards and Petitions Subcommittee 2017). For example, the EOO is subdivided into several patches if it spans large uninhabited regions so that errors or outliniers in observations at one time as opposed to the other time have less weight (IUCN Standards and Petitions Subcommittee 2017). However, the IUCN guidelines discourage the use of alpha hulls “because disjunctions and outlying occurrences accurately reflect the extent to which a large range size reduces the chance that the entire population of the taxon will be affected by a single threatening process. The risks are spread by the existence of outlying or disjunct occurrences irrespective of whether the EOO encompasses significant areas of unsuitable habitat” (IUCN Standards and Petitions Subcommittee 2017). Inappropriate exclusions thus might lead to underestimated EOOs (IUCN Standards and Petitions Subcommittee 2017). This is why the method of excluding unsuitable areas is also not recommended. For more accurate display of actual EOOs I will exemplary show the alpha hull method for 2018 with the package alphahull (Pateiro-López & Rodríguez-Casal 2016):

EOO.2018_alpha <- EOO.computing(orchids2018, method.range = "alpha.hull",export_shp = T, alpha=5)

par(mfrow=c(1,3))
sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Anacamptis pyramidalis\n(Alpha Hull) 2018")
sp::plot(EOO.2018_alpha[[1]][[2]], col="darkred", add=TRUE)  # first species (in alphabetical order)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Dactylorhiza viridis\n(Alpha Hull) 2018")
sp::plot(EOO.2018_alpha[[2]][[2]], col="darkred", add=TRUE)  # second species

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Himantoglossum hircinum\n(Alpha Hull) 2018")
sp::plot(EOO.2018_alpha[[3]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOOfor Limodorum abortivum\n(Alpha Hull) 2018")
sp::plot(EOO.2018_alpha[[4]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Neotinea ustulata\n(Alpha Hull) 2018")
sp::plot(EOO.2018_alpha[[5]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Ophrys apifera\n(Alpha Hull) 2018")
sp::plot(EOO.2018_alpha[[6]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Ophrys insectifera\n(Alpha Hull) 2018")
sp::plot(EOO.2018_alpha[[7]][[2]], col="darkred", add=TRUE)

sp::plot(ne_countries(country = 'germany'), # background map from rnaturalearth
         main="EOO Spiranthes spiralis\n(Alpha Hull) 2018")
sp::plot(EOO.2018_alpha[[8]][[2]], col="darkred", add=TRUE)

Computing Subpopulations

The approximate number of subpupulations are computed with cicular buffers (Resol_sub_pop gives buffer radius in km) around each observation and overlapping circles are merged to one subpopulation (Dauby et al. 2017).

par(mfrow=c(1,3))

spiralis1998 <- orchids1998 %>% filter(species == "Spiranthes spiralis")
sub.spiralis1998 <- subpop.comp(spiralis1998, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Spiranthes spiralis\nsubpopulations 1998")
sp::plot(sub.spiralis1998[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

spiralis2008 <- orchids2008 %>% filter(species == "Spiranthes spiralis")
sub.spiralis2008 <- subpop.comp(spiralis2008, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Spiranthes spiralis\nsubpopulations 2008")
sp::plot(sub.spiralis2008[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

spiralis2018 <- orchids2018 %>% filter(species == "Spiranthes spiralis")
sub.spiralis2018 <- subpop.comp(spiralis2018, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Spiranthes spiralis\nsubpopulations 2018")
sp::plot(sub.spiralis2018[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

anacamp1998 <- orchids1998 %>% filter(species == "Anacamptis pyramidalis") 
sub.anacamp1998 <- subpop.comp(anacamp1998, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Anacamptis pyramidalis\nsubpopulations 1998")
sp::plot(sub.anacamp1998[["subpop.poly"]], col="darkred", border = "grey",add=TRUE)

anacamp2008 <- orchids2008 %>% filter(species == "Anacamptis pyramidalis") 
sub.anacamp2008 <- subpop.comp(anacamp2008, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Anacamptis pyramidalis\nsubpopulations 2008")
sp::plot(sub.anacamp2008[["subpop.poly"]], col="darkred", border = "grey",add=TRUE)

anacamp2018 <- orchids2018 %>% filter(species == "Anacamptis pyramidalis") 
sub.anacamp2018 <- subpop.comp(anacamp2018, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Anacamptis pyramidalis\nsubpopulations 2018")
sp::plot(sub.anacamp2018[["subpop.poly"]], col="darkred", border = "grey",add=TRUE)

dacty1998 <- orchids1998 %>% filter(species == "Dactylorhiza viridis")
sub.dacty1998 <- subpop.comp(dacty1998, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Dactylorhiza viridis\nsubpopulations 1998")
sp::plot(sub.dacty1998[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

dacty2008 <- orchids2008 %>% filter(species == "Dactylorhiza viridis")
sub.dacty2008 <- subpop.comp(dacty2008, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Dactylorhiza viridis\nsubpopulations 2008")
sp::plot(sub.dacty2008[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

dacty2018 <- orchids2018 %>% filter(species == "Dactylorhiza viridis")
sub.dacty2018 <- subpop.comp(dacty2018, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Dactylorhiza viridis\nsubpopulations 2018")
sp::plot(sub.dacty2018[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

him1998 <- orchids1998 %>% filter(species == "Himantoglossum hircinum")
sub.him1998 <- subpop.comp(him1998, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Himantoglossum hircinum\nsubpopulations 1998")
sp::plot(sub.him1998[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

him2008 <- orchids2008 %>% filter(species == "Himantoglossum hircinum")
sub.him2008 <- subpop.comp(him2008, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Himantoglossum hircinum\nsubpopulations 2008")
sp::plot(sub.him2008[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

him2018 <- orchids2018 %>% filter(species == "Himantoglossum hircinum")
sub.him2018 <- subpop.comp(him2018, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Himantoglossum hircinum\nsubpopulations 2018")
sp::plot(sub.him2018[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

par(mfrow=c(1,2))
# lim1998 <- orchids1998 %>% filter(species == "Limodorum abortivum")
# sub.lim1998 <- subpop.comp(lim1998, Resol_sub_pop=10)  # no data for this year
# sp::plot(ne_countries(country = 'germany'), main="Limodorum abortivum\nsubpopulations 1998")
# sp::plot(sub.lim1998[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

lim2008 <- orchids2008 %>% filter(species == "Limodorum abortivum")
sub.lim2008 <- subpop.comp(lim2008, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Limodorum abortivum\nsubpopulations 2008")
sp::plot(sub.lim2008[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

lim2018 <- orchids2018 %>% filter(species == "Limodorum abortivum")
sub.lim2018 <- subpop.comp(lim2018, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Limodorum abortivum\nsubpopulations 2018")
sp::plot(sub.lim2018[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

par(mfrow=c(1,3))
neo1998 <- orchids1998 %>% filter(species == "Neotinea ustulata")
sub.neo1998 <- subpop.comp(neo1998, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Neotinea ustulata\nsubpopulations 1998")
sp::plot(sub.neo1998[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

neo2008 <- orchids2008 %>% filter(species == "Neotinea ustulata")
sub.neo2008 <- subpop.comp(neo2008, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Neotinea ustulata\nsubpopulations 2008")
sp::plot(sub.neo2008[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

neo2018 <- orchids2018 %>% filter(species == "Neotinea ustulata")
sub.neo2018 <- subpop.comp(neo2018, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Neotinea ustulata\nsubpopulations 2018")
sp::plot(sub.neo2018[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

api1998 <- orchids1998 %>% filter(species == "Ophrys apifera")
sub.api1998 <- subpop.comp(api1998, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Ophrys apifera\nsubpopulations 1998")
sp::plot(sub.api1998[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

api2008 <- orchids2008 %>% filter(species == "Ophrys apifera")
sub.api2008 <- subpop.comp(api2008, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Ophrys apifera\nsubpopulations 2008")
sp::plot(sub.api2008[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

api2018 <- orchids2018 %>% filter(species == "Ophrys apifera")
sub.api2018 <- subpop.comp(api2018, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Ophrys apifera\nsubpopulations 2018")
sp::plot(sub.api2018[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

insecti1998 <- orchids1998 %>% filter(species == "Ophrys insectifera")
sub.insecti1998 <- subpop.comp(insecti1998, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Ophrys insectifera\nsubpopulations 1998")
sp::plot(sub.insecti1998[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

insecti2008 <- orchids2008 %>% filter(species == "Ophrys insectifera")
sub.insecti2008 <- subpop.comp(insecti2008, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Ophrys insectifera\nsubpopulations 2008")
sp::plot(sub.insecti2008[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

insecti2018 <- orchids2018 %>% filter(species == "Ophrys insectifera")
sub.insecti2018 <- subpop.comp(insecti2018, Resol_sub_pop=10) 
sp::plot(ne_countries(country = 'germany'), main="Ophrys insectifera\nsubpopulations 2018")
sp::plot(sub.insecti2018[["subpop.poly"]], col="darkred", border = "grey", add=TRUE)

Spiranthes spiralis, Dactylorhiza viridis and Limordorum abortivum had only very few subpopulations.

Computing the preliminary IUCN assessment for multiple species

The output of the function IUCN.eval is a data frame that provides several paramaters and the preliminary IUCN category. Additionally, it gives a table with the output data in the working directory and a map for each species (Dauby 2018).

Beside the subcriterion B1 (EOO) this function calculates the subcriterion B2 (AOO, area of occupancy). This depends on occupied cells of a grid for a given resolution (default 2 km) (Dauby 2018). The argument Cell_size_AOO enables to modify the cell size, but 2 km are recommended by the IUCN guidelines. Moreover, the number of locations is estimated. This assessment normally requires contextual information about threats because it describes a distinct area where a threat can affect all individuals rapidly (IUCN 2012). In ConR the number of locations is the number of occupied cells of a grid for a given resolution (Dauby et al. 2017). By default, the size of grid cells is 10 km. This seems to be an appropriate estimation for orchid species representing a resolution at which subpopulations are equally affected by a certain threat. Another way is to consider protected areas where the threat is asumed to be lower (see second option below). A shapefile of protected areas in Germany can be downloaded from ProtectedPlanet.

The created maps show the final results of the IUCN assessment. The maps are saved in the working directory. I included them here as GIFs for vivid comparison. Results with protected areas are very similar but less easy to read, only an example is shown below.

# seperate runs for each year, two of them are commented because my computer cannot calculate them all together
# the function gives excel spreadsheets with the IUCN assessment information and an output map

# these calculations take some time!!

IUCN.eval(orchids1998, file_name="IUCN_all1998")
IUCN.eval(orchids2008, file_name="IUCN_all2008")
IUCN.eval(orchids2018, file_name="IUCN_all2018")


# considering protected areas
ger_pa <- rgdal::readOGR(".","WDPA_Feb2019_DEU-shapefile-polygons")
proj4string(ger_pa) # it has the right crs
colnames(ger_pa@data) # select unique ID -> WDPAID

IUCN.eval(orchids1998, file_name="IUCN_all1998PA", Cell_size_locations = 10, protec.areas = ger_pa, ID_shape_PA = "WDPAID", method_protected_area = "no_more_than_one")
IUCN.eval(orchids2008, file_name="IUCN_all2008PA", Cell_size_locations = 10, protec.areas = ger_pa, ID_shape_PA = "WDPAID", method_protected_area = "no_more_than_one")
IUCN.eval(orchids2018, file_name="IUCN_all2018PA", Cell_size_locations = 10, protec.areas = ger_pa, ID_shape_PA = "WDPAID", method_protected_area = "no_more_than_one")
# the following code is not necessary for anaylsis but makes the GIFs for more easy comparson. It depends on what is chosen for the file_name argument in the IUCN.eval function. To reproduce it, it is necessary to create a folder for every species in the working directory, name them according to the elements in spec_list and put the output maps of the IUCN.eval function for the three years accordingly in each of the folders. Then the GIFS are created and can be put into the R markdown document.


spec_list <- c("Anacamptis_pyramidalis", "Dactylorhiza_viridis", "Himantoglossum_hircinum", "Ophrys_apifera", "Neotinea_ustulata","Ophrys_insectifera", "Spiranthes_spiralis","Limodorum_abortivum") # same as folder names


# load each png file into working directory
for(s in seq_along(spec_list)){ # get number of list entry

  spec <- spec_list[s]

  for (i in c("1998", "2008", "2018")){
    fn <- paste0("IUCN_all", i, spec, ".png") # file name e.g. IUCN_all1998Anacamptis_pyramidalis.png (created by IUCN_eval function)

    png.dat <- image_read(paste0(spec,"\\",fn)) # give path and file name

    assign(paste0(spec, ".", i),png.dat) # png images as objects in environment
  }
}

rm(s, spec, i, fn, png.dat) # deletes all temporarily created objects in the environment


# creating the gif
# stack the objects for one location, and create animation


for (i in seq_along(spec_list)){  # get number of list entry
  
  for (y in c("1998", "2008", "2018")){
    objectname <- paste0(spec_list[i], ".", y) # e.g. Anacamptis_pyramidalis.1998
    if (y == "1998"){ 
      spec.objects.1 <- get(objectname)
      spec.objects.1 <- image_annotate(spec.objects.1, "1998", size = 50, color = "black", location = "+1800+150")
    } else if (y == "2008"){
      spec.objects.2 <- get(objectname)
      spec.objects.2 <- image_annotate(spec.objects.2, "2008", size = 50, color = "black", location = "+1800+150")
    } else {
      spec.objects.3 <- get(objectname)
      spec.objects.3 <- image_annotate(spec.objects.3, "2018", size = 50, color = "black", location = "+1800+150")
    }
  }

  spec.objects <- c(spec.objects.1, spec.objects.2, spec.objects.3) 
  rm(spec.objects.1, spec.objects.2, spec.objects.3) 
  img <- image_scale(spec.objects)# package calculates
  spec.objects.ani <- image_animate(img, fps = 0.5, dispose = "previous") # fps velocity
  image_write(spec.objects.ani, paste0(spec_list[i], "_animation.gif")) # output name e.g. Anacamptis_pyramidalis_animation.gif; takes quite a lot of calculation time!!
}

rm(i, y, spec.objects, objectname, img, spec.objects.ani) # deletes all temporarily created objects in the environment
# display GIFs
image_scale(image_read("Anacamptis_pyramidalis_animation.gif"))

image_scale(image_read("Dactylorhiza_viridis_animation.gif"))

image_scale(image_read("Himantoglossum_hircinum_animation.gif"))

image_scale(image_read("Ophrys_apifera_animation.gif"))

image_scale(image_read("Limodorum_abortivum_animation.gif"))

image_scale(image_read("Neotinea_ustulata_animation.gif"))

image_scale(image_read("Ophrys_insectifera_animation.gif"))

image_scale(image_read("Spiranthes_spiralis_animation.gif"))

Example of output map with protected areas for Anacamptis pyramidalis in 1998. Blue dots represent occurrences in protected areas.

The maps reveal interesting developments of the distribution of the species between 1998 and 2018. Anacamptis pyramidalis increased its locations from 4 to 29 and the number of occurrences from 10 to 56 (in 28 protected areas 2018 (3 in 1998)) so that the conservation assessment changed from EN to LC/NT. These are still very few observations but it shows a positive trend. Dactylorhiza viridis seems to be fragile because the conservation status changed from EN to LC/NT to VU indicating a need for monitoring and prioritized conservation efforts. It occupies more protected ares in 2008 and 2018 so that these might be important refuges. Himantoglossum hircinum increased its EOO and AOO distinctly and its locations from 3 to 73. Thus, the conservation status changed from EN to LC/NT. This species does not have priority in conservation efforts because it occupies 44 protected areas so that recent conservation efforts seem fruitful and the species might be relatively safe. Limodorum abortivum has not been recorded in 1998 and with only one occurrence in 2008 in a protected area and with 3 occurrences in 2 protected areas in 2018. It is assessed as CR in 2008 and EN in 2018. The species should be of high priority, but it has also a major area of distribution in the south of Europe which reduces the overall extinction risk. Neotinea ustulata shows a considerable increase from 1998 to 2008 and a slight decrease until 2018 being assessed as LC/NT. It occupied 23 protected areas in 2008 and 14 in 2018 so that the species should be monitored and reasons for losses in protected areas should be sought. Ophrys apifera increased its AOO considerably and occupies 20 more protected areas in 2018 than in 1998. The assessment changed from VU to LC/NT so that this species does not need current priority given the fact that it could be established in several protected areas. Ophrys insectifera shows similar trends. Spiranthes spiralis first increased its occurrences from 3 to 10 and the EOO from 10000 sqkm (2 protected areas) to 60000 sqkm (3 protected areas) but then shows a negative trend to an EOO of only 573 sqkm (2 protected areas). The assessment changed from EN to VU and again to EN. This species has a high priority because of ongoing decline. However, it is a good example for lacking data because I am aware of another occurrence in Germany which is not recorded. This can always conduct misleading results.

For better comparison of the output numbers, the created spreadsheet files (in the working directory, see above) are read in and a graph is produced to show differences between the years in the calculated parameters for every species.

a1 <- read_excel("IUCN_all1998PA.xlsx")
a2 <- read_excel("IUCN_all2008PA.xlsx")
a3 <- read_excel("IUCN_all2018PA.xlsx")
# tidy data
a1$year <- c(rep("1998",7))
a2$year <- c(rep("2008",8))
a3$year <- c(rep("2018",8))
a1$EOO <- a1$EOO/1000 # for better display in bar chart
a2$EOO <- a2$EOO/1000
a3$EOO <- a3$EOO/1000
alldata <- rbind(a1, a2, a3)

# for every single species:
spirdata <- alldata %>% filter(taxa == "Spiranthes spiralis") %>% gather("EOO","AOO", "Nbe_unique_occ.", "Nbe_subPop", "Nbe_loc", "Nbe_loc_PA", "Ratio_occ_within_PA",  key = "parameter", value = "value")

anadata <- alldata %>% filter(taxa == "Anacamptis pyramidalis") %>% gather("EOO","AOO", "Nbe_unique_occ.", "Nbe_subPop", "Nbe_loc", "Nbe_loc_PA", "Ratio_occ_within_PA",  key = "parameter", value = "value")

dactydata <- alldata %>% filter(taxa == "Dactylorhiza viridis") %>% gather("EOO","AOO", "Nbe_unique_occ.", "Nbe_subPop", "Nbe_loc", "Nbe_loc_PA", "Ratio_occ_within_PA",  key = "parameter", value = "value")

himdata <- alldata %>% filter(taxa == "Himantoglossum hircinum") %>% gather("EOO","AOO", "Nbe_unique_occ.", "Nbe_subPop", "Nbe_loc", "Nbe_loc_PA", "Ratio_occ_within_PA",  key = "parameter", value = "value")

neodata <- alldata %>% filter(taxa == "Neotinea ustulata") %>% gather("EOO","AOO", "Nbe_unique_occ.", "Nbe_subPop", "Nbe_loc", "Nbe_loc_PA", "Ratio_occ_within_PA",  key = "parameter", value = "value")

limodata <- alldata %>% filter(taxa == "Limodorum abortivum") %>% gather("EOO","AOO", "Nbe_unique_occ.", "Nbe_subPop", "Nbe_loc", "Nbe_loc_PA", "Ratio_occ_within_PA",  key = "parameter", value = "value")

apiferadata <- alldata %>% filter(taxa == "Ophrys apifera") %>% gather("EOO","AOO", "Nbe_unique_occ.", "Nbe_subPop", "Nbe_loc", "Nbe_loc_PA", "Ratio_occ_within_PA",  key = "parameter", value = "value")

insectidata <- alldata %>% filter(taxa == "Ophrys insectifera") %>% gather("EOO","AOO", "Nbe_unique_occ.", "Nbe_subPop", "Nbe_loc", "Nbe_loc_PA", "Ratio_occ_within_PA",  key = "parameter", value = "value")

spirdata$value <- as.numeric(spirdata$value)
anadata$value <- as.numeric(anadata$value)
dactydata$value <- as.numeric(dactydata$value)
himdata$value <- as.numeric(himdata$value)
neodata$value <- as.numeric(neodata$value)
limodata$value <- as.numeric(limodata$value)
apiferadata$value <- as.numeric(apiferadata$value)
insectidata$value <- as.numeric(insectidata$value)

# plot
mydatasets <- list(spirdata, anadata, dactydata, himdata, neodata, limodata, apiferadata, insectidata)
titles <- c("Spiranthes spiralis", "Anacamptis pyramidalis","Dactylorhiza viridis", "Himantoglossum hircinum","Neotinea ustulata", "Limodorum abortivum", "Ophrys apifera", "Ophrys insectifera")
t <- 1
for (i in mydatasets){
  p <- ggplot(i, aes(fill=year, y=value, x=parameter)) + 
  geom_bar(position="dodge", stat="identity")+
  scale_x_discrete(labels = c(expression(paste("EOO \n[1000 ", km^2,"]")), expression(paste("AOO [",km^2,"]")),'locations',expression(paste("locations \nin prot. areas")), expression(paste("sub- \npopulations")), "occurrences", expression(paste("ratio [%] of \noccurrences \nin prot. area"))))+
  theme_minimal()+
  theme(axis.text.x = element_text(vjust = -0.1),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())
  print(p + ggtitle(titles[t]))  
  t <- t+1
}

## Warning: Removed 1 rows containing missing values (geom_bar).

Conclusions

Overall, some of the orchid species show positive trends mostly due to more occurrences in protected areas. However, the species Dactylorhiza viridis, Limodorum abortivum and Spiranthes spiralis show fragile or negative developments indicating a need for prioritized conservation efforts. They are assessed in the categories VU and EN and thus can be described as threatened. The ConR package proves useful to get a personal overview of current species endangerment and occurrence developments over a time period. Accodung to Dauby et al. (2017) it is also suitable for official preliminary IUCN assessment.

References

Araújo M. B., Alagador D., Cabeza M., Nogués-Bravo D., Thuiller W. (2011): Climate change threatens European conservation areas. Ecology Letters 15(5): 484-492. DOI: https://doi.org/10.1111/j.1461-0248.2011.01610.x

Chamberlain S., Barve V., Mcglinn D., Oldoni D., Geffert L., Ram K. (2018): Package ‘rgbif’. Interface to the Global ‘Biodiversity’ Information Facility API, Version 1.1.0. URL: https://cran.r-project.org/web/packages/rgbif/rgbif.pdf (02.02.19)

Cousins S.A.O., & Eriksson O. (2008): After the hotspots are gone: Land use history and grassland plant species diversity in a strongly transformed agricultural landscape. Applied Vegetation Science, 11(3): 365-374. DOI: https://doi.org/10.3170/2008-7-18480

Dauby G., Stévart T., Droissart V., Cosiaux A., Deblauwe V., Simo-Droissart M., Sosef M.S.M., Lowry P.P. II, Schatz G.E., Gereau R.E., Couvreur T.L.P. (2017): ConR. An R package to assist large-scale multispecies preliminary conservation assessments using distribution data. Ecology and Evolution 7(24): 11292-11303. DOI: 10.1002/ece3.3704

Dauby G. (2018): How to use ConR (for beginners in R). URL: https://cran.r-project.org/web/packages/ConR/vignettes/my-vignette.html (02.02.19)

GBIF, Global Biodiversity Information Facility (2019): What is GBIF? URL: https://www.gbif.org/what-is-gbif (02.02.19)

Hornemann G., Michalski S. G., Durka W. (2012): Short-term fitness and long-term population trendsin the orchidAnacamptis morio. Plant Ecology 213:1583-1595. DOI: 10.1007/s11258-012-0113-6

IUCN, International Union for Conservation of Nature (2019): IUCN Red List of Threatened Species. URL: https://www.iucn.org/resources/conservation-tools/iucn-red-list-threatened-species#RL_importance (02.02.19)

IUCN, International Union for Conservation of Nature (2019a): The IUCN Red List Categories. URL: https://www.iucn.org/resources/conservation-tools/iucn-red-list-threatened-species#RL_categories (02.02.19)

IUCN Species Survival Commission (2012): IUCN Red List categories and criteria, version 3.1, second edition, Gland and Cambridge. URL: https://portals.iucn.org/library/sites/library/files/documents/RL-2001-001-2nd.pdf (02.02.19)

IUCN Standards and Petitions Subcommittee (2017): Guidelines for Using the IUCN Red List Categories and Criteria. Version 13. URL: http://cmsdocs.s3.amazonaws.com/RedListGuidelines.pdf (02.02.19)

Jaquemyn H., Brys R., Hermy M., Willemns J. H. (2005): Does nectar reward affect rarity and extinction probabilities oforchid species? An assessment using historical records fromBelgium and the Netherlands. Biological Conservation 121: 257-263. DOI: 10.1016/j.biocon.2004.05.002

Jaquemyn H., Brys R., Hermy M., Willemns J. H. (2007): Long-term dynamics and population viability in one of thelast populations of the endangeredSpiranthesspiralis(Orchidaceae) in the Netherlands. Biological Conservation 134: 14-21. DOI: 10.1016/j.biocon.2006.07.016

Kolanowska M., Kras M., Lipinska M., Mystkowska K., Szlachetko D. L., Naczk A. M. (2017): Global warming not so harmful for all plants - response of holomycotrophic orchid species for the future climate change. Scientific Reports 7. DOI: https://doi.org/10.1038/s41598-017-13088-7

Kull T., Selgis U., Villoslada Peciña M., Metsare M., Ilves A., Tali K., Sepp K., Kull K., Shefferson R. P. (2016): Factors influencing IUCN threat levels to orchids across Europe on the basis of national red lists. Ecology and Evolution 6(17): 6245-6265. DOI: https://doi.org/10.1002/ece3.2363

Pateiro-López B., Rodríguez-Casal A. (2016): Generalization of the Convex Hull of a Sample of Points in the Plane, Version 2.1. URL: https://cran.r-project.org/web/packages/alphahull/alphahull.pdf (02.02.19)

Pfeifer M., Wiegand K., Heinrich W., Jettschke G. (2006): Long-term demographic fluctuations in an orchid species driven by weather: implications for conservation planning. Journal of Applied Ecology 43(2): 313-324. DOI: https://doi.org/10.1111/j.1365-2664.2006.01148.x

Pykälä J., Luoto M., Heikkinen R. K., & Kontula T. (2005): Plant species richness and persistence of rare plants in abandoned semi-natural grasslands in northern Europe. Basic and Applied Ecology, 6(1): 25-33. DOI: https://doi.org/10.1016/j.baae.2004.10.002

Schatz G. E. (2002): Taxonomy and herbaria in service of plant conservation: Lessons from madagascar’s endemic families. Annals of the Missouri Botanical Garden, 89(2), 145-152. DOI: https://doi.org/10.2307/3298559

South, A. (2017): Package ‘rnaturalearth’. World Map Data from Natural Earth, Version 0.1.0. URL: https://cran.r-project.org/web/packages/rnaturalearth/rnaturalearth.pdf (02.02.19)

Strijker D. (2005). Marginal lands in Europe-causes of decline. Basic and Applied Ecology, 6(2): 99-106. DOI: https://doi.org/10.1016/j.baae.2005.01.001