Urban indicators for UK butterﬂies

Most people live in urban environments and there is a need to produce abundance indices to assist policy and management of urban greenspaces and gardens. While regional indices are produced, with the exception of birds, studies of the differences between urban and rural areas are rare. We explore these differences for UK butterﬂies, with the intention to describe changes that are relevant to people living in urban areas, in order to better connect people with nature in support of conservation, provide a measure relevant to human well-being, and assess the biodiversity status of the urban environment. Transects walked under the UK Butterﬂy Monitoring Scheme are classiﬁed as urban or rural, using a classiﬁcation for urban morphological zones. We use models from the Generalised Abundance Index family to produce urban and rural indices of relative abundance for UK butterﬂy species. Composite indices are constructed for various subsets of species. For univoltine and bivoltine species, where we are able to ﬁt phenomenological models, we estimate measures of phenology and identify urban/rural differences. Trends in relative abundance over the period 1995–2014 are more negative in urban areas compared to rural areas for 25 out of 28 species. For the composite indices, all trends are negative, and they are signiﬁcantly more negative for urban areas than for rural areas. Analysis of phenological parameters shows butterﬂies tend to emerge earlier in urban than in rural areas. In addition, some ﬂy longer in urban than in rural areas, whereas in other cases the opposite is the case, and hypotheses are proposed to account for these features. Investigating new urban/rural indicators has revealed national declines that are stronger for urban areas. For continued monitoring, there is a need for an urban butterﬂy indicator, and for this to be evaluated and reported annually. We explain how this may be interpreted, and the relevance for other monitoring schemes. The results of this paper, including the phenological ﬁndings, shed new light on the potentially deleterious effects of urbanisation and climate change, which require suitable monitoring and reporting to support policy and management, for example of urban greenspaces and gardens. © 2017 The Authors. by Ltd. This is an open access article under the CC BY license


Introduction
In recent decades the world has become increasingly urbanised, with over half of the global population now living in urban areas, a proportion which is predicted to increase to 66% by 2050 (United Nations, 2014). In 2009, 81.4% of the population in England was classified as living in an urban area, while in Wales, Scotland and Northern Ireland, 66.1, 81.6 and 63.6% of the population occupy urban areas, respectively (Pateman, 2011). Urbanisation has been identified as a cause of loss, degradation and fragmentation of habitats, and species richness is often reduced in highly urbanised areas, species rather than sites (Gregory et al., 2005), or for sites within specific city boundaries (Herrando et al., 2012). Van Dyck et al. (2009) found declines in abundance in urban areas for butterflies in the Netherlands and in Belgium. In the UK habitat-specific butterfly indicators have been developed for woodland and farmland (Brereton et al., 2011a).
Global biodiversity loss is well known and biodiversity indicators provide a means for monitoring changes and measuring performance in reaching biodiversity targets at national and global scales (Defra, 2016;Butchart et al., 2010). Urban biodiversity provides important cultural ecosystem services (Bolund and Hunhammar, 1999) and may contribute towards human well-being (Millennium Ecosystem Assessment, 2005). In this paper we derive and compare indicators for urban and rural areas of the UK, using butterfly monitoring data for illustration. Improved monitoring of species' abundance in urban areas may be used for assessing and refining conservation efforts, as well as engaging and educating the general public, the majority of whom live in urban areas (McKinney, 2002). Fox et al. (2014) attributed declines in occurrence of widespread moths in southern Britain to increases in urban and arable land cover in this area. Declines in the abundance of birds have similarly been found for southern and eastern England (Harrison et al., 2014;Massimino et al., 2015), as well as the extinction of some plant species (Preston, 2000). Unlike for UK birds, designations of sites as urban by the recorders are not available for butterflies, hence a suitable classification is required. Some studies have used land cover data at varying scales to define urban locations e.g. Botham et al. (2009), Evans et al. (2011. Butterflies respond sensitively and rapidly to changes in climate and habitat and may act as representatives for less well-monitored insect groups, with the exception of saproxylic species which depend upon dead or decaying wood (Thomas, 2005). Unlike most other insects groups, butterflies are well-documented, their taxonomy is understood, many species are easy to identify and there is a wealth of information on their ecology and life-histories (Thomas, 2005). Butterflies are also culturally important as demonstrated by their popularity amongst the general public and frequent appearances in art and literature . These attributes make butterflies potentially valuable biodiversity indicators, which has been recognized politically, with EU, UK and devolved governments adopting the status of butterfly populations as a measure of biodiversity and environmental health (Brereton et al., 2011a;van Swaay et al., 2016;Defra, 2016).
The pressures of urbanisation on butterflies have been identified in various studies undertaken in many countries, with these generally finding negative impacts on species richness and abundance (Blair, 1999;Jones and Leather, 2012;Ramírez-Restrepo and MacGregor-Fors, 2016).
Urban phenology has been widely studied, particularly for plants (Neil and Wu, 2006), and can be used to estimate urban heat island effects (Kaiser et al., 2016) and as an indicator for the effects of climate change (Jochner and Menzel, 2015). Surprisingly, Altermatt (2012) found delayed butterfly phenology in urban habitats, and suggested this could be due to dispersal among habitat types, and Diamond et al. (2014) also found delayed phenology due to the joint effects of urbanisation and climate change. By using generalised abundance index (GAI; Dennis et al., 2016a) models for analysis, we study parameters relevant to the phenology of individual broods, providing a new comparison of urban and rural phenology for UK butterflies.
In the following sections we describe the classification of urban areas used and the application of the GAI approach to produce urban and rural indices of abundance for a selection of UK butterfly species. Composite (multi-species) indices are then formed for groups of species, for comparison of urban and rural areas.
Estimates of parameters relating to phenology are also compared for urban and rural areas.

