Phenology and Robustness in plant-pollinator networks

Many metrics that describe the structure of mutualistic plant-pollinator networks have been found to be important for network stability and robustness. These metrics are impacted by a suite of variables, including species traits, species abundances, their spatial configuration, and their phylogenetic history. Here, we consider a specific trait, phenology, or the timing of life history events. We expect that timing and duration of activity of pollinators, or of flowering in plants, could greatly affect the structure of the networks in which they are embedded. Using plant-pollinator networks from 33 sites in southern British Columbia, Canada, we asked a) how phenological species traits, specifically timing of first appearance in the network and duration of activity in a network, were related to network structure, and b) how those traits affected network robustness to phenologically biased species loss. We found that long duration of activity increased connection within modules for both pollinators and plants and among modules for plants. We also found that date of first appearance was positively related to interaction strength asymmetry in plants but negatively related in pollinators. Networks were generally more robust to the loss of pollinators than plants, but robustness declined with loss of early-flying or long-duration pollinators. These pollinators tended to be among-module connectors. Our results show that changes in phenology have the potential to impact plant-pollinator networks, which may have conservation relevance in a time of changing climate.


Introduction
Species within communities form interaction networks, and many metrics that describe the structure of mutualistic plant-pollinator networks (e.g. interaction asymmetry) have been found to be important for the ability of networks to be resilient to perturbations . Many species attributes contribute to network structure, including traits like Tlower colour and feeding preferences (Vázquez et al. 2009, Valdovinos 2019, abundance (Vázquez et al. 2007, Valdovinos 2019, spatial in the interactions (Kaiser-Bunbury et al. 2010, Mello et al. 2011. Given the relationship between species phenological traits and network properties, and the relationship between network properties and network stability, we would expect that phenological traits themselves would be the mechanism that underlies the relationship between network structure and network stability. Indeed, Encinas-Viso et al. (2012) and Ramos-Jiliberto et al. (2018) studied the effect of phenological traits on the stability and robustness of networks.
Both of these studies used dynamical models to study this relationship. Encinas-Viso et al. (2012) found that as the length of the season in a network increases, diversity and resilience of the network also increases. Ramos-Jiliberto et al. (2018) went a step further and combined these dynamical models with empirical networks, and found that the loss of plants with earlier blooming dates and with longer active periods decreased pollinator persistence.
Here we use 33 mutualistic plant-pollinator interaction networks from Western Canada to ask how plant and pollinator phenology contribute to their network interaction structure. We focus on exploring four measures of network structure that are related to robustness: specialization, within-module degree, among-module connectivity, and interaction asymmetry. SpeciTically, we ask the following two questions: 1) How do date of Tirst appearance in a network, and length of activity during the season, affect individual species interaction patterns? We predict that species whose date of Tirst appearance is earlier, and that are active longer in the season, should be less specialized, have greater within-module degree, greater among-module connectivity, and have higher values of interaction asymmetry (they affect their partners more than the reverse). 2) How robust are networks to removal of species due to varying phenological "traits" (date of Tirst appearance early/late, duration of activity small/large)? Networks should be more robust to losing species whose date of Tirst appearance is later in the season, and species that are active during less of the season.

Study sites
A total of 33 mutualistic plant-pollinator networks were studied in British Columbia, Canada: oak savannah (12 networks), shrub-steppe (eight networks), and hedgerow restorations (13 networks). These three vegetation types comprised three different studies.
See Table A1 for site information, including latitude/longitude coordinates. For simplicity we use "pollinator" throughout this paper to refer to insects and hummingbirds observed visiting Tlowers and contacting reproductive organs, although their effectiveness at transfer of pollen has not been assessed.

Collection of mutualistic network data
Data were collected for two of three vegetation types using the plot method, and for the third using the transect method. Plots are generally more appropriate when the plant species in the community are very patchily distributed (Gibson et al. 2011), as they were in these regions. The plot method focuses on individual plant species, observing each plant species for an equal amount of time. For oak savannah sites (within the Coastal Douglas Fir biogeoclimatic zone), we collected data on species interactions in 1-ha plots at each of six sites in both 2009 and 2010. Each plot was surveyed about every 7-10 days, 10-12 times per season between late April and early July, the majority of the Tlowering period. Over the Tlowering period we attempted to visit sites morning, midday, and afternoon on different survey dates to reduce bias due to Tlight time differences among visiting insects. During each survey date, each plant species in Tlower was observed for a 10 min period by each of two surveyors, on haphazard walks throughout the plot. All Tlower visitors were collected and identiTied to the lowest taxonomic level possible. For the eight shrub-steppe sites (in the Bunchgrass biogeoclimatic zone), data were collected as for oak-savannah sites, but surveys were from the beginning of April through the end of July, 2010, for a total of 12 samples per site.
For hedgerow restoration sites, data were collected using a "transect" method, in which the plants along the transect were observed for a set amount of time, with time observed per plant species varying among species depending on their occurrence in the transect. Transects are more appropriate when plants are not clumped, but are widely scattered throughout a study site (Gibson et al. 2011), and in this case most of the restorations (within the coastal Western Hemlock biogeoclimatic zone) were linear, making transects efTicient. Sampling was equal across all sites, occurring approximately every 2 weeks, for a total of 9 samples between late April and the end of August, 2013. The transect was walked for 15 minutes by each of two observers during each sample date, and each site was again observed equally during morning, midday, and afternoon on different sample dates. All pollinators were collected for identiTication in the lab.

Species phenological variables
We collected the following phenological species variables for every pollinator and plant species: a) Tirst Julian day observed interacting in the network, and b) number of days observed in the network (last date observed -Tirst date observed). We treat each network as a replicate, and capture both variables for all species within the network. This means that there is no single value for Tirst day observed or total days observed for any particular species across networks. Phenological variables were calculated from the observation of the interactions. While we acknowledge this is an imperfect way to detect phenological variables, it is consistent across all 33 networks and has been previously used in the literature (Rasmussen et al. 2013). Because these two phenological variables could potentially be correlated (they are calculated from the same set of data), we calculated Pearson's correlation coefTicient for log10-transformed variables for each of plants and pollinators separately. We found that the two variables were weakly negatively correlated (pollinators: ρ = -0.5, P < 0.001, df = 1690; plants: ρ = -0.25, P < 0.001, df = 590), meaning that if a species has a very late First Julian day it cannot, by deTinition, be present many days in the network. In contrast if a species is has an early First Julian day, it can be present many days, or few days. The number of days and first Julian day were log10 transformed to improve assump;ons of normality.

