Thermally moderated firefly activity is delayed by precipitation extremes

The timing of events in the life history of temperate insects is most typically primarily cued by one of two drivers: photoperiod or temperature accumulation over the growing season. However, an insect's phenology can also be moderated by other drivers like rainfall or the phenology of its host plants. When multiple drivers of phenology interact, there is greater potential for phenological asynchronies to arise between an organism and those with which it interacts. We examined the phenological patterns of a highly seasonal group of fireflies (Photinus spp., predominantly P. pyralis) over a 12-year period (2004–2015) across 10 plant communities to determine whether interacting drivers could explain the variability observed in the adult flight activity density (i.e. mating season) of this species. We found that temperature accumulation was the primary driver of phenology, with activity peaks usually occurring at a temperature accumulation of approximately 800 degree days (base 10°C); however, our model found this peak varied by nearly 180 degree-day units among years. This variation could be explained by a quadratic relationship with the accumulation of precipitation in the growing season; in years with either high or low precipitation extremes at our study site, flight activity was delayed. More fireflies were captured in general in herbaceous plant communities with minimal soil disturbance (alfalfa and no-till field crop rotations), but only weak interactions occurred between within-season responses to climatic variables and plant community. The interaction we observed between temperature and precipitation accumulation suggests that, although climate warming has the potential to disrupt phenology of many organisms, changes to regional precipitation patterns can magnify these disruptions.

The timing of events in the life history of temperate insects is most typically primarily cued by one of two drivers: photoperiod or temperature accumulation over the growing season. However, an insect's phenology can also be moderated by other drivers like rainfall or the phenology of its host plants. When multiple drivers of phenology interact, there is greater potential for phenological asynchronies to arise between an organism and those with which it interacts. We examined the phenological patterns of a highly seasonal group of fireflies (Photinus spp., predominantly P. pyralis) over a 12-year period (2004-2015) across 10 plant communities to determine whether interacting drivers could explain the variability observed in the adult flight activity density (i.e. mating season) of this species. We found that temperature accumulation was the primary driver of phenology, with activity peaks usually occurring at a temperature accumulation of approximately 800 degree days (base 10°C); however, our model found this peak varied by nearly 180 degree-day units among years. This variation could be explained by a quadratic relationship with the accumulation of precipitation in the growing season; in years with either high or low precipitation extremes at our study site, flight activity was delayed. More fireflies were captured in general in herbaceous plant communities with minimal soil disturbance (alfalfa and no-till field crop rotations), but only weak interactions occurred between within-season responses to climatic variables and plant community. The interaction we observed between temperature and precipitation accumulation 2016 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
Much can be learned about biological systems by observation alone [1], and observational data are often captured incidentally as a result of human activity [2]. Incidental data can range from the very informal and uncontrolled (e.g. comments on a topic in a Web forum) to highly controlled and meticulously collected (e.g. unused data from scientific experiments). Indeed, research activities can produce systemic observational data of very high quality; for instance, insect trapping systems seldom only capture target taxa. This 'by-catch' can provide data that support investigations into entirely uninvestigated phenomena. In this study, we examine one such 'by-catch' dataset: a 12-year time series of firefly observations in southwestern Michigan, for their responses to environmental and habitat conditions.
Over 2000 species of firefly (Coleoptera: Lampyridae) have been identified across various temperate and tropical environments around the world [3]. As larvae, species within the family Lampyridae spend much of their time living underground feeding on earthworms, molluscs and other subterranean invertebrates [4]. As adults, most species abstain from feeding [5], with the exception of the species Photuris pennsylvanica, of which the female is a voracious predator of both conspecifics as well as other insects [5][6][7]. Few studies have been conducted on firefly conservation and broader-scale ecology in relation to changing environments and land uses, and little is known about how environmental parameters drive firefly life history. It has been demonstrated that the life history of at least one species of firefly is temperature-dependent; researchers found that P. pennsylvanica adult emergence could be artificially accelerated by exposing larvae to increased soil temperature [8]. However, much of the primary research on fireflies has focused on the bioluminescent properties of the firefly [9][10][11][12][13][14], while research describing basic population and community ecology of this important family is lacking.
In addition to the scientific importance of the Lampyridae for their bioluminescent properties and as model organisms for evolutionary investigations, fireflies are also among the most widely recognized and culturally valued insect families among non-scientists. Two US states have designated the firefly as their 'State Insect' [15]. Notably, fireflies also feature prominently in Japanese culture, where they have been designated as national natural treasures in many districts and have been used to generate support for biodiversity conservation efforts in Japanese agricultural regions [16][17][18]. They have also been touted as useful classroom tools for sparking student interest in biology [19]. Because of their popular appeal, it is unsurprising that public concern has grown about apparent declines of firefly populations from regions around the world where they occur [20].
Considering the paucity of ecological information about fireflies, their widespread popularity, ease with which adults are observed and concerns about their population viability, fireflies represent an ideal species for citizen science investigations. Citizen science efforts are currently underway seeking to gain information about the status, geographical distribution and phenology of fireflies [21][22][23], and peer-reviewed publications on fireflies have already been produced based on these volunteergenerated data [24,25]. The popularity of fireflies gives them great potential as a flagship and umbrella conservation species and potentially an indicator species of ecological degradation in agricultural regions [26]. However, to our knowledge, no long-term systematic study of firefly phenology and responses to environmental drivers has been published.
Phenology plays a significant role in regulating species abundance, distribution and biodiversity [27,28]. The timing of phenological events in insect life histories is strongly linked to climatic conditions [29][30][31], such as temperature and precipitation [27,32,33]. Changes in phenology can have communitywide consequences, and differential responses among various species within a community can lead to trophic mismatches [28,30]. For example, the timing of larval winter moth (Operophtera brumata) emergence was formerly largely synchronized with oak (Quercus robur) bud burst. Caterpillars that emerge too early lack a sufficient food source and will starve, while caterpillars that emerge too late will be exposed to older, poor-quality leaves, leading to negative physiological implications [34]. Increased spring temperature has resulted in changes in the timing of oak bud burst. However, the winter moth has yet to adapt to changing temperatures, which has led to disrupted synchrony between these two species [34]. Thus, phenological shifts can have both top-down and bottom-up consequences extending throughout multiple trophic levels. Long-term observations are important for understanding ecological trends and the merit of phenology as a predictor of ecological consequences. A long-term study on the Genji firefly (Luciola cruciate) in Japan found that population patterns changed in response to rainfall, potentially leading to early larval emergence and reduced foraging [35]. However, the ways in which climate change and other environmental events have impacted firefly species is less understood.
Developing a model for the emergence of adult fireflies is key to developing our understanding of firefly phenology, which can then be used to expand firefly conservation efforts, educational outreach and environmental research and to predict peak firefly display. In this study, we examine a 'by-catch' dataset documenting captures of fireflies at the Kellogg Biological Station over a 12-year period and place it in the context of other available data to gain insights into the long-term dynamics and phenology of this charismatic, but understudied taxon.

