Long‐term temporal patterns in flight activities of a migrant diurnal butterfly

Recent studies demonstrated that the Painted Lady (Vanessa cardui), a cosmopolitan diurnal butterfly performs long‐range migration between subtropical Africa and north‐western Europe, covered by individuals belonging to up to six generations. Here we analyze temporal patterns of complete annual migratory activity of the Painted Lady in Hungary, located in its Central European migratory route, almost completely unstudied before. To do so, we used field occurrence data collected between 2000 and 2019 and estimated temporal patterns in migratory activity by fitting kernel density functions on the daily mean number of individuals and observation frequency. The temporal distributions of kernel density estimates were analyzed as a function of time and key climatic predictors of the study area. We found that (i) the timing of spring arrivals has been advancing; (ii) the relative intensity of the first and last migratory peaks of the Painted Lady significantly increased during the past decades; and (iii) intensity of the last migratory peak is related to the mean temperature of the previous month, inferring that the migration is shifting to earlier dates and their volume of the migration has substantially intensified, evoking mutually nonexclusive, competing hypotheses. Our study indicates the strengthening migration activities of a southerly distributed, long‐distance migrant diurnal butterfly, most probably linked to the northward shift of wintering areas induced by warming trends of the southern parts of Europe. However, the complexity of the likely processes leading to changing migratory strategies calls up for further research in both breeding and wintering areas.


Introduction
During the past decades, evidence has been mounting that insects include a number of taxa which fully or partially migrate between breeding areas and wintering grounds with sufficient nutrient availability, even on large geographical or even intercontinental scales, involving regularity both in time and space (Rainey, 1963;Schaefer, 1969;Urquhart & Urquhart, 1978;Dingle, 1996;Hu et al., 2016). On global scales, among a number of migratory insects, a handful of iconic migrants have also been described, mostly consisting of large butterflies, dragonflies or major crop pests (Holland et al., 2006;Chapman et al., 2015). Moreover, other less conspicuous groups of insects do also migrate, such as hoverflies (Jauker & Wolters, 2008) and aphids (Tenhumberg & Poehling, 1995).
Recently, climate-induced phenological shifts in migration regimes have frequently been documented in a wide range of organisms, including vertebrates, plants and migratory insects (Cleland et al., 2007;Gordo, 2007;Taylor, 2008;Bell et al., 2015). For example, evidence is accumulating that both diurnal and nocturnal lepidopterans respond to current climatic trends by advancing first emergence dates and prolonging late autumn activities (Dell et al., 2005;Sparks et al., 2006;Végvári et al., 2015). These behavioral changes have been shown to affect migration and wintering strategies as well as population dynamics and extinction risks in a number of species (Bale & Hayward, 2010).
For example, several studies have been carried out in the migration strategies of the Monarch Butterfly (Danaus plexippus) (Zipkin et al., 2012;Zhan et al., 2014), conclusions of which highlight the complex relationship between climate and migratory performance and suggest that attempts to understand how Monarchs will be able to cope with predicted climate conditions will be challenging. This implies that current and predicted climatic processes might have a key relevance for the conservation of migratory lepidopterans (Lemoine & Böhning-Gaese, 2003;Oberhauser & Peterson, 2003). Among diurnal migratory lepidopterans, the Painted Lady (Vanessa cardui, Nymphalidae) has been proved to be an optimal model organism for insect migration research, as (i) it exhibits long-distance migratory movements between subtropical Africa and northern Europe, covered by up to six generations (Stefanescu et al., 2013;Talavera & Vila, 2016;Talavera et al., 2018); (ii) although the western flyway connecting NW-Europe and West-Africa has been thoroughly investigated, migratory routes further to the east are practically undescribed (Stefanescu et al., 2007;Stefanescu et al., 2013;Talavera & Vila, 2016;Talavera et al., 2018;Menchetti et al., 2019); (iii) as it is common and nonprotected, Painted Ladies can be easily collected and reared in artificial conditions; (iv) larval food include common plants (e.g., thistles Carduus and Cirsium spp. and nettles Urtica spp.) that are easy to be grown in laboratory studies (Stefanescu, 1997;Janz, 2005; own unpublished observations); (v) its cosmopolitan distribution covers all continents with the exception of South-America. Although the migratory movements of the western population of the Painted Lady have long been studied (Pollard et al., 1998), the effects of climatic variability on migratory strategies have not been investigated before.
Here we aim to analyze temporal patterns of the complete annual migratory activity in the Painted Lady throughout Hungary. To do so, we used field occurrence data collected by volunteers and all of the authors of this work, between 2000 and 2019; thus spanning 20 years, which has been shown to be long enough to detect climatic fingerprints on migration strategies (Végvári et al., 2015). We hypothesized that as a response to current climatic patterns, (i) the spring migration is advancing as a reaction to warming springs in North-Africa and South-Europe; (ii) the population size is growing, as the winter mortality decreases in this warm-adapted species, due to elevated temperatures; and (iii) the intensity of migration is related to mean temperature values and precipitation totals (Holland et al., 2006;Sparks et al., 2007).