UK Butterfly Monitoring Scheme (UKBMS)
The UK Butterfly Monitoring Scheme (UKBMS) began in 1976 with 34 sites, and has grown to a network of over 1400 transects monitored each year . Sites are self-selected but it is recommended that routes are chosen to include a fair representation of habitats present at the site, including areas that might become more suitable for butterflies in future. Should habitats change or degrade over time on sites, then that may contribute to declines in indices, however there is also an annual turnover of sites of approximately 13%, for example in 2015, 223 new sites were established .
Transects are typically 2-4 km long, where an observer counts butterflies under specified times of the day and weather conditions (Pollard and Yates, 1993). Counts are made weekly from the beginning of April until the end of September, the main season for butterfly activity in the UK. Roughly 30% of potential counts in the 26 week season are missed in practice, for example due to unsuitable weather conditions or recorder unavailability (Dennis et al., 2013).
The Wider Countryside Butterfly Survey (WCBS) has formed part of the UKBMS since 2009, and aims to contribute to improved population trend estimates for common and widespread butterfly species that are representative of the whole countryside (Brereton et al., 2011b). The scheme adopts the sampling framework used by the Breeding Bird Survey (BBS), which uses 1 km grid squares, and counts are made along two, evenly-spaced, parallel transects of 1 km length, which are surveyed on at least two separate days during July and August. In 2015 over 800 WCBS squares were sampled . In the following analysis we include data from both the WCBS and UKBMS.

Definition of urban
The Urban Morphological Zones data for the year 2000 (UMZ2000) provide a consistent definition of urban areas in all UK countries and across Europe and can be defined as "A set of urban areas lying less than 200m apart" (European Environment Agency, 2016a). Data were sourced from the European Environment Agency (European Environment Agency, 2016a), where further details on the definition and methodology can also be found, and are derived from Corine land cover data (European Environment Agency, 2016b). For the purposes of this study we use the urban classification for the UK, which is displayed in Fig. S1, and was extracted in QGIS (QGIS Development Team, 2016).

Data selection
Each UKBMS (and WCBS) transect was classified as either urban or rural according to the UMZ2000 classification, based on the grid reference of the middle of the transect. We compare the land cover composition for urban and rural coverage across the UK and among transects using a land cover map from 2007 (LCM2007; Morton et al., 2014). For each species monitored by the UKBMS, the total and minimum number of sites classified as urban and rural across years were identified. Species monitored at a minimum of 5 rural and urban sites each year and at least 20 across all years were retained for analysis. Clouded Yellow was additionally excluded since there were no positive counts on urban sites for 1995. This resulted in trends calculable for 15 and 28 species, based on sites monitored during 1980-2014 and 1995-2014, respectively, for which common Table 1 The 28 species used in the analysis for data for the period 1995-2014, together with the species codes used in the paper. The 15 species which were analysed for 1980-2014 are indicated by X, and flight periods (FP) are represented by U (univoltine), B (bivoltine) and M (multivoltine). and Latin names are given in Table 1. Mid-term (20 years) as well as long-term (35 years) trends were estimated to allow for more species to be included in the indicator in the former case.

