The influence of chalk grasslands on butterfly phenology and ecology

Abstract The influence of large‐scale variables such as climate change on phenology has received a great deal of research attention. However, local environmental factors also play a key role in determining the timing of species life cycles. Using the meadow brown butterfly Maniola jurtina as an example, we investigate how a specific habitat type, lowland calcareous grassland, can affect the timing of flight dates. Although protracted flight periods have previously been reported in populations on chalk grassland sites in the south of England, no attempt has yet been made to quantify this at a national level, or to assess links with population genetics and drought tolerance. Using data from 539 sites across the UK, these differences in phenology are quantified, and M. jurtina phenology is found to be strongly associated with both site geology and topography, independent of levels of abundance. Further investigation into aspects of M. jurtina ecology at a subset of sites finds no genetic structuring or drought tolerance associated with these same site conditions.


| INTRODUC TI ON
Changes in the phenology of Lepidoptera is a well-studied subject (Dell et al., 2005;Diamond et al., 2011;Hodgson et al., 2011;MacGregor et al., 2019;Roy et al., 2015;Roy & Sparks, 2000) in part because long-term monitoring data are available for a large number of species, allowing temporal changes in phenology to be measured (Roy & Sparks, 2000). Phenology, the annual timing of species life cycles, can be affected by a host of environmental factors such as temperature and daylength (Bale et al., 2002).
For example, increasing temperatures as a result of climate change have been shown to shift Lepidoptera phenology by advancing first flight dates (Roy & Sparks, 2000). Temperature varies spatially and temporally, resulting in changes in phenology based upon latitude. For example, butterfly flight periods have been shown to be shorter and begin later at northern latitudes (Brakefield, 1987).
However, equivalent changes to environmental factors occurring spatially and temporally can have different effects on phenology, with temporal changes in annual temperatures in fixed locations affecting phenology to a greater degree than differences in temperature spatially (Doi & Takahashi, 2008;Roy et al., 2015). This suggests that species may have some level of local adaptation as well as responding to fixed environmental conditions (Hodgson et al., 2011;Roy et al., 2015). cally heterogeneous (Diacon-Bolli et al., 2012;Mortimer et al., 1998) due to differences in vegetation structure and topography, resulting in variation in ground temperatures (Maclean et al., 2019). The resulting microclimates may allow individuals to persist in specific locations when surrounding areas of habitat are climatically unsuitable (Bennie et al., 2008;Suggitt et al., 2011), potentially broadening the flight period. Similarly, extreme warm temperatures in some microclimates may result in local drought conditions which are likely to affect larval development, for example, larvae of speckled wood (Pararge aegeria) reared on drought-stressed plants show longer development times and increased mortality rates (Gibbs et al., 2012(Gibbs et al., , 2018Talloen et al., 2004). Thus, longer development times for some individuals and climatically suitable patches may both contribute to the longer flight periods. Protracted flight periods are also observed in some M. jurtina populations in southern Europe (Haeler et al., 2014); however, this results from adult females entering a period of aestivation (Brakefield, 1984), which has hitherto not been observed in observations of UK populations. Additionally, as this appears to be controlled by geographic provenance and associated larval developmental conditions (Grill et al., 2013), it is unlikely this behavior is present in any UK population.
Although it seems likely that the variation in UK M. jurtina phenology results from differences in conditions that occur within chalk grasslands, the mechanisms that cause these responses are unclear.
A parsimonious explanation of the protracted flight period is that these sites contain more favorable habitat and therefore higher abundances, with the broad flight periods simply a result of the mathematical relationship between mean and variance (Taylor, 1961). If, however, the broad flight period of M. jurtina on chalk grasslands is not purely the result of high abundances, differences in the local site conditions and the ecology of populations at these sites may be affecting phenology.
The broader flight periods on chalk grasslands may be the result of genetic differences between populations, with some anecdotal suggestions of locally adapted races. Although we do not explicitly look at local adaptation here, we do investigate the potential for genetic structuring between populations, based upon the type of site that individuals are found in. Clear genetic clustering of individuals into chalk and nonchalk populations would suggest a high level of genetic differentiation, which may support the idea of locally adapted races as an explanation for the differences in flight periods.
To explore these possibilities, we examine the flight periods of M. jurtina in the UK at 539 sites differing in geology and topography and quantify the variability in phenology. We confirm that flight periods are protracted on chalk grasslands across a wide spatial scale, as previously reported at local sites (Goulson, 1993a). After controlling for abundance in our models, we then investigate levels of genetic diversity and differentiation, and drought tolerance at a subset of sites to determine whether differences in phenology are associated with genetic structuring of populations and whether there is evidence of increased drought tolerance from chalk sites that may influence the flight period length. Overall, we test the following:

| Long-term butterfly monitoring sites and landscape context
Abundance data from 539 long-term monitoring sites (1976 onward; if they had both relevant Natural England priority habitat map and digital elevation data (see below). UKBMS data are collected by volunteers using the "Pollard walk" method (Pollard & Yates, 1993). The UKBMS uses a two-step method (Dennis et al., 2013), using these data to fit generalized additive models which produce fitted weekly counts and an overall collated annual index of abundance at each site (Botham et al., 2020).
To quantify local site characteristics and capture the focal habitat within survey areas, we analyzed a 500-m radius buffer around the centroid of each of the 539 UKBMS sites, using data from the Natural England priority habitat maps (Natural England, 2019).
These maps capture a range of habitat characteristics, including lowland calcareous grassland (chalk grassland). Using a 50-m resolution digital elevation map (Morris & Flavin, 1990), topographic slope angles were estimated for the 539 UKBMS sites, using a systematic sampling of points at 50-m intervals within the 500-m radii of the site centroids, as described in Oliver et al. (2010). It should be noted that site steepness is positively correlated with increased variation in slope angles, that is, areas with steeper slopes are also more topographically variable (see Appendix 1, Tables A1 and A2).
For the population genetics analyses, distinct categories of sites were required. Sites were defined as either chalk or nonchalk sites based upon the presence of lowland calcareous grassland (Appendix 1, Table A3). The lowest percentage cover was 4.7% ("Dancersend" site). Although this represents a small percentage of the total site, it is worth noting that few sites across all UKBMS sites where lowland calcareous grassland is present (and associated with extended M. jurtina phenology from our monitoring data analysis) are dominated (>50% cover) by lowland calcareous grassland and that 25% of these sites (n = 70) have less than 4.3% cover. All of the chalk sites used in the analysis fall within the interquartile range of chalk cover across all UKBMS sites.

| Drought tolerance experiment
All drought experimentation was carried out following the methodology described in Gibbs et al. (2012). A summary of the methods is provided here. Potted host plants (Poa trivialis) were grown under standard conditions, with each plant watered via individual trays. Once mature, plants were randomly assigned to the treatment groups-drought-stressed or control. Control plants were watered daily from 20 days prior to larval hatching and then throughout the experiment. Plants were never oversaturated but watered enough to prevent soil drying and wilting. Drought-stressed plants received no water from 20 days prior to larval hatching and were then only watered every six days throughout the experiment. This treatment meant that green leaves were available at all stages of the experiment but ensured moderate drought stress occurred. At the end of the experiment, green leaves were still present on all plants. This ensured that food availability was not a factor limiting larval growth and survival. Rainwater was used in both treatments. F I G U R E 2 Fifteen sites around the Chiltern Hills from which Maniola jurtina samples were collected for genetic analysis. Large black circles signify where no chalk grassland was present within the 500 m radius of each site (n = 7), large white circles signify sites where chalk grassland occurred within a 500 m radius of the site centroid (n = 8; percentage cover 4.7%-21%). Circles with smaller dots at the centre were sites from which individuals were also collected for the drought experiment (n = 9). Main towns are marked with black diamonds.  Figure 2). Adults from these source populations were live-captured between the 21st of July and 4th of August, mated with individuals from the same population, and eggs were collected. In a common garden experiment, 12 newly hatched larvae from each source population were raised on three nondrought-stressed (control) host plants (four larvae, originating from the same source population, per plant) and 24 larvae were raised on six drought-stressed host plants (four larvae, originating from the same source population, per plant) under controlled conditions until eclosion, using the methods described in Gibbs et al. (2012).
A higher number of larvae were raised on drought-stressed plants due to an expected higher mortality rate (see Talloen et al., 2004), totaling 108 and 216 larvae on control and drought-stressed host plants, respectively. M. jurtina overwinter as small larvae, during which little growth occurs (Brakefield, 1984). As such, larvae were monitored at three time points: 49 days after the first larval hatch date (pre-overwintering), 162 days after hatching (post overwintering during larval growth), and 309 days after hatching (late larval growth and pupation phase). The number of larvae that survived until the third monitoring point was recorded. Individuals were monitored until they reached the pre-pupa stage, at which point they were removed.

| Phenology
To calculate butterfly flight periods, all weekly fitted count values for M. jurtina abundance were summed per UKBMS site in each year, and the day number of the recording period at which 10% of the total occurred was recorded as the flight period start date. The day at which 90% of the total occurred was recorded as the flight period and mean site aspect (cos((aspect × π)/180), such that 1 = due north, −1 = due south) were included as fixed effects and site and year as random effects. We included mean annual abundance as a covariate in these models because larger populations are likely to have a greater flight period range due to mathematical mean-variance relationship (Taylor, 1961). Northing was included in the model to account for the temperature gradient across the UK, with cooler average temperatures occurring at more northerly locations. This was necessary first because previous studies have shown that M. jurtina flight periods are shorter and begin later at northern latitudes (Brakefield, 1987) and second because temperature has been shown to affect M. jurtina phenology, with a predicted 4.7 and 5.4 days in advance to the first appearance and peak flight dates, respectively, per 1℃ increase (Roy & Sparks, 2000). Easting was included to account for longitudinal differences in site conditions, for example, differing levels of rainfall which can affect butterfly phenology (Roy et al., 2001). Site altitude and aspect were included to account for the effects these two factors might have on local temperatures. To reduce the range of magnitudes across the data, northing and easting were scaled by subtracting the mean from each value, followed by dividing by the standard deviation. Site and year were included as random effects to account for repeated measures at each site and variation in phenology between years, often associated with weather (Roy & Sparks, 2000).
All mixed-effects models were carried out in R (R Core Team, 2020) using the lmer function from the lme4 package (Bates et al., 2015). Model assumptions were checked using diagnostic plots for all mixed-effects models. Diagnostics from the initial model fits demonstrated that phenology at sites with very low abundances was much more variable, violating homoscedasticity. This is likely because at sites with very low abundances, there is increased detectability-related sampling error, increasing the uncertainty of the phenology estimate (McCarthy et al., 2013). To overcome this problem, all sites with an abundance index value of less than 20 were removed from the analysis.
where P is the phenology metric of interest (flight period start, mean, end day, or range), C is the percentage cover of chalk grassland per site, S is the mean slope angle per site, A is the site total abundance, N is the site northing, E is the site easting, H is the mean altitude per site, F is the mean aspect of each site, i is a random intercept for site, y is a random intercept for year, and ε indicates error term with zero mean and normal distribution.
All models were tested for spatial autocorrelation via Moran's I test. Residuals were extracted from each model and run against an inverse matrix of distance between sampling points using the Moran.I function from the ape package (Paradis et al., 2004).

| Drought tolerance
A generalized linear mixed-effects model was used to determine whether larval survival rates varied between sites in association with site characteristics. The model was fitted with a binomial error structure and with host plant drought treatment and percentage chalk cover (geology) as fixed effects with an interaction term, and population as a random intercept (Equation 2). The slope angle was not included due to a 0.8 Pearson's correlation with chalk cover.
where S is the larval survival rate, T is the treatment (drought/control), G is the geology of the origin site (percentage cover chalk grassland), p is a random intercept for the origin population of the larvae, and ε indicates error term with zero mean and normal distribution.
Population structure was estimated using STRUCTURE v.2.3.4 (Falush et al., 2007;Pritchard et al., 2000), using an admixture model and correlated allele frequencies with a 100,000 burn-in and 1,000,000 MCMC replications per chain. The potential number of genetic clusters (K) was tested from one to six, with 20 chains run per K. The likeliest K within the sample sets was estimated using the program STRUCTURE Harvester (Earl & VonHoldt, 2012) and visualized using CLUMPAK (Kopelman et al., 2015). Four separate where F is the pairwise F ST score between each pair of sites, G is the site geology (chalk/nonchalk), D is the Euclidian distance between sites, and ε indicates error term with zero mean and normal distribution.
As pairwise F ST values between sites are not independent, Mantel randomization tests with 999 permutations were conducted to assess whether the predictor variable (geology) was significant following the methodology described in Powney et al. (2012). The number of significantly different groupings within site type pairs was determined via a Tukey HSD test.

| Phenology
All phenology measures were significantly positively associated with differences in chalk cover (flight start date coefficient = 0.07, p = .009; mean date coefficient = 0.14, p < .001; end date coefficient = 0.19, p < .001; flight period range coefficient = 0.13, p < .001,  Table A4). Hence, average flight period dates were later on sites with greater levels of chalk cover or steeper slope angles, and average flight periods were longer on sites with greater levels of chalk cover or steeper slope angles. Northing and abundance were also significantly associated with all four measures of phenology, with two exceptions: (a) northing was not associated with flight period mean date and (b) mean local abundance was not associated with the flight period end date (Appendix 2, Table A4). Estimated model values for Equation (1) regarding abundance and northing can be found in Figures A1 and   A2. Aspect, altitude, and easting were not significantly associated with any measure of phenology; however, aspect and altitude were both marginally significantly associated with the flight period range

| Population genetics
All populations within the 15 sites in southern England displayed high levels of genetic diversity and low levels of genetic differentiation (Appendix 3, Tables A5-A9, Figure A3). In summary, no linkage disequilibrium occurred between any pair of loci (Table A5) (Table A7). All populations displayed a high level of heterozygosity, with high levels of allelic richness and infinite estimated effective population sizes (Table A8). Allelic richness was not significantly affected by site geology (p = .867), with a mean allelic richness of 8.2 for chalk sites and 8.3 for nonchalk sites ( Figure 7).
Pairwise F ST scores between pairs of sites were extremely low (mean = 0.002, variance = 0.00004), and none were significantly greater than zero (Table A9). However, when site pairs were grouped by geology (i.e., chalk and chalk, nonchalk and nonchalk, and chalk and nonchalk), combinations within site pairs had a significant effect (1) Figure 7), indicating evidence of weak population differentiation. The distance between sites had no effect on pairwise F ST (Table 1). No evidence of the population structure was found between these 15 populations. No population was found to be strongly genetically distinct from any other population, regardless of the number or allocation of sites included in the analysis ( Figure A3).