Painted Lady data
We collected occurrence data of adult Painted Ladies between 2000 and 2019 using (i) publicly available records of entomological websites, which collate data on a broad range of insects between late February and late November-thus covering the whole length of the migratory period of the easily identifiable study species in the study area-uploaded by citizen scientists the data of whom are regularly controlled by specialists (www.izeltlabuak.hu, lepketerkep.termeszet.org) and (ii) our observations recorded between 2014 and 2019, following the methodology of the Hungarian Butterfly Survey (www.lepkeszet.hu). The observations were collected using standard manual GPS devices or mobile phones, providing location coordinates with a resolution less than 10 meters. Applying these records, first we created a database including the following variables for each observation: date, time of the day (hours and minutes), geographical coordinates, and number of individuals counted along a transect. To control for observation effort, we grouped the complete set of observations by day and calculated the daily mean number of individuals. As the length of survey walks showed no temporal differences, we consider the daily mean value of observed individuals to be unaffected by sampling bias. Further, as we found no effect of weekend days on the number of individuals (linear regression with Poisson error distribution, F 1, 1945 = 1.131, P = 0.288, multiple R 2 = 0.001), we consider the dataset to be unaffected by temporal variance in observation activity.
To estimate potential spatial autocorrelation patterns in the records, we performed Moran's I test for the (i) mean and (ii) median number of individuals as well as for (iii) the observation frequency, defined as the number of species occurrences reported per day for each year. As a result, no spatial autocorrelation pattern emerged in the (i) mean number of individuals (Moran's I, observed = 0.010, expected = −0.0001, N = 1102, P = 0.268); (ii) median number of individuals (Moran's I, observed = 0.009, expected = −0.001, N = 1102, P = 0.310); and (iii) observation frequency (Moran's I, observed = 0.004, expected = −0.001, N = 1102, P = 0.619). This implies that the response variables are not affected by spatial autocorrelation.
To investigate the potential effects of the temporal distribution of sample days, we considered all sampling dates of the citizen science project for collecting butterfly data and analyzed the temporal patterns of the number of sampling days divided into 10-d periods. To do so, for each study year we fitted linear regressions on the number of sampling dates as a function of the number of 10-d periods during the flight season of the Painted Lady. This calculation showed that the number of survey dates were independent of time in all of the study years (linear regression, N = 16, P ≥ 0.227 for all cases). Therefore, we consider that the temporal distribution of the Painted Lady data was not influenced by temporal changes in sampling intensity in our dataset.
To test the effect of elevation on the temporal distribution of Painted Lady records, first we divided the dataset into data with low (< 300 m) and high (≥ 300 m) altitudes, as the elevation of 300 meters above sea level is supported by habitat distributions in the study region and also by the histogram of elevation data in our dataset (Mezősi, 2016). In the next step, we calculated the temporal distribution of the Painted Lady by fitting kernel density estimates and compared the distribution of the (i) Julian date and (ii) relative height of the 95 % quantile of the kernel density peaks by two-sample t-tests, as the kernel density curves are normally distributed. We found that both the 95 % Julian dates and relative heights were not correlated between high and low elevations (t-test, t = −0.055, df = 24.073, P = 0.619 for Julian dates and t = −1.652, df = 28.479, P = 0.109 for relative heights) indicating no difference in kernel density distributions between low and high elevations.

