The Effect of Ammonia Deposition on Butterfly Occupancy in Dutch Nature Reserves

(1)

The Effect of Ammonia Deposition on Butterfly

Occupancy in Dutch Nature Reserves

Bachelor Future Planet Studies

Author

H. (Hajna) J. Tijssen

(12191299)

Supervisors

Dr. S. H. (Henrik) Barmentlo

A. (Anne) G. Uilhoorn MSc.

Dr. R. (Roy) H. A. van Grunsven

Date

30 May 2021

(2)

Abstract

Excessive nitrogen (N) deposition in ecosystems is a cause for concern as it has been found to alter the flora and fauna composition of natural areas. Among the fauna, many butterfly species have been found to be sensitive to increases in N deposition. Changes in the presence of butterfly species within the Netherlands can be detected because an extensive monitoring system of butterflies has been developed. Furthermore, since 2005 the Measuring Ammonia in Nature Network (MAN) monitors the atmospheric concentrations of ammonia (NH3) as a proxy of dry NH3-deposition, which

is the prominent component of N deposition within the Netherlands. Therefore, this study aims to determine the extent of the effect of atmospheric NH3 concentrations on the occupancy of

butterflies within Dutch nature reserves from 2005 to 2018. Occupancy data was used, because this allows the probability of the presence of butterfly species to be studied. I used linear mixed effects models and a generalised linear mixed effects model, followed by a two-way analysis of variance test to discover potential effects of dry NH3-deposition on total butterfly occupancy and species-specific

occupancy. The results showed a small negative impact of NH3 on total occupancy, with the

magnitude and direction of the effect being strongly dependent on species, and therefore species-specific. Additionally, the total occupancy of highly N-preferring species was found to show a positive relation with NH3 concentrations. Furthermore, a differentiation was made between habitats of

different N-sensitivity, where there was found to be an effect of NH3 on sensitive habitats but not on

non-sensitive habitats. These results indicate that the occupancy of butterflies is influenced by dry NH3-deposition, however, that the magnitude and direction of this effect is strongly dependent on

the butterfly species, their relative sensitivity for N and the N-sensitivity of the habitat. Therefore, measures to reduce NH3 emissions close to nature reserves are essential in order to protect sensitive

butterfly species and habitats.

Keywords: Insects; Nitrogen deposition; The Netherlands; Measuring Ammonia in Nature

Network; Nitrogen sensitivity.

Het effect van ammoniakdepositie op de aanwezigheid van vlinders in natuurgebieden

Stikstofdepositie in natuurgebieden in Nederland brengt gevolgen met zich mee voor de

samenstelling van soortgemeenschappen. Omdat vlinders gevoelig zijn voor veranderingen in hun leefomgeving, kunnen ze worden bestudeerd voor het detecteren van gevolgen van

stikstofdepositie. Op de meeste locaties in Nederlands is de meest voorkomende verbinding van stikstof in de depositie ammoniak (NH3), wat voornamelijk afkomstig is van de landbouw. Ammoniak

kan door de werking van zwaartekracht (droog) en door neerslag (nat) deponeren, en in het algemeen is in Nederland het aandeel van droge depositie hoger. Hierom wordt er in dit onderzoek gekeken naar wat het effect van droge ammoniakdepositie is op de aanwezigheid van vlinders in natuurgebieden in Nederland. Hiervoor wordt de verspreidingsdata van de Vlinderstichting uit ‘occupancy-modellen’ vergeleken met luchtconcentraties van ammoniak gemeten door het Meetnet Ammoniak in Natuurgebieden (MAN) van 2005 tot 2018. De resultaten uit het onderzoek geven aan dat ammoniak een zwak negatief effect heeft op de verspreiding van vlindersoorten en dat dit effect erg soortafhankelijk is. Hierom werd er ook een onderscheid gemaakt in het effect van ammoniak op stikstofliefhebbers en stikstofgevoelige vlindersoorten, waarbij een toename in de verspreiding van stikstofliefhebbers met toenemende ammoniakconcentraties is gevonden en een relatieve afname in stikstofgevoelige vlindersoorten. Dit komt meest waarschijnlijk doordat planten die op rijke gronden groeien, planten die op schrale gronden groeien verdringen, waardoor de

waardplantensamenstelling van vlinders veranderd. Tot slot was er ook een afname gevonden in de verspreiding van vlinders in stikstofgevoelige leefgebieden, terwijl er bij niet gevoelige habitatten geen effect van ammoniak was gedetecteerd. Dit onderzoek belicht de noodzaak voor een reductie in ammoniakuitstoten, om diverse gemeenschappen aan vlinders te behouden in de Nederlandse natuur.

(3)

Introduction

Different anthropogenic activities have increased the amount of reactive nitrogen (Nr) in the

environment. Yearly, 210 Tg of non-reactive N2 from the atmosphere is converted to its reactive

forms (e.g. NH3, NOx, N2O), which is more than from all natural processes combined (Fowler et al.,

2013). N is mostly fixed during the Haber-Bosch synthesis of ammonia (NH3)for fertilizers and

feedstock, while NO and NO2 are mostly fixed during transport, industrial processes and biomass

combustion (Erisman et al., 2008; Fowler et al., 2013). From these anthropogenic fixation-processes Nr is emitted to the atmosphere: globally 40 Tg yr-1 of NH3 (mainly from agricultural sources) and 34

Tg yr-1_{of NO}

x (mainly from fossil fuel combustion; Fowler et al., 2013). Because this is more than

what can naturally be converted back to N2, this has resulted in an accumulation of Nr in the

environment (Erisman et al., 2013). Due to which there has been an increase in the deposition of Nr

to the Earth’s surface, including (semi-)natural ecosystems (Bobbink et al., 2010).

Through the additional nutrient input to ecosystems Nr deposition (hereafter: ‘N-deposition’) poses a

threat to biodiversity (Phoenix et al., 2006). The threshold at which N-deposition has an effect on different habitats, has been defined and estimated by the United Nations as the ‘critical loads’ (CL) for N-deposition (Sutton et al., 2008). Especially in temperate and northern Europe CL have been found to be exceeded, with the Netherlands having the highest percent (93% in 2020) of ecosystem area with CL in exceedance for biodiversity in Europe (Hettelingh et al., 2017). Because N is a limiting factor for plants, many terrestrial ecosystems have adapted to low N availability, creating diverse communities in these conditions (Bobbink et al., 1998). When Nr input increases, plants with high

maximum growth rates can outcompete slower growing species, which has led to the replacement of characteristic slow-growing nitrophobic (N-sensitive) plant communities by a few dominant fast-growing nitrophilic (N-preferring) species (Tilman, 1987; Sala et al., 2000; Bobbink et al., 2010). Furthermore, NH3 -deposition changes the nitrogen status of the soil, leading first to soil

eutrophication, whereafter acidification can take place when N is leached out of the soil in the form of NO3- (de Haan et al., 2008; Sutton et al., 2008). Even tough most of the research has concentrated

on the effect on plants, animal communities have also been found to be impacted by N-deposition (Adams, 2003). Although some direct toxic effects have been found on fauna, most of the impacts are through indirect effects, such as changes in the quantity and quality of food resources and host species, and cooler and moister microclimates (Nijssen et al., 2017).

One of the species groups which has been found to be heavily impacted by N-deposition is butterflies. In general, butterflies have been found to be in decline in Europe; the European

Grassland Butterfly indicator showed a 39% decline in abundance from 1990 to 2017 (van Swaay et al., 2019a). In the Netherlands in particular a large decline in butterflies has been found: a study spanning from 1890 to 2017 found an 84% decline for the average abundance of 71 species (van Strien et al., 2019). The intensification of farming, due to habitat fragmentation and pollutant emissions (including N), has been found to be most responsible for this decline in north-western Europe (van Swaay et al., 2019a). Additionally, in the Netherlands it has been found that N-deposition still exceeds the CL of many butterfly communities, and thus has been assigned specifically as one of the main causes for a decrease in butterfly diversity in the country (Wallis de Vries & Van Swaay, 2017; Feest et al., 2014).

