Phenological data

Since the 1950s, the first appearance dates for A. mellifera and P. rapae have been recorded by volunteer observers in hundreds of localities throughout Spain. These observers applied the standardized protocolsproposed by the Spanish Intituto Nacional de Meteorología (for more details about this scheme, see Gordo and Sanz (2006a, b)). Both species are widespread in Spain and are well-studied due to their relevance in agriculture. A. mellifera is the main pollinator for orchards and most fruit trees, while the larval stages of P. rapae are important pests of cabbage crops. Adults of P. rapae also may act as pollinators of some entomophilous crops.

Records dated between 1952 and 2004 from original files of the Intituto Nacional de Meteorología were collected and digitized. A total of 7,263 records from 737 localities were gathered (see Figure 1 for details for each species). Each date was transformed into a Julian day scale (1 = 1 January). For leap years, one day was added after February 28 to take the extra day into account. Before performing any analyses, temporal trends of phenological data were removed by the regression of appearance dates against the quadratic function of the year for each species (Gordo and Sanz 2006b). Residuals obtained from temporal regression models were used as a measure of the phenological variation independent of the year from which these dates were recorded (see also Gordo et al. 2007a, b). Temporally corrected dates were used thereafter.

The mean appearance date for all records from the same 100 km² (10 × 10 km) universal transverse mercator (UTM) cell was calculated for each species. Because some localities of the phenological network were in the same UTM cell, the final sample size available for calculations (i.e. different UTM cells) was smaller than the number of original localities (see Figure 1). Mean values for each UTM cell could be biased due to differences in the number of records. This possible dependence was tested by calculating Spearman rank correlations between mean values and number of records in each UTM cell (A. mellifera: rs = 0.034, t617 = 0.845, P = 0.395; P. rapae: rs = -0.011, t440 = -0.231, P = 0.817). Because the mean appearance date was not dependent on the number of records, all UTM cells with available records for both species were used.

Explanatory variables A. mellifera

A total of 47 explanatory variables were classified into five categories (spatial, topographical, climatic, vegetation productivity and land uses; see Table 1) and used to model the appearance dates of A. mellifera and P. rapae throughout Spain. Seven topographic variables were obtained from a digital elevation model (Clark Labs 2000) for each of the 100 km² UTM Iberian squares (n = 6063) using the IDRISI 32 Geographic Information System (Clark Labs 2001a). The mean, minimum and maximum altitude of each 100 km² UTM cell was calculated from all 1 km² pixels included in each 100 km² UTM cell. The altitude range, slope, aspect (the mean direction of the slope) and diversity of aspects were also obtained for each 100 km² UTM cell. Delayed appearance dates are expected in UTMs at higher elevations and/or with northern exposures (Scott and Epstein 1987; Weiss et al. 1988, 1993; Gutiérrez and Menéndez 1998; Cocu et al. 2005; Gordo and Sanz 2006b; Harrington et al. 2007).

Eighteen climatic variables were provided by the Instituto Nacional de Meteorología for each of the 100 km² UTM Iberian squares. Climatic variables were rainfall and minimum, mean, and maximum temperatures during each season (spring, summer, autumn and winter), together with the annual temperature range and an aridity index. Spring was defined as April, May and June, summer as July, August and September, autumn as October, November and December, and winter as January, February and March for all seasonal variables. The aridity index was calculated as:

Vegetation productivity also was evaluated as a possible explanatory variable for spatial patterns of the spring appearances of A. mellifera and P. rapae in the Iberian Peninsula. This variable was measured as the normalized vegetation difference index (NDVI). The NDVI is the normalized difference between red (0.55 – 0.68 µm) and infrared (0.73 – 1.1 µm) reflectance, as measured by the National Oceanic and Atmospheric Administration's polar orbiting satellite's advanced very high resolution radiometer sensor (Smith et al. 1997). The NDVI is determined by the degree of red wavelength absorption by chlorophyll, which is proportional to leaf chlorophyll density, as well as by the reflectance of near infrared radiation, which is proportional to green leaf density (Tucker et al. 1985). Therefore, the NDVI correlates well with variables such as green leaf biomass, leaf area index, total accumulated dry matter and annual net primary productivity (Nicholson et al. 1990). NDVI data were available from Clark Labs world images as monthly values from 1982 to 2000 at a spatial resolution of 0.1 degree (Clark Labs 2001b). A value of vegetation productivity for each 100 km² UTM of the Iberian Peninsula was calculated by averaging monthly images available for each season between 1982–2000. More productive areas were expected to support higher abundance and diversity of organisms of upper trophic levels, such as insects (Hawkins and Porter 2003; Bailey et al. 2004; Seto et al. 2004). Larger populations of insects may enhance an early appearance simply due to increased chances for detection of early individuals (Tryjanowski et al. 2005; Dennis et al. 2006). Alternatively, larger populations may indeed be related to an earlier appearance because of the greater genetic diversity and the consequent higher probabilities for early phenotypes.

Land use types were also included because features of the environment eventually determine the presence of insect species (Eyre et al. 2003, 2005b). Therefore, habitat availability can be considered as another surrogate for the abundance of populations at a local scale, which, in turn, can affect detection (see above). Furthermore, land use has been demonstrated to have effects on the phenology of other insect taxa, such as aphids (Cocu et al. 2005; Harrington et al. 2007). Consequently, it would be of interest to understand to what extent this evidence is applicable to different taxonomic groups, such as bees or butterflies. The distribution of 15 land use types for the Iberian Peninsula was obtained from Corine Land Cover 2000 at a 100 × 100 m resolution. The percentage of coverage in each category within each 100 km² UTM cell was calculated and used as 15 explanatory variables for the analyses (see Table 1). The heterogeneity of land use types within each UTM was summarized by the Shannon diversity index and included as another explanatory variable.

Finally, spatial variables were used to verify the existence of spatial gradients. They were defined as the central latitude and longitude of each UTM cell and were included in the analyses as a third degree polynomial (Legendre and Legendre 1998). The nine terms of the spatial polynomial can help to incorporate effects of other historical, biotic or environmental variables not otherwise taken into consideration (Legendre and Legendre 1998). Latitude and longitude were standardized (mean = 0, and standard deviation =1) as were topographic, climatic, and vegetation productivity variables in order to eliminate their measurement scale effects.

Table 1.

List of explanatory variables used in analyses.

Statistical analyses

Multiple regression models implemented in the General Regression Models module of STATISTICA (StatSoft 2001) were conducted to determine the relationship between response (appearance date of A. mellifera and P. rapae) and explanatory (spatial, topographic, climatic, vegetation productivity and land use) variables. Models were built in three sequential steps. First, the relationship between insect appearance and each explanatory variable was explored one-by-one. For each predictor, linear, quadratic or cubic relationships were sought. The functions whose terms were statistically significant (p < 0.05) were selected. Only those environmental variables that were significantly related to appearance dates were included in further analyses. In the second step, the modelling ability of each type of explanatory variable (i.e. spatial, topographic, climatic, vegetation productivity and land uses) was explored by including significant variables belonging to the same category in a single regression model. A procedure of backward stepwise selection was applied in all cases to obtain simplified models that included only significant variables. The model with spatial variables allowed assessment of the spatial structure of phenological data, while models with environmental variables (i.e., the rest of the predictors) constituted the rationale for a biological interpretation of the spatial structure of dates that was identified. In the third and final analytical step, the best explanatory model of appearance dates was sought. For this purpose, two different regression models were carried out. One was with environmental explanatory variables and the other included environmental and spatial variables. The nine terms of the third degree spatial polynomial included in the latter model seek to account for the potential effect of other non-considered variables that are spatially structured. A backward selection procedure was applied in both models to include only significant (p < 0.05) predictors. Predicted scores of the best final regression model were mapped and examined.