Flight activity
To assess temporal changes in migration phenology, first we computed two metrics of spring migration phenology, which controls for the nonindependency of the records collected on the same observation days: (i) earliest arrival dates, defined as the Julian date of the 5% percentile of the height of first migration wave; (ii) median arrival dates, defined as the Julian date of the 50% percentile of the height of the first migration wave, estimating the arrival time of the bulk of the population; as well as a temporal predictor of the last migratory peak, quantified as the Julian date of the 95% of the height of the last migratory peak, characteristic for late summer migration (Tryjanowski & Sparks, 2001;Végvári et al., 2015).
To identify separate peaks of flying activity which identify separate migratory waves, we applied kernel density estimation fitted on the (i) daily means of individual numbers and (ii) frequency of observations restricted for years with at least 15 individuals in total, using the default function available in the R 3.4.4 statistical programming environment (R Development Core Team, 2018). Kernel density estimates are considered as a 'smoothed' version of a histogram. The height of distinct modes assesses the relative intensity of migration waves. We used the default kernel bandwidth, which was derived from the raw data and which is scale-invariant (R Development Core Team, 2018). For example, in the resulting distributions, unimodal kernel density distribution indicates single migratory wave in a year, while bimodal distributions are proxies for two migratory waves in a given year. Consequently, we calculated the number of flight activity peaks by numerically computing the first derivatives of the curves and identifying local maxima by computationally finding zero slopes of the derivative (Végvári et al., 2015).

Climatic responsiveness
Temporal patterns in activity peaks have been analyzed by two approaches.
First, to investigate temporal patterns in the timing of activity peaks, we fitted linear regressions on the (a) Julian dates and (b) relative intensity (height) of the (i) 5% (FED), (ii) 50% (MED), and (iii) 95% (LED) percentiles of kernel density maxima fitted on the mean number of individuals per Julian date as a function of year. FED-shifts describe the modifications in the migration strategies of the earliest arriving part of the population, which is a robust predictor of population responses in fast changing climatic processes. In contrast, changes in MED represent the temporal shifts in the timing of arrival of the population bulk. Alternatively, changes in LED is a predictor of modification in wintering/hibernation strategies (Diamond et al., 2011;Végvári et al., 2015).
Next, we repeated kernel density analyses separated for the southern and the northern part of Hungary divided by the mean latitude of the latitudinal range of Hungary, as the southern part of the country is a candidate as a wintering area of Painted Ladies (Z. Varga, pers. comm.), which implies that the number of migratory waves can be different from that of the northern part, owing to differential dispersal and migration strategies.
In the second approach, we repeated these analyses for kernel density estimates fitted on the frequency of observations (Supp et al., 2015) To estimate relationships among flight activity and climatic parameters, we conducted two approaches. First, for each year, we calculated the precipitation totals and mean temperature averaged for Hungary and fitted linear regressions on the annual number of flight activity peaks as a function of climatic proxies. Second, for each activity peak, we calculated the precipitation totals and mean temperature for the (i) actual and (ii) previous months as well as for the (iii) previous year, averaged for the study area, and fitted multivariate linear regression on the relative height of the kernel density estimate as a function of Julian date of the activity peak as well as the precipitation sum and mean temperature of the (i) actual and (ii) previous month as well as of (iii) the previous year.

Results
Our dataset includes 1947 records in total, spanning 20 years (2000-2019, Fig. 1). The number of kernel density peaks ranged between 1 and 3 (Fig. 2) considering the whole of the country, and ranged between 1 and 5 in South-Hungary, and between 1 and 3 in the northern part of the country. In the whole of the country, the mean of annual medians of observation data was found to be 22 June (interquantile range: 15 May to 9 July, Table 1).

