DOI: 10.31038/GEMS.2022435

Abstract

Artefacts in the Williston Reservoir, British Columbia were collected and recorded by an archaeological firm over several years. While a large and extensive dataset of the locations of these artefacts was built up, several challenges to management and interpretation of the use of the landscape are presented. This study used the Variable Clumping Method (VCM) hierarchical clustering technique to detect clusters of artefacts for each object type. These detected clusters were then used for an object-type spatial correlation analysis, using Multivariate Local Indicators of Spatial Association (LISA). The results from LISA analysis included detection of significant global and local associations between several object-type pairs. For instance, strong spatial correlation was found between Scrapers and Points, which may suggest significant use areas, such as regularly used butchering sites or campsites. While it is not possible to draw a definitive conclusion of exactly what these relationships mean in terms of landscape use, they suggest a number of interesting hypotheses of possible uses of the area in the ancient time.

Keywords

Artefacts, Archaeological materials, Hierarchical clustering, Variable clumping method, LISA

Introduction

Archaeological data is intrinsically spatial in nature, and spatial analysis can often play a key role in interpreting the distribution of archaeological materials across a landscape. This approach is particularly useful in a situation where the artefacts are scattered across the site without stratigraphic relationships or known features to provide the context. In particular, spatial patterning among the artefacts and the behavioural interpretation across the site require rigorous intrasite spatial analysis and modelling so as to deduce the kind of activities that occurred at the site [1]. A range of intrasite spatial analytical techniques have been applied to date [2], but the treatment of the spatial patterns of the artefacts and other entities occasionally presents challenges at the point of interpretation as their spatial tendency is not always correctly recognised [3]. This is mainly due to the direct application of either (1) a non-spatial statistical approach whereby all observations are treated as spatially independent in nature and that spatial randomness is warranted as the base distribution behind the spatial arrangement of the artefacts, or (2) a spatial statistical approach with the understanding that spatial dependency exists among the intrasite spatial configuration of the artefacts but the spatial patterning is often sought with respect to a predetermined set of scales informed by the specific context of the site—in other words, the spatial statistical methods are usually applied with parameters set a priori, despite that contexts may not be always available. The spatial arrangement of the artefacts found or excavated at a site is likely to have spatial dependency among them, and the scale of activities (e.g. an individual engaging in a house task, or a community taking part in a social activity) will likely vary between the types of activities and the participants involved which are difficult to calibrate without stratigraphic narratives or known features that offer the context. The objective of this study is to identify spatial associations between different types of archaeological artefacts found across an excavation site and gain knowledge on the spatial configuration and the lifestyle of the ancient community that lived in this area. To overcome the challenges stated above, this study uses two types of spatial analytical methods that can extract the hierarchical spatial structure: (1) a local variant of a spatial autocorrelation method called Multivariate Local Indicators of Spatial Association (LISA) [4], and (2) a hierarchical cluster detection method called Variable Clumping Method (VCM) [5,6]. Clustering helps to simplify a large archaeological dataset to make spatial patterns easier to discern; but they need to be arranged at suitable scales with the recognition of the multi-scale across different levels of activities. Finding clusters in the distribution of the artefacts across multiple scales would provide more natural groupings for use in the subsequent analysis. Exploration of the spatial relationships between different types of artefacts through LISA could provide us with a clue to infer activities that took place in an ancient time and could greatly enhance our understanding on how the landscape was used at the time. The study focuses on artefacts retrieved from Williston Reservoir, British Columbia, Canada.

Literature Review

Spatial Analysis in Archaeology

Spatial analysis and GIS have been used in archaeology for years, as evidenced in a study by Kintigh and Ammerman (1982) [7] which highlights the usefulness of computer-based analysis, informed by expert interpretation, for identifying meaningful patterns in archaeological data. According to Carr (1984) [1], there are two levels of archaeological spatial analysis: inferential and operational. The goals of inferential analysis are to spatially delineate the activity areas as well as to identify the tool kits of artefact types. At the operational level, the focus is placed on identifying patterns in the spatial arrangement of different artefact types, including clustering, and relationships between different artefact types. In many cases, applying inferential analysis may not be suitable. Even where the spatial relationships between artefact types can be identified, Carr (1984) [1] cautions that these do not imply an activity area and that there may be a number of other possible explanations for the observed patterns. Dynamic processes, including human activities, occurring in both space and time, have contributed to the superimposition of materials to produce static patterns [8]. For instance, artefacts may have been found clustered due to intense human activity in that location in the past, but Wood and Johnson (1978) suggest it could also result from natural processes such as water transport and sorting, or recent human activity such as farming. Deriving the nature of these processes from the static patterns can therefore be extremely difficult. Nevertheless, spatial analysis has proven as an important step to understand the nature of the archaeological deposits. A number of methods for detecting patterns in the spatial distribution of archaeological deposits have been proposed in the past and applied by a number of studies [3]. For instance, Kintigh and Ammerman (1982) [7] used the k-means cluster detection technique to iteratively split and lump points in the dataset into a user-defined maximum number of clusters and, thereby, minimise the sum of the squares of the distances between each point and its cluster centre to find the optimal solution. More generally, Carr (1984) [1] refers to a number of methods used in intrasite archaeological contexts for determining presence of clustering as well as delimiting clusters in a point dataset. These include the Nearest Neighbour Analysis [9] used for determining presence of clustering; Whallon’s Radius Approach [10] for delimiting clusters using the frequency distribution of nearest neighbour distances between point observations; and conditional spatial patterning, a multi-scale method to test for spatial patterning in an archaeological dataset [8].Other applications of these methods vary in location, site type and the scale. For instance, in a cave site in Western Belize, artefacts scattered on different surfaces, including in niches and on ledges, were mapped and recorded in a GIS to analyse the spatial arrangement of the archaeological objects [11]. They adopted a clustering method, as some artefacts were found in multiple fragments in close proximity to one another and that the uses the individual object locations for analysis would have resulted in unequal weighting. A non-hierarchical method, k-means, was used to define natural groupings of objects in space, and the resultant clusters were analysed with regard to their association with cave features. Moyes’ (2002) [11] findings suggest that clustering and comparative analysis of the morphologic cave features could provide the context and the relative analytical units for use on a broader scale.

The examples above, however, are situations where the sites and the archaeological context are well defined. Many of these methods primarily deal with “more or less intact living surfaces” [12] that may not be applicable to other archaeological situations that lack an intact or a clear definition of the site. An example of such a situation is found in Australia; similar to the Williston Reservoir case, where the archaeology consists of open sites [13]. These open sites are extensive in area size, lack stratigraphy or clear boundaries, and contain few features [13]. In essence, it is a vast expanse of a surface scattered with primarily lithic artefacts. Site boundaries are difficult to determine in this type of landscape and, for the management purposes, arbitrary rules are often used to define boundaries that are not meaningful for analysis [13]. In these situations, a site-less analysis is required; i.e. the analysis is based on individual artefacts rather than sites [13]. As a result, the first step with open sites is to identify patterns in the spatial location of artefacts; after which, attempts can be made to interpret the patterns archaeologically. Holdaway et al. (1998) [13] address these concerns through geomorphological mapping of the landscape used for understanding the patterns of surface artefact density. This step determined that the depositional surfaces had lower artefact densities, which were excluded from the subsequent analysis to avoid bias in case undiscovered artefacts were buried under the sediments. Different buffer radii were used to look at different scales, and the resultant cluster patterns of individual artefact types were then used to identify assemblages and patterns between artefact types [13]. Other studies have also explored different methods for determining the associations between artefact types. These include a variety of statistical tests based on the number of artefact in partitioned units to determine patterns of aggregation and segregation between different artefact classes [14]. In this case, the definition of the unit shape and size may impact the results of the analysis. This is known as Modifiable Areal Unit Problem (MAUP) in spatial analysis [15,16]. Hietala and Stevens (1977) [14] therefore suggest changing the partition size to determine patterns for multiple partition options. Berry et al. (1980) [17] suggest the permutation test, which can use either grid-count values or point locations, and so can avoid the problems associated with defining a grid unit size. The permutation test uses as its test statistic the “average within-class distance” (Berry et al., 1980, p. 56) [17], and thus can test for significant associations between multiple classes. These methods detect global associations between artefact classes; i.e. they assess the overall density or the tendency of clustered-ness across the study area. However, as Premo (2004) [18] points out, being able to quantify the local spatial patterns (which enables us to identify the location, the extent and the intensity of each individual cluster) is important for the understanding and interpretation of archaeological material distributions. Methods for local spatial pattern detection used by Premo (2004) [18] include local Moran’s I and local G statistics, mainly for the purpose of detecting local spatial autocorrelation. These local statistics have a great potential for multiple archaeological applications, including that of multivariate analysis for identifying the association between artefact types or material types Premo (2004) [18].

Hierarchical Clustering Methods in Spatial Analysis

With the development of data collection technology such as Global Positioning Systems (GPS), remote sensing, and vast amount of spatial data is becoming increasingly available. Yet their interpretation is not always straightforward, and dealing with the data to extract meaningful information can be difficult at times. Exploratory spatial data analysis (ESDA) offers an important first step for deriving information from large spatial datasets [19]. Cluster analysis is a spatial extension of ESDA used across a broad range of disciplines, including crime analysis, disease analysis, and archaeology [6,7,20]. Cluster detection techniques for point datasets include techniques that simply determine the presence of clustering, such as nearest neighbour indices [21], as well as techniques that identify the individual points within a cluster [6]. Early clustering techniques include the Geographical Analysis Machine [22] for identifying hot spots or areas of high intensity; while more recent techniques are designed to identify which points belong to a particular cluster—some of these techniques can be classified further into partitioning, hierarchical, and graph-based techniques [21,23]. Partitioning methods, such as k-means, group all points in the dataset into a user-defined number of groups [21]. In addition to the disadvantage of having to specify a number of groups, which can lead to bias in the analysis, partitioning techniques are unable to identify cluster shapes that are not convex [24,25]. Hierarchical methods are typically either top-down or bottom-up; for top-down methods, points are either grouped into a single cluster, then that cluster split according to some function to create two clusters, those clusters further split, and so on [26]. Bottom-up methods, such as the nearest neighbour hierarchical clustering technique use some function to group the individual points into a number of clusters, then proceed to group the primary-level clusters into secondary-level clusters and so on, until there is a single cluster. One advantage of hierarchical methods is that the user does not specify how many clusters to generate; however, user-defined criterion are required to tell the software when to stop clustering, or to define the initial clustering criteria [21,24,26]. A limitation is that each level of clustering depends on the previous level. Graph-based techniques are those that compute a graph, where the points are vertices and edges are lines connecting pairs of points, with edge lengths representing the proximity of pairs of points [26]. The minimum spanning tree (MST) and Delaunay diagram are examples of graphs used in clustering algorithms. Some cluster methods falling into this category include AMOEBA [23], AUTOCLUST [26], and VCM [5]. AMOEBA and AUTOCLUST are similar techniques based on the Delaunay Diagram. AMOEBA uses the global mean and standard deviation of all edge lengths in the graph, compared to the local mean of all edges connected to a single point, to determine a tolerance value. Edge lengths exceeding this tolerance value are removed from the graph such that the remaining connected points form the clusters. The algorithm is then reiterated to generate sub-clusters from the primary clusters, and so on until no more edges are present in the graph, producing hierarchical clusters. AMOEBA detects clusters of different density, and also non-convex clusters [23]. AUTOCLUST is like AMOEBA, but it compares local mean and standard deviation of edge lengths for a point to the average local standard deviation of all points to determine the tolerance [26]. Like AMOEBA, it succeeds in identifying clusters of different density and arbitrary shape [26]. Other similar algorithms for finding clusters have been presented by various authors [25,27].