Explanatory variables were correlated to each other (i.e., multicollinearity) due to the cofluctuation imposed by spatial and environmental gradients of the Iberian Peninsula (e.g., northern and high altitude regions have cooler and moister climates). This fact hinders the estimation of the true relevance of predictors (Quinn and Keough 2002). A hierarchical partitioning of variance procedure (Birks 1996; MacNally 2000, 2002) was implemented to determine the relative importance of each type of explanatory variable. The relative importance can be estimated as the average effect of including each type of variable in all possible models built with the remaining types of variable. Therefore, 2^k functions must be constructed for k types of explanatory variables. In our case, k = 5 (spatial, topographic, climatic, vegetation productivity and land uses) and thus 2^k = 32 different models.

Residuals from all multiple regression models were examined to check for spatial autocorrelation. If residuals are spatially autocorrelated, one or several important spatially structured explanatory variables are left out of the models (Cliff and Ord 1981; Legendre and Legendre 1998). Moran's I autocorrelation coefficient with a Bonferroni-corrected significance level (Sawada 1999) was calculated against ten classes separated by a lag distance of 60 km (from 60 to 600 km).

Apis mellifera

Many localities with an early phenology (especially those with very early appearance dates, i.e., prior to mid-February) were located in the southern and coastal areas of Spain (Figure 1a). In contrast, late sites for bee phenology occur mainly in central Iberia (i.e. Northern Plateau) and some mountainous regions (e.g., Iberian System; see Figure 2). Nevertheless, this row pattern is blurred by variability in the smaller scale. There were large differences between neighbouring UTMs. This could explain why the assessment of the spatial structure of data (i.e., the spatial model) showed a moderate explanatory ability (Table 2). The only variable included in the spatial model was the cubic function of the latitude. In agreement both with the prediction and the visual inspection of raw data, the spring appearance of A. mellifera is later in northern areas (c. 3 days·degree of latitude^-1).

Figure 1.

Maps of the distribution of phenological records of the appearance dates of Apis mellifera. Each square represents a UTM with phenological records, and its colour depicts the mean date. Scale color bar in Julian days (1 = 1 January). The total number of records, localities, and UTM together with the mean value and the standard deviation (SD) of all records are also specified. A histogram with the distribution of mean dates per UTM is also shown (scale of x-axis in Julian days). High quality figures are available online.

The topographical model included only the minimum altitude, the cubic function of which was strongly fitted to appearance dates (Table 2). A. mellifera appears later in those localities that are at higher elevations (Figure 3a). The effect was especially marked in those UTMs with a minimum altitude above 300 m.

Table 2.

Best multiple regression models for appearance dates of the honey bee and small white.

Figure 2.

Topographical map of the Iberian Peninsula. Darkness is proportional to the elevation. Main geographical features cited in the text are also shown. Solid black lines indicate political borders. High quality figures are available online.

The climatic model was the best environmental model (Table 2), although just one variable — the spring minimum temperature — was included. The paucity of variables in this model is remarkable because all climate variables (except autumn rainfall, see Table 1) were each significantly related to A. mellifera appearance dates. UTMs with the earliest bee appearance were located in those areas with the warmest springs (Figure 3b).

The model obtained with vegetation productivity variables was the least explanatory (Table 2). The quadratic function of the spring NDVI was able to capture just 5.55% of the variability in appearance dates. The appearance of A. mellifera was earlier in regions that are less productive during the spring.

Up to six variables were significantly included in the land-use model, although they accounted for a low proportion of the variability in appearance dates overall (Table 2). The negative sign of all variables means that A. mellifera were recorded earlier in the UTMs with large coverage of urbanized (i.e. humanized) areas, vineyards, fruit trees, olive groves, croplands and wild vegetation.

Figure 3.

Scatterplots of insect appearance against the best environmentally-related variables. Each point represents the average date of first appearance of Apis mellifera or Pieris rapae in each sampled UTM and its corresponding value of a selected environmental variable (see x-axis). A. mellifera appears later in those UTM cells with a higher minimum altitude (a) and earlier in those with warmer spring temperatures (b). The first P. rapae are sighted in the southern (c) and driest UTMs (d). Solid lines depict the best polynomial fitted model. High quality figures are available online.