Mean number of individuals
Temporal trends We detected no temporal trend in the number of activity peaks, identified as the number of local maxima of kernel density estimates (linear regression, b = 0.018, df = 12, P = 0.579). Similarly, no temporal patterns emerged when considering the southern and the western part of the country separately (linear regression, b = 0.206, df = 6, P = 0.155 for South-Hungary and linear regression, b = 0.004, df = 11, P = 0.911 for North-Hungary). We found no correlation between the annual number of migration peaks of the separate regions (Pearson's correlation, r = 0.095, P = 0.824).

Climatic responsiveness
We detected significant advancement of temporal patterns in FED (linear regression, b = −2.788, df = 12, P = 0.017, Table 2, Fig. 3), but not in MED and LED dates (linear regression, N = 13, P ≥ 0.080 for both cases, Table 2). In contrast, temporal patterns in the intensity of activity showed significantly increasing trends in the timing of all migration metrics: (linear regression, df = 12, P ≤ 0.031 for all cases; Table 2; Fig. 4A-C).
Across the 13 years with representative sample size, we found no relationship between the annual number of flight activity peaks and annual precipitation totals or mean temperatures calculated either for the actual or the previous year (linear regression, N = 13, P ≥ 0.598 for all cases, Table 2).
In contrast, the relative intensity of the 95% percentile of the last flight peak was negatively associated to the Julian date (linear regression, b = −0.015, df = 24,  P = 0.047) and was positively related to the mean temperature of the previous month (linear regression, b = 0.111, df = 25, P = 0.008, Table 2). All other associations among 5%, 50%, and 95% percentiles of the kernel density estimates as well as monthly climatic parameters calculated for the actual or previous months proved to be nonsignificant (linear regressions, df = 25, P ≥ 0.055 for all cases, Table 2).

Frequency of observations
Similar to temporal patterns of FED-trends calculated for the mean number of individuals, we detected a significant temporal pattern in FED dates as a function of year (linear regression, b = −3.410, df = 12, P = 0.017), indicating the advancement of the arrival of the first individuals. In line with MED-and LED-trends for the mean number of individuals, no temporal patterns emerged for the timing of the bulk and last of the migrating waves (linear regressions, df = 12, P ≥ 0.075 for both cases).
In contrast to kernel density estimates fitted on the mean number of individuals, the relative intensity of the 5% percentile of the first migratory wave estimated by kernel density distribution fitted on observation frequency showed no temporal trend (b < 0.0001, df = 12, P = 0.459). Whereas the 50% percentile of the first migratory peak exhibited no temporal pattern across the years (b = 0.024, df = 12, P = 0.139), the volume of the last migratory peak, measured as the 95% percentile of the kernel density estimate increased significantly as a function of years (b = 0.161, df = 12, P = 0.043), similarly to kernel density estimates fitted on the mean number of individuals.
Across the 12 years with representative sample size, we found no relationship between the annual number of flight activity peaks and annual mean temperature (linear regression, b = −0.030, df = 11, P = 0.904), as well as between annual precipitation totals (linear regression, b = 0.001, df = 11, P = 0.587). Similarly, the annual number of migratory peak was independent of both annual temperature (linear regression, b = −0.417, df = 11, P = 0.086) and precipitation totals (linear regression, b = 0.001, df = 11, P = 0.648) of the previous year.
The relative intensity of the flight peak was neither related to the Julian date of the peak nor to any of the climatic proxies, calculated both for the actual and previous months (multivariate linear regression, df = 21, P ≥ 0.214 for all cases), similarly to the results of kernel density estimates fitted on the mean number of individuals.

Discussion
Our study provided the following key results. First, we consistently (i.e., supported by models fitted on kernel density estimates computed using both the mean number of individuals or observation frequency) detected the advancement of the arrival of the 5% percentile of the first migratory wave across the study years, which indicates the shift of the timing of the spring arrival of Painted Ladies to earlier dates. Second, the intensity of the first and last migratory peaks (which is a proxy of late summer migration), measured as the 5% and 95% percentile of the kernel density estimates increased significantly between 2000 and 2019. Third, the volume of the 95% of the kernel density estimate of the last peak during late summer migration was positively related to the mean temperature of the previous month, shown only by linear regressions fitted on the kernel density estimates computed using the mean number of individuals.
The advancement of the timing of the earliest arrivals is in line with a number studies that consistently show in a large taxonomic scale that even long-distance migrants change migration strategies by earlier arrivals in the breeding grounds which is considered an adaption to warmer temperatures (Jonzén et al., 2006, Sparks et al., 2007. Indeed, the earliest arriving individuals benefit by (i) reduced competition for resources and by (ii) avoiding predators which focus on the arrival of the bulk of the population, that is, which optimize for matching food peaks (Jonzén et al., 2007). In contrast, the timing of the median Julian dates of the first migratory wave and the 95% percentile of the last migratory wave did not change during the past two decades, implying that both the bulk of the population and those arriving during the late summer migration wave exhibits a conservative migration strategy, which we previously demonstrated in a larger taxonomic subset of a lepidopteran family (Végvári et al., 2015). Our finding that the relative intensity of both spring and late summer activity peaks of the Painted Lady significantly and consistently increased during the past decades infers that the volume of the migration waves of Painted Ladies has substantially intensified, independent of the timing of migration. To explain this pattern, we propose several mutually nonexclusive hypotheses. First, this finding is in line with a number of investigations demonstrating increasing population sizes in insects with southern distribution areas, which benefit from current warming trends (Vanhanen et al., 2007). Indeed, this pattern might be related to decreased mortality due to improved wintering conditions or owing to decreased distances between breeding and wintering sites, thus involving less predation risks and smaller likelihood of adverse climatic events. In line with this hypothesis, current field observations indicate that Painted Ladies might winter as far as both South-and North-Hungary (János Tóth, Zoltán Varga, unpublished data). Second, warming winters might increase the fitness of first emerging individuals, which thus might have higher breeding success that results again in increased population sizes arriving in Central Europe (Sparks et al., 2007). Third, bioclimatic studies indicate that climatic processes might induce substantial changes in migratory routes, which predict that during the study period we might have detected individuals from various wintering areas, the plausibility of which is also coupled with fast changes in precipitation patterns of the Sahel zone in Africa where Painted Ladies are known to winter (Biasutti, 2019). Further, several migratory animal species have changed migration strategies not only by advancing the onset of spring migrations, but also by shifting their wintering grounds northwards (Bókony et al., 2019). Indeed, current climatic scenarios consistently indicate the fast northward shift of the frost-free zones within Europe, supported by the novel overwintering of southerly distributed bird species, not experienced before the 1980s (Leito et al., 2015;Bókony et al., 2019).
Our finding that the volume of the last migratory wave was positively related to the mean temperature of the previous month supported by models developed using the mean number of individuals and not by observation frequency, indicates that increased temperatures improve survival rates, which parallels the preference of Painted Ladies for hotter and drier climates.
Parallel to our findings, the high complexity of relationships among climatic parameters and breeding phenology in Monarch Butterflies along their migratory route has already been documented by previous studies, which found that climate acts in conflicting ways during the spring and summer seasons, as various generations of the same migratory track experience various climatic trends (Zipkin et al., 2012). Similarly, a study on this migrant diurnal butterfly indicates increased climatic sensitivity by the predicted fast northward shift of the wintering areas (Batalden et al., 2014).
In line with our findings, previous studies anticipate the increase of the migratory population of the Painted Lady, eventually leading to its classification as agricultural pest, thus evoking further studies in its migratory behavior and food plant selection that might inform agricultural planning. For example, the caterpillars of migrant Painted Ladies can also feed on soybean and sunflower, if allowed by the timing of migration (Poston et al., 1977;Charlet et al., 1987), which has been shown to be related to current climatic processes (Hódar & Zamora, 2004).
The increasing migratory population size of the southerly distributed Painted Ladies parallels the findings of previous studies showing that the number of species of migratory lepidopterans in the south of the United Kingdom has been rising steadily, which is closely linked to rising temperatures in SW Europe. It is predicted that further climate warming within Europe will increase the number of migratory lepidopterans reaching the United Kingdom and the consequences of this influx also need urgent attention (Sparks et al., 2007).
Migratory animals have been shown to be especially vulnerable to current global climatic trends owing to their dependence on spatially distributed resources, which may be differentially influenced by climate (Lemoine & Böhning-Gaese, 2003). Thus, a number of populations of migratory lepidopterans (Oberhauser & Peterson, 2003) are considered as endangered by current climatic trends, highlighting the importance of longitudinal studies on the climatic responsiveness of model butterflies, that motivated us to analyze the temporal patterns of migratory activity in the Painted Lady throughout Hungary. Indeed, recent studies highlight the highly complex responsiveness of insects to current climatic trends, which is amplified by the interconnectedness of climatic processes acting on various (i) generations; (ii) stages of the life cycle; and (iii) sections of migratory routes (Forrest, 2016). These relationships provide exceptional challenges for predicting changes in phenology, breeding and wintering range as well as population dynamics and extinction risks in migratory insects.
In sum, our study indicates the strengthening migration activities of a southerly distributed, long-distance migrant diurnal butterfly, most probably linked to the northward shift of wintering areas induced by warming trends of the southern parts of Europe. However, the complexity of the likely processes leading to changing migratory strategies calls up for further research in both breeding and wintering areas.