Variable clumping method (VCM) is a hierarchical and graph-based method that uses a minimum spanning tree (MST) as the graph of edges [5]. It detects “clumps” of points at varying distances by iterating through the ordered (by length) set of edges in the MST and using each length in turn as the radius of circles centered on the points, so that a set of points within connected circles constitute a clump [5]. The use of the variable radii enables detection of clumps at multiple scales. VCM also uses Monte Carlo simulation to determine which clumps at each radius are significant, and only includes these significant clumps in the final set [5,6]. All these clustering techniques offer many advantages but they also have limitations. Selecting an appropriate clustering technique requires defining the needs of the particular analysis to determine which limitations are acceptable. While there are hierarchical clustering approaches that would serve the needs of this analysis (e.g. AMOEBA, AUTOCLUST), they are mostly unavailable as a ready-to-use software package. CrimeStat is an exception, as a freely available software package that includes a nearest neighbour hierarchical clustering routine [21]. However, preliminary tests with this routine revealed several undesirable qualities such as the need to input parameters of the minimum number of points to include in the cluster, as well as the threshold distance [21], but these parameters may be difficult to estimate in a case where no context or the extent is available. Also, after running the clustering routine on a sample dataset (a subset of the Williston artefacts dataset), some of the member points within a cluster turned out to be closer to points in another cluster, which violates the primary goal of clustering as described above. Finally, each subsequent cluster level only clustered the clusters, so that outliers from the primary-level clustering remained isolated in the secondary-level clustering and so on. These attributes made it unsuitable for this study. To overcome these challenges, this study adopts the VCM approach. While it does require a user-defined parameter, this is only needed for specifying the number of hierarchical levels to generate, which is constrained by the amount of complexity the analysis can handle (with more levels, the complexity increases). As detailed in the discussion below, it ensures that the distances between points within a cluster are minimised, through the use of natural breaks, and that any point within the cluster is closer to other points within the cluster than it is to points outside the cluster. Also, each cluster level is based on the full set of points, so outliers from the first cluster level are incorporated into higher-order clusters.

Methodology

Study Area

The study area of Williston Reservoir is located in northern British Columbia, Canada. It was created by the construction of the WAC Bennett Dam in the late 1960s, and is one of the largest reservoirs in the world [28]. Reservoir operation over the last 50 years has led to the exposure of primarily unvegetated large expanses of “beach”, made up of fine silts, clays, gravels, or sand. This area was heavily used by the First Nations peoples throughout history, and physical evidence of this use remains on the landscape in the form of exposed surface lithic artefacts and other cultural remains [29,30]. Millennia Research Limited, an archaeological consulting firm, has conducted archaeological surveys of the inundation zone of Williston Reservoir annually since 2008. The surveyed area consists of discrete beaches, which are confined and delineated by the natural boundaries, usually in the forms of large creeks or rivers. Artefacts found during these surveys were recorded using handheld Global Positioning Systems (GPS) units, with an estimated average accuracy of ±5 m. The reservoir is composed of three main “reaches”, the northernmost of which is Finlay Reach, where the majority of the archaeological work has been done. Several of the core beaches of Finlay Reach (those that are most densely populated with artefacts) form the study area for this project (Figure 1). The archaeologically surveyed portion of these beaches totals over 68 km², in which over 6,000 artefacts have been recorded for 2009-2011 alone. The environment of the study area and the characteristics of the archaeological remains present some difficulties for the management of the archaeological resource. In British Columbia, archaeological resources are protected by law, and the provincial Archaeology Branch maintains a registry of the archaeological sites [31]. Thus, for management purposes, a site definition is required. In this environment, with vast areas covered by scatters of artefacts without clearly defined features, defining the boundaries of the archaeological site can be challenging. Furthermore, the management decisions for archaeological resources are often dependent on their significance. As highlighted by Glassow (1977) [32], significance can be difficult to define. This is even more so if the unit of analysis is unclear because the extent and the boundary of the archaeological site is ambiguous. Ideally, the unit of analysis should be extracted from non-arbitrary grouping of associated archaeological materials, yet identifying the categories for such grouping can be very difficult, because of the sheer volume of the Williston Reservoir artefact dataset, which makes visual identification of patterns impractical. There are other challenges such as the lack of stratigraphic context, defined landforms, or geomorphological features; as well as the lack of intact archaeological features. In addition, the landscape is not static, which means that new artefacts are discovered each year in areas that had been previously surveyed. All of these factors pose challenges for interpretation, analysis and management of this vast archaeological resource.

Figure 1: Study area beaches and artefact locations

Hierarchical Clustering

VCM is a spatial analytical method for detecting statistically significant multi-level clumps [5]. Circles of a variable-sized radius are drawn around each observed point, and any connected set of circles represents a clump. Clumps have some key properties including the clump radius (r) and the clump size (k), i.e. the number of connected circles. The “variable” part of VCM comes from varying the radius, so that the clumps identified at varying radii, and thus multi-scale clumps, are identified. The clumping state at a specific radius C(r) is defined by the number of clumps of size k at radius r, or N(k|r), so that C(r) = (N(2|r), N(3|r), … , N(n|r)) [6]. Note that a clump of size one is not considered a proper clump [5]. A minimum spanning tree (MST) represents the distances between points, using the property of edge length (l), so that at radius r, any points connected by an edge with length l ≤ r form a clump. Because clumps will appear even in a random distribution of points, VCM also conducts a significance analysis to determine which clumps are significant [5,6]. A radius interval and maximum radius are specified to define a set of radii for the analysis. The set of clumping states of the observed points for this set of radii is then determined. Next, 10,000 simulations, using randomly generated point distributions and finding the MST for these point distributions, provide a frequency distribution of clumping states for the set of radii. From this distribution, using a significance level (α) of α = 0.05, the critical number for each clump size at each radius is determined from this frequency distribution. The null hypothesis (H₀) states that the number of clumps of observed clump size k at a radius r will be less than or equal to the critical value for k and r. The results from the observed data are then compared to these critical values, and H₀ is rejected where the observed number of clumps of size k at radius r is greater than the critical value [5,6]. A clustering method such as VCM is based on the MST of the points in the dataset. However, this study used a classification method to limit the number of distances at which to generate clusters. Individual line segments of the MST were classified by length using a natural breaks classification, to determine the distance thresholds at which each cluster level was defined. Thus, the number of classes chosen defines the number of levels in the resulting hierarchy. Figure 2 illustrates this process. In Figure 2a, the MST line segments have been classified into five classes. Note that the “Class 5” line is dashed – this is to indicate that this class is not used as a cluster level, as it would include all of the points in the dataset. Figure 2b shows the convex hulls of the resultant clusters. The process is cumulative, so that “Cluster level 2” includes all points connected by “Class 1” as well as “Class 2” line lengths.

Figure 2: An illustrative example of hierarchical clusters; (a) classified MST lines connecting artefact point locations, (b) hierarchical clusters shown as convex hulls surrounding the original artefact point locations.

The specific methodology used for the dataset in this study was to first separate the overall artefact dataset into individual datasets by beach. This was done so that the resulting clusters would reflect the distribution of artefacts on a particular beach. Some of the beaches are densely populated with artefacts, while others are more sparsely populated. As the beaches are well-defined by major landforms (creeks and rivers), it is sensible that they should be treated individually. An MST was generated for each dataset, and then classified using the Natural Breaks classification method, which is based on Jenks’ optimisation method. In essence, this clustering problem is very similar to the choropleth mapping problem discussed by Jenks (1967) [33] in that, while it may seem ideal to present every data value, or in the case of hierarchical clustering, every possible cluster level, a limited number of classes must be used in order to be able to understand and interpret the data. In that case, it is desirable that each class should contain very similar values, so that the within-class deviation is minimised.

Step 1: Select all of the primary-level cluster lines (i.e. lines that have a length less than or equal to the Class 1 breakpoint).

Step 2: Buffer the selected lines.

Step 3: Dissolve any overlapping buffer polygons.

Step 4: Spatially join the original points layer to the buffer polygons, so that they are assigned the buffer polygon ID, which becomes the “Cluster 1 ID”.

Step 5: Generate convex hulls, grouping by Cluster 1 ID.

Step 6: Repeat above steps for each cluster level, such that the each successively higher cluster level includes all line lengths smaller than the break value (i.e., each cluster includes points from lower-order clusters).

The resulting datasets included a point dataset with all of the artefacts assigned an ID value for each cluster level, as well as convex hulls of all of the clusters at each level, and a size one standard deviation ellipse for each cluster at each level.

Bivariate Autocorrelation Analysis for Object Type Relationships

Bivariate local Moran’s statistical analysis [34] is often used for describing the spatial correlation between the spatial distribution patterns of two variables; i.e. it is a local method that identifies the actual locations where significant spatial dependency was observed between the two variables. GeoDa software [35] was used to perform a multivariate LISA analysis on pairs of object types in order to determine if any statistically significant relationships exist between them locally or globally. Anselin (1995) [4] defines a LISA as any statistic that indicates the degree and significance of spatial clustering of similar values around each observation. In addition, the sum of LISA statistics for a set of observations must be proportional to a global statistic [4]. GeoDa uses the Moran’s I statistic for LISA analysis. The global Moran’s I value is an indicator of overall clustering within the dataset, while the local Moran’s I value indicates the locations of clusters [36]. A significant local association may occur that is not globally significant, or there may be patterns occurring locally that are opposite to the global trend [4]. When applied as a multivariate test of spatial correlation, the statistic compares values for one variable at a location to values for a second variable at neighbouring locations. Values are “standardised such that the mean is zero and standard deviation equals one” [33]. The standardised values at each observation location are compared to the spatially lagged, standardised values at neighbouring locations to produce a global multivariate Moran’s I [34]. The contributions of individual observations to this global value are also calculated to determine local multivariate Moran’s I statistics. These values can then be compared to the values expected under a scenario of complete spatial randomness in order to determine the significance of the relationship between the two variables, both globally and locally [34]. This is done by calculating the Moran’s I for a large number of randomised permutations, in which one of the variables is kept static, while the other is randomly reallocated amongst the observations. Running several thousand random permutations produces an indicator of how extreme, and therefore how significant, the observed values are [4]. Thus, within the GeoDa multivariate LISA analysis, it is possible to obtain both global and local indicators of significant spatial association between two different variables. For each pair of variables a pseudo-significance level for the calculated local Moran’s I statistic was determined using 9,999 randomised permutations. The results, including the local Moran’s I value for each cluster with neighbours, the spatial association type, and the significance p-value for the association, were saved to a table.

The outcome of the significance test of LISA classifies individual cluster location into four different categories:

High-High

If the LISA statistic is statistically significant, takes a positive value, and the standardised count/value of object type A is positive, then both the object type A and the object type B are significantly high.

Low-Low

If the LISA statistic is statistically significant, takes a positive value, and the standardised count/value of object type A is negative, then both the object type A and the object type B are significantly low.

Low-High

If the LISA statistic is statistically significant, takes a negative value, and the standardised count/value of object type A is negative, then the object type A is low but the object type B is high.

High-Low

If the LISA statistic is statistically significant, takes a negative values, and the standardised count/value of object type A is positive, then the object type A is high but the object type B is low.

Analysis