As a consequence of N-deposition, butterflies whose host plants favour nutrient-poor conditions have been found to generally decrease in abundance, while those whose host plants favour nutrient-rich environments increase (Weiss, 1999; Öckinger et al., 2006). This is most likely the case because the larvae of butterflies are mono- or oligophagous herbivores, and therefore strongly dependent on the local abundance of host plants, which are shifting to a few nitrophilic species due to N-deposition (Thomas et al., 2001; Lebigre et al., 2018). In addition, the quality of host plants has been found to decline by Nr, leading to higher digestion cost for butterflies (Audusseau et al., 2015; Lebigre et al.,

(5)

microclimate cooling - through enhanced mineralization and biomass accumulation of green plant material with high moisture content - which is unfavourable for the thermophilic butterflies (Wallis de Vries & van Swaay, 2006). Consequently, butterflies adapted to N-poor habitats and its

vegetation, such as heathlands and dunes, have been found to show an especially large decline in the Netherlands (Wallis de Vries & van Swaay, 2013; Feest et al., 2014).

Although the effect of Nr on butterflies has been studied in the Netherlands, the amount of

N-deposition has been based on modelled values (e.g. OPS-model: Feest et al., 2014; RIVM model: Wallis de Vries & van Swaay, 2013) and species-specific N optima (Wallis de Vries & van Swaay, 2017). Since 2005 the Measuring Ammonia in Nature network (MAN) provides measured values for atmospheric NH3 concentrations in protected areas in the Netherlands (RIVM, 2021a), allowing for a

more precise indicator to be used. As atmospheric NH3 concentrations are an indicator for the dry

deposition of NH3, the effect of dry NH3-deposition on butterflies can be investigated, which has not

yet been studied separately. Furthermore, in the Netherlands it has been approximated that 67% of the total N-deposition in 2016 was comprised of NH3-deposition, of which 65% in the form of dry

deposition (Wichink Kruit & van Pul, 2018). Thus, the dry deposition of NH3 is the dominant source of

N-deposition in the Netherlands and therefore an important factor to measure. Additionally, in contrast to NOx, a considerable extent of NH3 is deposited close to the source (40% within the first 5

km), at the beginning almost exclusively in the form of dry deposition (Bobbink et al., 2012;

Rijkswaterstaat, 2021). Thus, this allows the effect of emissions of Nr mostly within the Netherlands

to be studied.

Butterfly occupancy data is an indicator for the diversity of species within a specific area, because it shows the probability per species of its true presence or absence at a site. Occupancy modelling has allowed the correction of biases of opportunistic citizen data, such as detection and geographical bias, and for it to be considered together with standardised data of butterflies (van Strien et al., 2013). Consequently, by combining data from different sources butterfly occupancy data within the Netherlands has been able to become of high quality and coverage, which is why this data will be used for this study.

The aim of this research is to examine the effect of dry NH3-deposition on the occupancy of

butterflies. The focus will be on nature reserves, because especially semi-natural areas have been found to be affected by N-deposition (Bobbink et al., 1998). Therefore, the main research question is: ‘To what extent is butterfly occupancy affected by dry NH3-deposition in Dutch nature reserves?’. It is

expected that the occupancy of butterflies is negatively related with NH3 concentrations. Moreover,

because a difference in the response to N-deposition of N-sensitive and N-preferring butterflies has been found (Öckinger et al., 2006), a sub-question will be: ‘What is the difference between the extent of the effect of NH3 on the occupancy of N-sensitive and N-preferring butterflies?’. With the

hypothesis that N-sensitive butterflies show a decrease while N-preferring butterflies an increase with higher levels of NH3. Moreover, also a difference in N-sensitive and N-non-sensitive habitat

types will be made, resulting in the sub-question: ‘What is the difference between the extent of the effect of NH3-deposition on butterfly occupancy in N-sensitive and N-non-sensitive habitats?’. With

the expectation that in N-sensitive habitats there is a higher decline in the occupancy of butterflies.

Methods and Data

NH3 data

The NH3 data was retrieved from the MAN network, which has been developed by the Dutch

National Institute for Public Health and the Environment (RIVM) and provides monthly mean values of atmospheric NH3 concentrations (RIVM, 2021a). The MAN network has 301 measuring points

located mostly in Natura-2000 regions (Figure 1). NH3 is measured with Gradko diffusion tubes

(Figure 2), which are passive samplers for the long-term monitoring of atmospheric concentrations of gasses (Lolkema et al., 2015). Monthly measurements have been made since 2005, and for this study

(6)

measurements from March 2005 to January 2018 were used in µg m-3_{per month. Therefore, the}

study period corresponds with this timeframe (2005-2018).

Butterfly occupancy data

Butterfly occupancy data was used from the Dutch Butterfly Conservation (DBC). This occupancy data is based on butterfly data from the Dutch National Databank Flora and Fauna (NDFF; www.ndff.nl) and the Dutch Butterfly Monitoring Scheme (DBMS; CBS, 2021). The NDFF contains opportunistic presence-only data on butterfly species in the Netherlands, while the DBMS provides standardised data along transects of the presence and absence of butterfly species (De Vlinderstichting, 2020). The scheme started in 1990 and implements the method which was developed for the British Butterfly Monitoring Scheme (Pollard & Yates, 1993). The occupancy data was calculated using the method of van Strien et al. (2013) and is a mean occupancy data per year from 1990 onwards, with a 1-by-1 km resolution. The occupancy value ranges between 1 and 0, where 1 indicates that the species is certainly present and 0 that the species is not present. The occupancy of 53 species was studied in total per location (Table B1 in Appendix B), which corresponds to all butterfly species currently known to be present in the Netherlands (De Vlinderstichting, 2021b).

Combining the data

In order for the NH3 concentration and butterfly occupancy datasets to be combinable and

comparable, first several data cleaning and expanding steps were necessary (Figure 3). For example, the MAN dataset had to be converted from monthly to yearly values. Additionally, locations were added to the MAN dataset as RD New Amersfoort coordinates in the ArcGIS Pro software (ESRI, 2020), which were then used to download the corresponding grid cells of the occupancy data. This resulted in a total of 271 1-by-1 km grid cells, of which in 23 grids more than one MAN data point was present. After this, the datasets were formatted to be easily usable in R (version 3.6.0; R Core Team, 2019), where they were first explored separately by looking at, for example, their overall distribution and change over the years. Hereafter, the datasets were merged together based on the matching grid cell data, which allowed the data to be compared.

Figure 2. A Gradko diffusion tube, which is a passive

sampler used by the MAN network to measure NH3

concentrations. Air can enter at the bottom of the tube, where it is first filtered to decrease the change of disturbances of NH3 containing particles. At the top NH3

is transformed by an acid to NH4+. After a month, the

proportions are analysed in the laboratory, and the monthly NH3 concentrations are calculated.

Figure 1. Location of the MAN measurement points (red dots) and

(7)

Figure 3. Workflow diagram of the data cleaning, expanding and combining of the NH3 concentrations and butterfly

occupancy datasets.

Statistical analysis NH3 and occupancy

The effect of changes in atmospheric NH3 concentrations on the occupancy of butterflies was tested

using different statistical models on the combined dataset. First, a linear mixed effects model (LME) was fitted by the maximum log likelihood estimation on the combined dataset with NA values omitted. This was done in R (R Core Team, 2019) using the ‘lme’ function from the ‘nlme’ package (Pinheiro et al., 2021). For the model as the response variable yearly mean butterfly occupancy was used, while the predictor variables were NH3 concentration, species and year and as a random

variable location was included (model 1). The random variable was added in order to avoid pseudoreplication from the sampling of the same location repeatedly and to control unmeasured differences between sampling units (Walker, 2018). After the model, a type-II analysis of variance test (ANOVA) was run with the ‘anova.lme’ function from the ‘nlme’ package (Pinheiro et al., 2021). Furthermore, with the ‘ggpredict’ function from the ‘ggeffects’ package, the predicted values for the occupancy variable were computed with NH3 as the term (Lüdecke et al., 2020). After which the

coefficient of determination (R²) was calculated with the ‘r2’ function from the ‘performance’ package (Lüdecke, 2021).

(8)

Subsequently, a generalised linear mixed model (GLMM) with a binomial error distribution and Laplace approximation was performed using the same model (model 1). A GLMM was applied, because this allowed for a binomial distribution to be used, which was not possible for the LME, while the occupancy data did seem to show a quasibinomial distribution because of the high amount of 0 and 1 values present (Figure D1 in Appendix D). This was done with the ‘glmmTMB’ function from the ‘glmmTMB’ package (Magnusson et al., 2020). After the GLMM a type-II ANOVA was run with the ‘Anova’ function from the ‘car’ package (Fox et al., 2020), in order to calculate the

likelihood-ratio chi-squares and p-values for the different predictor variables and their interaction. This test is calculated according to the principle of marginality, which means that each term is tested after all others, with the use of a Wald test for GLMM’s (Fox et al., 2020).

Ellenberg nitrogen values

To the butterfly species present in the occupancy data Ellenberg Nitrogen Values (ENV’s) were added from Oostermijer & van Swaay (1998) and Feest et al. (2014). The ENV of butterfly species gives an indication of their relative sensitivity to N. The indicator value ranges from 1 to 8, where 1 indicates a very N-poor habitat preferring species, while 8 a N-rich habitat preferring species. From the 53 butterfly species studied, to 36 species an ENV could be added (Figure 4). A Tukey’s honest significant test was used to see whether the occupancy of species with an ENV value of 1 and 8 significantly differ from each other, with the ‘TukeyHSD’ function from the ‘stats’ package (R Core Team, 2019). Hereafter, a LME with model (1) was fit for the species with an ENV of 1 and 8 separately, in order to identify if there is a difference in the response of highly N-sensitive and N-preferring butterfly species their occupancy to NH3. After the LME, a type-II ANOVA test was applied.

Figure 4. Ellenberg Nitrogen Value (ENV) per butterfly species from Oostermijer & van Swaay (1998) and Feest et al.

(2014), with the species omitted for which no ENV’s were available. The ENV ranges between 1 and 8, where a low ENV indicates a species which prefers nutrient-poor situations, while a high ENV species which prefer nutrient-rich habitats.

Habitat nitrogen-sensitivity

To MAN locations Nitrogen Sensitivity Values (NSV’s) of the habitats present were added from the Aerius Monitor (RIVM, 2021b). Every MAN measurement point was located within the N-sensitivity map of the monitor, and in ArcGIS the NSV’s were added. The values added indicate the following:

1. Very N-sensitive: CL below 1400 mol ha-1_yr-1

2. N-sensitive: CL between 1400 and 2450 mol ha-1_yr-1

(9)

The monitor only considers Natura-2000 regions, therefore locations of MAN outside of these areas were set to not available (NA). After the NSV’s were added per location, this data was merged with the combined dataset of occupancy and NH3. Whereafter, a TukeyHSD test was used to determine

whether the mean occupancy of butterflies significantly differed between the habitats. This was followed by LME’s with model (1) separately for subsets of the data with different NSV’s, after which LME ANOVA tests were used.

Spatial autocorrelation

In order to test if there is a spatial correlation between the butterfly occupancy data or the NH3

concentrations, a Moran’s I test was conducted (Bivand & Wong, 2018). The data was first converted to a spatial dataframe, whereafter neighbours were calculated which were supplemented with spatial weights (Mendez, 2020). Consequently, with the ‘moran.test’ function from the ‘spdep’ package the test was run (Bivand et al., 2021). In order to also check for local spatial autocorrelation within the datasets, the ‘localmoran’ function was used separately for the occupancy and NH3 data as

well (Bivand et al., 2021). This data was then exported to a CSV and visualised on a map with scaled values in ArcGIS Pro (ESRI, 2020).

Results

NH3 concentrations (2005-2018)

The NH3 concentrations in the study area have a range of 0.37-20.04 µg m-3 over the study period

(2005-2018). From the violin plot per year (Figure 11) no clear increase or decrease is visible over the years. However, when comparing different years there are statistically significant differences in the mean NH3 concentrations (Table C1. in Appendix C). For example, the year 2016 shows a significantly

higher mean NH3 concentration compared to all the other years. Furthermore, the year 2017 has a

significantly higher mean than the other years, except for 2005 (p=.667). While the year 2018 shows a significantly lower mean NH3 concentration compared to all the other years (Figure C1 in Appendix

C), however, for this year only NH3 concentrations for the month January were used. When looking at

the data per location, in 55% of the locations an increase in NH3 concentrations was detected over

the years, while in 45% a decrease (Figure C1 in Appendix C).

Figure 5. Violin plot of the NH3concentrations (µg m-3) over the study period (2005-2018). The boxplots show the median

(horizontal line), mean (rhombus) and interquartile ranges (box) per year, which are accompanied with a probability density of the data per year smoothed by a kernel density estimator. Outliers (points) are data points further than 1.5* IQR (IQR=third quartile – first quartile) from the box. The first year (2005) is based on the months March to December, while the last year (2018) only on measurements from January.

(10)

Butterfly occupancy (2005-2018)

Only 13 butterfly species have an occupancy mean for all locations and years above 0.5, while the other 39 species have a lower mean occupancy (Figure 6). The most common species in the study area is the red admiral, while the least common species is the cranberry fritillary. The averaged occupancy for all locations does not show a strong increase or decrease over the years (Figure D2 in Appendix D). However, there are significant differences in the mean occupancy of butterflies between individual years present (Table D1 in Appendix D)

Figure 6. Butterfly occupancy boxplot for all studied butterfly species based on yearly values from 2005 to 2018 ordered by

mean values. For each species, the median (horizontal line), mean (rhombus), interquartile ranges (box) are shown. Outliers (points) are data points further than 1.5* IQR (IQR=third quartile – first quartile) from the box.

Effect of NH3 on mean butterfly occupancy

A slight negative linear relation was present between the yearly mean NH3 concentrations and yearly

mean butterfly occupancy. The slope of the relation between NH3 concentrations and occupancy was

-0.00087 (CI=-0.0015 - -0.0003), with an intercept of 0.27 (Figure 7). This effect of NH3 concentration

per year on occupancy was significantly dependent on the different species that were investigated (interaction of all three factors, F(52,136947)=1.59, p=.0043). However, separately NH3 did not show a

significant effect on the occupancy in this model (F(1,136947)=2.62, p=.1058), while species

(F(52,136947)=2006.05, p<.001) and year ((F(1,136947)=6.43, p=.0112) did (Table E1 in Appendix E). Of the

variance in occupancy, 41% could be explained by the fixed effects in the linear model, while 46% could be explained by the fixed and random effects combined (marginal R²=.411, conditional R²=.464).

(11)

When a binomial distribution was used, the NH3 concentrations still showed a slight negative effect

on the occupancy of butterflies (Est=-0.41, se=0.29). However, without considering the different

species and the year, the effect of NH3 on butterfly occupancy is not significant (X2(1, 137429) =0.167,

p=.68; Table 1). When the different butterfly species are considered as well, the effect of NH3 does

show a significant effect on the occupancy (X2

(52, 137429) =604.28, p=<.001). On the other hand, when

species, year as well as NH3 are considered, a non-significant effect was found (X2(676, 137429) =144.62,

p=1). Of the variance in occupancy 70% could be explained by the fixed effects in the binomial model and 75% by the whole model (marginal R²=.411, conditional R²=.751), which is more than the LME model could.

Table 1. Results of the type-II ANOVA for the GLMM with occupancy as response variable. The table shows the predictor

variables and their corresponding Chi-square value (Chisq), the degrees of freedom (Df), the p-value of the Chi-square (Pr(>Chisq) and its significance (Sign.).

Effect of NH3 on butterfly occupancy based on the N-sensitivity of species

The highly sensitive butterfly species have a significantly lower occupancy than the highly N-preferring butterfly species, with the overall means of 0.11 and 0.57 accordingly (Figure F1 & Table F1 in Appendix F). When looking at the effect of NH3 on the occupancy of highly N-preferring species,

a weak positive linear association is present with a slope of 0.0031 (CI=0.001-0.005; Figure 8). Furthermore, the joint effect of NH3, year and species on the occupancy of N-preferring butterflies is

significant (F(6,17853)=2.13, p=.046). However, there was no significant joint effect found for the

occupancy of highly N-sensitive butterflies (F(1,17853)=2.27, p=.131). Moreover, there was an effect of

species found on the occupancy of N-preferring (F(6,17853)=500.83, p<.0001) and N-sensitive

(F(6,17853)=771.73, p<.0001) butterflies (Table F2-3 in Appendix F).

Figure 7. Effect of NH3 on the occupancy of butterflies based on the LME used. The LME prediction (red line) with its

(12)

Effect of NH3 on butterfly occupancy based on habitat N-sensitivity

There is a significant difference between the mean occupancies of the different N-sensitive habitats (Figure G1 & Table G1 in Appendix G). The very N-sensitive habitats have the highest mean

occupancy (0.29), the limitedly N-sensitive the second highest (0.26) and the N-sensitive the lowest (0.23). All the different habitat types based on N-sensitivity show a small negative linear association between the occupancy of butterflies and the NH3 concentrations, with areas outside of the

Natura-2000 regions showing the largest slopes, and limitedly N-sensitive habitats the smallest (Table 2, Figure 9). However, there was no joint effect of NH3, year and species found on the occupancy of

butterflies in the limitedly N-sensitive habitats (F(52,40467)=1.3, p=.07) and the locations outside of the

Natura-2000 regions (F(52,18515)=1.3, p=.25;Table G2 in Appendix G). While for the very N-sensitive

(F(52,55857)=1.89, p=.0001) and N-sensitive habitats(F(52,21475)=1.63, p=.0026) a joint effect was found to

be present.

Table 2. The relation between NH3 and the occupancy of butterflies in habitats with different NSV values. The slopes were

predicted from the LME models with NH3 as the term. The lower and upper confidence interval of the prediction and the

intercept is also shown.

NSV Slope Lower CI Upper CI Intercept

1 -0.000981304 -0.001883 -7.975E-05 0.299691

2 -0.0012508 -0.001719 -0.0007821 0.244926

3 -0.000914132 -0.001948 0.0001196 0.267711

NA -0.001784583 -0.002154 -0.0014148 0.240071

Figure 8 A-B. The effect of NH3 on the butterfly occupancy of A) sensitive butterfly species with an ENV of 1 and B)

N-preferring species with and ENV of 8. The LME prediction (red line) with its corresponding confidence interval is shown on top of the NH3 and occupancy data (points).

(13)

Spatial autocorrelation of the data

Because the Moran’s I test under randomisation for NH3 is significant (p<2.2e-16) and the z-score is

positive (z=97.81), it can be rejected that NH3 is randomly distributed among the features in the

study area (Table H1 in Appendix H). Therefore, the spatial distribution of NH3 values in the dataset is

more spatially clustered than would be expected if underlying spatial processes were random. Because the calculated Moran’s I statistic is positive (I=0.48), the NH3 mean variable is positively

autocorrelated in the Netherlands. The Moran’s I test for occupancy is also significant (p<2.2e-16_{) with}

a positive z-score (z=19.35), due the which the null hypothesis that occupancy values of butterflies are randomly distributed in the study area can also be rejected (Table H1 in Appendix H). However, because the Moran’s I statistic for occupancy lies closer to 0 (I=0.07), the spatial distribution of this data is more random than the occupancy data, but still positively autocorrelated. When looking at the Moran’s I statistic per location, it is visible that the local Moran’s I is not only positively autocorrelated, but that local negative autocorrelation is also present (Figure 10).

Figure 9 A-D. The effect of NH3 on the butterfly occupancy of butterflies in different sensitive habitats. A) very