| D ISCUSS I ON
In this study, we quantified characteristics of M. jurtina flight periods with respect to geology and topography. We also determined whether differences in other aspects of ecology (population genetics and drought tolerance) were also associated with the same landscape attributes. We found significant, positive, associations between the phenology of M. jurtina and geology (chalk grassland) and topography (steepness of sites being a general proxy for topographical heterogeneity), that is, key flight dates are delayed with increasing chalk cover and slope angle. These associations remained after accounting for abundance, therefore, aspects of geology and topography are associated with phenology independent of mean local abundance. We found no strong evidence of genetic structuring of M. jurtina populations linked to geology, and only very weak evidence of genetic differentiation among populations. Finally, we found no effect of geology on larval survival (drought response).
Microclimatic heterogeneity may explain the longer flight periods on steeper (more topographically diverse) chalk grasslands.
Habitat and topographic diversity can allow species to persist in F I G U R E 4 Annual mean weekly fitted counts for M. jurtina between 1976-2015, across UKBMS sites split by chalk presence, where chalk sites are classified as sites with a greater than 0% cover of chalk within 500 m of the site centroid. Day number 1 = 1st April i.e. the start of the UKBMS recording window (1) for four measures of phenology for M. jurtina in relation to site mean slope angle areas of suitable microclimate when the surrounding climate is no longer favorable (Bennie et al., 2008), and habitat heterogeneity has been promoted as a method of improving species resilience under climate change (Crick et al., 2020). For example, south-facing chalk grassland hillsides were found to harbour populations of the warmth loving species the silver-spotted skipper (Hesperia comma), absent from other habitat types (Davies et al., 2006). However, increasing ambient temperatures at sites due to climate change has seen an expansion in the local distributions of this species (Lawson et al., 2014;Wilson et al., 2010). If microclimate heterogeneity alone causes the longer flight periods, we might expect to see a two-tailed expansion to the flight period on steep chalk sites, with suitable habitat patches available earlier as well as later in the year.