In the analysis, a two-tier cluster levels were used for testing the spatial patterns of artefacts at different scales. The primary-level clusters were combined into a single dataset, and the secondary-level clusters were combined into another dataset. The object types that were used in the analysis are listed in Table 1. Some object type categories had very few members to the extent that they would not sustain robust analysis and were therefore excluded from the analysis. One very large category of object types, Flake Debitage, was also excluded because these items were often recorded only cursorily, making the data for this object type unreliable for use in this analysis.

Table 1: List of object types used in the analysis, with number of clusters for each object type

GeoDa was used for generating a binary spatial weights matrix, in which a distance cut-off was applied to determine whether each cluster is considered a neighbour of another polygon. Clusters falling within the distance band are counted as neighbours in the analysis, and those falling outside the distance band are not considered [35]. GeoDa automatically calculates a distance that ensures that all clusters have at least one neighbour. However, the default distance turned out to be too large to provide archaeologically meaningful results, and it was adjusted through an exploratory process to calibrate it for localised analysis that also retains sufficient neighbours for most locations. Through this process, the primary-level clusters were set with the distance threshold of 250m, and the secondary-level clusters at 500 m. The resulting weights matrix for primary-level clusters had some isolated, neighbourless locations (Figure 3a). For the secondary-level clusters, every location had at least one neighbour (Figure 3b). This second level of analysis provides a more regional view of the relationships between object types, while the first level provides a more localised view.

Figure 3: Detected artefacts clusters (a) distribution of primary-level clusters with the omission of the neighbourless clusters, (b) distribution of secondary-level clusters.

Of the 20 different relationships between object types tested at the first cluster level, three global significant relationships were discovered, and significant local associations were present in all of the comparisons. Amongst the statistically significant global patterns which emerged from the analysis was a positive correlation between Alberta Points and Macroblades. The negative correlations observed between Impact Fractured Points and Scrapers, and Impact Fractured Points and Flake Tools were statistically significant. The secondary-level cluster analysis did not result in many global significant associations; however, a positive correlation between Scraper object types and Point object types was significant. The secondary-cluster level results provide a more regional summary of the relationships between object types. None of the significant primary-level global associations were present as global associations for secondary-level clusters. Table 2 summaries the significant local relationships for primary-level clusters by their spatial association types, and Table 3 offers the same for secondary-level clusters. The association types are perhaps the most informative with regards to the nature of the relationships, as the actual Moran’s I value is not necessarily easy to interpret on its own.

Table 2: Significant results from primary-level multivariate local Moran’s I analysis

Table 3: Significant results from secondary-level multivariate local Moran’s I analysis

Maps of selected variable pair results display some interesting trends (Figures 4-7). For example, Figure 4a shows the association between Flake Tool and Scraper object types for primary-level clusters. While there are many High-High associations (as seen in Table 2), most of these are grouped in two locations at the north end of the study area. However, when the secondary-level cluster map for this association is compared (Figure 4b), these groupings have disappeared. Instead, new locations of High-High associations appear at this different scale. For the association between Point and Scraper (Figure 5), a similar pattern shows, with most of the High-High associations for primary-level clusters (Figure 5a) grouped in two locations at the north end of the study area. Note that the southern of these two High-High groups for Point and Scraper is in the same location as one of the Flake Tool to Scraper primary-level cluster groups. Again, these significant associations disappear at the secondary-cluster level (Figure 5b). At the first cluster level, Battered Bifaces and Hammerstones do not tend to occur together (Figure 6a), except for a single High-High cluster. At the secondary-cluster level, this High-High association has disappeared, and a new High-High association location has appeared, towards the northern portion of the study area (Figure 6b). It is interesting to note that this also occurs in the same location as one of the High-High groupings seen in both Figures 4a and 5a. Elsewhere in the reservoir, the general non-association trend between Battered Bifaces and Hammerstones holds at the secondary-cluster level. Other trends of interest include the globally significant positive correlation between Alberta Points and Macroblades. When the local associations are viewed on the map (Figure 7a), all of the High-High associations are grouped in one area, towards the center of the study area. The primary-level global negative correlation between Impact Fractured Points and Flake Tools is also visible in the local association types, as there are no High-High associations (Figure 7b). At the secondary-cluster level, the global association between Scraper and Point object types is apparent as several groups of High-High associations (Figure 7c). However, there are also several significant High-Low locations, most notably the very large group of these two association types which occurs in the northern part of the study area.

Figure 4: Significant associations between Flake Tool and Scraper object types for (a) primary-level clusters, and (b) secondary-level clusters (insignificant associations and neighbourless clusters omitted).

Figure 5: Significant associations between Point and Scraper object types for (a) primary-level clusters, and (b) secondary-level clusters.

Figure 6: Significant associations between Battered Biface and Hammerstone object types for (a) primary-level clusters, and (b) secondary-level clusters.

Figure 7: Significant associations between (a) Alberta Point and Macroblade object types for primary-level clusters, (b) Impact Fractured Point and Flake Tool object types for primary-level clusters, and (c) Scraper and Point object types for secondary-level clusters.

Discussion

Key results of the LISA analysis include the four global significant associations, between Alberta Points and Macroblades, Impact Fractured Points and Flake Tools as well as Scrapers, and Scrapers to Points. A number of interesting local patterns were also discovered, including the suggested High-High groupings seen at the same location for multiple object-pairs (Figures 4-6), and scale differences in these local patterns as highlighted by the use of two different levels of clustering. The significant relationships may have a number of explanations. The negative correlations between Impact Fractured Points and Scrapers/Flake Tools suggest that the Impact Fractured Points are hunting losses, rather than retrieved from a kill and then discarded at a campsite or butchering site. If the latter were the case, they would be expected to be found more frequently in association with Flake Tools and Scrapers, both of which object types would typically be found at camp or butchering locations. Furthermore, the Impact Fractured Point to Scraper results for secondary-level clusters. None of the clusters have a significant local High-High relationship at this level and, while not globally significant, this seems to follow the trend of the primary-level cluster results, thus implying that this relationship may persist at this scale across the site. The global relationship between Scrapers and Points at the secondary-cluster level was composed of several distinct groupings of local High-High associations. These may suggest significant use areas, such as regularly used butchering sites or campsites. However, perhaps of equal or greater interest is the area with a large grouping of High-Low significant associations, which oppose the global trend of positive correlation, suggesting that something quite different may be occurring in this particular region. The Alberta Point to Macroblade global association could indicate that these items were used contemporaneously, thus providing a possible scenario of temporal patterning. Though this association is not significant at the secondary-cluster level, even at this scale there are no clusters with high numbers of Alberta Points and low numbers of Macroblades and one second level cluster retains a significant local High-High association, suggesting that Alberta Points are not typically found far from Macroblades. Other local relationships observed have indicated some interesting trends that may have a number of explanations. The groups of High-High clusters seen in the same locations for several of the variable pairs (Figures 4-7) suggest that these are not likely to be randomly strewn artefacts but may be an indicator of some important activity areas, such as camp locations, though further research would be required to confirm this hypothesis. Hierarchical clustering has proven to be an effective tool for comparing the change in clusters over time, as well as the discovery of associations between different artefact types across multiple scales. The object type analysis revealed significant patterns in the association between different artefact types at multiple scales, and while it is not possible to draw a definitive conclusion of exactly what these relationships mean in terms of landscape use, they suggest a number of interesting hypotheses of possible uses and provide direction for further studies.

References

1. Carr C (1984) The Nature of Organization of Intrasite Archaeological Records and Spatial Analytic Approaches to Their Investigation. Advances in Archaeological Method and Theory 7: 103-222.
2. Bevan A, Conolly J (2009). Modelling spatial heterogeneity and nonstationarity in artefact-rich landscapes. Journal of Archaeological Science 36: 956-964.
3. Gillings M, Hacıgüzeller P, Lock G (2020). Archaeological Spatial Analysis: A Methodological Guide, Routeledge.
4. Anselin L (1995) Local Indicators of Spatial Association—LISA. Geographical Analysis 27: 93-115.
5. Okabe A, Funamoto S (2000) An exploratory method for detecting multi-level clumps in the distribution of points – a computational too, VCM (variable clumping method). Jouranl of Geographical Systems 2: 111-120.
6. Shiode S, Shiode N (2009) Detection of multi-scale cluster in network space. International Journal of Geographic Information Science 23: 75-92.
7. Kintigh KW, Ammerman AJ (1982) Heuristic approaches to spatial analysis in archaeology. American Antiquity 47: 31-63.
8. Voorrips A, O’Shea JM (1987) Conditional spatial patterning: Beyond the nearest neighbour. American Antiquity, 52: 500-521.
9. Clark PJ, Evans FC (1954) Distance to nearest neighbor as a measure of spatial relationships in populations. Ecology 35: 445-453.
10. Whallon R (1972) A new approach to pottery typology. American Antiquity 37: 13-33.
11. Moyes H (2002) The use of GIS in the spatial analysis of an archaeological cave site. Journal of Cave and Karst Studies 64: 9-16.
12. Kintigh KW (1990) Intrasite spatial analysis: A commentary on major methods. In Mathematics and Information Science in Archaeology: A Flexible Framework, edited by Albertus Voorrips. Studies in Modern Archaeology 3: 165-200.
13. Holdaway S, Witter D, Fanning P, Musgrave R, Cochrane G, et al. (1998). New approaches to open site spatial archaeology in Sturt National Park, New South Wales, Australia. Archaeology in Oceania 33: 1-19.
14. Hietala HJ, Stevens DE (1977) Spatial Analysis: Multiple Procedures in Pattern Recognition Studies. American Antiquity 42: 539-559.
15. Openshaw S (1984) The Modifiable Areal Unit Problem. Concepts and Techniques in Modern Geography, No.38, Norwick: Geo Books.
16. Buzzelli M (2020) Modifiable areal unit problem. In Kobayashi, A. (ed). International Encyclopedia of Human Geography, 2^nd, Elsevier, pg: 169-173.
17. Berry KJ, Kvamme KL, Mielke PW Jr. (1980) A Permutation Technique for the Spatial Analysis of the Distribution of Artefacts into Classes. American Antiquity, 45(1), 55-59.
18. Premo LS (2004) Local spatial autocorrelation statistics quantify multi-scale patterns in distributional data: an example from the Maya Lowlands. Journal of Archaeological Science 31: 855-866.
19. Murray AT, Estivill-Castro V (1998) Cluster discovery techniques for exploratory spatial data analysis. International Journal of Geographical Information Science 12: 431-443.
20. Wheeler DC (2007) A comparison of spatial clustering and cluster detection techniques for childhood leukemia incidence in Ohio, 1996-2003. International Journal of Health Geographics 6. [crossref]
21. Levine N (2004) CrimeStat III: A Spatial Statistics Program for the Analysis of Crime Incident Locations (Version 3.3): Ned Levine & Associates, Houston, TX, and the National Institute of Justice, Washington, DC.
22. Openshaw S, Charlton M, Wymer C, Craft A (1987) A Mark 1 Geographical Analysis Machine for the automated analysis of point data sets. International Journal of Geographical Information Systems 1: 335-358.
23. Estivill-Castro V, Lee I (2000a) AMOEBA: Hierarchical clustering based on spatial proximity using Delaunay diagram. Paper presented at the 9th International Symposium on Spatial Data Handling, Beijing, China.
24. Qu Y, Lee K, Lee I (2011) Cluster polygonization and qualitative cluster reasoning overview. Interntational Journal of Advancements in Computing Technology 3.
25. Yang X, He L, Lu H (2009) A New Method for Non-spherical and Multi-density Clustering. Third International Symposium on Intelligent Information Technology Application, 2009, pg: 35-38.
26. Estivill-Castro V, Lee I (2000b) AUTOCLUST: Automatic Clustering via Boundary Extraction for Mining Massive Point-Data Sets. Paper presented at the 5th International Conference on Geocomputation, University of Greenwich, United Kingdom.
27. Malik M, Singh P, Sharma DAK (2011) A novel spatial clustering approach for outlier detection & cluster generation by probing various distance matrices & Delaunay triangulation. International Journal of Computer Science and Technology 2: 104-108.
28. BC Hydro (2012) Williston Reservoir Retrieved August 24, 2012.
29. Eldridge M, Brunsden J, Parker A, Eldridge R (2008) Permit 2008-0179 BC Hydro 2008 Williston Dust Abatement Project Archaeological Impact Assessment Final Report: Millennia Research Limited.
30. Eldridge M, Thiesson V, Parker A, Eldridge R (2010) BC Hydro 2010 Williston Dust Abatement Project Archaeological Impact Assessment Final Report: Millennia Research Limited.
31. The Province of British Columbia (2022) Ministry of Forests, Lands, and Natural Resource Operations: Archaeology.
32. Glassow MA (1977) Issues in evaluating the significance of archaeological resources. American Antiquity, 42: 413-420.
33. Jenks GF (1967) The data model concept in statistical mapping. International Yearbook of Cartography 7: 186-190.
34. Anselin L, Syabri I, Smirnov O (2002) Visualizing Multivariate Spatial Correlation with Dynamically Linked Windows. New Tools for Spatial Data Analysis: Proceedings of the Specialist Meeting.
35. Anselin L (2022) GeoDa (Version 1.20): GeoDa Center for Geospatial Analysis and Computation.
36. Anselin L, Syabri I, Kho Y (2006) GeoDa: An Introduction to Spatial Data Analysis. Geographical Analysis 38: 5-22.