N-sensitive habitat (NSV=1), B) N-N-sensitive habitat (NSV=2), C) Limitedly N-N-sensitive habitat (NSV=3), D) locations outside of the Natura-2000 areas (NSV=NA). The LME prediction with NH3 as the term its corresponding confidence interval is

shown (red line) on top of the NH3 and occupancy data (points).

C) D)

(14)

Discussion

Effect of NH3 on butterfly occupancy

The aim of this study was to investigate whether butterfly occupancy is affected by dry deposition of NH3. The results of this study did find an effect to be present of atmospheric NH3 concentrations on

the mean occupancy of butterfly species in Dutch nature reserves for the 2005 to 2018 period, when accounting for the different butterfly species. However, when looking at the effect of NH3 on the

occupancy of all the butterfly species combined this effect was weak. This is most likely because the relation between NH3-deposition and butterflies is strongly dependent on the butterfly species itself

and on whether it is susceptible to be impacted by NH3 or not. For example, when examining the

effect of NH3 on the occupancy of the red admiral butterfly, a highly positive relation is present

(Figure E1 in Appendix E), while for the silver-studded blue this relation is negative (Figure E2 in Appendix E). Although the abundance of butterflies has been found to be impacted in the

Netherlands by N-deposition (Wallis de Vries & van Swaay, 2013), the findings of this study show that atmospheric NH3 concentrations (as a proxy for dry deposition) also affect butterflies. This was