Network structural metrics
Before calculating network metrics, we normalized network matrices by dividing each cell value by the number of days the community was observed (zeros remain zeros).
We did this because there was unequal sampling effort across studies, even though there was equal effort across networks within studies. Normalizing resulted in non-integer values in some cases, but standardized individuals observed per unit time, therefore the network is weighted by the frequency of interactions for a given monitoring time.
We calculated four species-level network properties: (i) standardized specialization for each species (d′) following Blüthgen et al. (Blüthgen et al. 2006). We used this measure instead of species degree (number of other species the focal species interacts with), because degree is based on a binary matrix, and so does not utilize information on the frequency of interactions. Specialization (d′) calculates how strongly a species deviates from a random sampling of interaction partners available. We also calculated (ii) interaction strength asymmetry (ia), which measures the average mismatch between a focal species' effect on its interaction partners and the effect of the interaction partners on the focal species (Vazquez et al. 2007). In addition, for each network we identiTied modules -sets of species that are more connected to each other than to other species in the network (Olesen et al. 2007)-and for each species we calculated (iii) within-module degree (z, the standardized number of links per species within a module), and (iv) amongmodule connectivity (c, how well does the species connect different modules). We chose these because we were interested in how phenology affected network metrics, so needed to use metrics that were quantiTied at the species level (where phenology was varying). In addition, within-module degree (z) and among-module connectivity (c) characterize the roles that species play in a network, providing a rich way of understanding networks (Olesen et al. 2007, see also Figure 2).
We used the specieslevel function in the bipartite R package (Dormann 2011) to calculate d′, and ia. To calculate the species-level modularity metrics (z and c), we used a modularity-detecting algorithm, which maximized modularity using simulated annealing (SA) implemented in the command line function netcarto_cl in the C program Rgraph (Guimera andAmaral 2005a, 2005b). All other analyses were done with the programming language R (R core team 2018).

How do date of *irst appearance in a network, and length of activity during the season, affect individual species interaction patterns?
We tested for a relationship between the four species-level network structures (d′, ia, c, z) and phenological variables. For ia and z (as continuous values from 0 to inTinity) we used linear mixed effects models, with two phenology variables (date of Tirst appearance and days observed) as Tixed effects and network as a random effect. We also included a random effect for taxonomic group at the family level. This random effect was included as some groups will tend to have a longer duration in the season, for example, bumble bees have multiple generations throughout a season and will therefore be present for a longer period. Including these random effects allows different families to vary either in the intercept or the slopes. We compared Tive models for each predictor where we varied whether the family random effect was: (i) only for the intercepts (1|family), (ii) a random slope for the number of days but no covariance in between the intercept and slope (0+days| family), (iii) a random slope for the Tirst Julian day with no covariance between the intercept and slope (0+Tirst julian |family), (iv) a random slope for days with covariance between the intercept and slope (1 + days | family) and (v) a random slope for Tirst Julian day with covariance between the intercept and slope (1 | Tirst julian |family). For c and d′, variables that are proportions (values between 0 and 1), we used non-linear mixed effects models, with the same formula as above, but specifying a binomial distribution. We selected the best Titting model for the family random effect using AIC (Aho et al. 2014). Models were run separately for plants and pollinators, and for each network metric separately, for a total of eight models. The number of days and Tirst Julian day were log10 transformed to improve assumptions of normality and homoscedascity of residuals.