DOI: 10.31038/GEMS.2022434

Abstract

Impact of climate change is felt across a wide range of events on earth surface, including the geospatial movements by insects. This research aims to measure the impact of rising temperature on the Bombus (Bumble bee) population. The study firstly divides an 11-year accumulation of Bombus spotting data from Oxfordshire, United Kingdom into four categories that reflect the respective stages in the Bombus life cycle. Using regression analysis, the study investigates how the temperature has influenced the bee population at each stage. Results from the regression analysis revealed that the rise in temperature is affecting the population size of the Bombus population in all stages of their life cycle, particularly at an early stage where queen bees and worker bees emerge. Spatial analysis using the mean centre and the deviational ellipse also revealed that the queen bees are migrating generally towards north-west, and are becoming increasingly dispersed, expanding continuously in the direction of north-west to south-east.

Keywords

Bombus, Temperature, Phenology, Migration, Regression analysis

Introduction

A French botanist Charles Morren first defined the term phenology as “the science of theses sort of things” in 1833 [1]. The term has since evolved and, in modern research, it is used as a reference to the seasonal natural phenomena such as the life cycles of plant, animal and insect species in relation to seasonal changes. A main indicator of seasonal change is the change in temperature. Like many countries in the world, the United Kingdom is currently experiencing rise in the temperature as an effect of global warming with an accelerated pace in the past decade, and it is likely that an over-forty degrees Celsius temperature in the summer would become a norm within this century [2]. It would mean that the summer in the UK would not be dissimilar from that of the current state of Sahara Desert [3], something we did not previously anticipate in the region north of the temperate zone [4]. This situation makes the phenological life cycles particularly volatile. Bombus bumblebees are part of the Apidae (bee) family, who are mostly known for their furry bodies and large wings. They are mostly found in the temperate climates with seventy-nine species found in Europe [5]. They play a crucial role in the pollination of crops, often filling in the gaps for those that honeybees cannot cover [6]. They live in colonies which are led by one female referred to as the ‘queen.’ The queen is the only survivor from the previous season’s colony, and she emerges out of hibernation in early spring in search of food to begin the production of her nest. The locations of these nests vary with settlings occurring under thick vegetation in areas of grassland (O’Conner et al. 2014) or rodent holes and abandoned birds’ nests in urban environments [7]. Once a habitat is established, she lays her first set of eggs which would become the first batch of worker bees. These worker bees are responsible for the maintenance and sustenance of the nest by foraging throughout the season as well as rearing the constant new broods that emerge. The timing of a worker bees’ emergence is known to be crucial as this will indicate the queen’s physiology for the rest of the season [8]. Once late summer to early autumn reaches, the queen begins to give birth to both the male and female offspring, respectively. When the males have sexually matured, they permanently exit from the nest to seek female bees from other colonies—each batch of male bees exiting a nest is referred to as ‘drones’ [9]. Shortly after the drones depart, the females begin to exit to begin mating. It is said that the mating rate is heavily influenced by the temperature with the optimum occurring at 23 degrees [9] —an indication that mating can only occur in the late season. Once the mating season is over, all but the inseminated females survive, and they go on to seek shelter to hibernate over the winter period.

To this end, the seasonal temperatures have a profound impact on the life cycle of the Bombus population from triggering their release from hibernation to egg production and mating. As temperature in the UK during the summer has shown a steady increase over the last sixty years, its impact on the wider insect phenology is of great concern. This study investigates the effect of the temperature rise on the Bombus population, as they have a key role within the pollination community and, thereby, have serious implications on the local fauna and flora [10]. In particular, this study aims to unravel how the increasing seasonal temperatures in the UK has affected the Bombus phenology by delving into the main stages of their life cycle and investigating how temperature influences the size of Bombus population. The study also investigates how the seasonal temperature changes may affect the sphere of activities of the Bombus population at each stage of their life cycle by analysing their geographical trend overtime and movement.

Literature Review

Discussion surrounding the change in bee populations has been pursued for some time. Trends in the reduction and distribution of bee populations have been documented since the 1960’s with longitudinal and long-term analysis beginning to occur in the 1980’s [11]. Some of the first indicators of climate changes were studied by investigating the spatial distribution of certain Bombus species, with some disappearing from their native lands and others migrating to different areas (Williams 1982). The potential reasons for local declines have been explored since then, and several studies have reported the loss of agricultural sites [12], reduction in the availability of preferred flora [13], parasite invasion [14], the increased presence of electromagnetic radiation [15] and, above all, rise in the temperature. The UK has seen rise in local temperatures with the summer season being affected particularly heavily. Climate models predict that, without the mitigation of greenhouse gases, the possibilities of extreme heat (beyond 40 degrees Celsius) will be common in the foreseeable future [2]. The consensus amongst the literature points towards increasing temperatures negatively impacting the Bombus population, namely due to the deficiency in flight performance and the lack of colony productivity [16-18] reported that, due to the rising temperature in the UK, the start date of summer has advanced by an average of three days per decade since 1954. As of 2007, the summer start date stood as the 7^th of May: 18 days prior to that in 1954. This shift in the start date of summer implies that the exit from hibernation for Bombus species may have also advanced, as a queen’s emergence from dormancy is dictated by the temperature [19]. This assumption has been indeed recently proven by a group of researchers who studied the association between the temperature and the stages of the Bombus life cycle over a 35-year period in Central Europe [20]. They found that the rising temperatures are advancing the flight period by an average of 10-23 days. This piece of research had focused on four particular Bombus species and determined that, with an increase in temperature, both the emergence of queen bees and flight period had advanced. As mentioned earlier, the flight period begins when worker bees first exit the nest to begin foraging which continues throughout the entire summer. As the summer in the UK is advancing, this also indicates that the foraging period may not only be advancing but extending in length.

Kirbyshire and Bigg (2010) [18] also note a delay in the onset of autumn. The first frost date has generally been recorded later and later over a 50-year period, which suggests that winter is contracting on both ends. A more detailed study was also performed in the form of comparative analysis across 37 different climate models. The results showed that winter had contracted at an average rate of 2.1 days per year since 1952 and increased in temperature by 0.26 degrees every decade [21]. While findings by Guan and Yuping (2021) [21] is an aggregate across the entire northern hemisphere, it supports and develops on the study by Kirbyshire and Bigg (2010) [18]. As the Bombus mating period occurs between late summer and early autumn, the contraction of winter and the delay in the arrival of autumn may be also causing a delay in their mating period. Previous research on the association between the temperature and the Bombus lifecycle has reported varying findings. In many cases, they were studied under controlled temperatures where a set of bees were monitored over a period of time, and temperatures were artificially increased and decreased to study the optimum efficiency. These studies were aimed at identifying optimum temperatures to certain stages in a bee’s life cycle. For example, Kenna et al. (2021) [16] determined that the thermal performance of a bee is retained up to 27 degrees, suggesting that foraging activities during the summer months may be undermined if the temperature exceeds this limit. Another study showed that the queen exits hibernation between 5-9 degrees of temperature [19], suggesting that queens may emerge sooner if winter begins to become warmer. While research in the Bombus behaviour is rife, studies on the migration of or change in the spatial distribution of the Bombus population is limited. To date, most research papers have been in controlled environments, as it is difficult to monitor and assess the natural movement of bees. At the same time, there have been reports of a mass migration of bumble bees, mainly reported in the research domain of bird’s migration. The location of the new habitats of Bombus is unknown, but Fijen (2020) [22] suggests that they may be moving towards the north-east of their previous respective territories on the northern hemisphere. Research surrounding more local migrations has also been limited. One of the few exceptions was the study conducted by Williams (1982) [11] on the change of distribution in pre- and post-1960 Britain. They studied the distribution of all main British Bombus species and found that many species that were present pre-1960’s had either disappeared from the native area or had moved to neighbouring areas. They had also noted that the Bombus species have seen an overall reduction with those in more isolated areas becoming extinct. They concluded that this was due to the reduction in the favoured flora. However, in recent years, many of insect migration and related changes are being revisited and are considered to be triggered by the increase in temperature. Nevertheless, confirming the association between insect migration and temperature change is deemed quite complex, as it requires synthesis of knowledge from various domains such as meteorology, remote sensing and climatology [23]. This requirement for the combination of advanced knowledge from a variety of topics may be the reason why there is so little research studying the topic of temperature and insect migration, especially the correlation between bumble bee migration and temperature.

In a Bombus colony, the Queen bee is a lead figure and influences every stage of their life cycle. As mentioned earlier, her role as a Queen starts when she exits hibernation in spring. Makinson et al. (2019) [24] monitored the flight movement of Queen Bee’s in the county of Hertfordshire during a two-week period in 2015 to understand what dictates her flight and direction. They had inserted transponders into artificially hibernated Queens and tracked them using a harmonic radar system which recorded their GPS co-ordinates during their journey. Their results through random walk modelling (due to an expected unpredictability of queen’s flight movement) found that Queens spent more time resting than actively flying; on average, Queens fly 3 km to find their nesting site. They also note that there is no common pattern on their dispersal and, as Queens exit hibernation, they start moving at different angles randomly. The team also went on to analyse whether the Queens’ flight and distance were affected by wind direction, wind speed or solar hours; but concluded that none of them had any influence. This indicates that these climate features are not affecting the queen’s behaviour at exit, and that rise in the temperature may still be a key influencing factor.