Modelling approach
We use members of the generalised abundance index family (GAI; Dennis et al., 2016a) to estimate urban and rural indices of abundance for each species. To facilitate comparisons involving many species, only the Poisson GAI is used, although alternative distributions and goodness-of-fit have been explored in Dennis et al. (2016a). Briefly, each count y i,j for site i and visit j in a given year is treated as the realization of a Poisson random variable with expectation i,j = N i a i,j , where N i represents the relative abundance for the ith site, and a i,j denotes a function describing seasonal variation in counts.
An efficient concentrated likelihood approach was employed for obtaining maximum-likelihood estimates and is described in Dennis et al. (2016a). The number of parameters to estimate numerically is reduced appreciably by optimising a likelihood only with respect to the parameters associated with a i,j . The site parameters, N i , in a given year are then estimated by where T is the maximum number of visits to a site within a single year. The GAI approach is suitable for any number of visits to each site and can hence incorporate both UKBMS and WCBS data.
We use a phenomenological GAI to describe the annual flight period for univoltine and bivoltine species. Here a i,j in a given year is described by a mixture of B Normal probability density functions, corresponding to the number of broods, B, a species has per year, such that where w b , b and b correspond to the weight, mean and standard deviation, respectively, for the bth brood, and B b=1 w b = 1. In this paper we shall only consider the cases of B=1 and B=2. For simplicity, in the bivoltine case we assume a homoscedastic case, where = 1 = 2 , and denote w 1 = w, and w 2 = 1 − w, where w therefore represents the size of the first brood relative to the second brood. The standard deviation, , of the Normal distribution that is assumed for the flight period governs the length of the flight period: the larger is then the longer the flight period. Thus it is convenient to use the estimate of as a value which is effectively proportional to the estimated flight period length (for a given brood).
For species with more complex flight periods, a spline is required to describe seasonal variation in the counts. The spline GAI presented in Dennis et al. (2016a) does not include automatic selection of the level of smoothing. Hence to estimate seasonal variation, a i,j , we fit a simple generalised additive model (GAM) to the data using the mgcv package in R (Wood, 2006;R Core Team, 2016), which uses generalised cross validation to select an appropriate level of smoothing. The site effects, N i , are then estimated as for the phenomenological case, using the concentrated likelihood approach (Dennis et al., 2016a). Species' designated voltinism is given in Table 1. All analysis was performed in R (R Core Team, 2016), for which code is available in the Supplementary Material.
We fit a Poisson generalised linear model with year and site factors to the estimated site parameters N i , and use scaled predicted year effects as the index of abundance (Dennis et al., 2013). This accounts for variation in the sites sampled each year, which is particularly relevant in cases with small numbers of sites. Nonparametric bootstrapping with 500 replicates was used to produce confidence intervals for the indices of abundance, accounting for sources of uncertainty from all stages of modelling.
Multi-species (composite) indices of abundance were calculated from the geometric mean index across groups of species (Buckland et al., 2011). For each year, the composite index is estimated by the exponential of the average of the log abundance index for each species. Confidence intervals were derived by estimating a composite index for each bootstrap replicate. Composite indices were derived based on all species and for all resident species (i.e. excluding the two regular migrants Painted Lady and Red Admiral), as well as separately for univoltine, bivoltine and multivoltine species, and the three habitat specialist species (indicated in Table 1) with sufficient data to produce urban indices.
Trends in individual species and composite abundance were estimated by fitting simple linear regressions to the indices; associated errors were estimated by fitting simple linear regressions to each bootstrap replicate and the construction of suitable quantiles.
For the univoltine and bivoltine species, the phenomenological GAI provides estimates of additional parameters relating to phenology, namely b , the mean flight date for each brood, , representing flight period length, and for the bivoltine case, w, describing the size of the first brood relative to the second brood, from Eq. (1). We compare average annual estimates of the phenology parameters over 1995-2014 for urban and rural areas. Associated confidence intervals were derived from the bootstrapping approach. In models for bivoltine species we let 2 = 1 + d , where  Table S1 summarises the number of urban and rural sites monitored by the UKBMS and WCBS within each region considered, based on the UMZ2000 classification. Most urban transects are found in England (Fig. 1), which reflects the fact that the majority of urban areas in the UK are in England (Fig. S1).