How robust are networks to removal of species due to varying phenological "traits"?
We were interested in four scenarios with respect to network robustness: 1) species are removed from a network according to when they are Tirst active -i.e., the species that was Tirst active during the earliest date of the season is removed Tirst, and so on; 2) species are removed from a network according to the reverse order of activity -i.e., the species that was Tirst active on the latest date is removed Tirst, and so on; 3) species are removed from a network according to total duration of activity during the season -i.e., the species active the least number of days is removed Tirst, and so on; and 4) species are removed from a network according to reversed order of total duration of activity during the season -i.e., the species active the most number of days is removed Tirst, and so on. These can be referred to as FJ (Tirst julian), FJr (Tirst julian reverse), D (days), and Dr (days reverse), respectively for 1), 2), 3), and 4).
To measure network robustness we used the function second.extinct in the bipartite R package (Dormann et al. 2009), which works by removing a species, then counting how many species are subsequently lost due to the Tirst removal (i.e., secondary extinctions), and so on for each species removed. This is used to generate a curve of number of species going extinct in one group (plants or pollinators) based on extinctions in the other group.
The function robustness in the bipartite package measures the area under this curve to generate a statistic R, with the maximum value R=1 corresponding to a very robust network, whereas R=0 corresponds to a network that loses species quickly upon initial extinctions.
We ran eight separate simulations, four for each of plants and pollinators, with the four for each being one of the four scenarios outlined above. Each run for each network results in a single robustness (R) value. In addition, we calculated R for a set of 100 iterations randomly removing species from each network, rather than removing species based on their phenological traits. We compared the "observed" value of R for the distribution of each simulation run, and calculated a P-value based on whether the observed value fell outside of the 95% conTidence interval of the distribution of R from the randomly removed species set.
To compare results between sets of simulation runs (e.g., to compare FJ and FJr within plants), we used two-sample t-tests.

Do plants and pollinators differ in their species roles in the network?
A useful way to visualize modularity is plotting among-module connectivity (c) against within-module degree (z) following Olesen et al. (2007). In Figure 2,  We ranked models based on AICc, and identiTied top models based on a criteria of ΔAICc < 2.0 from the best model (Steel et al. 2013). On average, species level network metrics were positively related to phenology variables. For both plants and pollinators within-module degree (z) and among module connectivity (c) were positively related to the number of days in a network (plant and pollinator P < 0.001), such that plants and pollinators that were active longer in the season were more connected within the module and among modules (Table 1; Fig 2-3).
The degree of specialization (d') for plants and pollinators was negatively related with the number of days in a network (P = 0.008, P = 0.001 respectively), such that plants and pollinators that were active longer in the season were more generalized. For pollinators, the degree of specialization (d') was negatively related with the day of Tirst appearance (P = 0.035), such that pollinators that appeared later in the season were more generalized (Table 1; Figure 2,3,4).
Last, interaction asymmetry was positively related to the number of days active for both pollinators and plants (plants and pollinators P < 0.001), such that plants and pollinators that were active longer in the season affected the species they interacted with more than those species impacted them. The date of Tirst appearance was positively related with interaction asymmetry for plants (P = 0.019) and negatively related for pollinators (P < 0.001). Plants that appeared later in the season had a stronger effect on pollinators they interact with, while pollinators that appeared early the season had a stronger effect on the plants they interact with. On average pollinators were more strongly affected by the plants they interacted with ( Figure 3D and 3F, points below zero) while the plants affected pollinators more strongly ( Figure 4D and 4F, points above zero).