Pawlikowski et al. (2020) [20] studied the influence of the temperature on Bombus sightings. They recorded the mean dates of the emergence of each bee type and period: first queen, first worker, first male and the beginning of the main flight period, end of main flight period, main flight and duration of main flight and performed 8 separate regression models for each respective type and period. They discovered that there was an advance in bumble bee foraging period over a 35-year period. They also discovered that the most significant change was in the Bombus’s flight period which occurs between June and July, and that this was the part of the life cycle that was affected the most by the change in temperature. The above review suggests that, while the association between temperature change and Bombus population, their lifecycle and their geographical movement have been investigated by some, it remains largely understudied, especially with respect to their spatial distribution. This study aims to analyse how the change in temperature is affecting the phenology and the spatial distribution of the Bombus Bumble bees. Specifically, we will explore the following questions: (1) which part of the Bombus cycle is affected the most by the increasing temperatures? and (2) how does the spatial distribution of the Bombus bees change? We will address them by means of exploratory data analysis and regression modelling.

Methodology

Datasets

The study area is the entire extent of the County of Oxfordshire, United Kingdom. It is located in Central England and provides a reasonable representation on the UK’s temperature change (i.e. no extreme changes of temperature anticipated). Historic data from the weather stations are also easily accessible. The data on the average maximum temperature was taken from the UK Meteorological (Met) Office archives between the years of 2011 and 2021 (Table 1). It shows a gradual increase in temperature over the 11-year period. The study period of 2011 to 2021 was determined by the quality of data available for the bumble bee sightings.

Table 1: Average Maximum Temperature in Oxfordshire County between 2011 and 2021

Bombus sightings in the Oxfordshire region between 2011 and 2021 (Table 2) was sourced from the National Biodiversity Network. It records the location where a variety of insects and animals were spotted. Data was retrieved between the years of 2011 and 2021, as data preceding 2011 was not sufficiently reliable to sustain robust analysis. While these records do not offer an exhaustive account of every spotting in the county, it provides a reasonable-sized sample of the population for analysis. It holds a total of 4130 sightings recorded throughout the 11-year period ranging across 24 different species. The data fails to categorize the individual bees, i.e. queen or worker. We note that data used for the bee spotting’s may have a degree of ‘chance sighting’ which can perhaps reduce the credibility of findings obtained through its analysis. Much of the literature studying phenology have mentioned the same issue; i.e. difficulty in finding data and hence one of the reasons this area has remained unexplored. However, bees are spotted throughout the entire county and across the year, thus implying that the entire county is being roamed for sightings. Sighting in each season is also very much in proportion with what we would expect from the bumble bees, i.e. less sighting in the winter and magnitudes more frequent in the peak of summer, then a decline in the late summer and early autumn. For this reason, we assume that the data is sufficient for investigating the relative difference in the frequency of sightings between different seasons.

Table 2: Bee Spottings in Oxfordshire County between 2011 and 2021

Models

This study uses OLS regression analysis to assess the influence of seasonal temperature (Table 1) on the number of bees (Table 2). The data are grouped into four categories for separate modelling, reflecting the four different stages in the Bombus life cycle as mentioned earlier. Here are the four models:

Model A (The Queen Model)

This model takes bee sightings between 1 January and 20 April, which is when Queen Bees are expected to exit hibernation to forage and seek nesting sites. While they mainly exit between late March and early April, some sightings were recorded in earlier months.

Model B (The Worker Model)

Existing studies refer to this period as ‘foraging’ or ‘flight’. This model takes on bee sightings between 15 May and 30 July, which is when worker bees are most prevalent.

Model C (Drone and Mating Model)

This model takes bee sightings between the dates 15 August and 15 September which is when worker bees have mostly deceased. The male bees survive on average 30 days to search for a mate.

Model D (New Reign Model)

This model takes on sightings between the dates 20 September and 10 December. During this phase, the inseminated females forage and seek hibernation sites to overcome the winter.

Due to lack of information on the classification of individual bees (e.g. queen or worker), we assume that the bee sightings during the respective period when a specific type of bee is most active (e.g. a queen bee in a period specified for Model A) are of that type of bee.

Once the regression models are executed, the most significant model will undergo spatial analysis. In order to investigate the extent and direction of any geographical change of the spotting data over the 11-year period, the spatial mean centre and the directional ellipses analyses are carried out. The mean centre will give an indication as to which direction, how quickly and how far the bees are migrating or spreading. The directional ellipses will give an indication of how they are geographically distributed during each period.

Analysis

Influence of Temperature on Bombus Sightings

The result of regression analysis shows a positive correlation between temperature and Bombus sightings in all models (Figure 1), which confirms the relationship of a higher temperature and a more frequent Bombus spotting. However, the influence of temperature on the Bombus spotting was statistically significant for Model A (The queen model) and Model B (The worker model) only, with p<0.05 of t-statistic for the regression coefficient (Table 3). On the other hand, the influence of temperature in Models C and D is statistically insignificant (with p>0.05), which means that there is little influence of temperature on bee sightings. The degree of goodness-of-fit of the regression model is in the order of Model A, Model B, Model C and Model D. This is shown in the R-squared values (Table 3).

Figure 1: Results of the Regression Models: (a) Model A (The Queen model), (b) Model B (The worker model), (c) Model C (Drone & Mating Model), (d) Model D (New Reign Model).

Table 3: Key Statistics from the Regression Models

The result generally confirms findings from past research in that it is difficult to understand how temperature affects the final stages of the Bombus cycle, as the influence the end of a season has on both fauna and flora is a lot more difficult to measure than that of early season [18]. Sparks and Smithers (2002) [25] also noted a similar challenge and asserted that phenological responses to temperature change is more prevalent in the early seasons and less so later [25]. Although climatologists have stated that the onset of autumn is seemingly delayed [21], it has not been transposed sufficiently enough in the Bombus population according to the two late season models, namely, mating and new reign; which makes it difficult to delve deeper into these two life stages. The queen model’ yielded the highest goodness-of-fit value and provides the regression model of y = 4.114x – 24.669. According to this equation, Queen Bees only exit hibernation once the temperature exceeds 6 degrees. This is supported by Alford (1969) [19] who reported that hibernation exit usually occurs once temperatures reach 6-9 degrees. This equation is only credible when the average maximum temperature is between 5-23 degrees; as the model was specified and performed within this range. Once the temperature exceeds 30 degrees, bee spottings are likely to decline as intense temperatures and humidity can have a negative impact on Bee morphology, impacting their flight performance [16]. Therefore, the assumption that the rising temperature brings forward the release of Queen Bees from hibernation can be considered to hold to a certain extent. Once this degree is reached, the rise in the temperature may result in the reduction of the number of Queen Bees present. However, this is yet to be seen until we see British winters reaching temperatures beyond 30 degrees. As such, this cannot be determined – only predicted. The analysis here found that an increase in seasonal temperature has a potential to increase the emergence of Queen Bees from hibernation and increase the number of worker bees. Although, the periods of mating and new reigns also show some form of positive correlation, their degree is significantly smaller than the periods of the earlier seasons (confirmed by a statistical insignificance of the regression coefficient and the model itself).

Bombus’s Migration

This section investigates the geographical movements of Bombus over time by drawing the Mean Centre and the Deviational Ellipse for different periods. The analysis focuses on the queen bees, as the regression analysis in the previous section found out that, among the four types of bees, Queen Bees are the most sensitive to the temperature change, and therefore can be regarded as a representative of Bombus on the reaction to temperature change.

Figure 2 maps the change in the mean centres and the deviational ellipses during the four periods of life cycles between 2011 and 2021. The figure firstly shows that the mean centre has consistently moved farther to the north-west as the years have progressed, which indicates that Queen bees—and hence their colonies—have also moved in this direction. The mean centre has moved each year covering a total distance of 16.24 km north-west during this period (Table 4). Although we do not see a direct proportionality between the temperature and the distance travelled, we note that, in the years between 2014-2016 and 2017-2019, there was a temperature reduction of 0.14 degrees, and the distance moved within this time was at its smallest at 3.7 km. The years between 2011-2013 and 2014-2016 showed a movement of 7.88 km with the temperature increased 0.89 degrees. The years between 2017-2019 and 2020-2021 showed a movement of 4.66 km, in response to the temperature increase of 0.45 degrees. The deviational ellipses show how the distribution of the Queen Bees have generally become more dispersed over the years, which is reflected in the size of the directional ellipse. If the size of the ellipse in the first time period was set at 1.0, that for the second, third and the fourth is 1.3, 2.5 and 1.1, respectively. The orientations of the ellipses are between 113° to 129.9° throughout the years, meaning that the orientation of the dispersion constantly expands north-west to south-east, which broadly consistent with their movement direction.

Figure 2: Change in the mean centre and the deviational ellipse (2011-2022)

Table 4: Mean Centre and Deviational Ellipses results

Discussion

Findings from this study suggest that the increase in the temperature has clear influence on the Bombus phenology. Results from the regression analysis show that there is an impact of rising temperature on the emergence of Bombus in all stages of their life cycle, particularly the early two stages. The regression model performed most robustly for the Queen model, indicating that she is most vulnerable to the temperature changes. Therefore, it is logical to assume that as the UK climate warms and spring continues to advance, more and more queen bees will begin to emerge sooner. If her resources are abundant upon her release and she has the ability to forage off her desired flora, her colony may also emerge sooner. On the contrary, if her early emergence means her desired flora is not available, she will go on to produce a futile colony which will unlikely survive a season. Earlier emergence can also cause an imbalance of the rest of the cycle, meaning that the proceeding stages must not function at their optimum. This will cause the extinction and reduction of many Bombus species, which has already been recorded [26]. The result of the spatial analysis shows that bumble bees are gradually moving in the north-westerly direction, which coincides with findings in some of the literature [22]. In summary, the geographical patterns of the Queen bees are: (1) moving towards north-west, (2) becoming increasingly dispersed, and (3) their spread constantly expand in the directions of north-west to south-east. Overall, there seems to be some consistency in the movement of Queen’s suggesting that there is a common form of migration going on. Unlike Makinson et al. (2019) [24] who recorded that dispersion and flight of Queens are unpredictable; our study saw a consistent pattern with a clearly specified directionality. Oxfordshire has seen a general rise in its temperature during the winter to spring months over the past eleven years. Of course, temperature increases are not demonstrated as a definitive cause behind the bee migration, but this study shows that it could be a contributing factor to a Queen’s re-location. Increasing temperatures can also assist a Queen’s metabolic process and, therefore, she would be capable of flying farther afield to nest (Kenna et al. 2021) [16], the earlier she emerges from hibernation, the more time she is afforded to search for and building a nest, leading to a further distance being flown; whilst warmer winters are tampering with a queen’s thermal limits, possibly leading to flight confusion [27]. The theories for the Queen bee’s dispersal patterns are almost endless and perhaps open up other research areas.

While these findings are implied and may fall short of proving the causation for this shift to be the change in temperature, it is worth noting that much of the literature points toward warming climates as the reason for the general insect migrations. Past research has shown butterflies have chosen to move south and grasshoppers have migrated to neighbouring regions to accommodate their climate needs [28]. If Bombus continuously begin to move northbound, it may eventually create reduction in the Bombus population in the southern region; which in turn may have an adverse impact on the wildlife biodiversity due to pollination reduction. There is also an associated risk to the migration of the Bombus population, if they continue to migrate north. As bumble bees can only survive in temperate zones, their migratory locations are limited and, therefore, they may be forced into areas which cannot accommodate their needs, thus heading towards the death of a colony. On the other hand, Roff and Fairbairn (2007) [29] note that insects which migrate into areas that are inhospitable to their liking, could trigger a genetic variation to increase their chances of living in what may have been perhaps perceived as an extreme climate in the past. If another dataset with more details on Bombus sightings was to be identified, it would be interesting to see how other climate attributes may be potentially driving this migration such as rainfall and humidity. The result of the analysis would be also more accurate, if it contained classifications of the type of Bombus bee spotted; i.e. a queen or a worker, similar to that observed by Palowkowski et al. (2020). This would have enabled us to build a more clearly defined models based on the type of bees; as opposed to dividing the regression models using seasonal dates and assuming we assign the most likely bees in that respective period; and this forms another future aspiration. Further analysis on the Queen’s emergence from hibernation can be extended by repeating this project in other countries where Bombus bees are native, as this will give insights into whether the tendencies found in this study are a global phenomenon or a more localised tendency.