Defining urban transects
Unsurprisingly, a greater coverage of urban habitats is found for urban versus rural areas ( Table 2). The urban transects show lower average percentages of arable and horticultural land relative to rural transects, but still indicate relatively high percentages of improved grassland.
Each species was monitored on far more rural than urban transects (Table S2). Some species were monitored on only a minimal number of transects classified as urban, for example Grizzled Skipper and Dingy Skipper.

Comparing trends in abundance
Urban and rural indices of relative abundance for each species are presented in Figs. S2 and S3, where confidence intervals were omitted for clarity, and were significantly correlated for 25 out of 28 species for 1995-2014 and for all species 1980-2014 (Table S3). For 1995 onwards, the indices appear particularly similar for Brimstone, Painted Lady and Red Admiral, as well as Small and Large White, although the correlation coefficients were greater than 0.8 for 13 out of 28 species.
There was a significant overall correlation between population trends for 1995-2014 from rural and urban transects ( = 0.73, p < 0.001, Fig. 2), but not for 1980-2014 ( = 0.24, p > 0.05). For the period 1995-2014, population trends were mostly negative: for rural transects three species showed positive population changes, but only Orange-tip showed a significant increase, and for urban transects population trends for all 28 species were negative (Table 3). Significant declines (at the 5% level) were found for 17 and 18 species for urban and rural transects, respectively. The urban trend was more negative than the rural trend for 25 of the 28 species (Fig. 2), with a significant difference for 11 species. The urban trend was higher for Grizzled Skipper, Holly Blue and Small White, but not significantly (Table 3). For some species the confidence intervals for the urban trends are wide, which can be explained by fewer urban transects being sampled (Fig. S4). For 1980-2014, there was a smaller sample of 15 species, but most population trends were again negative (Table 4). Rural trends were negative with the exception of three species, of which Painted Lady and Red Admiral showed significant positive rural population trends yet significant negative urban trends. For 1980-2014 the urban trend was lower than the rural trend for 12 species, with significant differences for seven species.
Trends estimated from composite indices based on all species (Fig. 3) were 23.9 and 29.9% more negative in urban areas for 1995-2014 and 1980-2014, respectively (Table 5). These differences were significant at the 5% level. Composite trends based on subsets of species (Fig. S5) were also negative, and were generally significantly worse for urban compared to rural transects ( Table 5). The inclusion of the two regular migrants had a minimal effect on the composite indices. A comparison for habitat specialist species indicated urban declines were also greater for these species, although this was only based on three species.

Phenology
Annual estimates of the phenology parameters for urban and rural transects are shown in Figs. S6-S9, but here we focus on average values. Mean flight dates were earlier for urban versus rural transects for 16 of 20 univoltine and bivoltine species (Fig. 4a, based on the first brood for bivoltine species). Across all 20 species the mean flight date was approximately 2 days earlier, which increases to approximately 2.8 days when only the 16 species with earlier mean flight dates in urban areas were considered. Differences were generally greater for bivoltine versus univoltine species, with bivoltine species on average 3.2 days earlier in urban areas, however differences tended to be significant for univoltine rather than bivoltine species (Table 6). The largest difference, and the only significant difference (at the 5% level) among the bivoltine species, was for Table 3 Species trends in abundance for 1995-2014, with 95% confidence intervals, where values significantly different from zero at the 5% level are highlighted in bold.