expected because Wichink Kruit & van Pul (2018) found that dry NH3-deposition constituted about

36-45% of the total N-deposition from 2005 to 2016 in the Netherlands. The occupancy and abundance of species are usually positively related to each other, when the abundance declines, a decline in the occupied sites of the species is present, and vice versa (Gaston et al., 2000). However, changes in abundance are more common to be detected than changes in the area of occupancy of a species. Therefore, that the modelled occupancy was also found to be impacted by the dry

deposition of NH3 in Dutch nature reserves, is a substantial finding because this shows that even

distributions of species are affected by NH3-deposition, not just a few individuals within a population.

Effect of NH3 on relative N-sensitivity of butterflies

When distinguishing between highly N-sensitive and N-preferring butterfly species, a much lower occupancy of the sensitive butterflies is present. Furthermore, the occupancy of sensitive and N-preferring species is affected by NH3 concentrations, although this is again dependent on the

butterfly species. For the N-preferring species, the effect of NH on overall occupancy is positive,

Figure 10. The Local Moran’s I statistics for A) the NH3 dataset, B) the occupancy dataset. The differently scaled colours of

the locations indicate whether they are negatively or positively autocorrelated, and how random the values are compared to the neighbouring values (the closer to 0 the more random the distribution).

(15)

sensitive species was not significant, is probably because these species are much rarer, wherefor changes in their occupancy are more difficult to identify (Figure F3 in Appendix F). Feest et al. (2014), also found that butterfly populations are becoming more dominated by nitrophilic butterfly species, because an increase in the index of relative N-sensitivity of butterflies is present in the Netherlands, which the DBC also found to be present (van Swaay et al., 2019b). This observed shift in butterfly species has been found to be due to different indirect effects of N-deposition. The most important indirect effect being the shift in vegetation from nitrophobic to nitrophilic plant species as a

consequence of increased Nr inputs (Bobbink et al., 2010). This shift in vegetation changes the host-

and food-plant availability, which most butterfly species are strongly dependent on because their larvae are mono- or oligophagous herbivores (Thomas et al., 2001; Lebigre et al., 2018). Moreover, also the adult butterflies depend on the presence of particular plant species, mostly nectar rich flowers, because they are nectar-feeding insects, of which the composition and nectar quantity and quality has also been found to be impacted by an increase in Nr (Wallis de Vries et al., 2012; Stevens

et al., 2018).

Effect NH3 on butterfly occupancy of different N-sensitive habitats

When differentiating between the N-sensitivity of habitats, NH3 shows a significant negative effect on

the occupancy of butterflies in N-sensitive habitats, when accounting for the species and year, but not in the non-sensitive habitats. Other studies also found butterflies in N-sensitive habitats, especially heathlands, to be vulnerable to an increase in the Nr input through N-deposition in the

Netherlands (Feest et al., 2014; van Strien et al., 2013). In heathlands this is mostly because of a shift from dwarf shrubs (such as Calluna vulgaris) which many butterfly species prefer, to a few dominant grass species (Bobbink et al., 2010). Wallis de Vries & van Swaay (2013) also found a negative effect of N-deposition on butterfly diversity to be present in coastal dunes, Nardus grasslands and wet heathlands; which are all N-sensitive habitats. Therefore, a reduction in NH3-deposition in N-sensitive

habitats is important to protect the butterfly diversity in these habitats Spatial autocorrelation

For both the occupancy values and the NH3 concentrations a spatial autocorrelation is present. For

the NH3 concentrations this is most likely the case because the sources of NH3 also show a spatial

pattern. For example, stables are often present with a higher density in specific areas than in others within the Netherlands (van der Peet et al., 2018). Furthermore, when NH3 is emitted to the

atmosphere it is transported and deposited over large distances (Rijkswaterstaat, 2021), thus a large emission from the same source can influence several MAN measurement points. For the occupancy of butterfly species there is most likely a spatial autocorrelation present, because butterfly diversity is dependent on the presence of specific habitats (which are often larger than 1*1 km grid cells) and on temperature and rain patterns, because they are thermophilous organisms (Wallis de Vries & van Swaay, 2006), which also show spatial patterns within the Netherlands (Daniels et al., 2014).

However, within the statistical models used, the locations were considered as random effects with the purpose to quantify variation among units, and thus accounting for a part of the correlation in space (Bolker et al., 2009).

Possible factors influencing the findings

An important factor which could have influenced the findings of this study, is the presence of a lag-time in the response of butterfly species to elevated levels of NH3 concentrations in the atmosphere.

Butterfly species tend to respond quickly to changes in the quality of their habitat, because they have a relatively short life cycle (minimum of one generation each year) and are dependent on the

presence and quality of specific plant species (Thomas, 2005). However, the impact of N-deposition on the occupancy of plant species can take time, especially for slow growing species, which delays the shift in vegetation which indirectly affects butterflies. Furthermore, because there is also a lag-time in the recovery of an ecosystem after high Nr concentrations, in for example the plant and

(16)

when measured Nr concentrations are already lower than the CL. Further research could look into

whether there indeed is a lag-time in the response of butterfly occupancies to NH3 concentrations.