Data Source

Bombus Sightings between 2011 and 2021. Source: National Biodiversity Network. Available at: National Biodiversity Network (nbn.org.uk).

Oxfordshire Temperature between 2011 and 2021. Source: Met Office. Available at: National Meteorological Archive – Met Office.

References

1. Damarre G, Rutishasher T (2009) Origins of the Word “Phenology”. Eos, Transactions, American Geophysical Union 90.
2. Christidis N, McCarthy M, Stott PA (2020) The increasing likelihood of temperatures above 30 to 40°C in the United Kingdom. Nature Communications 11.
3. Cook K, Vizy E (2015) Detection and Analysis of an Amplified Warming of the Sahara Desert. Journal of Climate 28: 6560-6580.
4. Lauenroth WK, Epstein HE, Paruelo JM, Burke IC, Aguiar MR, et al. (2004) Potential effects of climate change on the temperate zones of North and South America. Revista Chilena de Historia Natural 77: 439-453.
5. Rasmot P, Guillaume G, Terzo M (2021) Bumblebees of Europe and neighbouring regions, NAP Editions.
6. Wahengbam J, Raut AM, Pal S, Banu N (2019) Role of Bumble Bee in Pollination. Annals of Biology 35: 290-295.
7. McFrederick QS, LeBuhn G (2006) Are urban parks refuges for bumble bees Bombus spp. (Hymenoptera: Apidae)? Biological Conservation 129: 372-382.
8. Sarro E, Sun P, Mauck K, Rodriguez-Arellano D, Yamanaka N, et al. (2021) An organizing feature of bumble bee life history: worker emergence promotes queen reproduction and survival in young nests. Conservation Physiology 9. [crossref]
9. Belsky JE, Camp AA, Lehmann DM (2020) The Importance of Males to Bumble Bee (Bombus Species) Nest Development and Colony Viability. Insects 11. [crossref]
10. Baur A, Strange JP, Koch JB (2019) Foraging Economics of the Hunt Bumble Bee, a Viable Pollinator for Commercial Agriculture. Environmental Entomology 48: 799-806.
11. Williams P (1982) The Districution and Decline of British Bumble Bees. Journal of Agricultural Research 21: 236-245.
12. Carvell C, Roy DB, Smart SM, Pywell RF, Preston C, Goulson D (2006) Declines in forage availability for bumblebees at a national scale. Biological Conservation 132: 481-489.
13. Bommarco R, Lundin O, Smith HG, Rundlöf M (2012) Drastic historic shifts in bumble-bee community composition in Sweden. Proceedings: Biological Sciences 279: 309-315.
14. Meeus I, Brown MJ, De Graaf DC, Smagghe G (2011) Effects of invasive parasites on bumble bee declines. Conservation Biology 25: 662-671. [crossref]
15. Lazaro A, Chroni A, Tscheulin T, Devalez J, Matsoukas C, et al. (2016) Electromagnetic radiation of mobile telecommunication antennas affects the abundance and composition of wild pollinators. Journal of Insect Conservation 20: 1-10.
16. Kenna D, Pawar S, Gill RJ (2021) Thermal flight performance reveals impact of warming on bumblebee foraging potential. Functional Ecology 35: 2508-2522.
17. Guirad M, Cariou B, Henrion M, Baird E, Gérard M (2021) Higher developmental temperature increases queen production and decreases worker body size in the bumblebee Bombus terrestris. Journal of Hymenptera Research 88: 39-49.
18. Kirbyshire AL, Bigg GR (2010) Is the onset of the English summer advancing? Climatic Change 100: 419-431.
19. Alford DV (1969) A Study of the Hibernation of Bumblebees (Hymenoptera:Bombidae) in Southern England. Journal of Animal Ecology 38: 149-170.
20. Pawlikowski T, Sparks TH, Olszewski P, Pawlikowski K, Rutkowski L, et al. (2020) Rising temperatures advance the main flight period of Bombus bumblebees in agricultural landscapes of the Central European Plain. Apidologie 51: 652-663.
21. Wang J, Guan Y, Wu L, Guan X, Cai W, et al. (2021) Changing lengths of the four seasons by global warming. Geophysical Research Letters 48.
22. Fijen TPM (2020) Mass-migrating bumblebees: An overlooked phenomenon with potential far-reaching implications for bumblebee conservation. Journal of Applied Ecology 58: 274-280.
23. Sujayanand GK, Karuppaiah V (2016) Aftermath of climate change on insect migration: A review. Agricultural Reviews 37: 221-227.
24. Makinson JC, Woodgate JL, Reynolds A, Capaldi EA, Perry CJ, et al. (2019) Harmonic radar tracking reveals random dispersal pattern of bumblebee (Bombus terrestris) queens after hibernation. Scientific Reports 9.
25. Sparks TH, Smithers RJ (2002) Is Spring getting earlier? Weather 57: 157-166.
26. Lopez O, Koch J, Morales M, Navarette D (2021) Reduction in the potential distribution of bumble bees (Apidae: Bombus) in Mesoamerica under different climate change scentarios: Conservation implications. Global Change Biology 27: 1772-1787. [crossref]
27. Oyen KJ, Giri S, Dillon ME (2016) Altitudinal variation in bumble bee (Bombus) critical thermal limits. Journal of thermal biology 59: 52-57. [crossref]
28. George M (2022) Migration of butterflies. Kerala Forest Research Institute 17: 844-847.
29. Roff D, Fairbairn DJ (2007) The Evolution and Genetics of Migration in Insects. BioScience 57: 155-164.

DOI: 10.31038/GEMS.2022433

Abstract

Green infrastructure refers to all manners of established green spaces that provides habitat for wildlife and ecosystem services for people. Areas with richer biodiversity is known to provide measurable benefits to the local wellbeing and the wider ecosystem, but it remains unclear if the size and the type of greenspace affect the extent of species richness. This study investigates the distribution of 92 species across green infrastructure in England to identify the variables which have the most effect on species richness. The results suggest that Forestry Commission Woodlands and Country Parks had the highest species richness, whilst Doorstep Greens and Village Greens returned the lowest species richness. ANOVA confirmed that the variation between the groups was significant. The main difference in the groups was the area, with Forestry Commission woods being on average 140 times larger than Doorstep Greens. When considering the greenspace with the highest species richness, habitats having highest mean area were woodland, managed grassland and waterfronts. These results confirm our intuition but also suggest that richer biodiversity can be only achieved in well preserved and managed woodlands and country parks occupying a sizable plot, and that regular patches of green infrastructure embedded within the urban areas do not offer high species richness.

Keywords

ANOVA, Biodiversity, Green infrastructure, Species richness, Woodlands

Introduction

Green Infrastructure (GI) is a relatively new term that is loosely defined as “a spatial structure providing benefits from nature to people [which] aims to enhance nature’s ability to deliver multiple valuable ecosystem goods and services such as clean air or water” [1]. More generally, it is a reference to established green spaces and new sites (of greenspace) that is considered to have a positive effect on our wellbeing and the wider society) [2], and they range from large areas such as agricultural land, forests, wetlands, woodlands and parks to individual street trees, private gardens, green roofs, rivers and transport corridors (UK Green Building Council 2015) [3]. Aside from the various social benefits and cultural services they offer to our society, GI has the potential to provide habitat and connectivity to wildlife. For instance, cities in the United States have seen an increase in the frequency and diversity of wildlife sightings, as more cities have invested in GI to improve their respective urban environment and achieve environmental resilience [4]. A systematic review of the literature [5] suggests that increase in biodiversity is generally considered to yield positive effect on human health. Also, GIs with greater diversity of avian species are known to evoke positive emotional responses from their visitors [6].

However, what triggers a greater range of biodiversity remains unclear. Lawton (1998) [7] once suggested that a coherent and resilient ecological network for biodiversity requires a series of larger habitat patches with good connectivity to other sites. In other words, larger areas usually hold more species present within them and the connectivity to other areas reduces the risk of localised extinctions due to changes in the local conditions. However, securing a large area of green space within and around the confined space of cities can be challenging. Also, it is unclear as to what types of green space are likely to offer high species richness. Most studies on GI and general greenspace have focused on a small sample size, or a short duration and were often unable to differentiate between greenspace type or quality. As a result, the association between the attributes of GI (the types of GI and the size in particular) and biodiversity remains understudied. Furthermore, with the predicted change in climate, a number of studies have investigated the ability of GI to alter the urban environment and alleviate the impact of climate change such as extreme temperature, rising sea levels, frequent extreme weather events and increase in the spread of diseases and pests [8-11]. Climate change is also considered to affect biodiversity with an increase in the intensity and frequency of droughts, storms and fires. For instance, the intense fires in Australia in 2019 and 2020 have increased the number of threatened species in GI and the vicinity of cities by 14% (The Royal Society 2021) [12].  Against this background, this study investigates which aspects of greenspaces would yield benefit for people and wildlife through maintaining and enhancing biodiversity.

Literature Review

Impact of Green Infrastructure on People

Green space is often considered to have a positive effect on our wider wellbeing, ranging from the reduction of post-surgery recovery time [13] to lowering cortisol levels and blood pressure [14]. However, this assertion is not conclusive, and some studies report a mixed or even negative impact of green space on our wider wellbeing. For instance, Shanahan et al. (2015) [15] systematically reviewed relevant studies and concluded that some studies reported reduced mortality and generally enhanced wellbeing, whilst some studies found no association between green space coverage and mortality, and one study even reported increase in mortality in relation to larger green space coverage. Similarly, Maas et al. (2009) [16] investigated the association between green space coverage and morbidity and other disease within a 1km window, and concluded that, the degree of urbanity of the study areas, rather than the area of green space, had stronger impact on the health outcomes. At the same time, Maas et al. (2009) also highlighted that having access to a greater green space coverage helped reduce anxiety related disorders. Others (Gong et al. 2014 [17], Coutts and Hahn, 2015) [18] report that, while accessibility of green space generally improves our health and wellbeing, the effect varies between different types of green space.

Impact of Green Infrastructure on Wildlife