Species
Urban   Brimstone, for which the mean flight date was on average 5.3 days earlier on urban transects (Fig. S6). The length of the flight period is represented by , which was described further in Section 2.4. For 15 out of 20 species the flight period is shorter in urban areas (Fig. 4b), and this subset included all 12 univoltine species, of which was significantly smaller for 9 species, and was on average 1.7 days shorter in urban areas. Green Hairstreak showed the greatest difference, with typically 5.9 days less (significantly) in urban areas. For five bivoltine species was larger than for most other species and was also greater for urban compared to rural transects, although not significantly at the 5% level.
For bivoltine species, as the difference in mean flight dates of the two broods increased, the difference between urban and rural transects showed a negative trend (Fig. 4c). Species such as Common Blue and Large White, with shorter time periods between broods, generally had a greater difference in mean flight dates of the two broods in urban areas, whereas Holly Blue and Brimstone, with a greater gap between broods, typically has a smaller difference between broods in urban versus rural areas. Comparisons between the weighting of the two broods, which reflects the size of the first brood relative to the second brood, between urban and rural transects tended to show a negative decline with an increasing mean value of the weighting (Fig. 4d), with the exception of Brimstone, which overwinters as an adult, and shows an increased relative size of the first brood in urban areas (Fig. S9).

Discussion
Indicators are used to measure and communicate progress in biodiversity conservation. We have illustrated the development of indicators for urban and rural areas for UK butterflies, and shown greater declines for urban populations at two time-scales, suggesting a need for regular monitoring of urban populations to provide a meaningful biodiversity indicator of the state of urban environments and the quality of life for people within them. The causes of differences between urban and rural areas are linked Table 5 Composite trends for varying subsets of n species, from (a) 1995-2014 and (b) 1980-2014, with 95% confidence intervals, where significant changes (at the 5% level) are highlighted in bold.  [1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013][2014] in weeks for urban and rural sites for each species, based on mean flight date, 1 (a), flight period length, measured by (b), and for bivoltine species, the difference in broods, d (c), and weight given to the first brood, w (d). The solid horizontal line indicates no difference between urban and rural. Univoltine and bivoltine species are indicated by blue triangles and black circles, respectively. Species codes represent species common names (Table 1). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 6
Difference in estimated phenology parameters (urban-rural), averaged over 1995-2014. Significant changes (at the 5% level) are highlighted in bold.

