Phenology in adult and larval Lepidoptera from structured and unstructured surveys across eastern North America

and


Introduction
Larval lepidopterans (caterpillars) play a critical role in forest ecosystems. As herbivores they may cause substantial defoliation with consequences for plants including the production of secondary metabolites, birds (Holmes et al. 1979, Holmes and Schultz 1988. The impact of caterpillars as both consumers and resources depends in part on their seasonal period of activity, or phenology, in relation to the timing of seasonal activity of adjacent trophic levels. A changing climate poses the risk that not all organisms will shift their periods of growth, activity, and breeding by the same amount, resulting in phenological asynchrony (Parmesan 2006, Thackeray et al. 2010, Renner and Zohner 2018, Samplonius et al. 2021. Phenological asynchrony between plants and caterpillars could have negative consequences for the caterpillars if they hatch prior to leaf emergence (Abarca and Lill 2015), or after chemical defenses in leaves have increased (van Asch andVisser 2007, Renner andZohner 2018). For foliage gleaning birds, asynchrony between peak caterpillar availability and brood rearing may result in lower fitness , and increasing asynchrony over time may be partly responsible for population declines (Both et al. 2006, Møller et al. 2008. While the phenology of other trophic levels in these tritrophic interactions has been well documented over broad spatial scales-bird phenology via the citizen science project eBird (Sullivan et al. 2009) and forest vegetation phenology with remote sensing (Zhao et al. 2009, Peng et al. 2017)-tracking caterpillar phenology over time at these scales remains a significant challenge due to lack of available data.
Although there are various regional to continental efforts to systematically monitor adult butterflies (e.g. North American Butterfly Association [NABA] counts (Taron and Ries 2015); eButterfly (Prudic et al. 2017)), few standardized monitoring efforts exist for caterpillars on woody vegetation. Those that do tend to be projects focused on single sites, such as at Hubbard Brook LTER in New Hampshire and Coweeta Hydrologic Laboratory in North Carolina. In recent years, the sampling methods employed by those long-term monitoring sites were adopted more broadly by the citizen science project Caterpillars Count! (Hurlbert et al. 2019), which has engaged participants in documenting the phenology of caterpillars (and other arthropods) at more than 100 sites throughout North America. These structured surveys are ideal for accurately representing the phenology of forest caterpillars, however, the sites are sparsely distributed across (mostly) eastern North America, sampling intensity at some sites is quite low, and most sites have only accumulated a few years of data so far. Thus, it is worth exploring the extent to which other unstructured larval and adult Lepidoptera datasets with stronger spatial and temporal coverage exhibit similar phenology and responses to temperature in space and time. For example, hundreds of thousands of opportunistic observations are available for a diverse array of Lepidoptera, including both larval and adult stages, from iNaturalist (iNaturalist.org), a platform for sharing and identifying primarily photo observations uploaded by users. While iNaturalist data have been used successfully to document the adult phenology of individual species or sets of species varying in traits (Li et al. 2019, Barve et al. 2020, Belitz et al. 2020, unknown observation processes coupled with potential biases in space, time, and taxonomic groups (Di Cecco et al. 2021) may introduce substantial uncertainty into estimates of caterpillar phenology at the community level. Compared to structured surveys in forest habitats, iNaturalist observations may underrepresent cryptic groups (e.g. Geometridae, Noctuidae, and Notodontidae) in favor of families that include large, charismatic species (e.g. Sphingidae, Nymphalidae, Papilionidae), and iNaturalist observations come disproportionately from areas of developed land cover (Di Cecco and Hurlbert 2022). However, adult observations are routinely identified to greater taxonomic resolution. Though butterflies represent only a fraction of Lepidopteran diversity, and primarily use open habitats, they are more highly represented in opportunistic observations as well as in documentation of life history traits which impact phenology.
Here, we compare both absolute and relative phenology metrics among three datasets --structured Caterpillars Count!-style surveys, unstructured opportunistic caterpillar observations from iNaturalist, and unstructured adult butterfly observations from iNaturalist. We partition the adult dataset, which is identified to species, by overwintering stage for comparison to phenology metrics from both caterpillar datasets, due to the importance of overwintering stage on flight phenology (Diamond et al. 2011). Specifically, we make the following predictions: 1) We predict that the phenology of all three datasets should exhibit similar spatial and interannual variation due to commonalities of Lepidopteran physiology, despite different taxonomic and habitat representation. 2) We predict that lags (i.e. differences in 10% date) between adult butterfly and caterpillar phenology from either data source will vary in the same order, based on expected adult phenology with species overwintering as pupae appearing first, followed by species overwintering as larvae, and then species overwintering as eggs. We expect these lags will be smaller at higher latitudes where the length of growing season is compressed, and in warmer years where phenology is accelerated for all groups. 3) We predict that phenological responses will be most similar between the two caterpillar datasets, rather than between either caterpillar dataset and the adult butterfly dataset. This expectation is based on presumed differences in taxonomic representation between the three datasets, with greater taxonomic overlap among the caterpillar datasets. Because the majority of forest (primarily moth) caterpillars overwinter as pupae (Wagner 2005), we also expect the butterflies that overwinter as pupae to be the adult subset with the most similar phenological patterns to either caterpillar dataset.

Study area and temperature data
We examined patterns in spring and summer phenology of Lepidoptera in North America east of 100° W and between 35-50° N. This region broadly encompasses the area of eastern deciduous forest with sufficient data density (Supplementary Material, Fig. S1). We estimated phenology across 2000-2020 from three types of Lepidoptera occurrence data: structured larval surveys, incidental larval records, and incidental adult records. To compare across datasets, we aggregated data to a common spatial scale by using hex grids at 70,000 km 2 resolution (approximately 300 km in diameter).
We obtained annual spring temperature data for each hex cell from Daymet (Thornton et al. 2019). Spring temperature in a given year consisted of averaged daily minimum and maximum temperatures at 1 km 2 within each hex cell for the months of April, May, and June. For each hex-year with sufficient Lepidoptera data, we also calculated interannual temperature anomalies, or the difference between spring temperature in a hex cell and the average spring temperature across years in that hex cell. The baseline period over which average temperature was calculated for a hex cell was the set of years for which lepidoptera phenology data were available, which varied by hex cell and lepidoptera dataset (see below). Positive temperature anomalies therefore indicate warmer than average years and negative values indicate colder than average years relative to the set of years for which phenology can also be evaluated.

Lepidoptera data sources
We aggregated structured larval surveys from a citizen science project Caterpillars Count! and a site-specific monitoring project. We first obtained records of caterpillars from structured surveys in the citizen science project Caterpillars Count! (Hurlbert et al. 2019), from 2011 to 2020. Caterpillars Count! observations come from systematic visual or beat sheet surveys of individual branches of woody vegetation, revisited multiple times over the spring and summer. Observers record all arthropods on the survey branch, identifying organisms to order and recording the length. Because we lack finer taxonomic level identifications, caterpillar phenology can only be examined in aggregate. We included data from 20 sites which met several data quality thresholds: at least 20 branches per site, with at least six weeks of data where at least 80% of total branches per site were surveyed on each survey date between April and July. For the analyses comparing phenology across years, we further restricted our dataset to the 6 sites which met those quality thresholds in at least two years. As a supplement to the Caterpillars Count! data, we included a long-term monitoring dataset using the same visual branch survey methodology from two sites near the Coweeta Hydrologic Laboratory in western North Carolina from 2010 -2018 (Lewis 2021). For simplicity, we refer to this combination of data from Coweeta and Caterpillars Count! as the "Caterpillars Count! dataset". Aggregated survey data provided 4,669 individual presence observations of caterpillars at 30 sampling sites across 14 hex cells.
Incidental records of caterpillar occurrences were obtained by exporting observations from iNaturalist that were part of the "Caterpillars of Eastern North America" project as of April 25, 2021 (https://www. inaturalist.org/projects/caterpillars-of-eastern-northamerica). All observations of Lepidoptera that have their life stage annotated as "Larva" on iNaturalist are added to this project. We restricted these observations to the years 2000-2020 and the prespecified eastern North American hex grids. Because Caterpillars Count! observations with photos are automatically shared with iNaturalist, these data were removed from the iNaturalist observations to ensure datasets were independent. We also removed observations with positional uncertainty greater than 50 km. In total, 195,011 observations of caterpillars remained in our iNaturalist caterpillar dataset. Although some of these records were identified to family, genus, or species, we considered caterpillar phenology in aggregate to facilitate comparisons with the Caterpillars Count! dataset. iNaturalist records can be biased by seasonal events such as the City Nature Challenge (Di Cecco et al. 2021), so unusually high numbers of observations on a specific date may reflect observer effort more than actual changes in seasonal abundance. We reduced the influence of dates with high effort by first determining the average number of observations for each overwintering stage in a particular hex grid and year on each calendar day, given at least one observation. Next, we detected dates with a higher than average number of observations in relation to a particular hex grid and year combination. For those dates, we thinned records to the average number of observations for that grid cell and year combination.
We generated a dataset of adult butterfly (Superfamily Papilionoidea) occurrence using 253,840 incidental records from iNaturalist via G B I F ( htt p s : / / w w w. g b i f. o rg /o c c u r re n c e / download/0006251-210914110416597). GBIF hosts research grade observations from iNaturalist, meaning observations are georeferenced, include photos, have a date, are not cultivated, and at least two-thirds of the community agree on the identification (Seltzer 2019). We removed iNaturalist observations that were included in the Caterpillars of Eastern North America iNaturalist project to ensure our dataset was primarily of adult butterfly occurrence records. Because butterfly observations were identified to species, observations were grouped according to the life stage in which they overwinter in the study region (Supplementary Material, Table S1). We were unable to assign overwintering stages to observations in either caterpillar dataset, because many caterpillar observations were not identified to species and because knowledge of overwintering stage is incomplete for those species that were identified.

Measuring phenology
Annual phenology metrics were estimated for each dataset using the ordinal day of presence records to fit Weibull distributions and calculate 10th percentile dates within hex cells. We fit Weibull distributions using the "phenesse" R package (Belitz et al. 2020) with 250 iterations to obtain each phenology metric.
Phenesse uses a parametric bootstrapping approach to generate bias-corrected phenological metrics for any percentile and has been demonstrated to produce accurate and unbiased estimates for simple seasonal abundance curves (Belitz et al. 2020). For iNaturalist caterpillars and adult lepidopteran records, we estimated 10th percentile dates at the hex scale. We generated phenometric estimates if each combination of hex-year (for caterpillars) or hex-yearoverwinter stage (for adults) had at least 10 unique days with at least one observation. For Caterpillars Count! records, we estimated 10th percentile dates at each site, and then averaged across sites within a hex cell because the dataset contained repeated observations throughout the spring and summer at each site within the cell.

Differences in phenology and phenological sensitivity between datasets
Using the 10th percentile date phenology metrics estimated for each dataset, we examined consistency in relative and absolute phenology across time and space between the three datasets. We first examined correlations in interannual anomalies in 10th percentile dates across hex-years between datasets, to see whether early and late years in one dataset corresponded to early and late years in another dataset. Interannual anomalies were calculated as the difference between a phenometric in a particular year from the mean phenometric across years within a given hex cell. For adult lepidopterans, we compared each overwintering group separately to the full set of iNaturalist or Caterpillars Count! larval phenometrics.
To quantify absolute differences in phenology between datasets, we calculated the difference between 10th percentile dates across datasets (again comparing each overwintering adult lepidopteran group to each larval dataset individually), to find the lag time between phenology metrics in days. Caterpillaradult lepidopteran lag times might be expected to be consistently positive or negative based on the life stage of the dataset and the timing of the sampling. If there is no systematic bias between the Caterpillars Count! and iNaturalist datasets, then the expected distribution of lag values across hex cells and years should be centered on zero.
To explain variation in absolute lags between 10th percentile date between datasets across hex cells (i) and years (y), we modeled the absolute value of the difference between 10th percentile dates across datasets (ΔP i,y ) as a function of latitude of each hex cell and interannual anomalies in spring temperature in that hex cell-year (T i,y ): We tested for and found no evidence for strong interactions between latitude and temperature anomalies, so only report the results of this simple additive model.
Finally, we examined the extent to which interannual anomalies in 10th percentile dates (reflecting early or late years) was linearly related to interannual anomalies in spring temperature (relatively cool or warm years) for each dataset, and for each overwintering stage separately within the adult dataset.
Across all hex-years with relevant data, the distribution of differences in 10 th percentile dates between iNaturalist caterpillars and Caterpillars Count! observations were not consistently biased towards positive or negative values, indicating similar estimates of phenology on average ( Fig. 2A). However, differences between datasets were often as great as 4 weeks in either direction depending on the hexyear. As expected, adult butterfly 10 th percentile dates tended to be earlier than those from either larval dataset (Fig. 2B, 2C), with the largest differences between iNaturalist caterpillars and adult butterflies that overwinter as pupae.
Lags in 10th percentile dates between datasets were generally smaller at higher latitudes and in warmer springs (Fig. 3), however, these effects were only strong (p < 0.05) for a few of the comparisons. The lag between iNaturalist caterpillars and adult butterflies overwintering as larvae varied as predicted with both latitude (p < 0.001; Supplementary Material, Table S2) and spring temperature anomaly (p < 0.001), while the lag between iNaturalist caterpillars and adult butterflies overwintering as pupae exhibited a strong temperature effect (p = 0.015). The lag between Caterpillars Count! caterpillars and adult butterflies overwintering as either larvae or pupae showed only weak effects, and the lag between Caterpillars Count! and adults overwintering as eggs was not modeled because there were only two hex cells in common between those groups.
Interannual anomalies in larval 10 th percentile dates from structured and incidental data were positively correlated (r = 0.56, p < 0.001), indicating agreement between the datasets in signals of early or late phenology within hex cells (Fig. 4A). Anomalies in iNaturalist caterpillar 10 th percentile dates were positively correlated with anomalies in adult butterflies (0.26 < r < 0.39 across overwintering stages and phenometrics; Fig. 4B), with butterflies overwintering as pupae having the strongest correlation. Anomalies in 10 th percentile dates for Caterpillars Count! were correlated with adult butterflies that overwinter as pupae (r = 0.61, p = 0.008), but not with those that overwinter as larvae or eggs (p > 0.05; Fig. 4C).
Among the 50 Lepidoptera families (Supplementary Material, Table S3) represented in this analysis, there was relatively low overlap in family composition between datasets. The greatest taxonomic overlap was between the two caterpillar datasets (Jaccard similarity = 0.20), while similarity was low (0.07 > Jaccard > 0.10) for each caterpillar dataset in relation to the adult butterfly dataset, which constitutes only 6 families.
Caterpillar phenology based on the Caterpillars Count! dataset tended to be earlier in hex-years with warmer than average spring temperatures, and later in hex-years with cooler than average spring temperatures, although this explained only 6% of the variance (estimate = -1.42 days/°C, R 2 = 0.06, p = 0.16, Fig. 5A). There was no relationship between spring temperature deviations and deviations in caterpillar phenology based on iNaturalist caterpillars (R 2 = 0.00, p = 0.89, Fig. 5B). Adult butterfly phenology also showed a negative, although noisy, relationship between deviations in butterfly phenology and deviations in spring temperature, with the strongest relationships for the subset of butterflies that overwinter as pupae (estimate = -2.69 days/°C, R 2 = 0.02, p = 0.02, Fig. 5C) and that overwinter as larvae (estimate = -3.04 days/°C, R 2 = 0.02, p < 0.01, Fig. 5C).

Discussion
Understanding how the phenology of forest caterpillars varies over time and space is critical for modeling the predicted impacts of climate change on caterpillars themselves, but also the cascading impacts on adjacent trophic levels such as forest trees and avian predators. Although estimates of caterpillar phenology based on standardized monitoring efforts are ideal, such data are dwarfed by the sheer number of opportunistic observations from projects like iNaturalist, including especially observations of often conspicuous and charismatic adult butterflies. We found that although phenology estimates between datasets rarely aligned in an absolute sense, the degree to which phenology could be considered early or late Figure 2. The distribution of differences in 10% date between datasets across all hex-years. A) shows the distribution of differences between the iNaturalist caterpillars and Caterpillars Count!, B) shows the differences between iNaturalist caterpillars and adult butterflies broken down by overwintering stage, and C) shows the differences between Caterpillars Count! and adult butterflies broken down by overwintering stage. within a given hex cell was generally consistent. This suggests some promise in the use of large opportunistic datasets for inferring changes in caterpillar phenology over time and over broad extents.
A strong positive correlation was found between caterpillar phenology anomalies in Caterpillars Count! and iNaturalist datasets. This correlation is impressive given the different biases in the habitat and land  cover sampled between the two datasets and the spatial resolution of the analysis, which aggregates observations across ~70,000 km 2 hex cells. Caterpillars Count! surveys are conducted on woody vegetation only, and the dataset represents repeated monitoring of the exact same branches over the course of the season. Relative to Caterpillars Count! sites, iNaturalist observations were more likely to come from areas with developed land cover (Di Cecco and Hurlbert 2022), with substantial variation in the specific locations within the hex cell being represented over time.
Other advantages of the Caterpillars Count! surveys are the known, and more or less constant, sampling effort per survey date, and the fact that participants are actively searching for caterpillars and other arthropods which results in the inclusion of even small or cryptic individuals. Because the iNaturalist dataset represents opportunistic observations, the caterpillars reported there tend to be larger and more conspicuous, often reflecting late instars or individuals on the ground during their pre-pupal wandering phase. For example, we examined 27,368 iNaturalist images of Danaus plexippus and Papilio polyxenes caterpillars and found that for both species, over 87% of observations were late (fourth or fifth) instars. While this bias should result in later estimated onset dates for iNaturalist relative to Caterpillars Count!, we found no such tendency across hex-years. It is likely that this bias is offset by the sheer volume of iNaturalist data spread across the entire calendar year which, all else being equal, results in earlier onset dates.
The correlation found between interannual deviations in phenology of the two caterpillar datasets was much stronger than a similar analysis conducted by Di Cecco and Hurlbert (2022). In addition to slightly different hex-years used in the analysis, here, we examined estimates of phenology onset based on the 10% date of a Weibull fit (Belitz et al. 2020) whereas the earlier study examined peak date and centroid date. Onset phenology (i.e. 10% date) should be less impacted by variation in summer conditions and may provide a more consistent signal compared to the other metrics that characterize modal or mean phenology. This may be especially true when considering the phenology of many species in aggregate as we do here. Onset date for the group will reflect the biological response of the subset of early season species (and first generations of multivoltine species), whereas peak, median, or centroid date can be equally affected by responses of the many mid-or late season species (and later generations of multivoltine species), and may thus be inherently noisier.
Despite the fact that the majority of forest caterpillars on Caterpillars Count! surveys were from moth families, we found similar interannual anomalies in adult phenology with iNaturalist butterflies. The correlation was strongest with the subset of adult butterflies that overwinter as pupae, which is the life stage at which most moth species in eastern North America overwinter as well (Wagner 2005). This subset of pupal overwinterers also exhibited the strongest correlation (r = 0.39, p < 0.001) with iNaturalist caterpillar phenology deviations. We expected this correlation to be even stronger, but perhaps even though both sets of observations come from iNaturalist, differences in observation process, organism conspicuousness, taxonomic representation, and habitat bias may exist between larval and adult subsets. Nevertheless, having completed most of their development, Lepidoptera overwintering as pupae are sensitive to temperatures over the narrowest window prior to emergence as adults compared to species with earlier overwintering stages. In contrast, species overwintering as eggs will be sensitive to spring temperatures over the longer period of larval plus pupal development. This finding highlights the importance of using ecologically relevant traits to refine predictions of phenology of Lepidoptera in particular, and insects more generally. While many of the most ecologically relevant traits (e.g., overwintering stage, voltinism, diet breadth) are well known for butterfly species, the vast majority of Lepidoptera are moths for which such information is lacking or incomplete. To better understand forest ecosystem dynamics, more work is needed to compile such information.
While we found a general correspondence in early versus late years between the three datasets, this interannual anomaly in phenology was poorly predicted by interannual anomaly in spring temperature, explaining only 2-5% of the variance. These results reflect a weaker link between phenology and climate than have been demonstrated in several other largescale analyses, highlighting the potential importance of analytical decisions related to characterizing both climate and phenological response. Di Cecco and Hurlbert (2022) found spring temperature deviations explained 36% of the variation in caterpillar centroid date using the same Caterpillars Count! dataset, although they found this centroid date was unrelated to temperature deviations for iNaturalist caterpillars. Caterpillars Count! surveys tend to be focused on the period of May through July while iNaturalist observations potentially span the entire year, which may be another reason for observed differences in sensitivity to a particular climatic window. Larsen et al. (in review) examined the impact of climate variables that incorporated both temperature and forest greenup on adult butterfly phenology, and found a strong signal of earlier onset with warmer springs and earlier green-up. However, a direct comparison of these results is difficult as the latter study used growing degree days rather than mean temperature in a more complex model involving many additional covariates.
Our results suggest that interannual anomalies in forest caterpillar phenology measured in structured surveys are reasonably well captured using opportunistic observations despite limited taxonomic overlap. However, there are a number of opportunities for improving this relationship. First, we need models that can produce accurate absolute rather than relative estimates of phenological timing if we are to address phenological mismatch between caterpillars and adjacent trophic levels. Ideally, a new modelling framework can be developed that facilitates the integration of diverse sources of structured and unstructured data in a way that mitigates rather than amplifies their inherent biases (Isaac et al. 2014, Isaac et al. 2020. Phenomenological models based on existing spatiotemporal variation and taking into account predicted lags as a function of latitude and climate are one obvious ingredient. A promising alternative for modeling the phenology of caterpillars from the phenology of more conspicuous adults is a mechanistic "caterpillar-cast" that hindcasts or forecasts (depending on overwintering stage) caterpillar presence based on species-specific growing degree day thresholds. Scaling this approach up to the entire caterpillar assemblage will require a more complete understanding of life history traits and developmental thresholds across a wider range of species and establishing which subset of data-rich adult butterfly species have larval phenologies most relevant to forest caterpillars. Finally, phenological mismatch with other trophic levels will depend as much on the total abundance or biomass of caterpillars as the seasonal timing (Durant et al. 2005, Shutt et al. 2019. Estimating abundance from opportunistic data remains a challenge (Di Cecco et al. 2021, Rapacciuolo et al. 2021, thus highlighting the need to continue and expand standardized monitoring efforts such as Caterpillars Count!

Author Contributions
GD, MB, EL, LR, RB and AH conceived of the project. RC and WL collected data associated with the Coweeta Hydrologic Laboratory. GD and MB conducted analyses. GD, MB, and AH wrote the first draft of the manuscript, and all authors contributed substantially to revisions.

Data Availability
Data used in this project and code to replicate analyses are archived in a Zenodo repository (https:// zenodo.org/record/5942449).

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Table S1. Overwintering stages of butterfly species in the study region. Table S2. Linear model results predicting lag dates between 10th percentile phenometrics across datasets. Table S3. Lepidoptera families included in analyses by dataset. Figure S1. Data availability for each dataset in years (time window from 2000-2020).