The direct and indirect effects of GI on humans is wide ranging and findings in the literature are also varied. The impact of GI on wildlife, on the other hand, is more straightforward in that it provides habitat connectivity in an increasingly urbanised world [19]. Several studies have attempted to assess which aspects of GI are important for different species, and their findings vary greatly, including size, connectivity/isolation, the management regime and type of habitats present within the GI [20-22]. Birds are among the most widely studied taxa in relation to GI, because they are relatively easy to survey and their life cycles are well understood [23]. For instance, Chamberlain et al. (2007) [24] investigated bird species richness in urban parks in London and found that site area was the most important factor related to bird species richness and that, in smaller sites, the presence of a high number of adjacent gardens caused the species richness to increase. In addition to size, those areas which held rough grass or waterbodies also returned higher species richness. The authors noted that the relationship between size and species richness simply reflects the fact that larger greenspaces tend to contain a higher number of habitats, which could attract a wider range of species. Negative associations were also recorded for the presence of buildings, roads and pavements across both the species groups and the seasons. In a similar study, Zorzal et al. (2020) [20] found that the taxonomic diversity of bird species recorded across six urban greenspaces in Brazil was positively associated with greenspace area. However, when analysing the bird diversity against a proximity index, the study found that there were no correlations with maximum noise or the degree of patch isolation. They also found that species richness had a positive correlation with the heterogeneity of the habitats within each greenspace. The study was limited to a relatively small sample size as it was constrained by the number of accessible urban greenspaces available. Other studies ranged from bird populations in Sweden [22], to bumble bee populations in relation to the urban features of San Francisco [25], amphibians in New Jersey, McCarthy and Lathrop (2011), and general insect diversity in green roofs in Nova Scotia [26], most of which found richer biodiversity in GI away from urban centres, and in GI with larger area.

Furthermore, a study in Slovenia by Groot et al. (2021) [27] found that out of the 2 most prevalent urban greenspace types in Ljubljana (urban forest and managed park) there was a lower abundance of bird species within parks than urban forest, but that species richness and species composition were not affected by greenspace type. Greenspace area significantly affected the bird abundance and greenspace type influenced which species was designated as indicative of the habitat. Although there were no significant differences in bird biodiversity between urban and peri-urban areas, Ljubljana holds relatively large areas of urban forest and the impact of urbanisation (e.g. disturbance, predation) are likely to be reduced. The link between biodiversity and the provision of ecosystem services has been explored further in a literature review by Sandifer et al. (2015) [5]. In their review, papers concentrating on the effects of biodiversity on ecosystem services were considered in order to provide a body of evidence on the effects of biodiversity on human health. The review highlighted that in the most part, biodiversity assists good human health and that in studies where relationships were not found, this could be attributed to a lack of suitable data, although there is the possibility of a confounding effect of socio-economic status on human health. Most of the studies investigated had inadequate sample sizes, study durations or measured controls, or were found to be lacking in rigorous statistics; objective data; unable to differentiate between greenspace type or quality; showing evidence of long-term effects; or suffering from the suite of health data collected being too narrow. Very few studies also attempted to identify the mechanisms through which the effects of biodiversity work.

Assessing the Quality of Green Infrastructure

Many of the studies highlighted above did not account for the quality of the GI in their analyses. Quality of GI can be measured in the field, which allows for detailed surveys of local areas, but this is time-intensive and cannot be easily extrapolated across larger areas. Another way to determine quality is through the use of GIS and spatial modelling where a set of spatial criteria are captured through remote sensing and/or a composite index was used for representing the regional environmental quality [28]. The quality of GI measured with respect to their benefit to humans as well as the wildlife can be analysed by assessing the current value of different Ecosystem Services within an area. One of these tools is the Outdoor Recreational Value tool (ORVal), which seeks to estimate the number of visits a greenspace might receive and the monetary equivalent of those visits. The ORVal tool estimates the value of a GI through the application of a recreational demand model which places importance on the cost in time of visiting a greenspace, and the quality of the recreational experiences at that greenspace [29]. The quality of recreational experiences is thus considered to be based on an extensive set of variables encompassing greenspace type, size, land cover, designations, points of interests and direct connectivity with other greenspaces. These studies suggest that GI has a generally positive effect on humans and the wildlife, but there is a gap in the literature where a systematic investigation into their effect and the contributing factors towards species richness is understudied at the larger, national scale.

Context and Data

As mentioned earlier, majority of the studies carried out on evaluating the impact of green space focus on a small sample size or a specific type of greenspace. This study will investigate the relationship between the recreational value of greenspace and the number of species they hold across England. England contains an estimated 27,000 public parks, and its National Planning Policy Framework (NPPF) mandates new developments to provide more GI—in the form of safe and accessible areas of greenspace with recreational, cultural and social facilities—as these areas are considered to help mitigate climate change effects and deliver wider benefits for nature [30]. Other areas of nature reserve and conservation areas are also protected from development. Despite that, the amount of green space in urban areas has dropped from 63% to 55% between 2001 and 2018 [31]. To mitigate the loss of sites, maintaining GI of high quality and the capacity to nurture biodiversity is essential. The quality of GI can be measured in a variety of different ways, including the proximity to urban areas and types of land cover, and it is difficult to extract a single attribute as a proxy to describe each GI. The Outdoor Recreational Value (ORVal) tool mentioned earlier is a tool developed by the University of Exeter in collaboration with DEFRA. Its primary purpose is to give an estimate of the quality for people of greenspaces across the United Kingdom both in monetary terms and as an estimated number of visits. At the same time, the tool also models the interactions between the habitats present, whether the site has any legal designations or points of interest, and whether there are any shared boundaries with other greenspaces; and it is these auxiliary data that this study will exploit. Habitat areas in ORVal are derived as a combination of the 25m-grid of the 2007 Land Cover Map (LCM), Ordnance Survey Master Map data, the Priority Habitat Inventory dataset from Natural England [32] and Open Street Map data.

The estimated greenspace visits are taken from the Monitor of Engagement with the Natural Environment (MENE) survey administered by Natural England. The estimated value of each greenspace is calculated through an opportunity cost model of recreational trip choice, whilst taking into account the socio-economic factors. While both values are estimates derived from the respective surveys and cost models, they have been calibrated through empirical studies to improve their accuracy [29]. Species data is taken from the National Biodiversity Network (NBN) which provides access to arguably the most comprehensive set of aggregated biodiversity data from multiple recording schemes across England and Wales. The NBN gateway limits the maximum number of downloads to 500,000 records per species and 10 million in total. To comply with these limits, the number of species included in this study was limited to amphibians, birds, butterflies, mammals (including bats) and reptiles listed in section 41 of the Natural Environment and Rural Communities (NERC) Act 2006; namely, those of principal importance for the purpose of conserving biodiversity and are considered most likely to be accurately recorded across the whole of England. Species groups such as beetles, moths and freshwater fish were excluded due to the lack of reliability in the measurement of their recording and/or because of their irrelevance to the terrestrial habitats under study.  Some species listed in section 41 are deemed sensitive, and to avoid reporting the specific locations of these species, the NBN gateway provides the location at a reduced geographical resolution. The exact location of these species records is stored internally and only available to individuals having enhanced access, with the data provided on the public interface being generalised to a lower resolution of 1km grid. To align all analyses, location data of the species records was renumerated to 1km grids. Further data processing and cleaning were conducted on QGIS platform and with R-Studio software. The ORVal website provides a GIS shapefile detailing the boundaries of the greenspaces. Greenspace types that may not be open to the public (e.g. golf course) were removed from the data set. This left 22,698 greenspaces in the analysis. Table 1 shows the number of each greenspace type after the data processing.

Table 1: Number of each greenspace type

Species records were plotted as points, and the points falling within each 1km grid square of England were aggregated to the respective grid so that the total number of records within each grid square was classified by the species’ groups. Additionally, the presence of each species was calculated within the grid squares; i.e. for each species group, the presence of species within that group was counted to give the number of unique species. The number of species records in each group is shown below (Table 2).

Table 2: Number of records for each species group

While different species and attribute data exhibited different patterns of distribution, a log(Y+1) transformation offered an overall best fit to bring each distribution closer to normal distribution and was therefore performed for all variables.

Analysis

The independent variables were identified through the literature as having potential correlations to species richness were plotted in scatter graphs (Figure 1). To understand the contributing factors for the biodiversity (unique number of species), this study conducted OLS regression. In other words, the OLS regression was used to determine which variables were significant in explaining the sum of unique species. Figure 2 shows the association between these variables and the number of species within greenspace. As a result, the following variables were used in the OLS regression model as the independent variables: area, wood, natural grass, managed grass, parking, urban percentage, rivers and canals, and the number of habitats. All variables were analysed using the transformed log(Y+1) data.

Figure 1: Scatterplots of the log of unique species against variables used in OLS regression

Figure 2: Variables affecting the number of species within greenspace

Figure 1 shows the scatterplots of the variables used in the OLS regression model against the unique species number. Whilst there is strong positive correlation between some variables, most have weak positive correlation even after transformation. The outputs of the scatterplots suggest that the analysis would be best achieved with non-parametric tests. However, due to the number of records used in the analysis, applying non-parametric tests to the entire dataset holds the risk of returning inaccurate results [33,34].

To explore the spatial concentration of these variables, the Local Indicator of Spatial Autocorrelation (LISA) was calculated in the form of local Moran’s I for each variable in the OLS regression (Figure 3). The maps clearly show that “area” size has very few hotspots or coldspots, with smaller low-low area clusters being generally around urban areas and larger high-high area clusters being generally clustered in rural areas. A large proportion of the greenspaces have no significant clustering of area, except for woodlands which show a clear contrast between the hotspots (or the areas with high concentration of woodlands surrounded by similarly high proportion of woodlands) and the remaining regions. In particular, counties with higher percentages of woodland area such as Surrey appear to be well represented in hotspot clusters in the woodland variable, with counties having typically lower woodland cover such as Suffolk and Somerset showing up as coldspots. There are areas such as Cornwall and the Lake District which have low overall woodland coverage but are both shown as having hotspot clusters, suggesting that in these counties the woodlands which are present occur in close proximity to each other rather than being dispersed across the wider county. Managed grassland habitat areas tend to be larger in urban areas whilst natural grassland habitat areas tend to be larger in rural greenspaces. Parking is again fairly sparse with large areas of no significant clustering. Whilst there appear to be some hotspots for parking in urban areas and some coldspots in rural areas, this does not hold true across the country.

Figure 3: Local Moran’s I clusters for each log variable: area; woodland; natural grassland; managed grassland; parking; rivers and canals; number of habitats; and urban percentage.

Rivers and canals have large significant hotspots in north-west and south-west with coldspots mainly concentrated in the south and south-east. The total number of habitats shows significant hotspots mainly in urban areas, however as man-made habitats were included in the calculation for number of habitats, there is perhaps a predisposition to urban areas with significant coldspots occurring in rural areas in a similar pattern to the managed grassland variable. The urban percentage clusters appear generally as expected, although there are comparatively few high-high clusters in urban areas. Table 3 shows the results for the OLS regression of unique species counts. It illustrates the highly significant relationship for all variables with a low standard error (≤0.5162). The variables explained 76 % of the variance of the species richness. Analysis of the t-values shows that woodland area has the strongest positive relationship with unique species number, and rivers and canals showing the weakest relationship if still significant. Urban percentage was shown to have a weak negative relationship.

Table 3: Results of OLS regression

To examine the relationship between greenspace type and species richness, a boxplot was produced (Figure 4). It highlights the difference in the number of unique species in each type of greenspace with Forestry Commission Woods and Country Parks having the highest species richness, whilst Doorstep Green and Village Green showing the lowest number of species. Woods and Nature had large variations in species richness, partly due to the frequent outliers present for these categories. An ANOVA test returned a significant result for the variation, where the sum of squares = 3184 and p = < 2.2e-16 with 1 DF.

Figure 4: Boxplot of log +1 unique species count and greenspace type

Further analysis of the means for the lowest and highest groupings are shown in Figure 5. Unsurprisingly, both country parks and forestry commission woods have higher areas of woodland cover than doorstep green does, although country parks also have more managed grassland. Parking, built habitats, and rivers and canals were also higher in country parks than the other greenspace types investigated. Whilst it was not plotted on the graph due to the large variation in size even at a log scale and the resultant skewing of the y-axis, Forestry Commission Woods were on average approximately 140 times larger than Doorstep Greens, and average area of Country Parks were approximately 50 times larger than Doorstep Greens.

