# Phenological shifts alter the seasonal structure of pollinator assemblages in Europe

François Duchenne, E. Thébault, D Michez, M Elias, M Drake, M Persson, J

S Rousseau-Piot, M. Pollet, P Vanormelingen, C Fontaine

► To cite this version:

François Duchenne, E. Thébault, D Michez, M Elias, M Drake, et al.. Phenological shifts alter the seasonal structure of pollinator assemblages in Europe. Nature Ecology & Evolution, 2019, 4 (1), pp.115-121. 10.1038/s41559-019-1062-4 . hal-02456569

HAL Id: hal-02456569

https://hal.science/hal-02456569

Submitted on 27 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Phenological shifts alter the seasonal structure of pollinator assemblages in

2 Europe

3

4 F. Duchenne*¹,², E. Thébault¹, D. Michez³, M. Elias⁴, M. Drake⁵, M. Persson⁶, J.S. Rousseau-Piot⁷, M.

5 Pollet⁸, P. Vanormelingen⁹ & C. Fontaine²

6 ¹Sorbonne Université, CNRS, Université Paris Est Créteil, INRA, IRD, Institute of Ecology and
7 Environmental Sciences of Paris, 75005 Paris, France, ²Centre d'Ecologie et des Sciences de la
8 Conservation, (CNRS, MNHN, Sorbonne Université), 75005 Paris, France, ³Laboratory of Zoology,
9 Research institute of Biosciences (University of Mons), 7000 Mons, Belgium, ⁴Institut de Systématique,
10 Evolution, Biodiversité (MNHN, CNRS, Sorbonne Université, EPHE, Université des Antilles), Paris
11 75005, France, ⁵Orchid House, Burridge, Axminster, Devon, UK, ⁶Department of Water Resources
12 Engineering (Lund University), 221 00 Lund, Sweden, ⁷Natagora, 1050 Bruxelles, Belgium, ⁸Research
13 Group Species Diversity (SPECDIV), Research Institute for Nature and Forest (INBO), Kliniekstraat
14 25, B-1070 Bruxelles, Belgium, ⁹Natuurpunt Studie, Coxiestraat 11, 2800 Mechelen, Belgium.

15

16 *: corresponding author, François Duchenne [francois.duchenne@mnhn.fr](mailto:francois.duchenne@mnhn.fr)
Pollinators play an important role in terrestrial ecosystems by providing key ecosystem functions and services to wild plants and crops, respectively. The sustainable provision of such ecosystem functions and services requires diverse pollinator communities over the season. Despite evidence that climate warming shifts pollinator phenology, a general assessment of these shifts and their consequences on pollinator assemblages is still lacking. By analyzing phenological shifts of over 2000 species, we show that on average the mean flight date of European pollinators shifted 6 days earlier over the last 60 years, while their flight period length decreased by 2 days. Our analysis further reveals that these shifts have likely altered the seasonal distribution of pollination function and services by decreasing overlap among pollinators' phenologies within European assemblages, except in the most northeastern part of Europe. Such changes are expected to decrease the functional redundancy and complementarity of pollinator assemblages and as such, might alter the performance of pollination function and services and their robustness to ongoing pollinator extinctions.

Numerous studies on plants, birds, amphibians and insects reveal that on average various phenological events - such as flowering or initiation of flight season – now take place earlier in the season than in the past decades because of climate warming¹. Despite this general trend, a substantial inter-specific variation is observed in these responses, spatially² (e.g. across latitudes), and temporally²,³ (e.g. spring *versus* summer species). This heterogeneity in species responses together with the fact that most studies focus on taxonomic rather than functional groups¹ challenges our ability to assess the consequences of phenological shifts for the functioning of communities and ecosystems across large spatial scales.

By modifying the set of species co-occurring in time, heterogeneity in phenological responses can induce mismatch among interacting species⁴, thereby affecting community structure and related functions. One key issue to our understanding of the impact of climate warming on ecological functions is thus to assess how phenological shifts combine themselves among the species assemblage involved in a given function. This requires to quantify the phenological responses of a large proportion of the species, not only in terms of mean flight date shifts but also of changes in phenology length, a
44 currently overlooked aspect of species responses⁵. The few studies that started to tackle this issue
45 revealed important changes in patterns of species temporal overlap in several local communities of
46 plants and amphibians, as a result of non-uniform phenological shifts⁴,⁶,⁷. However, these studies
47 remain restricted to a small set of functional or taxonomical groups and to a small set of local
48 communities.

49 Pollination is a key ecosystem function⁸,⁹ mainly performed by four insect orders in Europe: 50 Hymenoptera, Diptera, Lepidoptera and Coleoptera¹⁰. These flower visitors present a continuum of 51 pollination efficiency but the diversity within pollinator assemblage has been proved to increase 52 pollination performance¹¹. Current theoretical knowledge indicates that the level of heterogeneity in 53 phenological responses to climate warming among pollinators can strongly affect pollination 54 networks¹². However, the quantification of the phenological responses of pollinators to climate 55 warming is still limited, with studies focused on butterflies¹³,¹⁴ and to a lesser extent, on bees³ and 56 hoverflies¹⁵. A better understanding of the consequences of climate change on pollination thus requires 57 a much more complete assessment of changes in pollinator phenology, including more species and 58 changes in both timing and duration of the seasonal activities.

