High Arctic flowering phenology and plant–pollinator interactions in response to delayed snow melt and simulated warming

The projected alterations to climate in the High Arctic are likely to result in changes to the short growing season, particularly with varying predicted effects on winter snowfall, the timing of summer snowmelt and air temperatures. These changes are likely to affect the phenology of interacting species in a variety of ways, but few studies have investigated the effects of combined climate drivers on plant–pollinator interactions in the High Arctic. In this study, we alter the timing of flowering phenology using a field manipulation experiment in which snow depth is increased using snow fences and temperatures are enhanced by open-top chambers (OTCs). We used this experiment to quantify the combined effects of treatments on the flowering phenology of six dominant plant species (Dryas octopetala, Cassiope tetragona, Bistorta vivipara, Saxifraga oppositifolia, Stellaria crassipes and Pedicularis hirsuita), and to simulate differing responses to climate between plants and pollinators in a subset of plots. Flowers were counted regularly throughout the growing season of 2015, and insect visitors were caught on flowers during standardised observation sessions. As expected, deep snow plots had delayed snow melt timing and this in turn delayed the first and peak flowering dates of the plants and shortened the prefloration period overall. The OTCs counteracted the delay in first and peak flowering to some extent. There was no effect of treatment on length of flowering season, although for all variables there were species-specific responses. The insect flower–visitor community was species poor, and although evidence of disruption to phenological overlaps was not found, the results do highlight the vulnerability of the plant–pollinator network in this system with differing phenological shifts between insects and plants and reduced visitation rates to flowers in plots with deep snow.