The average relative explanatory capacity of each variable was notably lower than that obtained in previous models (Table 2). This could have been because of the high degree of collinearity among the explanatory variables. Nevertheless, the types of explanatory variables that showed the best fitted models also accounted for greater fractions of variability themselves. Interestingly, the variance hierarchical partitioning procedure revealed that the topographic model explained, on average, 2.4% more of the variability than the climatic model (Wilcoxon Matched Pairs test: Z₁₆ = 3.154, p = 0.002).

The complete environmental model included four variables and captured a bit more variability in appearance dates than the topographic or climatic models (Table 2). Minimum altitude and spring vegetation productivity had the same effect (i.e. sign) that was described above. The cubic function of the minimum altitude was the most relevant predictor of this model (partial R² = 22.81%). One climatic and one land use variable that were not included in their respective models were included in the final model. Although they did not show the best explanatory capacity within their groups of variables, they likely worked on a different part of the variability of dates. Consequently, summer precipitation and cover of dry farming were not overridden by altitude effects. According to the positive sign of summer precipitation and percentage of dry farming, localities with an early appearance of A. mellifera were related to regions with little precipitation during the summer and with a low proportion of dry farming. Spatial terms were not included in the complete model after the backward selection procedure, and, consequently, the model remained equal. Hence, no spatial structure remained in the phenological data after modelling with these environmental variables. This agrees with the absence of spatial autocorrelation in residuals of the complete model (Figure 4).

When predicted appearance dates from the complete model were mapped, the pattern of bee appearances in Spain was clearer than that offered by raw data (Figure 5). A. mellifera were observed for the first time on the southeastern Mediterranean coast of Spain. Other markedly early regions that were predicted included the southern regions and the Mediterranean coastal regions. Immediately later, A. mellifera appeared during the first 10 days of March in most of Spain. Finally, the late appearance of A. mellifera (from April 1st onwards) was expected in some central (Northern Plateau) and mountainous (Iberian System, Sierra Nevada or Pyrenees) regions of the Iberian Peninsula as a consequence of the high altitude.

Pieris rapae

A visual examination of raw data offered more evident spatial patterns for P. rapae appearances than observed for A. mellifera (Figure 6). Regarding the appearance of these insects, Spain can be unambiguously divided into two regions: one early area in the southern half of the Iberian and a late area in the northern half. Coastal areas both along the Mediterranean and along the Cantabrian seas also showed early appearance dates. Such patterns support a spatially well-structured P. rapae phenology. Indeed, all regression models for this species were better fitted (i.e., had a higher R²) than those obtained for A. mellifera (Table 2). This is especially striking because the explanatory variables included in the models were the same in most cases. Therefore, differences in the models' fitting between A. mellifera and P. rapae are likely not due to a different composition of the included explanatory variables. The spatial model showed a strong cubic latitudinal gradient (partial R² for latitude = 21.82%). The appearance of P. rapae was delayed almost one month from the southernmost point of Spain to the Cantabrian Mountains (Figure 3c). In the small latitudinal range from the Cantabrian Mountains to the Cantabrian coast, appearance dates tended to be earlier. Such a gradient was likely a result of the milder conditions due to the low altitude and to the proximity to the Atlantic Ocean.

Figure 4.

Spatial autocorrelation of model residuals. The isotropic correlogram represents the variation in the scores of Moran's I spatial autocorrelation statistic with the increasing separation distance between UTM cells, using a lag distance of 60 km and an active lag of 600 km. Spatial autocorrelation was absent at any distance in complete models (thick lines). High quality figures are available online.

Figure 5.

Maps of the predicted appearance dates for A. mellifera in Spain according to the best complete models. Scale color bar in Julian days (1 = 1 January).High quality figures are available online.