Data sources
Data were obtained through two publicly available datasets-a weather dataset that included daily maximum and minimum temperature and precipitation, as well as a dataset that focuses on ladybeetle observations, but also documents captures of the other insect species. Both datasets arise from Michigan State University's Kellogg Biological Stations (KBS), located in southwestern Michigan. The firefly abundance data were collected as a part of the KBS Long Term Ecological Research (LTER) Site within the Main Cropping System Experiment (MCSE) and forest sites starting in 2004. Fireflies were recorded to family alone; however, from spot-checks of the collected data, it appears that fireflies collected belonged to the genus Photinus, mainly Photinus pyralis, the big dipper firefly, although captures of other species cannot be excluded.
Within the MCSE, seven plant community treatments were established in 1989, ranging from a threeyear rotation of annual field crops (maize, soyabean and wheat) under four levels of management intensity (conventional, no-till, reduced input or biologically based), to perennial crops including alfalfa, poplar and early successional vegetation (i.e. abandoned agricultural fields maintained in an early successional state by yearly burnings; table 1 and figure 1). Each of these treatments is replicated six times across the MCSE site with each replicate consisting of a 1 ha plot. We also included three forest sites in our analysis; these sites were established in 1993 within 3 km of the MCSE site on KBS and represent one of three plant community treatments: conifer forest plantations, late-successional deciduous forest and successional forest arising on abandoned agricultural land (table 1 and figure 1). Forested treatment plots are also 1 ha in size but are replicated three times for each treatment.
Observations were taken on a weekly basis throughout the sampling season at five sampling stations within each replicate (both MCSE and forest sites). These insect abundance data are available publicly, online at http://lter.kbs.msu.edu/datatables/67. Insect abundance monitoring was done using unbaited two-sided yellow cardboard sticky cards (Pherocon, Zoecon, Palo Alto, CA, USA) suspended from a metal post within each sampling station, 1.2 m above the ground. Cards were deployed each week for a one-week exposure for the duration of the growing season. Sampling start and end dates varied each year depending upon planting date of the various crops; however, the length of sampling was fairly consistent (14 ± 1 weeks, on average, per year).
In addition to plant community treatment management information, we also included weather as an environmental factor to explain firefly abundance. These data were also obtained through a publicly available dataset, online at http://lter.kbs.msu.edu/datatables/7.