Introduction
The climate of Arctic regions is undergoing a significant period of change, and further changes in temperature and precipitation are expected (Thompson and Wallace 2001, Barber et al 2008, Comm. I.A.S. 2010). These relatively rapid alterations are likely to affect Arctic ecosystems via changes to an already short growing season (Cooper 2014). For example, it is anticipated that increasing winter and spring temperatures will reduce the proportion of precipitation falling as snow (Comm. However, studies from the High Arctic concerning the interaction of delayed snowmelt and increased temperature are lacking (see Aerts et al 2004, 2006, Legault and Cusa 2015. This study adopts a factorial experimental design to test the interactive effects of delayed snowmelt date and simulated warming on flower production in a High Arctic tundra system in Svalbard at 79°N. In addition, we investigate whether changes to flower timing impact insect visitation to dominant flower species using a subset of the plots. Whether future snow melt dates are earlier or later in the Arctic in a warmer future, the timing of key events in the life cycles of plants and animals, their phenology, are likely to be significantly affected Plant-pollinator phenological 'asynchrony' has been documented in some cases (e.g., Visser andHolleman 2001, Both et al 2006), and has been absent in others (Bartomeus et al 2011). As a result, it is unclear whether the phenomenon poses a critical threat, or whether it is actually a natural part of species interactions (Olesen et al 2008, Bartomeus et al 2013). However, some systems may be more at risk than others. Arctic plant-pollinator communities are less diverse than those in temperate regions, and this together with the short season length, is likely to contribute to instability among the interactions (Encinas-Viso et al 2012, Benadi et al 2014). In this study, we alter the timing of flower emergence using a field manipulation experiment in which snow packs are increased by fences that accumulate winter snow and temperatures are enhanced by open top chambers. While the insect pollinators were not manipulated by this design, we used a subset of plots from the experiment to simulate differing responses to climate between plants and pollinators. This approach has been adopted in temperate regions (Gezon et al 2016) and by using potted plants (Parsche et al 2011, Rafferty and Ives 2011), but our study represents the first time such a study has been conducted in the Arctic, and the first in depth study of insect pollinators on Svalbard. We hypothesised that plots with increased snow would exhibit delayed flowering of the six dominant species compared to control plots, but that this effect would be counteracted by increased temperatures. We further expected that flowering plant species with delayed phenology would receive fewer visits from pollinating insects due to a subsequent mismatch in the overlap of their populations. We expected that this would demonstrate the potential for phenological asynchrony in the High Arctic and highlight the vulnerability of such low diversity systems.

Field site and experimental design
The study area and experimental design has been described in detail elsewhere (Cooper et al 2011), but brief details are given here. The field site was situated in Adventdalen (78°10′N, 16°06′E), a large valley on Spitsbergen, Svalbard. Annual mean temperature is −6.7°C (1969-1990Svalbard airport;www.eklima. no) with July the warmest month (5.9°C). Annual mean precipitation is 190 mm, the majority of which falls as snow although some studies suggest that this may be an underestimation (Forland and Hanssen-Bauer 2003). Snow depth ranges from 3 m in hollows and gulleys to less than 20 cm on ridges (Cooper et al 2011). Other environmental information about the study site is detailed in Semenchuk et al (2016), Mallik et al (2011) andBlok et al (2015). Plots were established on two habitat types: (1) heath with rough stony soils, taller vegetation and dominated by plants such as Cassiope tetragona, Dryas octopetala, Salix polaris and Saxifraga oppositifolia, and (2) mesic meadow with shorter vegetation, flatter topography and dominated by Salix polaris, Luzula arcuate subsp. confusa and Alopecurus magellanicus. The study area was approximately 2.5 km×1.5 km, with two blocks of approximately 200 m×200 m established on each habitat type. Each block contained three snow fences and three control areas, although one fence from each habitat type was unusable during this study. This resulted in a total of 10 fences and 10 control areas. Fences were 1.5 m tall and 6.2 m long, aligned perpendicular to the prevailing winter winds and accumulated deeper snow behind them. Previous studies on this site have shown that the fences increase snow depth by a mean of 1 m compared to controls (Cooper et al 2011). Behind each fence and in each control, and where possible (table 1), six plots (0.75 m×0.75 m) were established, three containing C. tetragona and three containing D. octopetala as these species are expected to respond differently to treatments and do not always occur together. Three additional plots were established in a 'medium' snow depth area behind the fences. In addition, open-top chambers (OTCs) were used to simulate increased temperatures (Marion et al 1997). An OTC was placed over one medium snow plot, one control D. octopetala plot and one Deep snow D. octopetala plot at each fence. We were restricted in the number of OTCs available, so chose the D. octopetala plots as the medium snow plots tended to contain more D. octopetala than C. tetragona. OTCs were placed on the first snow free plot within each snow treatment.

Abiotic measurements
Soil surface temperatures were measured hourly by data loggers (Gemini Data Loggers, Tinytag, UK) placed in a subset of plots: at each OTC plot and three per fence in each level of snow treatment. Soil moisture was measured using a handheld moisture metre (ThetaProbe ML2x; Delta-T Devices, UK) whenever plant measurements were taken. The four corners of each plot were measured and a mean calculated. Snowmelt was determined daily between 20 May and 19 June, and the date that each plot became 50% snow free was recorded.

Flower phenology
The number of flowers of 6 focal species was counted every 3-4 days from the date of first flower to the date of the last flower. The focal species were D. octopetala, Bistorta vivipara, C. tetragona, Stellaria crassipes, Saxifraga oppositifolia and Pedicularis hirsuta. Flowers were considered to have opened when the flower buds had broken and anthers/stigmas were visible, and to have senesced when all petals had dropped off. Flower phenology was summarised by the timing in Julian days of first flower and peak flowering, and by the length of the flowering period (day of last flowering minus the day of first flowering) and the length of the prefloration period (day of first flowering period minus the day of snow melt). Peak flowering was calculated using the weighted mean date of occurrence (WMD), which is the arithmetic mean of all dates on which the species was observed, weighted by its abundance on each date.

Insect sampling
Due to limited time resources, insect sampling was restricted to the mesic meadow habitat and to a subset of plots (table 1). Sampling was carried out whenever the weather was suitable during the flowering season (wind calm, sunny days). This resulted in 31 days of insect sampling, from a total of 81 days of flowering season. Individual plots were observed between 10:00 h and 16:00 h for 10 min intervals each day. Whenever an insect visited a flower in the plot during this time, it was captured for later identification and the flower species it was visiting was recorded. The number of flowers available within the plot was also counted. The order of plot sampling was randomised each day. Insects were identified in the laboratory to species where possible. Good identification keys are not currently available for all Svalbard insect families however, so for the purposes of this study, identification to the lowest common level was used (Family level).

Statistics
All statistical analyses were carried out in the R programming environment (v 3.2.3, R Core team 2015) using linear mixed effect models (LME; lme4 packages; Bates et al 2015) or generalised additive mixed modelling (GAMM, mgcv package; Wood 2011) due to the unbalanced nature and nested structure of the experimental design. In all statistical tests, an initial maximal model was fitted with all fixed effects and no random effects. This model was then compared to models with different random effects structures using AIC. The model with the lowest AIC was chosen, and subsequent backward stepwise selection was used to eliminate insignificant fixed effects.
Snow melt out dates in Julian days were analysed using a LME with Snow treatment (3 levels) and habitat type (2 levels) as fixed factors, and fence number nested within block as random effects. The OTC treatment was left out of this analysis because OTCs were placed on the first available snow free plot behind fences. Soil surface temperatures and soil moisture readings were averaged over time and analysed in the same way.
The four flowering phenology variables were analysed using LMEs with snow and OTC treatments as fixed effects, and plot number nested within snow treatment, nested within fence number, nested within Block as random effects, and an additional random effect for plant species. Treatment effects were also analysed on each plant species separately, using similar models but without the random effect for species. Insect and plant phenological synchrony was tested using LMEs with the WMD of each species on each plot as response variable, and treatment, trophic level and the interaction between each treatment and trophic level as explanatory variables. We expected that if delayed flowering was disrupting plant-pollinator synchrony, the difference in peak occurrence between treatments would be greater for plants than for insects. Here the nesting structure detailed above and species identity were random factors. WMD for insects was calculated based on number of individuals divided by the number of flowers in the plot.
An alternative indicator of phenological synchrony is the correlation between insect visits and number of flowers (Kudo 2014), with a strong relationship expected to result from phenologically matched species. The number of insect visits (log(n+1)) was therefore modelled against the number of flowers per plot (log(n+1)), snow and temperature treatments, Week (continuous), time of day and maximum temperature and rainfall on day of sampling, with 2-way interactions included between the treatments and number of flowers, and week and number of flowers. To account for non-linear effects, a GAMM was used, initially fitting a cubic shrinkage smoothing term for each variable. Insignificant smoothing functions were removed and replaced as linear fixed effects. Only a smoothing term for time of day was included in the final model. The random terms in the model include a random intercept for plant species and plot number, and a random slope for day of sampling (continuous). An auto-regressive correlation structure was also included to account for temporal pseudoreplication.
Finally, to test for differences between treatments in insect visitation, a series of LMEs were implemented using the total number of insect visits, the rate of insect visitation and the number of days of visitation by insects as response variables, the two treatments and their interaction as fixed effects and the nested structure described above as random effects. Response variables were log transformed where required to normalise residuals.

Abiotic measurements
The snow fence treatment and habitat type had a significant effect on the plot melt out dates (table S1). As with previous years, plots within the meadow became snow free earlier than those within the heath, and control plots were free of snow before medium snow plots, which were in turn free of snow before deep snow plots. OTCs significantly increased the soil surface temperature of warmed plots by approximately 1°C (±0.3 s.e.; table S2, figure S1). OTCs also significantly reduced soil moisture (table S3; figure  S3), but there was only a weak significant difference between the control and deep snow treatments when the data for June only were analysed (table S3).  Species varied in their overall numerical response to treatments, with C. tetragona and B. vivipara producing more flowers in deep snow treatments, and D. octopetala producing more in control snow plots. In addition, flowers of D. octopetala, S. crassipes and B. vivipara were more numerous in OTC plots compared to ambient plots.

Flower phenology
When all species were analysed together (table 2), snow treatments significantly delayed the first and peak flowering dates of the focal community of flowers, and the deep snow treatment exhibited a significantly shorter prefloration period than that of control, ambient plots. The warming treatment advanced both the first and peak flowering dates compared to control, ambient plots, but did not impact the prefloration period. The weak snow x OTC effect on peak flowering date suggests that the warming treatment counteracted some of the delay in peak flowering. However, due to the highly unbalanced nature of the presence of focal species in the plots this result should be treated with caution. There were no effects of the two treatments on the length of the flowering season overall (results not shown).
When analysed separately, the focal plant species showed differing responses to the experimental treatments (figure 3). Snow depth treatment significantly delayed the date of first and peak flowering and shortened the prefloration period in all species except for S. oppositifolia. The snow treatment only had a significant negative impact on the length of the flowering season for C. tetragona and B. vivipara, however this was increased for these two species under simulated warming (weak effects). Plot warming by OTCs had the effect of advancing the first and peak flowering dates of D. octopetala, P. hirsuta and S. crassipes (weak effect), and the peak flowering of S. oppositifolia. The only significant interaction between the treatments was for the length of the B. vivipara blooming period  where the deep snow and warming treatments combined to increase the length of the flowering season.
Insect and flower phenological synchrony A total of 263 insects belonging to 8 families were captured visiting the six plant species. The majority of insects were from the Order Diptera (72%), with the remaining from a single family within the Hymenoptera, the Ichneumonidae (parasitoid wasps). Within the Diptera, the majority of individuals were tiny darkwinged fungus gnats (Sciaridae), followed by the dagger flies (Empididae), non-biting midges (Chironomidae) and house flies (Muscidae). Most observed flower visits were made to D. octopetala flowers (139 visits mostly from Empididae and Sciaridae), followed by S. crassipes (69 visits, mainly by Ichneumonidae) and B. vivipara (54 visits, mainly by Sciaridae and Ichneumonidae). Only 3 flower visits were recorded from P. hirsuita and 5 from C. tetragona, so these species, together with Saxifraga oppositifolia which was not recorded within the mesic meadow plots, were excluded from plant-insect analyses. The only specialist insect family was the Empididae, which solely visited D. octopetala flowers. The numbers of the main insect visitors to flowers of these species over the season are shown in figure 4. These graphs clearly show that the insects generally begin visiting flowers at or around the beginning of the flowering season, and peak at times within the main flowering period, regardless of treatment.
The results of the analysis between insect and flower WMD are shown in table 3. In all three cases, there was a significant effect of snow treatment, suggesting that the peak occurrence of both insects and plants was delayed in the deep snow plots. The lack of trophic level effect suggests that within control plots, the difference in peak occurrence for the two groups does not differ, indicating phenological matching. However, the significant snow x guild interaction implies that the delay in peak occurrence in response  to the deep snow treatment was greater for plants than for insects ( figure 5).
The results of the GAMM for the number of pollinator visits analysed against the number of flowers are shown in table 4. Overall, the number of visits was strongly linked to the number of flowers present in the plot. This relationship changed over time and within the different snow treatments: the overall relationship between insect visits and flowers was stronger in the deep snow treatment than in controls. However, within the deep snow plots, flowering and insect visits occurred for approximately one week less than in control plots and the change from positive to negative relationship occurred at different times ( figure 6).
Total insect visits were significantly greater in OTC plots compared to ambient controls (table 5), but there was no effect of snow treatment. The same result was found for the number of days on which pollination was observed. However, the rate of visitation was significantly lower in deep snow treatment.

Flower phenology
In this study we provide further evidence that High Arctic plants alter their phenology in species-specific ways in response to climate factors, supporting Table 3. Effects of snow-fence treatment on peak occurrence of insects and flowers in Julian days (189=8 July 2015). The effects of OTC and Guild were removed as they were insignificant.    vivipara, flowering phenology may be more sensitive to temperature or day length in order to benefit from stable, peak summer temperatures (Molau 1997). Our results partly support these expectations: increased snow depth delayed the first and peak flowering dates of five of our study species. We did not get this result for the early flowering S. oppositifolia possibly due to generally low numbers of this species in our plots. Additionally, simulated plot warming advanced the phenology of three early-and one late-flowering species. The implications of these changes are that most of our focal species can track changes in environmental conditions, but the delay in snow melt is likely to be most problematic for later-flowering species as the compressed, shorter growing season may compromise seed ripening (Thorhallsdottir 1998, Cooper et al 2011, Semenchuk et al in press (this issue of ERL)).
Plant responses can also be interpreted by the classification into periodic or aperiodic species (sensu Sørensen 1941). Periodic species, are those that exhibit a defined growth period controlled by genetic constraints, whereas aperiodic species often show plastic responses to variable environmental conditions. Of our focal species, B. vivipara (Starr et al 2000, Rumpf et al 2014) and C. tetragona (Rosa et al 2015) have been previously classified as periodic species. However, these two were the only species to decrease flowering season length under deep snow treatments, and to increase it under simulated warming. Thus, although periodic species have fixed vegetative growth periods and may suffer from competitive disadvantages in a warming world due to an inability to take advantage of extended growth seasons (Lechowicz 1995), the plasticity of their flowering period may compensate when temperatures increase, or exacerbate the situation if snow melt is delayed in isolation. Conversely, the lack of plasticity of aperiodic species' flowering periods may be beneficial in a climate change scenario, because floral duration is positively linked to reproductive success   competitive advantages if they can maintain a normal length of flowering season.
The prefloration period of all species except S. oppositifolia was shortened under the deep snow treatment, and this has been demonstrated elsewhere (Bienau et al 2015) and can be stronger for late-flowering species (Van der Wal et al 2000). It is thought that when snow melts later into the summer, plants emerge to warmer temperatures and can therefore flower quicker (Kawai and Kudo 2011). This may ensure that plants synchronise flowering across snow gradients (Bienau et al 2015) which in turn is likely to be important for attracting insect pollinators during optimal environmental conditions. However, the lack of a warming effect on the prefloration period for all species does not support this and indicates that there are limits to the ability of species to accelerate flowering.
Given the patterns demonstrated at the plot-level in this study, a number of landscape-scale implications may be envisioned. As previously hypothesised, there are likely to be 'winners' and 'losers' at this scale in the event of widespread snowmelt delays, potentially resulting in at least short-term community changes (Legault and Cusa 2015, Wahren et al 2005) and general declines in productivity (Wipf and Rixen 2010). Early-flowering and aperiodic species may be least detrimentally affected by the shortening of the growing season, while late-bloomers and periodic species will probably suffer from such changes. However, the shortened growing season may lead to more synchronised flowering between species, and subsequent increased competition for pollinators (Molau 1997, Wipf andRixen 2010). Furthermore, to some extent these effects will be offset by changes in temperature and summer precipitation, phenotypic adaptation

Insect phenology
Mixed responses to treatments were also found for insect visits in this system, and while the results are not clear cut with respect to phenological synchrony, they do highlight a number of aspects of the Svalbard plant-pollinator community that could be sensitive to future change. For example, the peak dates of occurrence of insects and flowers are not sufficient to show phenological mismatches in a delayed snow melt scenario because the insect visit data are not independent from flower data (Forrest 2015). However, the data do show that if insects and flowers differ in their emergence times, peak occurrence is likely to differ more with delayed snowmelt implying a reduction in the overlap between plant and insect populations. Such timing shifts, although not constituting phenological asynchrony, could be enough to result in major ecological changes (Fabina et al 2010). On the other hand, the insect visitation data (figure 4) do suggest that a certain amount of asynchrony is within the normal pattern of plant-pollinator interactions in this system: some flower visitation occurred outside the mean flowering period for all three of the main species. This may indicate that sufficient floral resources exist for insects to maintain populations across the range of flowering times seen here.
Plot-level delayed flowering may have also extended the season for flower visiting insects in the surrounding landscape. This is supported by the results of the GAMM analysis, where the positive relationship between insect visits and number of flowers weakens at around week 31 in control plots, but is maintained in the deep snow plots. This suggests that insects were switching to the latter, perhaps because of greater apparency of the delayed floral resources compared to senescing flowers in ambient plots and the wider landscape. Variable flowering timing is to be expected in a landscape characterised by ridges and hollows, and strong-flying insects are likely to be attracted to flower patches as they become available. In a landscape-scale deeper snow scenario however, insect emergence would also be delayed somewhat, and in a shortened season fewer weeks of interactions can be accommodated ( figure 6). However, this finding is not supported by the overall number of days of insect-flower interactions.
In general, the total number of visits did not differ between snow treatments, but the rate of visitation was significantly reduced. Thus, the same number of visits were being paid to fewer flowers in control plots (data not shown). This suggests either higher visitation efficiency (insects are better at finding flowers) or higher numbers of available pollinators (likelihood of capture is greater) in control conditions. In either case, the probability of visitation is greater in control conditions which may be because of a slight difference in the insect visitation community. For example more large strong-flying Diptera are likely to be available early in the season, while late-emerging parasitoid wasps preferring S. crassipes would have been available during peak or late flowering in deep plots. There were insufficient data to test this, but previous studies have demonstrated negative impacts of differing pollinator community assemblages on the effectiveness of pollinators and plant reproduction (Rafferty and Ives 2012, Gezon et al 2016). In a landscape-wide delayed snow melt scenario, the delayed emergence of early fliers may bring them into closer temporal synchrony with late-fliers, increasing interspecific competition for resources. However, further study is required to confirm the likelihood of this.
We found little effect of the warming treatment on insect visits, and this is probably due to the low number of replicates and because the OTC treatment did not significantly advance flowering phenology of the three main species. The increased total number of visits and number of days in OTC plots may have been an artefact of the treatment therefore, with warmer temperatures inside the chambers attracting insects and temporarily trapping them inside. Nevertheless, the combined effect of air temperature and snow melt date on insects is not well understood (Hoye and Forchhammer 2008b) and further studies examining the mechanisms driving insect emergence and the end of insect activity are needed to make realistic predictions. While insect emergence is likely to be closely related to snow melt timing (Hoye and Forchhammer 2008a), the air temperatures when insects do emerge will impact their activity, abundance, life span and the timing of end of activity (Hoye and Forchhammer 2008b, Iler et al 2013, Scaven and Rafferty 2013). All of these aspects of insect biology will combine to determine how aligned the populations are with the flowering plants they visit. Furthermore, with this experimental design we have only manipulated one side of the plant-pollinator mutualism. We have not tested whether different groups of insects differ in their responses to snowmelt delay and how this affects phenological synchrony (Rafferty and Ives 2011). We also only considered plant responses at the plot-level: while we delayed the phenology of flowering in deep snow plots, flowering and pollination in the surrounding landscape continued as normal and pollinator populations may have been declining when the deep snow plots reached peak flowering. However, in a landscape-scale delayed snowmelt scenario the availability of pollinators to flowers is likely to be greater assuming relative synchrony of emergence, so our plot-level visitation rates may have been underestimated and should be viewed with caution.

Conclusions
Despite the interesting trends demonstrated by this study, the findings should be viewed with an understanding of the limitations of this kind of study. While plot level studies provide important insights into possible future impacts of changing environmental factors, the interpretation of results to landscape scales is difficult. The wide variation in spring snow depths in Adventdalen alone cover the variability in our experimentally increased snow depths, and the plasticity of plant flowering phenology has probably evolved in response to spatial and temporal variability. Therefore, although we have demonstrated delays in flowering phenology, these delays may well be within the ordinary range of plastic responses to environmental conditions. The sampling of insect flower visitors at the flowers themselves is also a limited method because the two datasets cannot be independent of each other. For example, insects may have emerged long before flowering, underestimating our observed effects for example. Nevertheless we have demonstrated that delays in snow melt timing can significantly delay the timing of flowering of several Arctic tundra species, and the differential responses of our focal species suggest further impacts on community composition. Furthermore, our findings suggest that some form of differential response to climate change of flowering-visiting insects is likely to lead to altered patterns in visitation for key species. In such a low diversity ecosystem with short season lengths, this is likely to have connotations for both the insect and plant communities and highlights an area of research that requires further study.