How robust are networks to removal of species due to varying phenological "traits"?
Robustness values (R) varied from about 0.40 (less robust) to about 0.92 (more robust; Figure 5). Interestingly, values of R were on average higher for pollinator removals than plant removals, suggesting that networks are more robust to removal of pollinators than to removal of plants. For plants, when species were removed Tirst if they appeared early ( Figure 5A) in the network 18 of 33 (55%) were signiTicantly different from random species removal. In contrast, when species were removed Tirst if they appeared late in the network, only 9 of 33 (27%) networks were signiTicantly different from random species removal. Overall, network robustness did not differ with removal of plant species appearing early (mean ± 1 s.e.; 0.62 ± 0.01) vs. appearing late (0.59 ± 0.01; Welch's 2 sample t-test, P = 0.136; n=33 for all comparisons). When species were removed Tirst based on the shortest number of days in a network ( Figure 4B) 29 of 33 (88%) were signiTicantly different from random species removal, whereas only 9 of 33 (27%) networks were signiTicantly different from random species removal when species were removed starting with those active the most number of days. Overall, network robustness was greater with removal of species active shortest (0.67 ± 0.01) vs. active longest (0.56 ± 0.01; P < 0.001).
For pollinators, when species were removed Tirst when they appeared early ( Figure   5C) in the network 0 of 33 were signiTicantly different from random species removal, whereas 29 of 33 (88%) networks were signiTicantly different from random species removal with removal starting with species that appeared late in the network. Overall, network robustness differed for those with removal of species appearing early (0.65 ± 0.01) vs. appearing late (0.85 ± 0.01; P < 0.001). When species were removed Tirst based on the shortest number of days in a network ( Figure 5D) 31 of 33 (94%) were signiTicantly different from random species removal, whereas only 1 of 33 (3%) networks was signiTicantly different from random species removal with removal starting with species that were active the most number of days. Overall, network robustness was greater for those with removal of species active shortest (0.86 ± 0.01) vs. active longest (0.67 ± 0.01; P < 0.001).

Discussion
We asked how species phenological traits (date of appearance in a community, and duration of activity in that community) inTluenced resulting network structure of 33 plant-pollinator mutualistic networks across three habitat types in British Columbia, Canada.
Importantly, we showed in one of the few empirical examinations that phenology, a species attribute that can be sensitive to climate change (Molnár et al. 2012, Bartomeus et al. 2013a), can be an important predictor of structure. This has implications for how mutualistic plant-pollinator networks may respond to climate change via changes in species phenologies.
Species within a module are linked more tightly together than they are to species in other modules. Identifying which species are in these modules can help us understand the mechanisms that structure these networks (Olesen et al 2007). Overall we found that plants were network and module hubs, while pollinators were not. Network hubs -those species that had high within and among network connectivity-and module hubs -species that had high within module but low among module connectivity-tended to be species with radially symmetrical Tlowers and "easy access" Tloral rewards. Morphological adaptations that mean Tlowers can provide food for a diverse array of visiting insect species may increase connectivity within and among modules (e.g. Asteraceae, Asparagaceae, Amaryllidaceae in our study; Wolfe andKrstolic 1999, Sargent 2004).
Pollinators were neither module hubs nor network hubs. This is partly explained by ecological traits of plants and pollinators. For example, pollinator species that were connectors tended to be species that were eusocial (Bombus) or likely social (Lasioglossum (Dialictus) pacatum) which have long activity times as they have multiple generations per season (Michener 2000). In addition, this result can also be partly explained by a statistical effect, because (as in most network research) the number of pollinators and plants differed greatly. Across all studies, the number of pollinators was three times that of plants.
Therefore any given plant is more likely to interact with multiple sets of pollinators than any given pollinator.

How do date of Jirst appearance in a network, and length of activity during the season, affect individual species interaction patterns?
Our goal was not only to describe species' network properties (such as modularity), but also, whether phenology was the determinant of those species network properties.
Phenology was an important predictor of species-level network structures for both plants and pollinators in this study. This is consistent with Tindings from numerical simulation studies (Encinas-Viso et al. 2012) that found that phenology could be extremely important in structuring networks. For plants we found that species active longer in the season were more connected both within their module and among modules. This suggests that species are more likely to interact with many different partners the longer they remain active. In addition, plant species that were active longer in the season were also less specialized in their interactions, and had a stronger effect on pollination partners than partners had on them. This suggests that plant species that are active longer during the season are important structurally in the network as they interact with multiple species, exert strong effects on interaction partners and are generalists (Guimarães et al. 2007, González et al. 2010). Some of the species that were active the longest were Cornus stolonifera (Family: Cornaceae), Symphoricarpus albus (Caprifoliaceae) and Ranunculus repens (Ranunculaceae). Ranunculus is a weedy forb while Symphoricarpus and Cornus are a long lasting shrubs. These species additionally have radially symmetrical Tlowers allowing many types of pollinators to access the Tlowers (Lovett-Doust et al. 1990, Wolfe and Krstolic 1999, Stewart-Wade et al. 2002, Sargent 2004).
For pollinators, species that were active longer in the season had more connections among and within modules, were less specialized and had a higher interaction strength asymmetry. These results are consistent with the results we found for plants. Some of the species that were active the longest in the season were Bombus centrals (Apidae), Sphaerophoria weemsi (Syrphidae), and Lasioglossum pacatum (Halictidae). These tend to be highly generalized in Tloral visit patterns (Laverty and Plowright 1988). Bombus for example have multiple generations of workers per season (Michener 2000) and although individual workers may specialize on particular Tloral resources, the species as a whole uses diverse plants over an extended Tlight period. Hover Tlies (Syrphidae) also are mutivoltine and use a diverse array of generalized Tlowers. In contrast, solitary bees in our region are active for only a few weeks, making it unsurprising that they are not within or among module connectors. Similarly, we found that pollinators were less specialized as Tirst Julian day increased. Therefore pollinators that emerged later in the season were more generalist, which in our data would include wasps and some hover Tlies.
We also found that plant species that emerged later in the season (i.e. their Tirst julian day was higher) had a higher interaction strength asymmetry. Therefore species that emerged later in the season affected their interacting partners more strongly than their partners affected them. Some of the plant species that emerged later in the season were Galeopsis tetrahit (Family: Lamiaceae), Lythrum salicaria (Lythraceae) and Polygonum persicaria (Polygonaceae). Unlike plants, we found that for pollinators the interaction strength asymmetry decreased as species emerged later in the season (i.e. their Tirst julian day was higher). Therefore species that emerged later in the season were more strongly affected by interacting partners than the reverse, and species that emerged earlier in the season affected interacting partners more strongly than partners impacted them. Some of the pollinator species that emerged earlier in the season were solitary bees Andrena nigrihirta, Andrena merriami, Andrena sladeni, Andrena trizonata and Andrena porterae (Andrenidae). In two of the ecosystems studied, oak-savannah and shrub-steppe, we observed that at the beginning of the season, many plant species bloom in high density, but temperatures are not yet reliably warm enough for insect activity, so plant reproduction may be pollinator limited (Schemske et al. 1978, Kudo andIda 2013). By the end of the season the amount of food and nutrients available for the pollinators is lower, but their populations have grown; it may be that late emerging plant species are therefore extremely important resources for pollinators (Mattila andOtis 2007, Garbuzov andRatnieks 2014).
On average, pollinators were more affected by plants (values below 0 on Tigure 3F and E), while plants affected the pollinators more (values above 0 on Tigure 4F and 4E). Overall these results suggest that resource limitation may shift along the season from pollinator limitation to Tlower limitation, a hypothesis that could be tested.

Robustness
We observed that networks were on average more robust to removal of pollinators than to removal of plants. This makes sense because plants more often link together the plant-pollinator network (hubs organize around plants), whereas pollinators are less often important hubs (see Figure 2). This pattern was also seen in another study -Olesen et al. (2007) found that plant species were more often module hubs and network hubs than pollinators (see Figure 2 in Olesen et al.). The networks we used have 3x as many Pollinators that appeared early in the season were more important for network robustness than the pollinators that appeared late in the season. These results are consistent with our results for plants and pollinators relating interaction-strength asymmetry to the Tirst julian day of appearance. Again, at the beginning of the season, the system may be more pollinator limited while at the end of the season pollinators may be resource limited (Schemske et al. 1978 (Bartomeus et al. 2011, Kudo andIda 2013).
Increasing the temporal mismatch in the spring can further increases the pollen limitation experienced by plants early in the season (Schemske et al. 1978). This mismatch can also reduce the robustness of pollinator networks as shown in our study.

Conclusion
Our results show that across a large sample of 33 networks, species phenology can be an important predictor of network structure. In particular, the number of days a species is active in a network predicted how connected they are for both plants and pollinators. We also found that plants tended to be more important in the network at the end of the season, while pollinators were important at the beginning of the season. Future work should build on the work presented here by exploring how experimental or natural changes in phenological variables, like time of Tirst appearance or duration of activity in a community, inTluence network structure.  Figure 1. Visualization of the phenology of species in each network, for all 33 networks.
The top set of horizontal lines in each panel are pollinators, and the bottom set are plants.
Note that each panel follows the same x-axis.    last day minus Tirst, in number of days). Black symbols were signiTicantly more robust than the null model (random species removal), while grey symbols did not differ from the null model. Note: networks in each of the four panels are in the same order as Figure 1; scales in all four panels are the same; drop lines connect points that belong to the same network in each panel.