Phenology determines the robustness of plant–pollinator networks

Plant–pollinator systems are essential for ecosystem functioning, which calls for an understanding of the determinants of their robustness to environmental threats. Previous studies considering such robustness have focused mostly on species’ connectivity properties, particularly their degree. We hypothesized that species’ phenological attributes are at least as important as degree as determinants of network robustness. To test this, we combined dynamic modeling, computer simulation and analysis of data from 12 plant–pollinator networks with detailed information of topology of interactions as well as species’ phenology of plant flowering and pollinator emergence. We found that phenological attributes are strong determinants of network robustness, a result consistent across the networks studied. Plant species persistence was most sensitive to increased larval mortality of pollinators that start earlier or finish later in the season. Pollinator persistence was especially sensitive to decreased visitation rates and increased larval mortality of specialists. Our findings suggest that seasonality of climatic events and anthropic impacts such as the release of pollutants is critical for the future integrity of terrestrial biodiversity.


Supplementary Information S1: Model details
The model we developed governs the biomass dynamics of a complex interaction system of flowering plants and their insect pollinators. Each plant species is composed of four state-variables: (1) Adult plants, able to produce flowers and seeds, (2) flowers, that can be visited by pollinators to be fertilized, (3) floral rewards, which are used by pollinators as food, and (4) seeds, produced by adult plants and assumed to be obligate dormants.
For simplicity, all flower-visiting insects are considered as pollinators. They are structured into two developmental states: (1) Adult insects, able to visit flowers, consume floral rewards and produce larvae, and (2) Larvae, a single juvenile stage assumed to be the dormant state of pollinators. For all states, biomass density (g/m 2 ) is used as the unit.
Transition from a dormant state, either seed or larva, to the adult state occur when two conditions are simultaneously met: (a) a given latency period has finished, implying that the dormants reached a minimum age, and (b) environmental conditions are favorable for seed germination or insect recruitment, which is dictated by time-dependent functions f G (t) and f R (t) respectively (see Fig.4).
Adult plants can produce flowers in a given period of the year, and flowers produce floral rewards for pollinators. The phenological flowering dynamics is governed by the time-dependent function f F (t). During flowering, adult insects visit flowers and consume floral rewards. As a consequence of visits, some flowers are fertilized, getting unavailable for further interactions, pollinators obtain food that will be converted into offspring and plants reproduce giving rise to dormant seeds.
The mathematical model proposed is a system of integro-differential equations on partial derivatives. Phenophases of germination, flowering and insect recruitment are described by time-dependent functions, named here "phenology functions" ( f (t)), that describe the within-year time evolution of the level of environmental favorability for each of the phenological events. All the state-variable values depend on chronological time (t). For dormants (seeds and larvae), state-variable values depend on age (u) too.
Variation on seeds biomass is given by : where the first term on right side is an advection term, which describes loss of seeds of age u due to aging. The second term is seed mortality and the last term represents biomass loss due to recruitment from dormant to adult stage. This term is composed of parameter s i , representing the rate at which seeds transit to adult plants, environmental favorability for germination f G i (t), and the maturation function (g(u)) that describes the biological readiness for germination as a function of seed age: with u 0 being an age threshold above which seeds are mature and ready to germinate. Seed production is given by: and is proportional to the visitation rate (j i j ) to flowers of plant species i (F i ) by pollinators of specie j (A j ). This function was modeled as a Beddington-DeAngelis-like functional response: where the set R( j) are plant species visited by pollinator j and the set C(i) are pollinator species that visit plant i.
Biomass dynamics of adults plants are given by: Plant biomass growth, due to seed germination, is limited by intra-and inter-specific competition, which is regulated by parameters a i and a k . The last term is mortality rate of plants.
Biomass dynamics of flowers are described by: where the first term represents limited flower growth sustained by plant biomass, and favored by the phenology function f F i . The second term represent time-dependent mortality, composed by a basal mortality rate µ F i and a seasonally-driven extra mortality that is maximal, with value r F i outside the flowering phenology period dictated by f (t). The third term is a removal of fertile flower biomass due to visitation-driven fertilization.
The dynamic equation for floral resources is: which has a term of limited growth sustained by flower biomass, and a term of biomass loss due to consumption by pollinators. Biomass growth rate of insects is given by: This equation has a structure equivalent to the equation of seeds (1). Production rate of larvae is described by the boundary condition: Production of larvae depends on consumption of floral resources by adult insects, and it is limited by intra-and inter-specific competition, governed by w i and w j .
Finally, biomass density growth rate of adults insects is: The population growth rate of insects is governed by recruitment, in an equivalent form than for plants (eqn. 5). Insect mortality has, similar to eqn. 6, a basal mortality term µ A j and a seasonal mortality with maximum r A j outside the recruitment phase dictated by f R j . All parameter values were drawn at random from uniform probability distributions, centered at estimated mean values (Table ) ±25%. Initial estimation of mean parameter values was based on biological plausibility and in some cases guided by literature data available for some species. After that, mean values were adjusted as to avoid very low species persistence.