Data preprocessing
All analyses performed are available as an R script at https://github.com/cbahlai/lampyrid/. Analyses presented in this manuscript were run in R v. 3.3.1 'Bug in Your Hair' (R Development Core Team 2016). Firefly data were extracted from the database held at the KBS data archive and combined with relevant agronomic data (which are encoded in plot and treatment numbers in the main database) and are hosted at figshare at https://ndownloader.figshare.com/files/3686040. Data were subjected to quality control manipulations to remove misspellings in variable names that had occurred with data entry. Observations with missing values for firefly counts were excluded from analysis. Because subsample data were zero-biased, we used reshape2 [36] to sum within date, within plot observations, and created an additional variable to account for sampling effort (which was usually consistent at five traps per plot per sampling period, but on occasion traps were lost or damaged).  Weather data (daily maximum and minimum temperatures were reported in degrees Celsius and daily precipitation in millimetres) were downloaded directly from the Kellogg Biological Station Data archive (http://lter.kbs.msu.edu/datatables/7.csv). To overcome errors in calculations requiring accumulated annual weather data caused by rare missing data points (most often occurring during winter, in periods of extreme cold leading to equipment malfunction), we created a function to replace missing values in the temperature data with the value that was observed for that variable from the day before the missing observation.
We created a dummy variable for 'start day' to enable the user to test the sensitivity of our conclusions to varying our within-year start of accumulation of environmental conditions. We empirically determined (by the Akaike information criterion (AIC)) that 1 March (start = 60) provided the best compromise between capturing early growing season weather variation and negating brief variation in winter conditions; however, the selection of the precise day did not dramatically influence the overall trends in the results unless changed by more than 15 days.
We then created a function to calculate daily degree-day accumulation and season-long degree-day accumulation based on Allen's [37] double sine function, using our daily maximum and minimum temperature data. We created a dummy variable for our minimum development threshold to facilitate sensitivity analysis, but set it to a default value of 10°C. We did not use a maximum development threshold in our calculation, assuming that temperatures exceeding its hypothetical value (often more than 30°C for temperate insects) were relatively rare. Accumulations were calculated from the start day variable, as described above. We also created functions that calculated the accumulation of precipitation over the sampling week, the accumulation of precipitation over the growing season, from the start date and the number of rainy days in a sampling period. Weather data were merged with firefly data to facilitate subsequent analyses.

Data analysis
We used ggplot2 [38] to visualize trends in captures of fireflies by plant community treatment over years. We then conducted a multivariate analysis to determine whether firefly plant community use patterns changed within or between years, and what environmental factors were associated with plant community use patterns. To accomplish this, data were cast as a date-by-treatment matrix at two resolutions (weekly observations and yearly observations), transformed using the Wisconsin standardization, and Bray-Curtis differences were subjected to non-metric multi-dimensional scaling (NMDS) in vegan [39]. Environmental parameters were fit to the NMDS plots using envfit to determine whether patterns were influenced by weather.
To examine patterns in firefly captures over time, and the interactions of these captures with environmental variables, we visualized trends in capture data by sampling week and degree-day accumulation. Noting that degree-day accumulation was associated with the clearest patterns in firefly captures (see Results), with some variation owing to plant community, we built a generalized linear model (GLM) with a negative binomial structure to explain these patterns. The model included the degree-day accumulation in linear and quadratic forms as continuous variables, year and plant community treatment as factors, and trapping effort as an offset variable (to account for lost or compromised traps). Model structure was determined empirically by AIC. After fitting the model, we used the resultant regression parameters to generate predicted values, so we could visually compare the performance of the model to the raw data.
Because the model found year-to-year variation in the activity peak that was not explained by degree-day accumulation, we extracted the activity peaks from each year as predicted by the model to a new data frame, and matched these data to other relevant environmental variables in the weather matrix (week the peak occurred in, precipitation variables corresponding to that week). We visualized the relationship between activity peak and other variables, and then constructed a generalized linear model for a quadratic relationship between the activity peak, by degree days, and the precipitation accumulation at the activity peak.
For all frequentist analyses, a significance level of α = 0.05 was used.