Figure 5: Mean habitat areas by greenspace type

Discussion

Figure 3 showed that the contributing factors towards high species count (hence, high biodiversity) tend to be spatially aggregated into known species diversity hotspots such as the New Forest and Jurassic Coast, Ainsdale NNR, the Cambridge fens and the Norfolk and Suffolk coasts. The presence of hotspots mainly outside of urban area validates the use of the total unique species to determine species richness rather than using species abundance. The presence of coldspot clusters also confirms results consistent with findings from previous studies whereby species richness decreases with increasing urbanicity [25,26].  The decision to aggregate species counts by 1km² may have affected the clustering, and a smaller grid could have given a more detailed picture of the species distributions. The large blocks of hotspots and coldspots may have been an artefact of the large data set, and analysis using smaller grid units may have resulted in more nuanced distribution of clusters.

The locations of the high-high clusters suggest that the best greenspaces for both people and wildlife mainly occur on the edges of suburbs around large urban centres, or in larger urban parks, although this may again simply be highlighting an issue with the extraction method for the species data. Analysis of the variables which make up each of the cluster grouping suggests that those greenspaces which are composed of woodlands with managed grassland are more likely to benefit wildlife. However, clusters of high species-low value greenspaces had a high mean woodland area, without a high mean managed grassland area. The large impact of woodlands on high biodiversity may also be an artefact of the species data which was used. Bird records accounted for more than 60% of the records and studies investigating the effects of landcover on bird diversity have reported a strong relationship between species richness and woodland cover [22,35]. The results of the OLS models confirm that the woodland area of a greenspace is the most important factor when predicting the number of unique species present within that greenspace. This finding reiterates the suggestion above that the heavy skew towards bird records within the original species data, and the strong relationship between bird biodiversity and woodland area, may be over-emphasising the relationship.

Further analysis of the grouped species data would be necessary to identify if the same trends are found in other species groups. The relatively low result for area in the regression model was surprising, with the estimate being ranked 5^th out of the 8 variables modelled. This may be due in part to the negative multicollinearity with both natural grassland and rivers and canals. It may also be that since there are many factors that could affect each species differently, if the species groups were modelled separately then the size of each greenspace may have shown more significance in explaining species richness for generally lower mobility groups such as butterflies. The relatively high effect of parking suggests that even though attempts were made to reduce recorder bias, bias may still exist within the data as it gives the impression that an increase in parking area resulted in higher unique species counts. The variance in species richness caused by these variables is unlikely to be geographical in origin, as the geographic variation cannot be discerned at this spatial scale. The exceptions are the managed grassland and natural grassland variables, both of which broadly follow the respective distributions. A potential alternative method would have been to geographically subset the data in order to provide the local Moran’s I for each area. Also, splitting the data down into distinct geographical regions such as the south-east, south-west and so on may have provided a method to investigate how the variables behave with an adaptive bandwidth. Indeed, the effect of greenspace type on biodiversity appears to reflect their size and ruralness. Both country parks and forestry commission woods tend to be large areas of greenspace, which may explain the positive relationship between these and the number of unique species. Country parks have a set criterion for designation including a minimum area, facilities and accessibility whilst being a predominantly semi-natural landscape. Country parks should be over 10ha in size and as an increase in area typically results in an increase in the number of habitat types and thus available habitat niches, as well as reducing population isolation and disturbance [36]. As Forestry Commission Woods are the largest supply of sustainably managed timber in the UK, and commercial forestry operations require large amounts of land, Forestry Commission Woodlands are on average the largest greenspace type within the analysis.

Conclusion

This study found that green space outside the urbanised areas has a high level of species diversity than those within urban areas. Variables that contributed to high species richness were area, woodland cover, grassland cover (both managed and natural), rivers and canals, parking, number of habitats and percentage urban cover. Greenspace types with large areas (e.g. country parks and forestry commission woods) corresponded to a higher number of unique species present than those that are usually small in size (e.g. village greens and doorstep greens). Where size is not dictated by greenspace type, such as in the case of woods and commons, these result in a lower average unique species with a large number of outliers being present in the data. Potential avenues for further investigation into this topic include analysis of different species groups and comparing their findings with those from this study; and identifying similar greenspaces with similar recording effort across the country in order to test whether the differences shown in species richness during this study are geographical or an artefact of over- and under-recording. Despite these challenges, this study confirmed that richer biodiversity can be generally only achieved in well preserved and managed woodlands and country parks occupying a sizable plot, and that regular patches of green infrastructure embedded within the urban areas do not offer high species richness. This has policy implications in that, while we can explore the positive effect of having a varying extent and types of green infrastructure within urban areas, we cannot expect high biodiversity in relation to such green infrastructure.

References

1. EU (2013) Green Infrastructure (GI) — Enhancing Europe’s Natural Capital.
2. Natural England (2009) Green Infrastructure Guidance, York.
3. UK Green Building Council (2015) Demystifying Green Infrastructure, London.
4. Hunold C (2019) ‘Green Infrastructure and Urban Wildlife: Toward a Politics of Sight’, Humanimalia 11: 89-109.
5. Sandifer, P. A., Sutton-Grier, A. E. and Ward, B. P. (2015) Exploring connections among nature, biodiversity, ecosystem services, and human health and well-being: Opportunities to enhance health and biodiversity conservation. Ecosystem Services, Elsevier, 12: 1-15.
6. Cameron RWF, Brindley P, Mears M, McEwan K, Ferguson F, et al. (2020) ‘Where the wild things are! Do urban green spaces with greater avian biodiversity promote more positive emotions in humans?’, Urban Ecosystems 23: 301-317.
7. Lawton JH (1998) Daily, G. C. (Ed.). 1997. Nature’s services. Societal dependence on natural ecosystems. Island Press, Washington, DC. 392 pp. ISBN 1-55963-475-8 (hbk), 1 55963 476 6 (soft cover).’ Animal Conservation, Cambridge University Press, vol. 01, no. 01, p. S1367943098221123.
8. Arntz, W., Strobel, A., Moreira, E., Mark, F., Knust, R., Jacob, U., Brey, T., Barrera-Oro, E. and Mintenbeck K (2012) Impact of climate change on fishes in complex Antarctic ecosystems. Advances in Ecological Research 46: 351-426.
9. Forzieri G, Feyen L, Russo S, Vousdoukas M, Alfieri L, et al. (2016) ‘Multi-hazard assessment in Europe under climate change’, Climatic Change 137: 105-119.
10. Emilsson T, Sang ÅO (2017) ‘Impacts of Climate Change on Urban Areas and Nature-Based Solutions for Adaptation’, ln: Kabisch N, Korn H, Stadler J, and Bonn A (eds) Nature-Based Solutions to Climate Change Adaptation in Urban Areas: Linkages between Science, Policy and Practice, Cham, Springer International Publishing, pg: 15-27 [Online].
11. IPCC Secretariat (2021) Scientific review of the impact of climate change on plant pests – A global challenge to prevent and mitigate plant pest risks in agriculture, forestry and ecosystems., Rome.
12. The Royal Society (2021) Climate Chang eand Biodiversity – Interlinkages and Policy Options, London.
13. Ulrich RS (1984) View Through a Window May Influence Recovery from Surgery, Science, American Association for the Advancement of Science 224: 420-421.
14. Grazuleviciene R, Vencloviene J, Kubilius R, Grizas V, Dedele A, et al. (2015) ‘The Effect of Park and Urban Environments on Coronary Artery Disease Patients: A Randomized Trial’, Castelnuovo, G. (ed), BioMed Research International, Hindawi Publishing Corporation.
15. Shanahan DF, Lin BB, Bush R, Gaston KJ, Dean JH, et al. (2015) Toward improved public health outcomes from urban nature. American Journal of Public Health 105: 470-477. [crossref]
16. Maas J, Verheij R, de Vries S, Spreeuwenberg P, Schellevis F et al. (2009) ‘Morbidity is related to a green living environment’, Journal of Epidemiology and Community Health 63: 967-973. [crossref]
17. Gong, Y., Gallacher, J., Palmer, S. and Fone, D. (2014) Neighbourhood green space, physical function and participation in physical activities among elderly men: The Caerphilly Prospective study. The International Journal of Behavioral Nutrition and Physical Activity 11. [crossref]
18. Coutts C, Hahn M (2015) Green Infrastructure, Ecosystem Services, and Human Health. International Journal of Environmental Research and Public Health 12: 9768-9798. [crossref]
19. Zhang S, Ramírez F (2019) Assessing and mapping ecosystem services to support urban green infrastructure: The case of Barcelona, Spain. Cities 92.
20. Zorzal R, Diniz P, Oliveira R, Duca C (2020) Drivers of avian diversity in urban greenspaces in the Atlantic Forest. Urban Forestry & Urban Greening 59.
    
21. Holtmann L, Philipp Becke, Fartmann T (2017) Effects of habitat and landscape quality on amphibian assemblages of urban stormwater ponds. Urban Ecosystems 20.
    
22. Mörtberg U, Wallentinus HG (2000) Red-listed forest bird species in an urban environment — assessment of green space corridors. Landscape and Urban Planning, Elsevier 50: 215-226.
23. Leveau LM, Isla FI (2021) Predicting bird species presence in urban areas with NDVI: An assessment within and between cities’, Urban Forestry & Urban Greening, Urban & Fischer 63.
24. Chamberlain DE, Gough S, Vaughan H, Vickery JA, Appleton GF (2007) Determinants of bird species richness in public green spaces. Bird Study, Taylor & Francis 54: 87-97.
25. McFrederick QS, LeBuhn G (2006) Are urban parks refuges for bumble bees Bombus spp. (Hymenoptera: Apidae)? Biological Conservation 129: 372-382.
26. MacIvor JS, Ksiazek-Mikenas K (2015) Invertebrates on Green Roofs. Pg ; 333–355.
27. de Groot M, Flajšman K, Mihelič T, Vilhar U, Simončič P et al. (2021) Green space area and type affect bird communities in a South-eastern European city. Urban Forestry & Urban Greening, 63.
28. Rahman Md R, Shi ZH, ad Chongfa C (2014) Assessing regional environmental quality by integrated use of remote sensing, GIS, and spatial multi-criteria evaluation for prioritization of environmental restoration’, Environmental Monitoring and Assessment 186: 6993-7009.
29. Day B, Smith G (2017) The ORVal Recreation Demand Model: Extension Project, Exeter.
30. Communities and Local Government Committee (2012) The National Planning Policy Framework: Eighth Report, London.
31. Public Health England (2020) Improving Access to Greenspace: A new review for 2020, London.
32. Day B, Smith G (2016) Outdoor Recreational Valuation (ORVal) Data Set Construction, Exeter.
33. Skovlund E, Fenstad G (2001) Should we always choose a nonparametric test when comparing two apparently nonnormal distributions? Journal of Clinical Epidemiology 54: 86-92.
34. Fagerland MW (2012) t-tests, non-parametric tests, and large studies—a paradox of statistical practice? BMC Medical Research Methodology 12.
35. Jokimäki J (1999) Occurrence of breeding bird species in urban parks: Effects of park structure and broad-scale variables. Urban Ecosystems 3: 21-34.
36. Zhao JM, Zhou LZ (2018) Area, isolation, disturbance and age effects on species richness of summer waterbirds in post-mining subsidence lakes, Anhui, China. Avian Research 9.