Species
1 d w to a number of environmental pressures. These include increased pollution; influences on habitat availability, for example habitat loss and fragmentation due to development; reductions in garden size; and loss of brownfield land (McKinney, 2002;Matteson and Langellotto, 2010;Muratet and Fontaine, 2015). Differences in the declines of butterflies among urban and rural areas may be better understood by species' traits (Lizée et al., 2011;Diamond et al., 2014;De Palma et al., 2015), for example responses may be influenced my mobility and habitat preference (Olivier et al., 2016). Indices of abundance for Painted Lady and Red Admiral showed high correlation in urban and rural areas, which might be expected for highly mobile migrant species. There were also similarities in the indices for Brimstone, Small White and Large White which are mobile generalist species. Greater differences between urban and rural areas were found for many species, some of which might be more limited in habitat preferences and less typically found in urban areas such as gardens, for example Brown Argus, Wall and Small Skipper.
Confidence intervals for the urban trends were relatively wide for some species due to small sample sizes, hence the true significance of the greater declines in urban areas is likely to be stronger than predicted here. The greater sampling of rural compared with urban areas is further demonstrated in the size of the confidence intervals for the composite indices. Additional monitoring effort could be targeted at improving urban coverage in the UKBMS and WCBS, which could also potentially allow for more species to contribute to the urban indicators. However for some species the small sample sizes may be due to limits in a species range or habitat with respect to urban areas.
Monitoring scheme data are typically collected by amateur experts, requiring a high level of commitment and identification skills. These data are valuable for assessing changes in species' abundance, but more typical citizen-science data collected by the general public may also be valuable for urban ecology (Bates et al., 2015;Wei et al., 2016). For UK butterflies, opportunistic data sources, such as the Butterflies for the New Millennium, could be explored to produce urban and rural indicators of occurrence using occupancy modelling (Kéry et al., 2010;Fox et al., 2015;Dennis et al., 2015). This could produce urban indicators for more species, although species misidentification may be an issue for opportunistic data (Miller et al., 2011). Multi-species occupancy models could be explored, as adopted by Mata et al. (2014). Data from the Big Butterfly Count, a mass-participation citizen science project, could be valuable for studying urban populations of common butterfly species. Due to the nature of the scheme, the data provide good coverage of urban areas, particularly gardens (Dennis et al., 2016c), which have an important role for urban wildlife including butterflies (Goddard et al., 2010;Bergerot et al., 2011;Fontaine et al., 2016).
Using the phenomenological GAI allowed for new analyses of parameters relating to phenology, demonstrating the potential of an urban indicator to provide a measure of climate change. Mean flight dates were earlier in urban areas for most species, in contrast with the findings of Altermatt (2012) and Diamond et al. (2014). Advanced urban phenology could be explained by urban heat island effects (Kaiser et al., 2016), or could be biased by the southern distribution of many of the urban transects. The phenology of butterflies is known to be delayed in cooler parts of their range (Roy and Asher, 2003), although varying with a lesser extent than betweenyear climatic variation , and this aspect could be investigated by additionally accounting for this spatial variation, for example by incorporating covariates (Dennis et al., 2016a). The difference between the mean flight date was smallest for Gatekeeper (when the difference in flight period length was also minimal), for which emergence is thought to be synchronised across the UK (Roy and Asher, 2003). In contrast the greatest difference in mean emergence in urban and rural areas was for Brimstone, which is supported by findings that species which overwinter as adults show greater phenological advancements (Diamond et al., 2011), although the UKBMS season may also miss the beginning of this species' flight period.
The univoltine species generally emerged earlier in urban areas and the flight periods were shorter; the latter feature could be due to shorter periods of emergence and/or shorter life expectancies. For bivoltine species, such as the three Pieris species which showed increases in flight period length in urban areas, warmer temperatures could lead to additional broods or an increased likelihood of early and late records (Altermatt, 2010), however producing an extra generation may be detrimental to some species, for example the Wall butterfly (Van Dyck et al., 2015). Phenomenological models for multivoltine species with more than two broods could be tested, to allow for further study of phenology, or more typical phenology summaries could be used for these species.
We compared mean estimates of phenological parameters (across years) for urban and rural areas, but a wider study could also compare urban and rural trends in phenology. Phenological variation for insects is known to be complex (Forrest, 2016), and further analysis of the relevant parameters from a GAI could lead to new insights.
In this paper phenomenological and spline GAIs were used, however an alternative is to use a stopover model (Matechou et al., 2014;Dennis et al., 2016a). This mechanistic description of the counts allows for the separation of emergence and survival, which could produce more fine-tuned urban-rural comparisons of phenology and novel study of species' life-expectancies in urban areas. Furthermore, the use of dynamic models (Dennis et al., 2016b) would produce estimates of productivity for each brood for urban and rural areas. Matechou et al. (2014) also demonstrate how, if available, a temperature covariate may be used to account for varying detection probability, however this requires a sufficient number of sites and counts, and variation in detection is generally assumed to be random and minimised by the standardised sampling conditions (van Swaay et al., 2008). Urban and rural trends may also be better studied by building the trend into the GLM stage of the GAI rather than posthoc estimation. Alternatively smoothed indices may be estimated (Fewster et al., 2000;Knape, 2016), and change points in the indicator may be identified (Fewster et al., 2000).
We used a single urban classification for the year 2000 to represent the entire study period, but the UMZ classification is also available for 1990 and 2006, so the effects of changing urbanisation could be studied. The classification of transects was based on the grid reference of the middle of each transect, but with suitable mapping the classification over the length of each transect could be verified. For the purpose of developing an urban indicator, we used a binary urban-rural classification, but analysis along an urban-rural gradient (see for example Blair, 1999) could be more revealing of the impact of urbanisation on declines.
Urban indicators could be developed for particular regions or cities of interest of relevance to conservation and policy, for example an indicator for London, where many of the urban transects were located, could be developed. With greater or alternative datasets, UK country-level urban indicators could be derived. Localised measures of urbanisation could be beneficial for widening public awareness and engagement and the need for conservation in urban areas.
New urban indicators could also be produced for various other taxa, or multi-taxa indicators (Hayhow et al., 2016;Butchart et al., 2010). A European indicator exists for grassland butterfly species (van Swaay et al., 2016), as well as for other taxa (Gregory et al., 2005) and given the geographical scope of the urban classification used here, European urban indicators could be developed to aid the assessment of the effects of urbanisation on a larger scale. Implementation and regular reporting of urban indicators are vital for monitoring population changes under the increasing pressures of urbanisation and anthropogenic climate change in order to support future policy and management.