Results
Over   is not a strong driver of firefly emergence. However, plotting firefly numbers instead against degreeday accumulation dramatically reduced the asynchrony of emergence peaks and indicated that a single activity peak occurred in each year (figure 5b). Thus, degree-day accumulation appears to be a better predictor of firefly populations than week of year or associated variables. Our model for firefly activity incorporating degree-day accumulation, plant community treatment and year performed well at predicting the timing of the activity peaks (figure 6), accounting for more than 40% of the variation in the raw data. However, model selection favoured the inclusion of a year term as a factor, suggesting that another factor in addition to degree-day accumulation was varying from year to year and impacting firefly activity. Activity peaks varied from year to year by nearly 180 degreeday units, varying from 720 ± 38 DD in 2004 to 898 ± 55 DD in 2012 ( figure 7). However, we found the year-to-year variation was well explained by precipitation accumulation: a quadratic relationship occurs between degree days at peak emergence and precipitation accumulation (pseudo-R 2 = 0.456, p = 0.026; figure 8).

Discussion
The greatest proportion of fireflies was captured in alfalfa and no-till plant communities (figure 2), indicating that areas with moderate soil disturbance and primarily herbaceous plant communities favoured firefly emergence. This result was unexpected; because fireflies spend much of their life cycle in the soil, it might be expected that plots with little soil disturbance (coniferous, deciduous and   . Number of firefly adults captured, as predicted by GLM and, as observed, by observation number. Predicted values were generated using GLM accounting for variability owing to plant community treatment degree-day accumulation, and year as a factor variable. Details of GLM can be found in the data analysis section in Material and methods. successional forests) would foster the greatest populations of fireflies. However, these plots produced capture rates similar to those observed in the intensively managed and tilled conventional row crop plots. Our result contrasted with observations of another genus of fireflies in Malaysia (Pteroptyx), where researchers found that plant canopy structure was the most important determinant of abundance [40]. Also surprising was the relatively low capture rate in early successional plots, which are primarily herbaceous, with a no-till management regime. Thus, the yearly burnings may play a role in suppressing firefly populations in these plots. An alternative explanation for these variations in captures could be   differences in trapping efficiencies between plant communities. However, if this were the case, we would expect trapping efficiencies in the three other row crop treatments (conventional, organic and reduced input management) not to differ appreciably from that of the no-till row crop plant community. When plotted over sample years (figure 3), captures of fireflies by treatment seem to suggest an intriguing cyclical dynamic, with alternating peaks and troughs in captures on an approximately 6-year cycle. Our time series only spans 12 years, meaning more data will be required to elucidate this pattern and its drivers. Similarly, analysis of plant community use patterns was inconclusive (figure 4). At the weekly resolution there was a trend away from woody treatments over the growing season (figure 4b; i.e. with both increasing week and increasing degree-day accumulation). Although this pattern was not strong, it could result from fireflies overwintering in forest habitats and then moving to lower-canopy herbaceous habitats for mating displays. We observed very similar performance of both degree day and week, probably owing to autocorrelation between the two variables that cannot be resolved at the sampling resolution used over the course of the study.
The degree-day model GLM suggested that activity peaks occurred at a degree-day accumulation of approximately 800 DD, accumulated from March 1. A model for Photinus carolinus in the Great Smokey mountains found peak display occurred at approximately 1100 degree days (using a base of 10°C and the same start date as our model) [25]. The difference between heat units required for peak activity observed between this and our study may be a result of species or locality differences (i.e. more southern firefly populations are probably adapted to warmer spring conditions). Similarly, differences in methodology for calculating degree days may account for some of these differences. However, both studies support the observation that degree-day accumulation is the dominant cue governing the activity patterns of temperate fireflies.
Although both photoperiod and degree-day accumulation can both play a role in the phenology of insects, our results suggest that degree-day accumulation is the dominant driver of firefly flight activity. The model was unable to account for between-trap variation within a single sampling day (figure 6), though it was able to capture the overall trends in activity quite well, using only degree-day accumulation, plant community treatment and year as predictors. Nevertheless, degreeday accumulation was not the sole driver in within-season variability. Our model found year-to-year variability in activity peaks that could not be explained by degree-day accumulation alone. We found that this variation in activity peak by degree-day accumulation had a quadratic relationship with precipitation, indicating that both drought and heavy rainfall in the time period leading up to their activity peak can delay the peak (figure 8). Assuming an approximate 20°C daily average temperature at this site in late June and early July, this could translate to a 10-or-more day change in activity peak owing to precipitation extremes in any given year. Yet, there are several alternate explanations for the pattern we observed, and some patterns detected may have been driven by statistical outliers. For example, very high rainfall in 2015 at the site strongly influences our conclusion that a quadratic relationship exists between degree-day accumulation and precipitation accumulation in explaining firefly activity peaks. Indeed, if observations from 2015 had not been included in our analysis, we would probably have concluded that degree-day accumulation had a negative, linear relationship with precipitation accumulation; that is, increasing rainfall would cause fireflies to emerge earlier, given constant temperatures. This result would align with previous work showing that firefly abundance in Japan is generally negatively correlated with rainfall [35]. However, considered within the context of firefly biology, it seems unlikely that the relationship between these parameters would be linear throughout the range of possible precipitation values, as soil-dwelling larvae of non-aquatic firefly species and/or their prey are probably adversely affected by abnormally waterlogged soils.
As the sampling at our study site continues, we will watch rainy years with particular interest to determine whether population data collected in these years support or refute this pattern, or if an alternate driver can explain more of the variation. Indeed, firefly activity may have been driven by factors not considered in this study. Although using a start date of 1 March was favoured in our analysis (i.e. the AIC of the model using this start date was minimized), when the start day was changed in a sensitivity analysis, the relationship between degree-day accumulation and precipitation in firefly activity changed or disappeared. This result could suggest that alternate drivers not accounted for in this study may be driving aspects of firefly activity. Factors such as winter snow cover and variations in winter temperature are known to affect the phenology of temperate insects [41], and thus these factors should be considered in subsequent work.
In this study, we have clearly demonstrated a taxon whose phenology varies in response to multiple drivers. Species with phenological responses to multiple drivers are not rare [42]. Yet ecological interactions among species with multiple drivers of phenology may be complex and unpredictable [43,44], potentially leading to dire consequences in a changing environment [45]. Our study examined the phenological responses to environmental conditions of adult fireflies; however, data on larvae or sex of the adults were unavailable. Adult Photinus fireflies are non-feeding [14], so shifts in their activity are unlikely to have direct consequences through phenological asynchronies. Shifts in adult activity probably correspond to shifts in development or activity among larvae, potentially leading to asynchronies between larvae-prey populations at this critical development time period. Resources acquired during the predaceous larval stage are important in determining mating success among adult fireflies: males provide an energetically costly nuptial gift to the female in the form of a spermatophore [46]. If sex differences in phenological responses to environmental conditions exist, asynchronies between males and females may additionally reduce mating success and fecundity [47]. Male fireflies were always observed earlier than females in Elkmont, TN, USA. In fact, in that system, females were often found during or after peak emergence of males and thus this should be an area of emphasis in future study [25]. Additionally, phenological shifts in fireflies may lead to consequences at other trophic levels. For example, generalist ground-dwelling predators like firefly larvae and other predaceous beetles are known to have dramatic effects on the establishment of agricultural pests early in the growing season [48]. Similarly, although distasteful and avoided by many predators, some birds, lizards and frogs are known to feed on adult fireflies [49], thus shifts in firefly activity may have dietary consequences for animals at higher trophic levels.