Moreover, the results in the present study could be distorted by a shifting baseline syndrome

(Clavero, 2014). For example, from 1890 to 2000 butterfly abundance in the Netherlands is estimated to have declined by 65% (van Strien et al., 2013). Additionally, Wallis de Vries and van Swaay (2017) found that most of the increase in the N-indicator in the Netherlands was already present before 1990. Therefore, there might be a weaker effect of NH3 in the studied period because most of the

sensitive butterfly species already disappeared from Dutch nature reserves before 2005. An important assumption to be aware of when looking at the findings of this study, is that

atmospheric NH3 concentrations are only an indicator for the dry deposition of NH3 and it was not

considered at what rate the vegetation and soil present take up NH3. For example, the height of the

vegetation strongly influences the dry deposition rate of NH3, with higher vegetation taking up more

NH3 because of the greater aerodynamic roughness (Bobbink et al., 2012). However, in order to

measure the dry deposition of NH3, expensive and labour-intensive techniques are needed, due to

which the creation of a detailed spatial network as the MAN would not be possible at an acceptable cost (Noordijk et al., 2020). However, the MAN measurements are calibrated with six high accuracy instruments which measure dry-deposition, and models are used to calculate the proportion of the NH3 taken up by the vegetation (Noordijk et al., 2020). Furthermore, the atmospheric NH3

concentrations are indicative for relative differences in dry deposition.

An important further step of the research is to use GLMM’s with binomial distributions for the calculation of the effect of NH3 on the occupancy of N-sensitive species and habitats as well, because

of the quasibinomial distribution of the occupancy data due to the many 1 (present) and (0) absent values. Moreover, for the parameter estimation, instead of the Laplace approximation, the Markov chain Monte Carlo could be used, which is more accurate, but was not used because it is technically more challenging and has a slow computation time (Bolker et al., 2009). Additionally, model order reduction could be applied, to reduce the variance taken up by other variables than NH3 and to

reduce the eigenvalues which were quite large for the models (Rommes, 2007).

The results of this study indicate that further action within the Netherlands to limit the effect of NH3

on nature reserves is necessary. Habitat restoration measures could be applied, but most importantly reductions in NH3 emissions close to protected areas with N-sensitive butterflies and N-sensitive

habitats should urgently be taken.

Conclusions

This research aimed to identify if the occupancy of butterflies is impacted by the dry deposition of NH3. Through the use of atmospheric NH3 concentrations and occupancy data in Dutch nature

reserves from 2005 to 2018, I found a significant negative effect of NH3 on butterfly occupancy. This

effect was found to be highly variable per species, with N-preferring species actually showing a positive relation, most likely due to a shift in the vegetation composition to nitrophilic species. Additionally, habitats which were found by earlier research to be N-sensitive, showed a decline in butterfly species as a consequence of NH3-deposition. Therefore, this study illustrates that dry NH3

-deposition is an important driver for the decline of butterfly species in Dutch nature reserves. This could indicate that the diversity of other insect groups, such as bees or ants, are also in decline, which are all species with extremely important ecological functions. For this reason, a reduction in NH3 emissions, especially close to sensitive species and habitats, is urgently needed.

Acknowledgements

First of all, I would like to thank Henrik Barmentlo for making this thesis project possible, because it was his idea to create this project in the first place, and he has helped a lot throughout the whole

(17)

project. Additionally, I would like to thank Roy van Grunsven for the help with accessing the butterfly occupancy data and his useful comments for the project. Furthermore, Emiel van Loon for his comments about the statistical analysis. Last but not least, I would like to thank the reviewers for their constructive criticism of the proposal and the final thesis (Márti Domokos, Sanna Wessel, Yishi Hu and Pieter Beaufort).

References

Adams, M. B. (2003). Ecological issues related to N deposition to natural ecosystems: research needs.

Environment International, 29(2-3), 189-199. https://doi.org/10.1016/S0160-4120(02)00179-4 Audusseau, H., Kolb, G., & Janz, N. (2015). Plant fertilization interacts with life history: Variation in

stoichiometry and performance in nettle-feeding butterflies. PloS one, 10(5), e0124616. https://doi.org/10.1371/journal.pone.0124616

Bivand, R. S., & Wong, D. W. (2018). Comparing implementations of global and local indicators of spatial

association. Test, 27(3), 716-748. https://doi.org/10.1007/s11749-018-0599-x

Bivand, R. S., Altman, M., Anselin, L., Assunção, R., Berke, O., Bernat, A., … & Westerholt, R. (2021). Package ‘spdep’ – Spatial Dependence: Weighting Schemes, Statistics. R Package, version, 1.1-7. https://github.com/r-spatial/spdep

Bobbink, R., Bal, D., Van Dobben, H. F., Jansen, A. J. M., Nijssen, M., Siepel, H., ... & De Vries, W. (2012). The effects of nitrogen deposition on the structure and functioning of ecosystems. Ecol.

Applications, 39-79. Retrieved on 11 May 2021, from

https://ec.europa.eu/environment/nature/natura2000/platform/documents/part-i-chapter_2_nov-2012_2013-09-10_en.pdf

Bobbink, R., Hicks, K., Galloway, J., Spranger, T., Alkemade, R., Ashmore, M., ... & De Vries, W. (2010). Global assessment of nitrogen deposition effects on terrestrial plant diversity: a synthesis.

Ecological applications, 20(1), 30-59. https://doi.org/10.1890/08-1140.1

Bobbink, R., Hornung, M., & Roelofs, J. G. (1998). The effects of air‐borne nitrogen pollutants on species diversity in natural and semi‐natural European vegetation. Journal of ecology, 86(5), 717-738.

https://doi.org/10.1046/j.1365-2745.1998.8650717.x

Bolker, B. M., Brooks, M. E., Clark, C. J., Geange, S. W., Poulsen, J. R., Stevens, M. H. H., & White, J. S. S. (2009). Generalized linear mixed models: a practical guide for ecology and evolution. Trends

in ecology & evolution, 24(3), 127-135. https://doi.org/10.1016/j.tree.2008.10.008

CBS, (2021). Dutch Butterfly Monitoring Scheme. Centraal Bureau voor de Statistiek. Retrieved on 15

March 2021, from

https://www.cbs.nl/en-gb/news/2019/13/over-80-decline-in-butterflies-since-late-1800s/dutch-butterfly-monitoring-scheme

Clavero, M. (2014). Shifting baselines and the conservation of non‐native species. Conservation Biology,

28(5), 1434-1436. https://doi.org/10.1111/cobi.12266

Daniels, E. E., Lenderink, G., Hutjes, R. W. A., & Holtslag, A. A. M. (2014). Spatial precipitation patterns and trends in The Netherlands during 1951–2009. International journal of climatology, 34(6),

1773-1784. https://doi.org/10.1002/joc.3800

de Haan, B.J., Kros, J., Bobbink, R., van Jaarsveld, J.A., De Vries, W. & Noordijk, H. (2008). Ammoniak in Nederland. Rapport Planbureau voor de leefomgeving, 500125003, Bilthoven. Retrieved on 15

May 2021, from https://www.pbl.nl/sites/default/files/downloads/500125003.pdf

De Vlinderstichting, (2020). Meetnet Vlinders. Retrieved on 12 March 2021, from

https://www.vlinderstichting.nl/wat-wij-doen/meetnetten/meetnet-vlinders/

De Vlinderstichting, (2021a). Gallery of butterfly images [Pictures]. Retrieved on 5 May 2021, from https://www.vlinderstichting.nl/vlinders/overzicht-vlinders/

De Vlinderstichting, (2021b). De Vlinderstichting website. Retrieved on 5 March 2021, from https://www.vlinderstichting.nl/

Erisman, J. W., Galloway, J. N., Seitzinger, S., Bleeker, A., Dise, N. B., Petrescu, A. R., ... & de Vries, W. (2013). Consequences of human modification of the global nitrogen cycle. Philosophical

Transactions of the Royal Society B: Biological Sciences, 368(1621), 20130116.

https://doi.org/10.1098/rstb.2013.0116

Erisman, J. W., Sutton, M. A., Galloway, J., Klimont, Z., & Winiwarter, W. (2008). How a century of ammonia synthesis changed the world. Nature Geoscience, 1(10), 636-639.

https://doi.org/10.1038/ngeo325

ESRI, (2020). ArcGIS Pro software. Version 2.5.0. www.arcgis.com

Feest, A., van Swaay, C., & van Hinsberg, A. (2014). Nitrogen deposition and the reduction of butterfly biodiversity quality in the Netherlands. Ecological Indicators, 39, 115-119.

(18)

