Modelling asynchrony in phenology considering a dynamic representation of meteorological variables

Study area and species

Ireland is situated in the north-west of Europe in the North-East Atlantic basin, presenting a temperate maritime climate (McCarthy, Gleeson & Walsh, 2015). Future annual temperatures (mean, maximum, minimum) and rainfall (annual, seasonal) are expected to increase in both the short-term (2040) and long-term (2100) (O’Brien & Nolan, 2023; O’Brien et al., 2024; Wang et al., 2024), which will have implications for the timing of phenological events. The landscape is fragmented with a patchy configuration of grassland pastures, agricultural areas, and forest, with the predominant types of woodland being broad-leaved, bog, coniferous, mixed, and transitional woodland-shrub (Carlier et al., 2021).

Changes in phenology in the Republic of Ireland (hereafter Ireland) are well documented, particularly advanced spring phenology for different species groups (O’Neill et al., 2012; Donnelly, Salamin & Jones, 2006; Donnelly, Yu & Liu, 2015). For example, earlier time of emergence and extended flight season of moths (O’Neill et al., 2012), advancements and/or delays in time of arrival and departure of migrant birds as well as extended length of stay (Donnelly, Geyer & Yu, 2015) and earlier start and duration of growing season in trees (Donnelly, Salamin & Jones, 2006). Consequently, we have selected Ireland as a case study to apply our models for phenological research while emphasizing these methodologies have a broad scope and can be applied to phenological systems elsewhere.

We explored phenological synchrony/asynchrony for three different trophic levels: woodland vegetation, insects (lepidoptera; butterflies and moths), and birds. Species selection for insects and birds was performed to select common and migrant species, respectively, while screening for species with a high number of observations and preferably those that are interlinked through the trophic network. Table 1 shows the common names and codes used for animal species and woodland vegetation habitats. Indicator species for the woodland habitats (canopy, shrub, and field layer) are listed by Perrin et al. (2008).

Species data