Conclusion
Fireflies are a charismatic and important taxon with ties to trophic function, economic importance and culture. Although empirical evidence of specific declines of Photinus fireflies has not been clearly demonstrated in longitudinal studies, naturalists and citizen scientists perceive a decline in their number [21], leading to interest in their conservation. Our study has offered new insight to support conservation efforts and to direct future research. Photinus pyralis appears to thrive in habitats with moderate soil disturbance. Thus, efforts to foster no-till and perennial agricultural systems [50,51] will probably benefit the species. Climate warming may advance the activity of fireflies to progressively earlier in the growing season, but other extremes of climate in the form of precipitation may introduce unpredictable elements to this, and add the possibility of inducing asynchrony with other systems.
The availability of long-term observational data, made freely accessible to the public, was an essential factor in the discoveries made in this study. Although the study that provided these data was not initiated with this purpose in mind, we were able to empirically demonstrate and disentangle the effect of multiple drivers on firefly phenology simply because we had the statistical power to do so. Although species that respond to multiple, interacting environmental drivers are relatively common, data supporting investigations of this kind are rare [52]. We therefore encourage all practising ecologists to curate their species observation data and make them publicly available, to foster long-term, broad-scale investigations in the future [53][54][55].
Data accessibility. All data and analysis code produced by this study are publicly available. Lampyrid abundance data are available at https://ndownloader.figshare.com/files/3686040. Weather station data are available at http://lter. kbs.msu.edu/datatables/7.csv. These data are automatically downloaded when the R script file at https://github. com/cbahlai/lampyrid/ is run.