Fowler, D., Coyle, M., Skiba, U., Sutton, M. A., Cape, J. N., Reis, S., ... & Voss, M. (2013). The global nitrogen cycle in the twenty-first century. Philosophical Transactions of the Royal Society B:

Biological Sciences, 368(1621), 20130164. https://doi.org/10.1098/rstb.2013.0164

Fox, J., Weisber, S., Price, B., Adler, D., Bates, D., Baud-Bovy, G., …. & Zeileis, A. (2020). Package ‘car’

– Companion to Applied Regression. R Package, Version, 3.0-10.

https://CRAN.R-project.org/package=car

Gaston, K. J., Blackburn, T. M., Greenwood, J. J., Gregory, R. D., Quinn, R. M., & Lawton, J. H. (2000). Abundance–occupancy relationships. Journal of Applied Ecology, 37, 39-59.

https://doi.org/10.1046/j.1365-2664.2000.00485.x

Hettelingh, J. P., Posch, M., & Slootweg, J. (2017). European critical loads: database, biodiversity and ecosystems at risk: CCE Final Report 2017. Coordination Centre for Effects , Bilthoven,

FoosNetherlands, RIVM Report 2017-0155. https://doi.org/10.21945/RIVM-2017-0155

Lebigre, C., Vanderbeken, C., Turlure, C., & Schtickzelle, N. (2018). Host plant nitrogen enrichment has both positive and negative effects on the larval growth of a specialist butterfly. Ecological

Entomology, 43(4), 494-505. https://doi.org/10.1111/een.12525

Lolkema, D. E., Noordijk, H., Stolk, A. P., Hoogerbrugge, R., van Zanten, M. C., & van Pul, W. A. J. (2015). The measuring ammonia in nature (MAN) network in the Netherlands. Biogeosciences,

12(16), 5133-5142. https://doi.org/10.5194/bg-12-5133-2015

Lüdecke, D. (2021). Package ‘performance’ – Assessment of Regression Models Performance. R Package,

Version, 0.7.2. https://easystats.github.io/performance/

Lüdecke, D., Aust, F., Crawley, S., & Ben-Shachar, M. (2020). Package ‘ggeffects’-Create Tidy Data Frames of Marginal Effects for “ggplot” from Model Outputs. R Package, Version, 0.8.0.

https://strengejacke.github.io/ggeffect

Magnusson, A., Skaug, H., Nielsen, A., Berg, C., Kristensen, K., Maechler, M., ... & Brooks, M. M. (2020). Package ‘glmmTMB’ – Generalized Linear Mixed Models using Template Model builder. R

Package, Version 1.0.2.0. https://github.com/glmmTMB/glmmTMB

Mendez C. (2020). Spatial autocorrelation analysis in R. R Studio/RPubs. Retrieved on 24 April 2021,

from https://rpubs.com/quarcs-lab/spatial-autocorrelation

Nijssen, M. E., Wallis de Vries, M. F., & Siepel, H. (2017). Pathways for the effects of increased nitrogen deposition on fauna. Biological Conservation, 212, 423-431.

https://doi.org/10.1016/j.biocon.2017.02.022

Noordijk, H., Braam, M., Rutledge-Jonker, S., Hoogerbrugge, R., Stolk, A. P., & van Pul, W. A. J. (2020). Performance of the MAN ammonia monitoring network in the Netherlands. Atmospheric

Environment, 228, 117400. https://doi.org/10.1016/j.atmosenv.2020.117400

Öckinger, E., Hammarstedt, O., Nilsson, S. G., & Smith, H. G. (2006). The relationship between local extinctions of grassland butterflies and increased soil nitrogen levels. Biological Conservation,

128(4), 564-573. https://doi.org/10.1016/j.biocon.2005.10.024

Oostermeijer, J. G. B., & van Swaay, C. A. M. (1998). The relationship between butterflies and environmental indicator values: a tool for conservation in a changing landscape. Biological

conservation, 86(3), 271-280. https://doi.org/10.1016/S0006-3207(98)00040-8

Phoenix, G. K., Hicks, W. K., Cinderby, S., Kuylenstierna, J. C., Stock, W. D., Dentener, F. J., ... & Ineson, P. (2006). Atmospheric nitrogen deposition in world biodiversity hotspots: the need for a greater global perspective in assessing N deposition impacts. Global Change Biology, 12(3),

470-476. https://doi.org/10.1111/j.1365-2486.2006.01104.x

Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., Heisterkamp, S., Van Willigen, B., & Maintainer, R. (2021). Package ‘nlme’ - Linear and nonlinear mixed effects models. R Package, version,

3.1-152. https://svn.r-project.org/R-packages/trunk/nlme/

Pollard, E., & Yates, T. J. (1994). Monitoring butterflies for ecology and conservation: the British butterfly

monitoring scheme (1st ed.). Chapman & Hall.

R Core Team, (2019). R: A language and environment for statistical computing. R Foundation for

Statistical Computing, Vienna, Austria, Version 3.6.0. https://www.r-project.org/

Rijkswaterstaat, (2021). Ammoniak en veehouderij. Rijkswaterstraat, InfoMil. Retrieved on 16 March

2021, from https://www.infomil.nl/onderwerpen/ruimte/omgevingsthema/ammo/

RIVM, (2021a). Meetresultaten MAN. Ministry of Health, Welfare and Sport, RIVM. Retrieved on 15 March

2021, from https://man.rivm.nl/

RIVM, (2021b). Aerius Monitor. Ministry of Health, Welfare and Sport, RIVM. Retrieved on 25 March

2021, from https://monitor.aerius.nl/

Rommes, J. (2007). Methods for eigenvalue problems with applications in model order reduction

(Doctoral dissertation). University of Utrecht. http://dspace.library.uu.nl/handle/1874/21787

Sala, O. E., Chapin, F. S., Armesto, J. J., Berlow, E., Bloomfield, J., Dirzo, R., ... & Wall, D. H. (2000). Global biodiversity scenarios for the year 2100. Science, 287(5459), 1770-1774.

(19)

Stevens, C. J. (2016). How long do ecosystems take to recover from atmospheric nitrogen deposition?.

Biological Conservation, 200, 160-167. https://doi.org/10.1016/j.biocon.2016.06.005 Stevens, C. J., David, T. I., & Storkey, J. (2018). Atmospheric nitrogen deposition in terrestrial

ecosystems: Its impact on plant communities and consequences across trophic levels. Functional

ecology, 32(7), 1757-1769. https://doi.org/10.1111/1365-2435.13063

Sutton, M. A., Erisman, J. W., Dentener, F., & Möller, D. (2008). Ammonia in the environment: From ancient times to the present. Environmental Pollution, 156(3), 583-604.

https://doi.org/10.1016/j.envpol.2008.03.013

Thomas, J. A. (2005). Monitoring change in the abundance and distribution of insects using butterflies and other indicator groups. Philosophical Transactions of the Royal Society B: Biological

Sciences, 360(1454), 339-357. https://doi.org/10.1098/rstb.2004.1585

Thomas, J. A., Bourn, N. A. D., Clarke, R. T., Stewart, K. E., Simcox, D. J., Pearman, G. S., ... & Goodger, B. (2001). The quality and isolation of habitat patches both determine where

butterflies persist in fragmented landscapes. Proceedings of the Royal Society of London. Series

B: Biological Sciences, 268(1478), 1791-1796. https://doi.org/10.1098/rspb.2001.1693 Tilman, D. (1987). Secondary succession and the pattern of plant dominance along experimental

nitrogen gradients. Ecological monographs, 57(3), 189-214. https://doi.org/10.2307/2937080

van der Peet, G., Leenstra, F., Vermeij I., Bondt, N., Puister, L., van Os, J., (2018). Feiten en cijfers over de Nederlandse veehouderijsectoren 2018. Wageningen Livestock Research, Rapport, 1134.

Retrieved on 26 May 2021, from https://edepot.wur.nl/464128

van Strien, A. J., van Swaay, C. A., & Termaat, T. (2013). Opportunistic citizen science data of animal species produce reliable estimates of distribution trends if analysed with occupancy models.

Journal of Applied Ecology, 50(6), 1450-1458. https://doi.org/10.1111/1365-2664.12158 van Strien, A. J., van Swaay, C. A., van Strien-van Liempt, W. T., Poot, M. J., & Wallis de Vries, M. F.

(2019). Over a century of data reveal more than 80% decline in butterflies in the Netherlands.

Biological Conservation, 234, 116-122. https://doi.org/10.1016/j.biocon.2019.03.023 Van Swaay, C. A. M. (2019). Basisrapport Rode Lijst Dagvlinders 2019 volgens Nederlandse en