59 We took advantage of recent developments of large biodiversity databases and museum collections 60 and we compiled a database of over 19 million records of flower visitor occurrences (Supplementary 61 Table 1), spanning the period 1960-2016. This database includes 2023 European species from the 4 62 main insect orders of pollinators: Hymenoptera, Diptera, Lepidoptera and Coleoptera (Extended Data 63 Fig. 1). Numerous species exhibit distinct modes in their phenology, either because they are 64 multivoltine (i.e. multiple generations per year) or because the phenology differs between sexes or 65 social casts. Since different modes from a species can shift in a different direction, we studied each 66 mode separately, leading to 2248 phenology modes (see Methods). For each phenology mode, we 67 estimated changes in mean flight date and flight period length over the years by modeling the mean 68 and variance of collection dates (see Methods). Similarly to previous studies working with historical 69 records³, due to the lack of long-term standardized monitoring for many flower visitor taxa and at 70 large spatial scale, our analysis relies on opportunistic data. However such datasets have been shown
71 to give estimates of phenological shifts quantitatively consistent with those based on standardized
72 monitoring data¹⁵,¹⁶.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Mean flight date (MFD) shifts of European flower visitors between 1960 and 2016. (a) Phylogeny of studied species and MFD shifts (n=2248). The bars around the phylogeny tips are proportional to the MFD shifts and colored in blue and red for phenological advancement and delay, respectively. Values below -0.5 and above 0.5 days/year are truncated to preserve readability. Histograms show MFD shifts for all studied species of Coleoptera (b, red, n=194), Diptera (c, blue, n=305), Hymenoptera (d, light green, n=322) and Lepidoptera (e, magenta, n=1427). Full bars represent number of species with values significantly distinct from zero whereas open bars correspond to the number of species with a value non-significantly distinct from zero. MFD shifts shown here are predicted for the averaged latitude, longitude and altitude of each species' records.</figcaption></figure>

73

74

75

76

77

78

79

80

81

82

83
## 84 Results

We find that the mean flight date changes on average at a rate of -0.104±0.004 days/year (mean ± SE) implying that European pollinators are flying on average 5.8 days earlier in 2016 than in 1960, a value consistent with previous estimations on bees³ and butterflies¹³. Climate warming appears as a likely cause as mean flight date shift mainly occurred after 1980, following the temperature increase (Supplementary Method 1, Extended Data Fig. 2). Considering flight period length, we find that on average the standard deviation of collecting dates decreases slightly with time, at a rate of -0.016±0.003 days/year (mean ± SE), which corresponds to a decrease of 1.8 days of the flight period length over the last 56 years. This reduced flight period length might be due to a reduced genetic variability on phenology caused by a directional selection on phenology advancement. Indeed, we know that a directional selection on a phenotypic trait can reduce the variance of this trait¹⁷, and the significant positive Pearson correlation between the changes over time of mean flight date and flight period length (r = 0.09, t<sub>df=2246</sub>=3.89, p-value = 1e-4) can suggest such a mechanism. However, whether these changes are adaptive or not, and the mechanisms underlying these responses (adaptation vs phenotypic plasticity), remain unknown.

Despite these overall trends, we observe a substantial heterogeneity among species in the response of mean flight date and flight period length (Fig. 1, Extended Data Fig. 3). 13% of the phenologies studied exhibit a significantly delayed mean flight date whereas 30% do not exhibit any significant shift (Supplementary Table 2). Such heterogeneity is even more striking for flight period length changes, where 27% of the phenologies studied are significantly lengthened and 43% are unchanged (Supplementary Table 2). If an increase of winter temperature is known to advance species phenology by reducing the development time¹⁸, some species also react in an opposite way¹⁸,¹⁹, which might explain observed variations in mean flight date shifts. Turning to the heterogeneity in flight period length responses, a temperature increase can either reduce or increase flight period length, for example by reducing insect lifespan²⁰ or by increasing the number of generations within years¹⁴.