The minimum altitude also played an important role in the spatial variability of P. rapae phenology (Table 2). Those UTMs with higher minimum altitudes showed later appearance dates. The effect of minimum altitude was especially pronounced at higher altitudes, as was stressed by the quadratic relationship (see Figure 3a).

The climatic model was also the best model for P. rapae, although in this case, the explanatory capacity was markedly higher than the topographic model (Table 2). Temperature was the best explanatory variable. The linear relationship of appearance dates against mean spring temperatures accounted for 41.18% of the variability of P. rapae phenology. P. rapae is recorded earlier in areas that are warmer during the spring. The quadratic function of the summer rainfall was also included in the climatic model. This variable notably modeled appearance dates (R² = 26.2%; Figure 3d). The earliest appearances were linked to the driest areas. The early appearance in the moistest UTMs is linked to the latitudinal gradient. The moistest areas of the Iberian Peninsula are located on the Cantabrian coast, which has unusually mild conditions due to its proximity to the sea. The strong explanatory capacity of temperature and precipitation models did not add up as a result of the strong collinearity of both variables. Consequently, the full climatic model fit slightly more than previous models with a single climatic variable (Table 2). Nevertheless, the mean spring temperature was the most important variable in the climatic model (partial R² = 38.50%).

Vegetation productivity variables produced the worst explanatory model (Table 2). Spring patterns of productivity were related to appearance dates in a quadratic form. Early dates were linked to regions with low NDVI values during the spring.

Figure 6.

Maps of the distribution of phenological records of the appearance dates of Pieris rapae. Each square represents a UTM with phenological records, and its colour depicts the mean date. Scale color bar in Julian days (1 = 1 January). The total number of records, localities, and UTM together with the mean value and the standard deviation (SD) of all records are also specified. A histogram with the distribution of mean dates per UTM is also shown (scale of x-axis in Julian days). High quality figures are available online.

The land use model included a large number of variables (Table 2) and showed a moderate explanatory capacity of P. rapae phenology (Table 2). The first butterflies are recorded in areas with large coverage of olive groves and fruit trees, but with a low coverage of dry farming, coniferous forests, moorlands and transitional areas from shrublands to forests and grasslands. This result concurs with the spatial meaning of previous models. Fruit trees and olive groves are located mainly in southern Spain and on the Mediterranean coast. However, dry farming dominates landscapes from the Northern Plateau, and coniferous forests, moorlands and grasslands are typical for mountainous regions with high elevations, low temperatures and moist climates.

The average percentage of variability that accounted for each type of variable (see relative R²_adj in Table 2) was much lower, considering the high degree of collinearity among explanatory variables. Nevertheless, the relevance order was maintained. Topographic and climatic variables accounted for greater fractions of variability by themselves. Interestingly, true differences in the explanatory capacity of topographical and climatic variables were not too large, according to the relative R² values. On average, vegetation productivity and land use captured a negligible part of the variability.

The best complete model was built with climatic variables, minimum altitude and the coverage of dry farming (Table 2). The addition of minimum altitude and coverage of dry farming to climatic variables minimally increased the explanatory capacity previously achieved by the climatic model. The inclusion of spatial terms improved explanatory capacity of the model. This was due to the replacement of climatic variables with latitude. Residuals from both complete models (with and without spatial variables) did not show significant autocorrelation scores at any lag distance (Figure 4). Therefore, all spatially structured variations were included in P. rapae models.

Figure 7.

Maps of the predicted appearance dates for P. rapae in Spain according to the best complete models. Scale color bar in Julian days (1 = 1 January). High quality figures are available online.

Predicted date scores from the best complete model showed plainly differentiated phenological regions (Figure 5b). The earliest appearances of P. rapae were in the southern coast of Spain on both the Atlantic and the Mediterranean sides. However, the latest appearance dates were in the inner regions of the Northern Plateau and the Iberian System, as well as in the mountainous regions of the Pyrenees and Sierra Nevada. Intermediate regions also show a clear latitudinal gradient with later appearance dates in the Cantabrian coast and the Ebro valley than the dates in the Levant coast and southwestern Spain.