F I G U R E 5 Estimated model values from Equation
In contrast, we found that all measures of phenology were positively associated with chalk cover, including start date. This means that sites with more chalk have later start dates, creating a long, single-tailed extension to the flight period later into the season. Drought conditions have been shown to lead to lengthened larval development times and later emergence dates, in species such as the speckled wood, as a result of physiological stress (Gibbs et al., 2012).
In habitats with heterogeneous microclimates, such as hilly chalk grasslands, certain microhabitats (e.g., with thinner soils on southfacing slopes) may lead to host plants becoming particularly droughtstressed. This would result in a certain proportion of individuals at a site with delayed emergences and a more protracted flight period overall but one that is single-tailed. One point to note is that the fixed effects in our models account for relatively little variation within the data (7%-15%), and the majority of variation (46%-77%) is explained by the random effects for site and year. This is unsurprising as year captures weather effects, which are known to have a large effect on butterfly phenology (Mills et al., 2017;Roy et al., 2001), although there may be differences between sites that are not captured in our relatively coarse scale topographic descriptors (e.g., local vegetation and microclimatic factors that further mediate phenological responses; Davies et al., 2006;Hindle et al., 2015).
A limitation of this study is that UKBMS data do not fully en-   (Baxter et al., 2017;Goulson, 1993b;Habel et al., 2009;Schmitt et al., 2005;Thomson, 1987). Despite being statistically significant, the differences in genetic differentiation between site types (as indicated by pairwise F ST scores) are extremely low, being below the 0.05% threshold typically viewed as indicative of genetic differentiation (Freeland, 2011). This suggests that populations on chalk sites are marginally more genetically distinct from populations on other chalk sites than those from populations in the surrounding environment. Additionally, no population structuring was found via any combination of sites, possibly due to the dispersal ability of M. jurtina (Schneider et al., 2003) and the ubiquity of its host plants.
Therefore, it appears that all populations included in the study belong to a single, large population, with properties similar to the one at panmixia with random mating. Very low levels of differentiation are present, although insufficient to have any great effect on population structuring. The suggestion that populations of M. jurtina on chalk grasslands form a distinct genetic race is not supported; in fact, the opposite is found, with populations on chalk sites being more distinct from each other, although these levels of differentiation are very low. Therefore, we conclude that the differential phenology associated with geology and topography found in this study is unlikely to be explained by differentially adapted host races. However, it should be noted that due to the high correlation found between chalk percentage cover and site steepness, we cannot determine the effect of topography with this experimental setup. Therefore, caution in interpretation is required as our other analyses have shown that site topography can have an effect on aspects of M. jurtina ecology.
Contrary to our expectations, we found no association between the percentage of chalk cover from source sites and larval survival when exposed to drought conditions. However, these results should be interpreted with caution owing to the relatively small sample size and spatial scale of the analysis, and the fact that slope could not be included in the drought models, despite affecting phenology.
Additionally, in wild situations, larvae would be able to move from plant to plant, ensuring that a sufficient quantity of food could be consumed. In the experimental setup, larvae were constrained to single pots containing food plants and therefore unable to move to fresh sites, even though sufficient green plant material was available throughout the experiment and remained at the end to ensure that food quantity was not a limiting factor in larval growth. Our results suggest that although drought conditions reduce larval survival rates, the effects do not appear to be mitigated by local adaptation specific to chalk sites.
In conclusion, we found butterfly phenology varied at the national scale associated with geology and topography. We found neither evidence of genetic structuring of populations based upon these site conditions nor any differences in drought tolerance.
Future research may benefit from a detailed analysis of other ecological factors influencing phenology such as host plant distribution and quality at different sites. This may allow a greater understanding of why phenology is affected by both chalk percentage cover and site topography. Additionally, factors affecting the potential for local adaptation could also be investigated, for example, slope aspect, microclimate, vegetation cover, and habitat management (Bennie et al., 2006;Brakefield, 1987;van Noordwijk et al., 2012). Such studies will become increasingly important for understanding and predicting species responses to a rapidly changing climate.

CO N FLI C T O F I NTE R E S T
We declare no conflicts of interest. Site characteristics and correlations Note: Data compiled from Natural England priority habitat maps and a 50 m resolution digital elevation map (Morris & Flavin, 1990). Slope Angle = Degrees from horizontal, such that 0 = flat, 90 = vertical. Aspect (East) = Mean Eastness of aspect in landscape around site (Eastness = sin((aspect × PI)/180), such that 1 = due East, −1 = due West). Aspect (North) = Mean Northness of aspect in landscape around site (Northness = cos((aspect × PI)/180), such that 1 = due North, −1 = due South). Altitude = Mean height above sea level (m).