Symbol Description
Value Unit µ S Basal seed mortality 1.3 · 10 2 week 1 µ P Basal plant mortality 2.8 · 10 3 week 1 µ F Basal flower mortality 8.8 · 10 1 week 1 µ A Basal adult insect mortality 4.4 · 10 2 week 1 µ L Basal larva mortality 1.0 · 10 2 week 1 r F Maximal winter extra mortality of flowers 10.0 week 1 r A Maximal winter extra mortality of insects 10.0 week 1 r F Production rate of flowers by unit plant biomass 2.0 week 1 r N Production rate of floral resources by unit flower biomass 4.2 week 1 n Conversion rate from seeds to plants 1.0 · 10 3 dl y Conversion rate from larvae to adult insect 4 · 10 1 dl s Germination rate of mature seeds 1.1 · 10 1 week 1 c Recruitment rate of mature larvae 1.0 week 1 k Flower carrying capacity per unit plant biomass 2 · 10 2 dl q Floral resources carrying capacity per unit flower biomass 1 · 10 1 dl t Maximum visitation rate per unit insect biomass 1 · 10 6 g 1 m 2 week 1 d Resource consumption per visit 1.5 · 10 2 dl g Flower biomass loss per visit 7 · 10 3 dl p Maximum seed production per visit 1.2 · 10 1 dl h Conversion efficiency of floral resources to larvae production 4.1 · 10 1 dl b Half-saturation parameter of visitation rate function 1.0 dl e Half-saturation parameter of per-visit consumption function 2.5 · 10 2 dl a Coefficient of competition among flowers for insect visitation 1.0 · 10 3 g 1 m 2 b Coefficient of competition among insects for visiting flowers 1 · 10 2 g 1 m 2 a ii Coefficient of intra-specific competition among plants for space 1 · 10 3 g 1 m 2 a i j Coefficient of inter-specific competition among plants for space 1 · 10 4 g 1 m 2 w ii Coefficient of intra-specific competition among insect larvae for space 5 · 10 2 g 1 m 2 w i j Coefficient of inter-specific competition among insect larvae for space 5 · 10 1 g 1 m 2 Symbol Description Value Unit Extinction threshold P i (0) Adult plants 9.9 · 10 1 g m 2 1 · 10 3 Floral resources 3 · 10 2 g m 2 0 A i (0) Adult insects 3 · 10 3 g m 2 1 · 10 8 S i (0, u) Seeds of age u 6.8 g m 2 week 1 1 · 10 5 L i (0, u) Insect larvae of age u 7 · 10 4 g m 2 week 1 1 · 10 8 Supplementary Information S3: Simulation results for all tested plant-pollinator networks  Figure S1: Simulation results for the Zackenberg (1996) network. Left/right plot shows decreases in plant/pollinator species persistence (mean ± 95% C.I.) after exerting different pressures, relative to species persistence without exerting any pressure (see text for details). Experimental pressures, applied one by one, were increased plant competition a, decreased visitation rate t, increased pollinator mortality µ A , larvae mortality µ L , plant mortality µ P and seed mortality µ S . Pressure intensities were assigned to species according to their attributes: R: random, D: degree. Ps: phenophase start, Pe, phenophase end, Pd: phenophase duration. In blue, pressure magnitude increased on species of increasing attribute values (low to high D, earlier to later Ps or Pe, and small to large Pd). In orange, pressure magnitude increased on species of decreased attribute values. The smaller bars are shown in front for better view.  Figure S5: Results of non-metric multidimensional scaling (NMDS) analysis for the 12 networks under study, belonging to four locations: Zackenberg (blue), Nahuel Huapi (red), Athens (violet) and Villavicencio (green).

Supplementary Information S5: Sensitivity analysis
Here we report results of simulations for the most recent networks at each of the four locations. These simulations were run with the same mean values of parameters and initial conditions, but with a random variation of 10%, 50%and75% around the mean. Results shown in Supplementary Information S3 were obtained with a variation of 25%.  Table 2. Mean (variance) of corrected impacts on plant species persistence of pressures exerted on plant competition (a), visitation rate (t), adult pollinator mortality (µ A ), larval mortality (µ L ), plant mortality (µ P ) and seed mortality (µ S ) of species, according to their degree (D), phenophase start (Ps), phenophase end (Pe), and phenophase duration (Pd). Pressure magnitude increased on species of higher attribute values (e.g. high D, later Ps or Pe, and large Pd). Only positive mean values were considered a t µ A µ L µ P µ S D -0.007 (