IUCN-criteria. Rapport VS2019.001, De Vlinderstichting, Wageningen. Retrieved on 15 March 2021,

from https://assets.vlinderstichting.nl/docs/ea504174-725d-4b69-9401-ba364f4ed3de.pdf

van Swaay, C. A. M., Dennis, E. B., Schmucki, R., Sevilleja, C., Balalaikins, M., … & Roy, D. B. (2019a). The EU butterfly indicator for grassland species: 1990-2017: Technical report. Butterfly

Conservation Europe & ABLE/eBMS. Retrieved on 4 March 2021, from www.butterfly-monitoring.net

van Swaay, C. A. M., Wallis de Vries, M. F., van Grunsven, R. H. A. (2019b). Stikstofindicator Vlinders en

Libellen. De Vlinderstichting, Rapport VS2019.030. Retrieved on 5 March 2021, from

https://assets.vlinderstichting.nl/docs/3791766b-a428-42a5-87ed-dcd92a5ef186.pdf Walker, J. A. (2018). Elements of Statistical Modeling for Experimental Biology. Retrieved on 25 May

2021, from https://www.middleprofessor.com/files/applied-biostatistics_bookdown/_book/

Wallis de Vries, M. F., & van Swaay, C. A. (2006). Global warming and excess nitrogen may induce butterfly decline by microclimatic cooling. Global Change Biology, 12(9), 1620-1626.

https://doi.org/10.1111/j.1365-2486.2006.01202.x

Wallis de Vries, M. F., & van Swaay, C. A. (2013). Effects of local variation in nitrogen deposition on butterfly trends in The Netherlands. In Proceedings of the Netherlands Entomological Society

meeting (Vol. 24, pp. 25-33). Retrieved on 2 March 2021, from

https://www.nev.nl/pages/publicaties/proceedings/nummers/24/25-33.pdf

Wallis de Vries, M. F., & van Swaay, C. A. (2017). A nitrogen index to track changes in butterfly species assemblages under nitrogen deposition. Biological Conservation, 212, 448-453.

https://doi.org/10.1016/j.biocon.2016.11.029

Wallis de Vries, M. F., Van Swaay, C. A., & Plate, C. L. (2012). Changes in nectar supply: a possible cause of widespread butterfly decline. Current Zoology, 58(3), 384-391.

https://doi.org/10.1093/czoolo/58.3.384

Weiss, S. B. (1999). Cars, cows, and checkerspot butterflies: nitrogen deposition and management of nutrient‐poor grasslands for a threatened species. Conservation Biology, 13(6), 1476-1486. https://doi.org/10.1046/j.1523-1739.1999.98468.x

Wichink Kruit, R. J., & van Pul, W. A. J. (2018). Ontwikkeling in Stikstofdepositie. RIVM Briefrapport,

0117, 2018. Retrieved on 19 May 2021, from

(20)

Appendices

Appendix A. Data repository

The data is saved on the University of Amsterdam server by the supervisor Henrik Barmentlo. Once further analysis with the data is performed and an official publication is published, a link to the data will follow.

Appendix B. Methods

Table B1. List of the butterfly species studied. With their name in English, Dutch, their Latin family name and their subfamily. The family and subfamily of the species were added based on the Dutch species name from the DBC website (De Vlinderstichting, 2021b).

Species English name Species Dutch Name Family Subfamily

Speckled wood bont zandoogje Nymphalidae Pararge aegeria

Wood white boswitje Pieridae Leptidae sinapis

Brimstone citroenvlinder Pieridae Gonepteryx rhamni

Purple hairstreak eikenpage Lycaenidae Favonius quercus

Pale clouded yellow gele luzernevlinder Pieridae Colias hyale

Alcon blue gentiaanblauwtje Lycaenidae Phengaris alcon

Large skipper groot dikkopje Hesperiidae Ochlodes sylvanus

Large copper grote vuurvlinder Lycaenidae Lycaena dispar

Purple emperor grote weerschijnvlinder Nymphalidae Apatura iris

Silver-washed fritillary keizersmantel Nymphalidae Argynnis paphia

Tree grayling kleine heivlinder Nymphalidae Hipparchia statilinus

Old World swallowtail koninginnenpage Popilionidae Papilio machaon

Clouded yellow oranje luzernevlinder Pieridae Colias crocea

Camberwell beauty rouwmantel Nymphalidae Nymphalis antiopa

Brown hairstreak sleedoornpage Lycaenidae Thecla betulae

Cranberry blue veenbesblauwtje Lycaenidae Agriades optilete

Glanville fritillary veldparelmoervlinder Nymphalidae Melitaea cinxia

Red admiral atalanta Nymphalidae Vanessa atalanta

Holly blue boomblauwtje Lycaenidae Celastrina argiolus

Comma gehakkelde aurelia Nymphalidae Polygonia c-album

Large white groot koolwitje Pieridae Pieris brassicae

Green-veined white klein geaderd witje Pieridae Pieris napi

Small white klein koolwitje Pieridae Pieris rapae

Small tortoiseshell kleine vos Nymphalidae Aglais urticae

Map landkaartje Nymphalidae Araschnia levana

Wall brown argusvlinder Nymphalidae Lasiommata megera

European peacock dagpauwoog Nymphalidae Aglais io

Painted lady distelvlinder Nymphalidae Vanessa cardui

Orange tip oranjetipje Pieridae Anthocharis cardamines

Brown argus bruin blauwtje Lycaenidae Aricia agestis

Common blue butterfly icarusblauwtje Lycaenidae Polyommatus icarus

Queen of Spain fritillary kleine parelmoervlinder Nymphalidae Issoria lathonia

Essex skipper zwartsprietdikkopje Hesperiidae Thymelicus lineola

Meadow brown bruin zandoogje Nymphalidae Maniola jurtina

(21)

Niobe fritillary duinparelmoervlinder Nymphalidae Fabriciana niobe

Small heath hooibeestje Nymphalidae Coenonympha pamphilus

Ringlet koevinkje Nymphalidae Aphantopus hyperantus

Grizzled skipper aardbeivlinder Hesperiiidae Pyrgus malvae

Small skipper geelsprietdikkopje Hesperiidae Thymelicus sylvestris

Dark green fritillary grote parelmoervlinder Nymphalidae Argynnis aglaja

Small copper kleine vuurvlinder Lycaenidae Lycaena phlaeas

Gatekeeper oranje zandoogje Nymphalidae Pyronia tithonus

Sooty copper bruine vuurvlinder Lycaenidae Lycaena tityrus

Grayling heivlinder Nymphalidae Hipparchia semele

Chequered skipper bont dikkopje Hesperiidae Carterocephalus palaemon

Ilex hairstreak bruine eikenpage Lycaenidae Satyrium ilicis

Green hairstreak groentje Lycaenidae Callophrys rubi

Silver-studded blue heideblauwtje Lycaenidae Plebeius argus

Silver-spotted skipper kommavlinder Hesperiidae Hesperia comma

Cranberry fritillary veenbesparelmoervlinder Nymphalidae Boloria aquilonaris

Large heath veenhooibeestje Nymphalidae Coenonympha tullia

Figure B1. Map of the MAN measurement points which were used (green points) and omitted (red) form the dataset, due to double measurement points in the grid cells of the occupancy data. Natura 2000-regions are also shown (white outline).

(22)

Appendix C. Ammonia

Table C1. Results of Tukey’s HSD test for ammonia concentrations over the years. The years which were compared are shown with their corresponding difference (diff), the lower (lwr) and upper (upr) 95% confidence interval of the difference and the corresponding p-value (p adj).

Compared

years diff lwr upr p adj

Compared

years diff lwr upr p adj