We further show that this heterogeneity in phenological responses is related to the evolutionary history of species as shown by the strong phylogenetic signal in mean flight date shifts (Pagel's λ =
0.75, p-value < 0.001) and in flight period length changes (Pagel's λ = 0.82, p-value < 0.001). This phylogenetic signal is related to strong differences among orders in these phenological shifts, Diptera

<figure><img src="image_3.png" /></figure>

and Coleoptera advancing their mean flight date more than Lepidoptera and Hymenoptera while Coleoptera decrease their flight period length more than other orders (Supplementary Table 3). However the phylogenetic signal remains significant within orders for mean flight date shifts (Supplementary Table 3). Such phylogenetic signal indicates that species traits underlying phenological responses are conserved across the phylogeny.

**Figure 2:** Spatial and seasonal heterogeneity in phenological shifts among species. Mean flight date (MFD) shifts (top panels) and changes in flight period length (FPL, bottom panels) against species MFD (a, d) and averaged latitude (b, e) and longitude (c, f) of species records. The horizontal grey lines show the 0 value whereas red lines are PGLS predictions. Estimates and standard errors are shown in Supplementary Table 4.

Mean flight date and flight period length responses also demonstrate spatial and seasonal heterogeneity among species. Species with southern and western distribution areas show a stronger mean flight date advancement than species with northern and eastern distribution areas (Fig. 2b-c, Supplementary Table 4), matching previous results on European plants². We also find that species with northern and western distribution areas experience a smaller decrease in flight period length than
species with southern and eastern distribution areas (Fig. 2e-f, Supplementary Table 4). In addition we find a seasonal pattern where spring species experience a significantly greater advancement than summer/autumn species (Fig. 2a), consistently with previous results on American bees³ and European plants². Regarding flight period length, we find that earlier species shorten their flight period more than later species (Fig. 2a-2d, Supplementary Table 4). Such differences could be explained by the fact that summer/autumn and northern species might rely more on photoperiod, a determining factor of insect phenology²¹, than spring and southern species. Such patterns have been shown for plants²²,²³, but studies on this point for insects are missing.

We further show that the mean flight date shifts vary within species in a way that echoes the patterns found at the inter-specific level. Indeed, we detect a significant positive interaction between latitude and year effect for 29% of species, indicating that southern populations experience a stronger shift of their mean flight date towards earlier dates than northern populations (Supplementary Table 2). By contrast, no longitudinal pattern was found. The seasonal pattern of stronger advancement earlier in the season is also found at the intra-specific level. Amongst the 190 species with multimodal phenology and sufficient data to study them, 59% have their first mode advancing significantly faster than their second mode while the opposite pattern occurs only in 10.5% of the species (Extended Data Fig. 4).
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Changes in within- and among-orders average overlaps in phenology between 1980 and 2016 across Europe. (a) Average phenology over all species in 1980 (solid lines) and in 2016 (dashed lines) for one grid cell (centroid = 55,0) by orders: Coleoptera (red), Diptera (blue), Hymenoptera (light green) and Lepidoptera (magenta). The average phenology is calculated by averaging all probability density functions (Gaussians representing phenologies) over all species of each orders, assuming identical species abundances. Observed changes in the average overlap among phenologies between 1980 and 2016, within orders (b) and among orders (c). Uncolored cells are underprospected. Number of species by order across Europe is shown in Extended Data Fig. 6.</figcaption></figure>
156

To assess the consequences of these phenological shifts for the seasonal structure of pollinator assemblages across space, we analyzed changes in the phenological overlap of species co-occurring within locations of 5°×5° grid cells in Europe, between 1980 and 2016. We used the linear models for mean flight date and flight period length to predict the phenologies of each species for each grid cell predicted for both years (Fig. 3a). Considering that all phenologies are unimodal, we modelled them by Gaussian density distribution, in order to calculate the pairwise phenological overlap among all pairs of pollinators present in a grid cell (see Methods). We averaged these measures among pollinators belonging either to the same or to different insect orders (see Methods).

157

158

159

160

161

162

163

164

First, we show that species co-occurrence in time increases towards the beginning of the season and then abruptly decreases in the second half of the season (Extended Data Fig. 5), consistently with the average advancement of pollinator mean flight date. This indicates that the advance of mean flight dates have likely shifted the pollination function and services earlier in the season. Secondly, assuming no changes in abundance/distribution of species, we show that both within- and among-orders average overlaps in phenology have decreased within the last 36 years in most parts of Europe, except in the extreme northern part (Fig. 3b-c. The observed increase of the overlap among phenologies in the northern Europe is likely due to the fact that there, in contrast to other regions, the average mean flight date shift is almost null whereas the flight period length slightly increases (Fig. 2). Sufficient data on long-term dynamics of plant-pollinator networks is currently missing to fully assess the consequences of such changes in the seasonal structure of pollinator assemblages on pollination function. However, the within order and among order overlaps should be related to temporal redundancy and complementarity within pollinator assemblages, respectively. Indeed, the pervasive phylogenetic signal within pollination networks indicates that related pollinators tend to visit the same plants²⁴,²⁵. This implies that species with overlapping phenologies and belonging to the same insect order should visit the same set of co-flowering plant species and as such, belong to the same pollinator functional group. On the contrary, species with overlapping phenologies but from different insect orders are expected to provide complementary pollination function, by visiting different sets of co-flowering plant species.

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182
183

As such, the observed decrease in the overlap within insect orders, by lowering the temporal
redundancy among pollinators, might decrease the robustness of plant-pollinator interaction networks
to pollinator extinction²⁶. A decrease in the overlap may also have beneficial effects for pollinators by
decreasing competition for nectar and pollen resources, but such competition release might in turn
restrict pollinator visits to the most profitable plant species following optimal foraging theory
predictions²⁷. Turning to the observed decrease in phenology overlap among pollinator orders, it
suggests a decrease in temporal complementarity within pollinator assemblages, thereby weakening
the pollination function delivered to plant communities²⁸. This result echoes theoretical findings on
pollination networks showing that the more phenologies are scattered over the season, the more
community diversity decreases²⁹.

184

185

186

187

188

189

190

191

192

193

Discussion

194

195

196

197

198

199

200

201

202

203

204

205

206

Our results show that flower visitor responses to climate warming depend on their evolutionary
history, geographical location and seasonal earliness. This high variation in species phenological
responses is expected to drive heterogeneity in the consequences of climate warming on pollination
function across Europe and across the season. For most parts of Europe, the observed modifications of
the seasonal structure of pollinator assemblages are expected to have negative consequences on
pollination, while in northeastern Europe they might have positive effects on pollination as they result
in an increased phenology overlap, both within and among pollinator orders (Fig. 3). Moreover, in
most parts of Europe, observed changes are expected to have a positive effect on pollination
performance and robustness early in the season but a negative effect from the middle to the end of the
pollination season (Extended Data Fig. 5). Thus, our results highlight the importance to assess
responses at large spatial and temporal scales and to include many species, in order to capture the high
spatial and seasonal heterogeneity in the consequences of climate change on pollinator assemblages
and related function.

207

208

209

Climate warming is recognized as a major threat to biodiversity. Our results suggest that climate
warming, by reducing pollinator co-occurrence in time within seasons has had negative effect on the
delivery of pollination function as well as on its resistance to further perturbations, in most part of
Europe. Such findings raise the question of potential interactive effects between climate warming and other pressures related to global change such as agricultural intensification³⁰,³¹, which could amplify expected negative effects on pollination. In addition to its effect on species phenology, climate warming is expected to affect the spatial distribution³² and the abundance³¹ of flower visitors, and so are other drivers of global change. How such effects combine with those observed in this study remain currently unknown. This stresses the need to explore multiple responses of species to multiple drivers of global change in order to assess potential synergistic effects among species responses to global change drivers over large scale.
218 **Material and Methods**

219 **Constructing the database on flower visitor phenologies**

220 **Assembling data on flower visitor occurrences in time and space**

221 European flower visitors mainly belong to four insect orders – Coleoptera, Diptera, Hymenoptera and Lepidoptera¹⁰. We first looked for occurrence data (i.e. sighting at a given date and location) of 222 species that belong to these insect orders and that are defined as floricolous in scientific or grey 223 literature. We restricted our search to European species listed in Fauna Europaea³³. Data are from 15 224 distinctive sources, summarized in the Supplementary Table 1, with a high proportion from the Global 225 Biodiversity Information Facility (GBIF). After the removal of duplicates (same species, date and 226 locality), the database initially included about 30 million of occurrences between 34° and 72° of 227 latitude North and between -15° and 32° of longitude. 228

229 **Modelling multimodal phenologies and removing larval records**

230 Numerous species exhibit distinct modes in their phenology, either because they are multivoltine 231 (i.e. multiple generations per year) or because the phenology differs between sexes or social casts. 232 Since different modes in the same species are temporally distant, they might not respond to the same 233 environmental cues. As a consequence, each mode might potentially shift in a different direction and 234 should thus be studied separately. Additionally, larvae might be easier to spot than adults for 235 Lepidoptera and some Coleoptera. So, a substantial proportion of records may actually be larvae, 236 which are not floricolous and should be removed from the analysis. To split the occurrences of 237 multimodal imago phenology into distinct modes as well as to identify larval occurrences, we 238 developed the following method. 239

239 The first step of the method was to detect multimodality. Since phenologies vary spatially, 240 multimodality can be the product of sampling in different localities. In order to take this spatial 241 variation into account, for each species separately, we fitted the following linear mixed-effects model 242 accounting for spatial variables on the Julian day of records. 243

$$
243 \qquad Y_{ik} = \mu + \rho_1 \times \text{latitude}_k + \tau \times \text{longitude}_k + \theta \times \text{altitude}_k + \varphi_i + E_{ik} \quad (1)
$$
Where $Y_{ik}$ is the Julian day of the observation $k$ of the year $i$, $\mu$ is the grand mean (intercept), $\rho_1$ and $\tau_1$ are latitude and longitude effects respectively, while $\theta$ is an altitude effect. $\phi_i$ is a random year effect (factor) and $E_{ik}$ is the error term (independent and identically distributed, following $N(0, \sigma^2)$). The residuals of this model thus represent the collection dates once spatial and altitudinal variations have been removed. To detect multimodality in the distribution of these residuals, we smoothed the distribution with the R function *density*, using the value 1.3 for the *adjust* parameter and counted the number of local maximums (*nbmax*) which reaches 7% of the highest mode. We used this cut-off in order to remove small peaks on the edges of the phenology and we defined the value of the threshold after a visual inspection of phenologies. Several modes were initially detected for 494 species. For each of these species, we checked in scientific and grey literature whether a multimodal phenology was expected. In 208 cases, there was no strong biological support of existing multimodal phenology and we thus considered these species had one single mode. After this step, 288 remaining species showed a multimodal phenology ($nbmax > 1$). We applied the second step only for these species.

The second step of the method was to attribute each record to a specific mode. To do so we used clustering Gaussian mixture-models implemented in the *mclust* R package³⁴, considering a number of gaussians in 1 to nbmax. This clustering model allow us to initialize the attribution of each record to a given mode. Using the classification given by these clustering models, we run linear mixed-effects models, similar to the one described in equation (1) but with the addition of a mode effect ($\beta_j$):

$$Y_{ijk} = \mu + \rho_1 \times \text{latitude}_k + \tau \times \text{longitude}_k + \theta \times \text{altitude}_k + \beta_j + \phi_i + E_{ijk} \quad (2)$$

We kept the number of modes that minimize the Bayesian Information Criterion (BIC) of this linear mixed-effects model. We then manually changed the mode of poorly predicted points. If the change improved the likelihood of this mixed-effects model we retained it and continued this process iteratively. We stopped the process when changing the mode of poorly predicted points did not further improve the likelihood of the model. The R script of the full method is available here: https://github.com/f-duchenne/Flower-visitors-phenology. Although the mode effect ($\beta_j$) is independent from spatial variables and altitude in equation (2), our method still allows to take into account spatial and altitudinal variation in the number of modes (Extended Data Fig. 7). We
confronted the relevance of detected modes regarding what we know on the biology of species. We
found that our method distribute records among modes in a highly consistent way. Some examples can
be seen in Extended Data Fig. 7. We identified 19 species for which we had a mode corresponding to
larval phenology, and we removed the corresponding records. Overall, this analysis lead to 2473
unimodal phenologies from 2179 species.

276 *Database after selection process*

Following the separation of distinct phenological modes for each species and the removal of larval
records, we selected phenologies (or phenological modes) with at least 400 records during the period
1960-2016 and with at least 40 records from the period 1960-1980, to be able to study phenological
shifts between early and more recent period. We removed species (n=30) with phenology peaking
during winter by excluding species with a mean flight date before 60 or after 306 Julian days.
Studying the phenology of such species raises methodological questions that we will not address here.
We also removed records with imprecise localization (above 1km²) except for small countries
(Luxembourg, Belgium, Switzerland, Netherlands, Denmark, Lichtenstein, Monaco, Andorra and
Kosovo). Thus, our dataset includes some records with imprecise localization (above 1km²) but they
represent less than 0.1% of the final dataset. This selection process lead to 19 845 792 occurrence
records with 2248 phenologies for 2023 species (Supplementary Table 1). The repartition of records
among insect orders and throughout the study period is presented in Extended Data Fig. 1.
Supplementary Table 1 indicates the amount of data coming from the various data sources.

290 **Analyses of species phenological shifts over time**

291 *Estimating species phenological shifts*

To estimate changes in both the mean flight date (MFD) and the flight period length (FPL), we
modeled jointly the mean and the variance of collection dates using the *dispmod* R package³⁵ which
performs two nested linear models, one for the mean and one for the variance. Due to computational
limits it was not possible to use one model including the whole dataset, modelling both, MFD shifts
and FPL changes, and modelling spatial effects properly for each species. Thus, we studied each
species and phenology mode separately. For each phenology modes, the model for the mean collection date was:

$$ 
\begin{aligned}
Y_k = & \mu + (\pi + \alpha \times \text{latitude}_k + \delta \times \text{longitude}_k) \times \text{year}_k + (\rho_1 + \gamma_1 \times \text{longitude}_k) \times \text{latitude}_k + \\
      & (\rho_2 + \gamma_2 \times \text{longitude}_k^2) \times \text{latitude}_k^2 + (\rho_3 + \gamma_3 \times \text{longitude}_k^3) \times \text{latitude}_k^3 + \tau_1 \times \\
      & \text{longitude}_k + \tau_2 \times \text{longitude}_k^2 + \tau_3 \times \text{longitude}_k^3 + \theta \times \text{altitude}_k + E_k
\end{aligned} 
\quad (3) $$

$Y_k$ is the Julian day of the observation $k$, $\mu$ is the grand mean (intercept), $\pi$ is the time effect on the mean collection date as well as on its variation across latitude ($\alpha$) and longitude ($\delta$). $\rho_1, \rho_2, \rho_3$ and $\tau_1, \tau_2, \tau_3$ are linear, quadratic and cubic effects for latitude and longitude, respectively, $\gamma_1, \gamma_2, \gamma_3$ are spatial interactions terms, $\theta$ is an altitude effect and $E_k$ is the error terms (independent and identically distributed, following $N(0, \sigma^2)$).

The joint model for variance of collection date was:

$$ \log(\sigma^2) = \mu_v + \rho_v \times \text{latitude}_k + \tau_v \times \text{longitude}_k + \theta_v \times \text{altitude}_k + \pi_v \times \text{year}_k \quad (4) $$

Where $\sigma^2$ is the variance of the collection date, $\mu_v$ is a constant term, $\rho_v, \tau_v, \theta_v$ and $\pi_v$ are latitude, longitude, altitude and year effects respectively. We performed model simplification based on the Akaike Information Criterion (AIC), first on the model for the mean collection date, removing only polynomials effect of latitude and longitude ($\gamma_1, \gamma_2, \gamma_3, \rho_2, \rho_3, \tau_2, \tau_3$) and interactions between spatial variables and time effect ($\alpha$ and $\delta$), and second on the model for the variance in collection date.

The MFD shifts presented in the paper are $\pi + \alpha \times \text{latitude} + \delta \times \text{longitude}$ from equation (3) where 'latitude' and 'longitude' are averaged latitude and longitude of the species records respectively. The FPL changes are the $\pi_v$ from equation (4) for each species.

## Phylogenetic analysis

In order to get a phylogeny of all the studied species we combined several published phylogenies. We used the phylogeny from Rainford et al.³⁶ as the backbone to which we added some available and recent phylogenies to get a phylogeny at the genus level for Papilionoideae³⁷, Vespidae³⁸ and Apoidea³⁹. For all other families, genus (as defined by the GBIF taxonomy) were inserted on a
polytomy positioned midway between the family origin and the tip. Then species from each genus were placed on a polytomy positioned midway between the genus origin and the tip. Such method does not allow a good estimation of the recent evolutionary history but because there is no phylogeny of insects at the species or genus level, it is the only way to include all species responses and take in account for intra-family heterogeneity. Moreover, because these polytomies were not too old relative to the entire phylogeny, it should not affect strongly our results. Because they are not present in our phylogeny three families of Diptera (Heleomyzidae, Limoniidae and Pediciidae) and two Lepidoptera species (Sphrageidus similis, Lymantriidae, and Agria desoptilete, Lycaenidae) were excluded from phylogenetic analysis.

We estimated phylogenetic signal in phenological shifts using Pagel's $\lambda^{40}$ implemented in the phylosignal R package$^{41}$, because it is much more robust to polytomies than Blomberg's K$^{42}$.

### Links between phenological traits and phenological shifts

To test if the seasonal precocity and the spatial distribution of species were linked to phenological shifts, we used the following phylogenetic generalized least squares model (PGLS) implemented in the caper R package$^{43}$ controlling for the Pagel's lambda of the residuals at the maximum likelihood:

$$PS_z = \mu + (\alpha) \times MFD_z + (\beta) \times \text{latitude}_z + (\delta) \times \text{longitude}_z + E_z \quad (5)$$

Where $PS_z$ is the phenological shift (i.e. MFD shift or FPL change) of the species z, $\mu$ is the grand mean (intercept), $\alpha$ is the effect of the mean flight date calculated with recent records (from 2000), $\beta$ is the effect of the average latitude of records, $\rho$ is the effect of the average longitude of records and $E_k$ is an error term following $N(0,\sigma^2)$.

### Analyses of the seasonal structure of pollinator assemblages

#### Predicting species phenology in different locations and years

To assess the effect of phenological shifts at the scale of the full pollinator assemblages, we calculated changes in the overlap among phenologies. Because phenological shifts depend on location, we discretized the studied area in cells of 5°×5°. This size was chosen in order to smooth the differences in sampling effort among localities. To ensure a representative pollinator assemblage, we
348

only included grid cells with at least 3 insect orders with 20 species with at least 30 records each. The
remaining cells were considered as under-prospected. Thus species are considered present in a grid-
cell if it has at least 30 records between 1960 and 2016. By doing so, we assume that the composition
of species assemblages are the same in 1980 and in 2016, which allows to study only the effect of
phenological shifts on seasonal structure. We considered that all species have the same abundance, and
a circular Gaussian phenology. We used wrapped circular normal distributions instead of a classical
Gaussian distribution in order to take phenologies that span winter into account. We estimated the
mean and the standard deviation of these Gaussians for the years 1980 and 2016 and for each grid cell,
using the predictions of the linear models used to estimate phenological shifts, described in equations
(3) and (4).

349

350

351

352

353

354

355

356

357

358

Calculation of phenological overlaps within assemblages

359

For each sufficiently prospected grid cell we calculated pairwise overlap among pollinator
phenologies present in the given grid cell. We considered that all species have the same abundance,
and a circular Gaussian phenology. The overlap between two phenologies is the integral of the
minimum of both gaussians. We calculated two overlap measures for each grid cell: the first one
focusing on the overlap within insect orders and the other one among insect orders. To give equal
weight to each insect order, and thus avoid over-representation of Lepidoptera, we first calculated the
mean overlap by insect order, or by pair of insect orders, respectively for the overlap within and
among orders. Second, we averaged these mean values per grid cell. Finally, in order to have more
robust values, we repeated this overlap calculation after shifting segmentation of the latitude and the
longitude by 1.25°, 2.5° and by 3.75°. Then we averaged values obtained by 1.25°×1.25° grid cells for
both measures, overlap within and among orders.

360

361

362

363

364

365

366

367

368

369

370

In order to study the seasonal dynamic of overlap changes, we calculated a proxy of the
phenological overlaps day by day in 1980 and in 2016 for each grid cell (Extended Data Fig. 5). We
do not use exactly the same calculation of overlap as previously for computational reasons. To
simplify the calculation method, we aggregated predicted phenologies at the order level to get a
density distribution by order, henceforth called order phenologies, as presented in Fig. 3a. Then we

371

372

373

374
calculated the pairwise overlap among order phenologies day by day for both years, 1980 and 2016, 375 and for each grid cell. We also evaluated the day-by-day density value for each order phenologies for 376 both years, 1980 and 2016, and for each grid cell. This density value is a proxy of the phenological 377 overlap within order, because we assume that every species has the same constant abundance. Then we 378 calculate the daily changes of these both indexes between 1980 and 2016 (Extended Data Fig. 5). We 379 did so for one grid pattern only (i.e. without sliding windows). 380

381 **Authors contribution**

F.D., C.F. and E.T. conceived the project. F.D. assembled dataset and performed the statistical 382 analysis. F.D., C.F., E.T., D.M. and M.E. interpreted the results and M.D., M.Pe., J.S.P., M.Po. & P.V. 383 provided data, biological expertise on studied species. F.D. wrote the paper with contributions from all 384 authors. 385

386 **Data availability**

The final dataset analyzed in this paper is available here: https://zenodo.org/record/3480120. 387

388 **Code availability**

Codes used to extract data from the GBIF, to separate modes of multimodal phenologies and to 389 estimate phenological shifts are available here: https://github.com/f-duchenne/Flower-visitors- 390 phenology. 391

392 **Competing interests**

The authors declare no competing interests. 393

394 **Acknowledgements**

We especially thank Nicolas Loeuille and Tom Van Dooren for their comments on this work, 395 Emmanuelle Porcher, Adrien Perrard, Elsa Teulière, Emile Faure, Elise Kerdoncuff, Benoit Perez and 396 Théophile Olivier for fruitful discussions. We also thank all naturalists who provided data to complete 397 our database, especially Lionel Casset, Fred Chevaillot, Laurent Guilbaud (INRA Avignon), Jean- 398 Laurent Hentz (Nature du Gard) and Guillaume Lemoine, as well as Christophe Daugeron, Rodolphe 399 Rougerie, Claire Villemant, Emmanuel Delfosse, Jérôme Barbut and Olivier Montreuil who helped us 400
401 to access insect collections of the French National Natural History Museum. This project was funded
402 by the Ministère de la Transition Ecologique et Solidaire as part of the project "What is the sensitivity
403 of pollinators to global warming in France" (convention n° SJ 3-17) lead by Christophe Daugeron and
404 Colin Fontaine, and by the Institut de la Transition Ecologique, Sorbonne Université, as part of the
405 project Yapludsaison.

406 .

407 References

408 1. Parmesan, C. Influences of species, latitudes and methodologies on estimates of phenological
response to global warming. *Glob. Change Biol.* **13**, 1860–1872 (2007).

410 2. Menzel, A. et al. European phenological response to climate change matches the warming pattern.
*Glob. Change Biol.* **12**, 1969–1976 (2006).

412 3. Bartomeus, I. et al. Climate-associated phenological advances in bee pollinators and bee-pollinated plants. *Proc. Natl. Acad. Sci.* **108**, 20645–20649 (2011).

414 4. Theobald, E. J., Breckheimer, I. & HilleRisLambers, J. Climate drives phenological reassembly of
a mountain wildflower meadow community. *Ecology* **98**, 2799–2812 (2017).

416 5. Parmesan, C. Ecological and Evolutionary Responses to Recent Climate Change. *Annu. Rev. Ecol.*
*Evol. Syst.* **37**, 637–669 (2006).

418 6. Diez, J. M. et al. Forecasting phenology: from species variability to community patterns. *Ecol.*
*Lett.* **15**, 545–553 (2012).

420 7. Carter, S. K., Saenz, D. & Rudolf, V. H. W. Shifts in phenological distributions reshape
interaction potential in natural communities. *Ecol. Lett.* **21**, 1143–1151 (2018).

422 8. Ashman, T.-L. et al. Pollen Limitation of Plant Reproduction: Ecological and Evolutionary
Causes and Consequences. *Ecology* **85**, 2408–2421 (2004).

424 9. Klein, A.-M. et al. Importance of pollinators in changing landscapes for world crops. *Proc. R.*
*Soc. Lond. B Biol. Sci.* **274**, 303–313 (2007).

426 10. Willemstein, S. C. *An Evolutionary Basis for Pollination Ecology*. (Brill Archive, 1987).
427

11. Garibaldi, L. A. et al. Wild Pollinators Enhance Fruit Set of Crops Regardless of Honey Bee Abundance. *Science* **339**, 1608–1611 (2013).

428

429

12. Memmott, J., Craze, P. G., Waser, N. M. & Price, M. V. Global warming and the disruption of plant–pollinator interactions. *Ecol. Lett.* **10**, 710–717 (2007).

430

431

13. Roy, D. B. & Sparks, T. H. Phenology of British butterflies and climate change. *Glob. Change Biol.* **6**, 407–416 (2000).

432

433

14. Altermatt, F. Tell me what you eat and I’ll tell you when you fly: diet can predict phenological changes in response to climate change. *Ecol. Lett.* **13**, 1475–1484 (2010).

434

435

15. Hassall, C., Owen, J. & Gilbert, F. Phenological shifts in hoverflies (Diptera: Syrphidae): linking measurement and mechanism. *Ecography* **40**, 853–863 (2017).

436

437

16. Robbirt, K. M., Davy, A. J., Hutchings, M. J. & Roberts, D. L. Validation of biological collections as a source of phenological data for use in climate change studies: a case study with the orchid Ophrys sphegodes. *J. Ecol.* 235–241

438

439

440

17. Walsh, B. & Blows, M. W. Abundant Genetic Variation + Strong Selection = Multivariate Genetic Constraints: A Geometric View of Adaptation. *Annu. Rev. Ecol. Evol. Syst.* **40**, 41–59 (2009).

441

442

443

18. Fründ, J., Zieger, S. L. & Tscharntke, T. Response diversity of wild bees to overwintering temperatures. *Oecologia* **173**, 1639–1648 (2013).

444

445

19. Stefanescu, C., Peñuelas, J. & Filella, I. Effects of climatic change on the phenology of butterflies in the northwest Mediterranean Basin. *Glob. Change Biol.* **9**, 1494–1506 (2003).

446

447

20. Miquel, J., Lundgren, P. R., Bensch, K. G. & Atlan, H. Effects of temperature on the life span, vitality and fine structure of Drosophila melanogaster. *Mech. Ageing Dev.* **5**, 347–370 (1976).

448

449

21. Bale, J. S. & Hayward, S. *Insect overwintering in a climate change*. **213**, (2010).

450

22. Way, D. A. & Montgomery, R. A. Photoperiod constraints on tree phenology, performance and migration in a warming world. *Plant Cell Environ.* **38**, 1725–1736 (2015).

451

452

23. Saikkonen, K. et al. Climate change-driven species’ range shifts filtered by photoperiodism. *Nat. Clim. Change* **2**, 239–242 (2012).

453
454

24. Rezende, E. L., Lavabre, J. E., Guimarães, P. R., Jordano, P. & Bascompte, J. Non-random coextinctions in phylogenetically structured mutualistic networks. *Nature* **448**, 925–928 (2007).

455

456

25. Ibanez, S., Arène, F. & Lavergne, S. How phylogeny shapes the taxonomic and functional structure of plant-insect networks. *Oecologia* **180**, 989–1000 (2016).

457

458

26. Memmott, J., Waser, N. M. & Price, M. V. Tolerance of pollination networks to species extinctions. *Proc. R. Soc. Lond. B Biol. Sci.* **271**, 2605–2611 (2004).

459

460

27. Fontaine, C., Collin, C. L. & Dajoz, I. Generalist foraging of pollinators: diet expansion at high density. *J. Ecol.* **96**, 1002–1010 (2008).

461

462

28. Fontaine, C., Dajoz, I., Meriguet, J. & Loreau, M. Functional Diversity of Plant-Pollinator Interaction Webs Enhances the Persistence of Plant Communities. *PLOS Biol.* **4**, e1 (2005).

463

464

29. Encinas-Viso, F., Revilla, T. A. & Etienne, R. S. Phenology drives mutualistic network structure and diversity: Effects on the dynamics of mutualistic networks. *Ecol. Lett.* **15**, 198–208 (2012).

465

466

30. Potts, S. G. et al. Global pollinator declines: trends, impacts and drivers. *Trends Ecol. Evol.* **25**, 345–353 (2010).

467

468

31. González-Varo, J. P. et al. Combined effects of global change pressures on animal-mediated pollination. *Trends Ecol. Evol.* **28**, 524–530 (2013).

469

470

32. Devictor, V. et al. Differences in the climatic debts of birds and butterflies at a continental scale. *Nat. Clim. Change* **2**, 121–124 (2012).

471

472

33. Jong, Y. de et al. Fauna Europaea – all European animal species on the web. *Biodivers. Data J.* **2**, e4034 (2014).

473

474

34. Scrucca, L., Fop, M., Murphy, T. B. & Raftery, A. E. mclust 5: clustering, classification and density estimation using Gaussian finite mixture models. *R J.* **8**, 205–233 (2016).

475

476

35. Scrucca, L. *dispmod: Modelling Dispersion in GLM*. R package version 1.2. (2018).

477

36. Rainford, J. L., Hofreiter, M., Nicholson, D. B. & Mayhew, P. J. Phylogenetic Distribution of Extant Richness Suggests Metamorphosis Is a Key Innovation Driving Diversification in Insects. *PLOS ONE* **9**, e109085 (2014).

478

479
480

37. Chazot, N. et al. The Trials and Tribulations of Priors and Posteriors in Bayesian Timing of Divergence Analyses: the Age of Butterflies Revisited. *bioRxiv* 259184 (2018). doi:10.1101/259184

481

482

483

38. Bank, S. et al. Transcriptome and target DNA enrichment sequence data provide new insights into the phylogeny of vespid wasps (Hymenoptera: Aculeata: Vespidae). *Mol. Phylogenet. Evol.* **116**, 213–226 (2017).

484

485

486

39. Hedtke, S. M., Patiny, S. & Danforth, B. N. The bee tree of life: a supermatrix approach to apoid phylogeny and biogeography. *BMC Evol. Biol.* **13**, 138 (2013).

487

488

40. Pagel, M. Inferring the historical patterns of biological evolution. *Nature* **401**, 877–884 (1999).

489

41. Keck, F., Rimet, F., Bouchez, A. & Franc, A. phylosignal: an R package to measure, test, and explore the phylogenetic signal. *Ecol. Evol.* **6**, 2774–2780 (2016).

490

491

42. Molina-Venegas, R. & Rodríguez, M. Á. Revisiting phylogenetic signal; strong or negligible impacts of polytomies and branch length information? *BMC Evol. Biol.* **17**, 53 (2017).

492

493

43. Orme, D. et al. *caper: Comparative Analyses of Phylogenetics and Evolution in R.* (2013).

494