Start of season (hereafter, “green-up”, “GU”) dates were calculated from national scale Landsat data. These data were obtained in a pre-processed form from the United States Geological Survey (USGS) EROS Science Processing Architecture On Demand Interface (ESPA) (https://espa.cr.usgs.gov/). Spatially and temporally explicit data on insects (first flight) and birds (date of arrival) were collated from citizen science projects. Butterfly data were obtained from The Irish Butterfly Monitoring Scheme from the National Biodiversity Data Centre (https://www.biodiversityireland.ie), and moth data from MothsIreland (http://www.mothsireland.com), and bird data from eBird (2017) (Sullivan et al., 2009) and BirdTrack (https://www.bto.org/our-science/projects/birdtrack).

Meteorological data

Three meteorological drivers were selected for inclusion as explanatory variables in the statistical models: maximum temperature (TMAX), minimum temperature (TMIN), and total precipitation (TPPT). These variables were calculated using hourly values obtained from Met Éireann’s ReAnalysis (MÉRA, Whelan, Gleeson & Hanley, 2018; Gleeson, Whelan & Hanley, 2017) data set. This high resolution regional reanalysis spans an area covering Ireland, the UK and part of northern France on a 2.5 km horizontal grid from 1981 to August 2019. Daily precipitation totals and maximum temperatures were computed using the sum and maximum value respectively between the period 9 UTC to 9 UTC the following day. Daily minimum temperatures were computed using the minimum temperature for the same period, but the value is assigned to the following day. We utilised a high-resolution (30 m) digital elevation model (DEM) dataset to apply a correction to the temperatures (daily TMAX was too low and daily TMIN was too high on average), to account for a mismatch between the orography in the HARMONIE-AROME (Bengtsson et al., 2017) model used for MÉRA and the actual orography (Gleeson & Whelan, 2020). The DEM dataset was firstly projected onto the MÉRA grid using the nearest neighbour method in Climate Data Operators (CDO). The height differences between the MÉRA orography and those of the DEM were used along with the international standard atmosphere lapse rate of 6.49 K/km, as defined by the International Civil Aviation Organization (ICAO), to apply temperature corrections based on the height differences between the model and DEM orographies.

Data processing

Green-up dates in the form of Julian day were extracted from daily Normalised Difference Vegetation Index (NDVI) values derived from Landsat imagery using the half-amplitude method (Misra et al., 2018; White et al., 2009). Pre-processed Landsat time series NDVI data masked for cloud cover for the years 2007–2018 (See Zenodo for detailed methodological steps). The National Survey of Native Woodlands 2003–2008 (Perrin et al., 2008) and Ancient Long-Established Woodlands (Perrin & Daly, 2010) survey-based polygons were intersected to maximise the information available within each dataset obtained from the National Parks and Wildlife Service, Ireland (NPWS, 2012; NPWS, 2019). The resulting intersected polygons consisted of areas that have been pristine for the last two centuries, have 100% coverage of trees and information on the dominant species. These polygons were used to extract mean values of NDVI per raster layer and aggregated to monthly maximum values to reduce the frequency of missing observations due to the previous pre-processing (Fensholt & Proud, 2012). Residual missing observations in the NDVI time series of each polygon were then filled using mean annual values, also referred to as filling with climatology (Kandasamy et al., 2013). Each NDVI time series was subsequently smoothed to remove spikes in the time series data and interpolated to daily values using a LOESS function (Hufkens et al., 2019; Yuan et al., 2021).

First flight (FF) and date of arrival (DA) were calculated for insects and birds, respectively. To ensure sufficient species records to allow a robust sample size for use in the statistical models and to overcome uncertainty associated with whether the first observation of species sightings from the citizen science datasets was reflective of observer effort, we re-scaled all data to a spatial resolution of 10 km × 10 km. We required a minimum of 3 species sightings per 10 km grid to be included in the statistical analysis, with a minimum of 20 grids required per species and per year (see Zenodo). The number of grids represents the sample size, herein referred to as N. This resulted in us merging the eBird and BirdTrack datasets to ensure sufficient data of bird arrival. We took the first (i.e., earliest) date of the species sighting as the date of the phenological event. We should note that the date in which each phenological event is registered for the first time in a grid might differ within grids and within years, and, consequently, the meteorological conditions in which the event took place might be quite different to those experienced by the same species taking place in another grid.

Synchrony-asynchrony index

To evaluate phenological asynchrony/synchrony between different trophic levels we used the following interaction indices: (1)${\displaystyle A{L}_{\left(a,b\right)}=\left(D{A}_{ai}-F{F}_{bj}\right)}$ (2)${\displaystyle A{L}_{\left(a,c\right)}=\left(D{A}_{ai}-G{U}_{ck}\right)}$ (3)${\displaystyle A{L}_{\left(b,c\right)}=\left(F{F}_{bj}-G{U}_{ck}\right)}$

where a, b, and c refer to the trophic levels incorporated in the analysis, respectively, which are level 3 or secondary consumers (i.e., birds), level 2 or primary consumers (i.e., insects) and level 1 or producers (i.e., woodland vegetation). DA refers to the date of arrival of birds, FF refers to the date of first flight of insects, and GU refers to green up of vegetation. i, j, and k refer to the different species incorporated in the index for trophic levels a, b, and c, respectively.

When the phenological event for the higher trophic level takes place after the lower levels, this index results in a positive number. This infers that the higher trophic level is either exhibiting synchrony or delayed synchrony with the lower levels. Conversely, when it takes place before, this index results in a negative number. This infers that the higher trophic level is exhibiting asynchrony, or in other words birds are arriving before the insects take flight.

For ease of interpretation, AL values were standardised between a consistent maximum and minimum, with the maximum and minimum values for each group of species within a trophic level. Below is an example for trophic levels a and b. (4)${\displaystyle A{L}_{\left(a,b\right)}^{\prime }=\frac{\left(A{L}_{\left(a,b\right)}-min\left(A{L}_{\left(a,b\right)}\right)\right)}{\left(max\left(A{L}_{\left(a,b\right)}\right)-min\left(A{L}_{\left(a,b\right)}\right)\right)}}$

Values close to 0 reflect asynchrony, 0.5 reflect perfect synchrony, and 1 reflect delayed synchrony. Delayed synchrony means that species are undertaking their phenological events in the anticipated order, but the delay between them may still result in trophic cascades.

To assess interaction across the multiple trophic levels, we integrated Eqs. (1) and (3) to estimate asynchrony among three trophic network levels using a qualitative integrated (a)synchrony index (QAI) Eq. (5): (5)${\displaystyle QAI=A{L}_{\left(a,b\right)}^{\prime }|A{L}_{\left(b,c\right)}^{\prime }}$

We took the median value for AL’_(a,b) and AL’_(b,c) and substituted the numeric values by “A” (asynchrony) when AL’ was lower than 0.5 or “S” (synchrony) when AL’ was equal or greater than 0.5. We then put together those qualitative categories as shown in Eq. (5) with the aim to easily identify whether synchrony was preserved over trophic levels, by comparing bird-insect — insect-vegetation. We carried out a total of ninety combinations for the QAI through 3 trophic levels. Qualitative combinations of QAI were: AS, AA, SS or SA. For example, SS was reported when green up occurred before first flight, which subsequently occurred before date of arrival, while AS was reported where green-up occurred before first flight, but first flight occurred after date of arrival.

Interaction index of climate window movement

To evaluate changes in the sensitivity to meteorological variables and synchrony-asynchrony, we developed the Interaction Index of Climate Window Movement (IICWM). This novel approach to quantifying synchrony-asynchrony includes a meteorological dimension building on our AL index. For the purposes of this index, the use of climate and meteorology can be considered synonymous, as the approach could be extended from daily resolution weather data to decadal climate data. The meteorological dimension is the critical time window for which a particular meteorological event resulted in the ‘best’ statistical model of the phenological event calculated using a relative sliding time window (SWR) analysis from Climwin (Bailey & Van De Pol, 2016), in R (R Core Team, 2014). Climwin fits all potential time windows within a selected range for each phenological record and employs nested for-loops varying the start and end time of these windows. In our case we wanted to test both the effect of spring phenology, but also if our models could capture autumn signals. Therefore, we chose a temporal range of between 0 and 180 days prior to the event. The meteorological dataset is then subset, only taking values that match the tested window, providing a “window open” value and “window close” value. We ran the relative sliding time window analysis for each species and year, across the study period. For a full methodological account of the sliding time window analysis, please see Zenodo.

We then combined the synchrony-asynchrony index calculated through Eq. (1) with the variation in the date of the opening critical time window for TMIN, TMAX, and TPPT, for each bird-insect interaction. (6)${\displaystyle I\text{ICWM}=\left(\frac{A{L}_{\left(a,b\right)}}{WO}\right)}$ where WO refers to window open, calculated using Eq. (7). We used window open as this parameter represents the initial date in which the model captures a signal from the meteorological cue. (7)${\displaystyle WO=\left(WO\left[TMINa\right]-WO\left[TMINb\right]\right)+\left(WO\left[TMAXa\right]-WO\left[TMAXb\right]\right)+\left(WO\left[TPPTa\right]-WO\left[TPPTb\right]\right)}$ with the value a combination of the three environmental variables, TMIN, TMAX and TPPT run iteratively for trophic levels a and b. Using the non-standardised AL_(a,b) and WO values for each particular year and species combination we identified three possible categorisations (+,-,*) and two scenarios:

• Positive (+) IICWM represents a bird-insect combination that is synchronous, and situations where bird time windows take place further in the past and/or insect windows take place closer to first flight (scenario A).

• Asterisk (*) IICWM represents a bird-insect asynchronous combination together with bird time windows taking place closer to DA or insect time windows taking place further back before first flight (scenario B).

• Negative (-) IICWM will be returned in two cases: (1) when a particular bird-insect combination is asynchronous in combination with scenario A, or (2) when a synchronous combination occurred together with scenario B.

Phenological events from 2008–2018

There was high variability in the timing of first flight, date of arrival, and green-up between 2008–2018 for all groups of species (Fig. 1). Five out of nine moth species took their first flight between days 120–140 (Fig. 1A), while first flight among butterflies was more variable (Fig. 1B). Four species of butterfly initiated first flight before day 100, while six species showed similar median first flight values (Fig. 1B). In the case of birds, barn swallow, and northern wheatear showed an early spring arrival while greater whitethroat and sedge warbler showed arrivals later in the spring. Three species, willow warbler, sedge warbler, and greater whitethroat presented an early median date of arrival (Fig. 1C). Vegetation green-up dates showed the greatest variation, yet median dates were similar for four out of seven vegetation groups (Fig. 1D).

When phenological events were analysed per species (Zenodo), certain trends emerged, but the consistent result was varying phenology. For example, green-up dates showed great variation during the period 2008–2018, with four vegetation groups (BW, MBC, OAH, and WAA) advancing their phenology from 2008, two groups (MBW and OBH) delaying their phenology, while one group (WOA) remained constant. Five moth species (BM, DA, ET, HD, and WE) had a later first flight in 2018 than in 2008, while three (CMC, FS, and SS) showed an earlier first flight across the decade. A similar result was observed in butterflies, with five species presenting a later first flight in 2018. Arrival of migrant birds also fluctuated, but date of arrival generally advanced over the decade for all species, except for the greater whitethroat.

Synchrony-asynchrony indices of the trophic levels

The indices of synchrony-asynchrony between all birds and insects (AL_(3,2)), insects and woodland vegetation (AL_(2,1)), and birds and woodland vegetation (AL_(3,1)) were performed for 251 species combinations. Median values across 2008–2018 showed that 77 species combinations were asynchronous (AL’ < 0.5) and 174 were synchronous (AL’ > 0.5). The most asynchronous interactions included barn swallow, northern wheatear, sedge warbler, greater whitethroat, willow warbler, common marbled carpet, and early thorn.

A network of potential species interactions at the national level for all 30 species was developed using the median AL’ scores from 2008–2018 (Figs. 2–3). Species that were synchronous (Fig. 2) are potentially more likely to interact, but thinner lines represent delayed synchrony. Species that were asynchronous (Fig. 3) show those that are occurring in the ‘wrong’ order, with thicker lines showing a longer time discrepancy. Birds reported the greatest asynchrony (Fig. 3) with several delayed relationships (i.e., thick lines) reported. For example, barn swallow showed generalised asynchrony with AL’_(3,2) ranging between 0.41–0.12 with all but two insect species, common marbled carpet (0.7) and early thorn (0.63). The common marbled carpet showed asynchrony with all vegetation groups with AL’_(2,1) ranging from 0.47–0.43. Early thorn was asynchronous with three vegetation groups, bog woodland, mixed broadleaved/conifer woodland, and wet pedunculate oak-ash woodland, with AL’_(2,1) ranging between 0.489–0.49.

Interaction index of climate window movement

Meteorological drivers & time windows of phenological events

Our climwin models illustrated a marked interannual variation in the timing of phenological events. The opening and closing of relevant temporal windows were estimated to influence the three meteorological variables at group and species level (Zenodo). For example, in the case of woodland vegetation, all three meteorological variables showed similar estimated influences for all seven types of vegetation, and Pc indicated statistical significance across drivers and windows (Pc < 0.5), as shown in Table 2 (see Zenodo for a definition of Pc). Data for all time windows across all thirty species and their Pc estimates are shown in Zenodo. Model coefficients showed a negative relationship with increased TMAX, TMIN, and TPPT and green-up dates, indicating that an increase of 1 degree in temperature (either maximum or minimum) and one mm³ of TPPT over the identified time window is expected to advance green-up between 0.25 and 10 days.

Table 2:

Climwin results for OAH vegetation group showing yearly time windows in which green-up was most influenced by each one of our meteorological variables.

Window open reflects the number of days prior to GU in which a driver started to influence GU, and takes place further back in time, while window close (also shown in days prior GU) always takes place after window open, closer to the phenological event. GU represents the Julian day of the year in which green-up took place at national level (in ordinal date). N means sample size, number of grids in which GU was measured. Asterisks denote statistical significance showed by climwin Pc values.

Oak-ash-hazel woodland			TMIN		TMAX		TPPT
Year	N	GU	W Open	W Close	W Open	W Close	W Open	W Close
2008	129	64	168	147***	179	108***	89	56***
2009	129	38	165	118***	179	121***	82	42***
2010	129	57	175	141***	142	89***	147	96***
2011	129	42	175	140***	151	78***	60	1***
2012	129	41	179	144***	72	19***	90	57***
2013	129	54	15	0***	169	152***	89	49***
2014	129	37	177	148***	175	97***	56	6***
2015	129	37	46	5***	77	9***	61	33***
2016	129	35	178	144***	179	118***	76	12***
2017	129	32	180	159***	60.5	27***	29	1***
2018	129	45	37	7***	143	123***	92	14***

DOI: 10.7717/peerj.18653/table-2

Notes:

*P-value < 0.5.

**P-value < 0.1.

***P-value < 0.001.

In the case of the time windows of butterflies and moths, TMIN and TMAX were often the meteorological variables that showed higher influence on first flight, while TPPT time windows had less estimated influence, sometimes no better than random (see Zenodo). Coefficients for all three meteorological variables fit for butterfly and moth phenology showed high variation among species and years, with both significant positive and negative effects on first flight (Zenodo). Time windows that influence bird date of arrival and model significance was highly variable (Zenodo), with no clear trends in model coefficients for all time windows across species (Zenodo).

The period that the meteorological time windows were open fluctuated across species and groups between 2008–2018 (Figs. 4 and 5). Woodland vegetation time windows presented the longest time lapse of critical influence of each meteorological driver across all four groups, with the longest windows corresponding with TMIN, TPPT and TMAX, respectively (Figs. 5D–5F). First flight of butterflies and moths showed a similar pattern in the way each meteorological driver influenced their time windows (Figs. 4A–4F). TMIN had a greater influence around 0-70 days prior to first flight, with the median window open (∼50 days prior) and window close (∼25 days prior) being particularly low (Figs. 4A, 4D). Among migrant bird species, the length of time windows varied (Fig. 5), with TMIN median values for most species (Fig. 5A) ranging from 55 to 25 days prior to date of arrival. For a detailed description of the time window results broken down into species across the 10-year period, see Zenodo.