2006-2005 -0.759 -0.91414 -0.60386 0 2010-2009 -0.49206 -0.61117 -0.37294 0 2007-2005 -0.83078 -0.9836 -0.67796 0 2011-2009 0.10373 -0.01341 0.220867 0.149115 2008-2005 -0.84709 -1.00018 -0.694 0 2012-2009 -0.03188 -0.14814 0.084385 0.999734 2009-2005 -0.51037 -0.66319 -0.35755 0 2013-2009 -0.3348 -0.45079 -0.21881 0 2010-2005 -1.00242 -1.14264 -0.86221 0 2014-2009 0.249328 0.136379 0.362277 0 2011-2005 -0.40664 -0.54518 -0.2681 0 2015-2009 -0.39522 -0.50752 -0.28291 0 2012-2005 -0.54225 -0.68005 -0.40445 0 2016-2009 0.761662 0.649169 0.874156 0 2013-2005 -0.84517 -0.98273 -0.7076 0 2017-2009 0.596557 0.484375 0.708739 0 2014-2005 -0.26104 -0.39605 -0.12603 0 2018-2009 -1.24973 -1.36282 -1.13665 0 2015-2005 -0.90559 -1.04006 -0.77111 0 2011-2010 0.595787 0.495646 0.695927 0 2016-2005 0.251294 0.116661 0.385927 0 2012-2010 0.460178 0.36106 0.559295 0 2017-2005 0.086189 -0.04818 0.220562 0.665361 2013-2010 0.157257 0.058463 0.25605 8.4E-06 2018-2005 -1.7601 -1.89523 -1.62498 0 2014-2010 0.741385 0.646177 0.836593 0 2007-2006 -0.07178 -0.20815 0.064589 0.891659 2015-2010 0.096839 0.002396 0.191282 0.038085 2008-2006 -0.08809 -0.22477 0.048588 0.657875 2016-2010 1.253719 1.159052 1.348386 0 2009-2006 0.248631 0.11226 0.385002 1E-07 2017-2010 1.088614 0.994318 1.18291 0 2010-2006 -0.24343 -0.36551 -0.12134 0 2018-2010 -0.75768 -0.85305 -0.66231 0 2011-2006 0.352361 0.232209 0.472513 0 2012-2011 -0.13561 -0.23234 -0.03888 0.000223 2012-2006 0.216752 0.097452 0.336052 1E-07 2013-2011 -0.43853 -0.53493 -0.34213 0 2013-2006 -0.08617 -0.2052 0.032862 0.459678 2014-2011 0.145598 0.052878 0.238318 1.24E-05 2014-2006 0.497959 0.381886 0.614031 0 2015-2011 -0.49895 -0.59088 -0.40701 0 2015-2006 -0.14659 -0.26203 -0.03114 0.001687 2016-2011 0.657933 0.565769 0.750096 0 2016-2006 1.010293 0.894664 1.125922 0 2017-2011 0.492827 0.401044 0.58461 0 2017-2006 0.845188 0.729862 0.960514 0 2018-2011 -1.35346 -1.44635 -1.26058 0 2018-2006 -1.0011 -1.11731 -0.8849 0 2013-2012 -0.30292 -0.39826 -0.20759 0 2008-2007 -0.01631 -0.15034 0.117728 1 2014-2012 0.281207 0.189593 0.372821 0 2009-2007 0.320413 0.186691 0.454136 0 2015-2012 -0.36334 -0.45416 -0.27252 0 2010-2007 -0.17164 -0.29076 -0.05253 0.000118 2016-2012 0.793541 0.70249 0.884592 0 2011-2007 0.424143 0.307006 0.541281 0 2017-2012 0.628436 0.537771 0.719101 0 2012-2007 0.288534 0.172271 0.404798 0 2018-2012 -1.21786 -1.30964 -1.12608 0 2013-2007 -0.01439 -0.13037 0.101601 1 2014-2013 0.584128 0.492865 0.675391 0 2014-2007 0.569741 0.456792 0.682691 0 2015-2013 -0.06042 -0.15088 0.030046 0.600094 2015-2007 -0.0748 -0.18711 0.0375 0.604548 2016-2013 1.096463 1.005765 1.187161 0 2016-2007 1.082076 0.969582 1.194569 0 2017-2013 0.931357 0.841046 1.021668 0 2017-2007 0.91697 0.804789 1.029152 0 2018-2013 -0.91493 -1.00636 -0.82351 0 2018-2007 -0.92932 -1.04241 -0.81624 0 2015-2014 -0.64455 -0.73108 -0.55801 0 2009-2008 0.33672 0.202685 0.470755 0 2016-2014 0.512335 0.425556 0.599113 0 2010-2008 -0.15534 -0.2748 -0.03587 0.001079 2017-2014 0.347229 0.260855 0.433603 0 2011-2008 0.44045 0.322956 0.557943 0 2018-2014 -1.49906 -1.58661 -1.41152 0 2012-2008 0.304841 0.188218 0.421463 0 2016-2015 1.15688 1.070942 1.242819 0 2013-2008 0.00192 -0.11443 0.118267 1 2017-2015 0.991775 0.906245 1.077305 0

(23)

2015-2008 -0.0585 -0.17117 0.054178 0.90135 2017-2016 -0.16511 -0.25088 -0.07933 0

2016-2008 1.098382 0.985518 1.211247 0 2018-2016 -2.0114 -2.09835 -1.92444 0

2017-2008 0.933277 0.820723 1.04583 0 2018-2017 -1.84629 -1.93284 -1.75974 0

2018-2008 -0.9130 -1.0264 -0.7995 0

(24)

(25)

Figure 11. Histogram of yearly NH3 mean (µg m-3) values per location from 2005-2018. The histogram of the dataset

indicates a right-skewed distribution, which is also substantiated by the order: mode > median > mean (4.22 > 4.31 > 4.51 µg m-3_).

Appendix D. Butterfly occupancy

Figure D1. Histogram of the mean butterfly occupancy per year for the study period (2005-2018). Because most of the

(26)

Figure D2. Violin chart of the mean butterfly occupancy per location per year from 2005 to 2018. For all years, the median

value (horizontal line) is very close to 0, while the mean value (rhombus) lies between 0.25 and 0.3. For all of the years a quasibinomial distribution is present.

Table D1. Results of Tukey’s HSD test for the butterfly occupancy data between years. The table shows which years are

compared to each other (Comp years), what the difference is between these years (Diff), what the lower and upper confidence interval is (Lwr CI, Upr CI) and the corresponding p-value (p (adj)).

Comp years Diff Lwr CI Upr CI p (adj) Comp years Diff Lwr CI Upr CI p (adj)

2006-2005 0.0156 1.49E-06 0.031 0.0499 2015-2006 -0.009 -2.46E-02 0.007 0.8144

2007-2005 -0.007 -2.31E-02 0.008 0.9453 2015-2007 0.0141 -1.48E-03 0.03 0.1255

2007-2006 -0.023 -3.87E-02 -0.01 6E-05 2015-2008 0.0137 -1.96E-03 0.029 0.1635

2008-2005 -0.007 -2.26E-02 0.009 0.9677 2015-2009 -0.014 -2.97E-02 0.002 0.1291

2008-2006 -0.023 -3.82E-02 -0.01 0.0001 2015-2010 -0.004 -1.91E-02 0.012 1

2008-2007 0.0005 -1.51E-02 0.016 1 2015-2011 0.0107 -4.95E-03 0.026 0.5612

2009-2005 0.0208 5.14E-03 0.036 0.0007 2015-2012 0.0052 -1.04E-02 0.021 0.998

2009-2006 0.0051 -1.05E-02 0.021 0.9982 2015-2013 -0.018 -3.36E-02 -0 0.0086

2009-2007 0.0282 1.26E-02 0.044 1E-07 2015-2014 -0.006 -2.14E-02 0.01 0.9937

2009-2008 0.0278 1.21E-02 0.043 2E-07 2016-2005 0.0088 -6.77E-03 0.024 0.8274

2010-2005 0.0102 -5.44E-03 0.026 0.6408 2016-2006 -0.007 -2.24E-02 0.009 0.9752 2010-2006 -0.005 -2.11E-02 0.01 0.9967 2016-2007 0.0163 7.02E-04 0.032 0.0307 2010-2007 0.0177 2.03E-03 0.033 0.0111 2016-2008 0.0158 2.21E-04 0.031 0.043 2010-2008 0.0172 1.55E-03 0.033 0.0163 2016-2009 -0.012 -2.75E-02 0.004 0.3682 2010-2009 -0.011 -2.62E-02 0.005 0.576 2016-2010 -0.001 -1.70E-02 0.014 1 2011-2005 -0.004 -1.96E-02 0.012 0.9999 2016-2011 0.0129 -2.77E-03 0.028 0.2452 2011-2006 -0.02 -3.53E-02 -0 0.002 2016-2012 0.0074 -8.25E-03 0.023 0.9512 2011-2007 0.0035 -1.22E-02 0.019 1 2016-2013 -0.016 -3.14E-02 -0 0.0446 2011-2008 0.003 -1.26E-02 0.019 1 2016-2014 -0.004 -1.93E-02 0.012 1

2011-2009 -0.025 -4.04E-02 -0.01 9E-06 2016-2015 0.0022 -1.34E-02 0.018 1

2011-2010 -0.014 -2.98E-02 0.001 0.1225 2017-2005 0.0074 -8.22E-03 0.023 0.9495

The Effect of Ammonia Deposition on Butterfly Occupancy in Dutch Nature Reserves