Abstract
Bark beetles (Coleoptera: Curculionidae: Scolytinae) are mortality agents to multiple tree species throughout North America. Understanding spatiotemporal dynamics of these insects can assist management, prediction of outbreaks, and development of “real time” assessments of forest susceptibility incorporating insect population data. Here, dispersal of Douglas-fir beetle (Dendroctonus pseudotsugae Hopk.) is estimated over four regions within Colorado and Wyoming from 1994 to 2010. Infestations mapped from aerial insect surveys are utilized as a proxy variable for Douglas-fir beetle (DFB) activity and analyzed via a novel GIS technique that co-locates infestations from adjacent years quantifying distances between them. Dispersal distances of DFB infestations were modeled with a cumulative Gaussian function and expressed as a standard dispersal distance (SDD), the distance at which 68% of infestations dispersed in a given flight season. Average values of SDD ranged from under 1 kilometer for the region of northwestern Colorado to over 2.5 kilometers for infestations in Wyoming. A statistically significant relationship was detected between SDD and infestation area in the parent year, suggesting that host depletion and density-dependent factors may influence dispersal. Findings can potentially provide insight for managers—namely, likelihood of DFB infestation increase for locations within two to five kilometers of an existing infestation.
1. Introduction
The Douglas-fir beetle, Dendroctonus pseudotsugae Hopkins, (DFB hereafter) is a major mortality agent of Douglas-fir, Pseudotsuga menziesii (Mirbel) Franco, across the Western United States [1, 2]. As a native insect and a natural disturbance agent, it is always present as an endemic influence, playing an important ecological role by killing diseased or otherwise stressed trees. DFB exhibits one generation per year and attacks new hosts every year during its dispersal flight from early spring and through the summer depending on geographic location. The insect overwinters primarily in the adult stage and as larvae inside the host tree [2]. Population levels increase periodically, resulting in widespread tree mortality [3] which can impact management resource objectives and ecosystem services. These eruptive populations usually develop after other disturbances such as fire [4], windstorms [5], or defoliation [6] events which provide an abundance of stressed trees that the insect can exploit. Once stressed trees are no longer a suitable resource, populations can disperse into surrounding stands, where the insect prefers to infest areas with an abundance of large-diameter trees growing in dense environments and exhibiting reduced growth [7, 8].
Land managers and forest health protection specialists often use models based on current stand conditions to estimate the likelihood of infestation development in a stand or the potential tree mortality levels when populations increase [8, 9]. These models do not incorporate data about active bark beetle populations and are hampered by limited knowledge about DFB dispersal at the landscape level. Accurate ways to characterize and predict movement of Douglas-fir beetle populations across landscapes are needed to develop spatially referenced predictive models that can identify forests stands more likely to be infested due to the presence of nearby beetle populations. Such information can then be presented in the form of landscape-level assessments of the potential impact of epidemic populations. This holds the promise of directing management activities at susceptible high value areas in a timely manner. Shore and Safranyik [10] present such information for mountain pine beetle, Dendroctonus ponderosae Hopkins, in British Columbia.
Powers et al. [11] point out that infestation occurrence and severity depend on several factors involving multiple scales of space and time. In a multiscale study of DFB involving area wide and local spatial scales, the authors found that on larger spatial scales strong winds and prolonged drought preceded epidemics, and on more local scales, greater likelihoods of infestation were linked to drier southern facing stands, sites situated at lower elevation, and areas with greater abundance of host.
At all scales, outbreaks depend on the population dynamics of DFB. Without an abundant number of beetles, outbreaks will not occur. In many cases, what causes a beetle population to increase suddenly can be unknown and difficult to determine. However, beetle populations are a function of fecundity, mortality, immigration, and emigration, with immigration and emigration being two forms of dispersal.
Among available data sources that can be used to describe the presence of the DFB, the most prolific is the USDA Forest Service Aerial Insect and Disease Survey (IDS). Obtained from low altitude flights and conducted by regional forest health specialists, the data set provides a comprehensive and spatially explicit representation of the locations, sizes, intensities, and the geographic shapes of various forms of forest damage, including tree mortality from DFB attacks [12, 13]. Using the above data source, Dodds et al. [14] characterized patch-level characteristics and the spatial relationships between DFB infestations in northern Idaho, observing that when overall beetle populations were clearly at eruptive levels, frequencies of nearest neighbor distances between infestations in two adjacent years showed a distinct concentration at lower distances with a gradual dropoff as distances increased, becoming near zero at between 2000 and 2500 meters. The present work builds on the study by Dodds et al. [14] by examining dispersal distances of DFB. We use the locations of new infestations as a surrogate for beetle dispersal. Our objective is to examine and quantify spatial dispersal of DFB using aerially detected new infestations at the landscape scale in parts of Colorado and Wyoming.
2. Materials and Methods
2.1. Site
Our analysis focused on four areas in Colorado and Wyoming where DFB populations frequently occur—the mountainous regions of the Colorado eastern slope (ECO), southwestern Colorado (SWCO), northwestern Colorado (NWCO), and Wyoming (WY) (Figure 1). The Colorado eastern slope region (approximately 37–41°N, 104–106.5°W), extending from the Continental Divide down to the urban corridor of the Colorado Front Range, includes the Roosevelt National Forest, Rocky Mountain National Park, an eastern portion of the Arapaho National Forest, and extends south to include the Pike-San Isabel National Forest. The northwestern Colorado region (approximately 39–41°N, 105.67–109°W) includes the White River and Routt National Forests as well as the western portion of the Arapaho National Forest. To the southwest, this region (approximately 37–39.4°N, 105–109°W) comprises the Gunnison, Rio Grande, Grand Mesa, San Juan, and Uncompahgre National Forests. Lastly, the Wyoming Region (approximately 42.25–45°N, 106–111°W) comprises the Bighorn and Shoshone National Forests as well as portions of the Medicine Bow and Bridger-Teton National Forests.
2.2. Aerial Survey Data Source and Background
Insect and disease survey (IDS) data is gathered by the USDA Forest Service Forest Health Protection group and state forestry agencies cooperators and is generally collected by low-flying aircrafts with aerial observers marking infestation locations in paper maps. More sophisticated methods using maps in touch screen devices linked to global positioning systems are also available but not used as frequently. Data is stored as polygons in one geographic information system shape file for each year. Observers are trained to recognize the characteristic “signature” of disturbances such as the crown fading of bark beetle-killed trees [12, 13]. Since Douglas-fir trees fade in coloration a year after attack by DFB, IDS data for one year actually represents infestations initiated the year before. Observers record delineations of affected areas by an agent, including relative size, shape, and location with the best accuracy possible. An estimate of the number of trees affected is also recorded. Each polygon shows an area of tree damage, mortality, or both, either by bark beetles, defoliators, pathogens, or abiotic factors, and each polygon has an entry in an associated attribute table that stores damage agent information.
Observations from moving aircrafts in complex terrain are inherently difficult, and even though observers are highly skilled, there is unavoidable spatial error of position, size, and shape of affected areas by the different agents. IDS data presents a picture of forest conditions at a landscape scale [15]. Johnson and Ross [16] examined the spatial accuracy of an IDS flight by using ground plots to confirm aerial observations. They reported that, when allowing tolerances in positional error of 0 m, 50 m, and 500 m, the overall accuracy of polygon location was 61%, 68%, and 79%, respectively. The authors concluded that these levels of error make the data suitable for coarse-scale landscape analysis but not for fine scales. Therefore, we assumed here that these issues of spatial accuracy produce negligible consequences with analyses designed to detect landscape-level phenomena on the scale of multiple kilometers, which is the scale at which we present our study. For future work, Wulder et al. [17] indicate that the spatial accuracy issues have improved over time with the advent of spatial mapping systems and may be improved in the future by the further use of real-time imagery during the mapping process.
2.3. Characterizing Douglas-Fir Beetle Dispersal
Quantification of dispersal requires the postulation of causational dependence between infestations in a given year and infestations nearby in a subsequent year, more specifically between all infestations in a given previous year (heretofore called “parent infestations”) and all infestations in a corresponding upcoming year (heretofore called “child infestations”). This becomes a reasonable assumption when we consider that the insect is univoltine [2]. Dispersal is estimated by quantifying distances between such polygons in adjacent years. Using a 50 m grid we rasterized infestation polygons from all years and created for each polygon what in GIS terms are referred to as “allocation zones” [18], in which each 1994 infestation was then surrounded by a larger zone containing all raster grid locations having the given infestation as the nearest parent infestation (Figure 2). This was then overlaid with similarly rasterized child infestations from the following year, and specific parent-child causational dependence was assumed based upon which parent allocation zone each child infestation was located—that is, which previous-year parent infestation was closest to each child infestation (Figure 3). This procedure was performed on all pairs of consecutive infestation years beginning with 1994-1995 through 2010-2011, and independently for each of the four regional study areas.
2.4. Quantifying Dispersal Distances
In light of the complex geometric shapes inherent to many observed DFB infestations across the landscape, a systematic sampling technique was used to quantify distances of dispersal. Straight lines from each raster midpoint in a parent infestation were connected with all possible raster midpoints in an associated child infestation in the upcoming year (Figure 3). Point-to-point connections were not allowed to cross allocation zone boundaries, since each child infestation was assumed to only be associated with the nearest parent infestation. This resulted in a histogram of sampled distances for each region and pair of consecutive years (e.g., Figure 4). These histograms were further normalized by dividing each histogram bin value by the maximum number of possible grid cells that can be located at the given distance from a single cell, a value that roughly increases by the square of the distance but can exhibit some error because of the finite raster resolution.
With the exception of the restriction of a child infestation to only a single-parent infestation, the above method approximates the more known methodology of calculating a spatial autocorrelation [19]. The autocorrelation coefficient at a given spatial distance is defined as follows: where is the autocovariance function: and is the variance function: For our purposes, these functions are greatly simplified by the fact that we are dealing with binary data describing “detection” or “no detection”: Then the equation for becomes where defines the number of instances of positive detection taking place at a given distance . Similarly, becomes Thus, which is the number of sampled distances between parent and child infestations occurring at a distance divided by the total possible sampled distances at the same distance (note also that spatial autocovariance is here assessed only for detections, not the absence of detections, since the latter is not considered relevant and so is excluded from (5)).
2.5. Evidence for Causal Relationships
To be able to examine distances of DFB infestations between years, it was important to confirm that the distribution of new infestation across the landscape was independent from the locations of forest stands which contained the host tree and that there was indeed a relationship between the location of parent and child infestations. To confirm this, we examine if the histograms of sampled distances possessed a statistically significant departure from results expected from random spatial infestation occurrence confined exclusively within a spatially aggregated host as infestations in our study sites are restricted to locations of Douglas-fir stands or mixed stands composed of Douglas-fir and ponderosa pine (Pinus ponderosa P. & C. Lawson). Hence, histograms like the one shown in Figure 4 were compared against a null hypothesis of infestations occurring randomly throughout Douglas-fir stands. For each parent infestation year, this was done by repeating the analysis after replacing the observed child infestation layer with a rasterized layer showing locations of all nearby stands containing Douglas-fir. This procedure gave equivalent histogram results for infestations systematically distributed throughout all possible host locations, which for this analysis provides results equivalent to those of a random distribution throughout host. Rejection of the null hypothesis would support the assumption of causality between parent and child infestations.
Normalized cumulative versions of both histograms were compared using the Kolmogorov-Smirnov test, which compares the maximum differences between two cumulative histogram curves. Since the null hypothesis histogram uses allocation zones from previous year infestations and, thus is potentially compromised by observational error inherent in the IDS data, the two-sample test was used, which is more conservative than the one-sample test that assumes one distribution randomly samples from the other [20]. This analysis was performed for two locations from which vegetative data locations were obtained—Pike-San Isabel National Forest and White River National Forest.
2.6. Modeling Dispersal Distances
To describe the complex histogram distributions with a single statistic, cumulative versions of each histogram distribution were modeled using a cumulative Gaussian function [21]: where the parameter is the estimated value of a standard dispersal distance (SDD) for each flight season, from which beetle dispersal for a given region and year is collectively described with a single parameter, a parameter which will be compared with coincident values of infestation area in the parent year in an effort to investigate a relationship via linear regression [22].
2.7. Sensitivity Analyses
Optimal usage of any spatially explicit data source requires an understanding of the degree to which any spatial uncertainties in the data may or may not affect the outcome of any analytical results. A sensitivity analysis was performed by modeling the histogram results with varying degrees of positional error. This was done by modeling each raw histogram bin value as a Gaussian distribution, where negative tails were made into positive distances by taking the absolute value (Figure 5).
The procedure resulted in two adjusted histogram functions that resulted from using Gaussian standard errors of 540 and 764 meters, respectively. As the spatial errors of aerial detection survey polygons indicate locations that are reliably accurate to the landscape scale of a watershed catchment (i.e., within 2 to 3 kilometers) [16], then the above standard errors are representative of the actual spatial uncertainties associated with these data. The effect observed is one of smoothing the histogram, but not one of significantly changing the histogram’s overall shape, or its estimate of SDD. Sensitivity analyses were also conducted to test our methodology with spatial resolutions other than 50 meters. From these alternate resolutions, the resulting histograms produced identical values of SDD.
All GIS procedures were performed using the ESRI ArcGIS package versions 8.3, 9.0, 9.1, 9.3, and 10.0 [23]. Model fits were estimated using PROC NLIN from the SAS statistical package v9.1 [24]. Additional statistical analyses were performed and graphs generated using the statistical package [25].
3. Results
Standard dispersal distances are shown in Table 1 with infestation areas of their parent years (e.g., infestation area in 1997 is shown with the SDD value for 1997-1998). Regional averages for SDD ranged from about 800 meters (NWCO) to about 2,600 meters (Wyoming) with annual values ranging from 204 m to over 10 km in a Wyoming season that also displayed the highest value of overall infestation area (Table 1). Results were omitted where curve fits were not estimable.
After what appears to be a momentary episodic increase in infestation area in 1994, populations on the Colorado eastern slope reached consistently elevated levels of over 1000 ha in 1998. This was accompanied by a more long-term upward trend, with total infested area reaching 6,000 to over 7,000 hectares in 2003 and 2004, followed then by a decline and a resurgence in 2010. The largest proportional increases in total infestation area occurred at four times—in 1998, 2003, 2006, and 2009. These moments approximately coincide with the largest values of SDD occurring in 1998, 2001, 2005, and 2008. There is one infestation area peak in 1994 that for lack of data cannot be described in terms of annual growth and may be purely episodic. It is, however, approximately coincident with a similar peak SDD value in 1995.
Results for southwestern Colorado appear to show two different sets of dynamics across two separate time frames. From 1994 through about 2001 infestation areas were relatively low and values of SDD more variable across years, even momentarily reaching over 4 km in 1995 and dipping to less than 500 m just two years later. For the remaining years of 2002 through 2010, infestation areas were consistently higher than those in the previous timeframe, and values of SDD were more consistent in the 1–2.5 km range, even displaying an apparent semiannual oscillation between the top and bottom of this range in the years from 2001 through 2007.
Infestation areas in northwestern Colorado were the smallest of all four study areas. Populations in the mid-1990s were small enough that no dispersal analysis was possible. In addition, when values of SDD were calculable, such values were generally smaller than those of other study areas. Infestation area appears to have increased in stages from 1999 through 2007, where two episodically high values of SDD in 2001 and 2004 preceded two separate quantum increases in infestation area in 2003 and again in 2006. The high infestation area in 1996 occurs in isolation from any other coincident phenomena.
Wyoming showed the highest values of total infestation area as well as values of SDD, with the results of 2005 showing over 60,000 ha of infestation area and having an SDD value of over 10 km. Infestation area from 1998 through 2005 showed two periods of apparently exponential growth, from 1998 through 2002 and again from 2003 through 2005. Values of SDD in the same time period show a similar pair of consistent increases from 1999 through 2011 and again from 2002 through 2005.
Rejection of the null hypothesis comparing histogram results against a null hypothesis of random infestation occurrence confined to locations of Douglas-fir host was typical, indicating that the spatial aggregation of Douglas-fir host stands cannot fully explain the location of infestations. Figure 6 illustrates an example of the results for the Kolmogorov-Smirnov test that compared cumulative histogram results with an alternate cumulative histogram computed by replacing child infestations with the entire host area, effectively showing what the results would be if the child infestations were to occur either completely throughout the host, or randomly spaced throughout the host. This latter histogram is displayed with a 95% confidence region, where any departure from this region by the dashed line would indicate a rejection of the null hypothesis. The histograms for actual child infestations almost invariably showed more aggregation than the host; therefore, host distribution solely cannot explain the spatial location of infestations. Although there are no means of directly connecting our primary histogram results to beetle flight, these comparisons against a null hypothesis imply beetle dispersal as the driver of the results. Two exceptions to this are the 2004-2005 and 2005-2006 flight seasons in the Pike and San Isabel National Forests (not shown), where the results were not statistically distinct from results expected from infestations spread completely throughout all Douglas-fir host, which is reasonable considering the magnitude and spread of infestation in those years, which likely affected most of the susceptible stands in the landscape.
Values of SDD demonstrated a significant difference among regions (), with average regional values ranging from about 1 to 2.5 kilometers. In addition, a statistically significant relationship was detected () between values of SDD in one year and total infestation area in the following year. In light of the process of dispersal, this relationship is intuitively reasonable, since larger dispersal distances allow beetles to have access to larger areas to infest in the following year.
Lastly, a regression analysis (Table 2, Figure 7) was performed to relate SDD and infestation area in the current year. The association was detected as statistically significant (, ), where and , respectively, refer to SDD and infestation area in a given year, and the second, third, and fourth variables are binary indicators of which region is under consideration. The various trend lines show an increase in SDD with increases in infestation area. Values in the above binary variables are consistent with statistically significant differences in SDD across regions, with the Wyoming region having the largest SDD values and northwestern Colorado having the smallest.
4. Discussion
Years of increasing infested area were oftentimes approximately coincidental with years of elevated values of SDD. The instances of increased dispersal after high infestations may be associated with local host depletion and the presence of distant favorable hosts, thereby fostering increased dispersal distances. In addition, density-dependent factors influencing insect population dynamics may also be indicated as triggering increased dispersal distance.
Dodds et al. [14] estimated DFB infestations to have average adjacent year nearest neighbor distances mostly residing in the 2–5 kilometer range. Our results corroborate this finding as well as that of the theoretical dispersal studies of Byers [26], who modeled dispersal of Ips typographus (L.) and indicated that beetles can disperse quickly over large areas but are influenced by factors such as wind, tree density, and angle of turn by the insect. Laboratory studies have shown that DFB attached to rotating mills flew continuously up to nine hours at velocities of 30 to 40 meters per minute [27], which would translate to dispersal distances potentially exceeding 20 km. The present analysis suggests that this distance, although apparently biologically possible, does not represent what commonly occurs in nature.
Dispersal distance is partially a function of resistance to movement as progeny of the parent population moves into new host environments. Although specific factors that drive beetle dispersal and establishment of new infestations in the landscape are largely outside of the scope of this study, Barclay et al. [28] demonstrated the seminal role that host connectivity can play in the landscape dispersal dynamics of the closely related mountain pine beetle. The work of Negron [8] would suggest the development of new infestations in areas of high tree density, high percentage of Douglas-fir basal area, and reduced host growth rates. At larger scales, Powers et al. [11] would stress the importance of aspect, elevation, windthrow, presence of drought, and spatial aggregation of old conifers. Resistance sources include scarcity of host, diversity of host, and distance between suitable hosts. According to Kausrud et al. [29], “mean connectivity depends not only on landscape connectivity but also on factors like dispersal mortality, fat content of the insect, and attraction arrestant (windfalls and pheromone sources) density.” The need for further study on this is clearly indicated.
The present work focuses on a representation of beetle dispersal that gives an indication of dispersal distance, but not dispersal direction. Directional anisotropies in dispersal patterns would likely be the result of wind patterns and stand structure [26], among other factors requiring subsequent analysis. Of additional consideration here would be the relationship of dispersal to differences in habitat connectivity or the ability to traverse landscape topography [30] and related aspects of habitat connectivity and fragmentation [31]. Coulson and Witter [32] divided the process of dispersal into two categories, active and passive, with active dispersal taking the form of flying or walking and passive dispersal involving windborne transport. Although the authors above classify bark beetle dispersal as active, comparing directional anisotropies in wind patterns to those of beetle dispersal may elucidate the degree that such dispersal is a mixture of processes inherent to both categories.
Ecological processes operate at a range of spatial scales, and the range of spatial scales inherent to a given analysis can have an impact on which processes are detected by that analysis [33]. In reference to the findings of Powers et al. [11] aspect is likely a particularly strong influence in Colorado as Douglas-fir occurs mostly in north facing slopes which provide more moisture. Elevation would also likely be a strong factor as Douglas-fir in Colorado has a well-defined elevation range in which it occurs. At the patch scale (tree and stand level), growth rates have also been correlated with DFB attacks [8]. In addition, the spatial distribution of DFB-caused mortality in stands with suitable diameter classes have been found to aggregate strongly around intrastand neighborhoods of high basal area [34].
Although the use of variograms in our simplified binary form (4) presents some novelty, the use of variograms to model spatially explicit ecological data is well documented [35–37], and curves used to represent the results tend to take the form of exponential, Gaussian, Whittle-Matern, or spherical. Examples of such methods being applied specifically to dispersal, however, are limited. In the context of dispersal of seeds, graphs of the form of Figure 4 would be referred to as “seed shadows” (e.g., Holbrook and Smith [38]). Further, when applying curve fits to variograms describing dispersal, the use of a Gaussian form has some biological backing, since the dispersal process can be thought of (and modeled as) a sequence of a large number of random smaller-scale movements [39], and the Central Limit Theorem in statistics indicates that distributions resulting from sums of many independent random processes theoretically approach a Gaussian distribution [40].
Schröder and Seppelt [41] suggest four foundational core questions to any model-based analysis of pattern-process interactions. These questions involve (1) the detection and quantification of spatiotemporal patterns, (2) the analysis of such patterns for underlying ecological processes, (3) the understanding of the scale dependencies of such analyses and processes, and (4) the development of methodologies for simulation of such processes. To these four, we suggest adding a fifth core question, namely, (5) the influence on management decision making. The present analysis addresses the first question by quantifying spatial patterns at the landscape level and demonstrates that spatial distribution of infestations is statistically distinct from those that would arise from spatially random infestation occurrence. Powers et al. [11] and Dodds et al. [14] begin the explanation of underlying processes of insect dispersal (question 2), yet the problem requires further investigation. The third question is addressed by clearly defining the spatial scales inherent to the analysis—namely, those of landscape dispersal. With underlying ecological processes not addressed in the present work, the question of simulating such processes (question 4) is similarly left to future analysis. Lastly, question five, management decision making, can be addressed by providing this information to forest health specialists who can monitor potential outbreak development by examining the locations of current infestations in light of these dispersal distances. The information can also be of use in planning documents by forest managers such as Forest Plans. Ultimately, combining DFB dispersal knowledge with stand susceptibility models can provide a “real-time” assessment of stands at risk of infestation as presented for mountain pine beetle in Canada [10].
5. Conclusions
Our four study areas showed four different populations of DFB which together explained a wide array of population dynamics. The results among these four study areas, though quantitatively different, retain some apparent qualitative similarities. A model is presented which, although lacking in predictive capability (), indicates that increases in total infestation area are associated with increases in SDD. In the infestation risk model of Shore and Safranyik [10], threats of mountain pine beetle infestation show a dependence on the proximity of any neighboring infestations, with threats reaching lowest values when the nearest infestations are more than 4 km away. With average values of SDD found to be in the range of 1–2.5 km, the results of this study can provide more quantified estimation of preventive management costs, by indicating that such management, to be successful, would require management and monitoring efforts to span at least 3–5 kilometers from existing beetle outbreaks (assuming an appropriate cutoff to be two standard dispersal distances which would theoretically encompass 95% of all beetle dispersal). Eventually such estimates of dispersal can be incorporated into spatially explicit landscape models to more effectively predict the risk of future infestations [31]. Lastly, the above technique can be very feasibly applied to other beetle species whose presence is described by aerial detection surveys or similar data sources such as aerial photography or other remote sensing applications.
Acknowledgments
Sincere thanks are extended to Rudy King at the USDA-FS RMRS facility in Fort Collins, CO, USA, for providing ongoing statistical consultation for this project. This project was supported by funding from the USDA Forest Service National Fire Plan Project 01.RMS.C4, Rocky Mountain Research Station, Fort Collins, CO, USA.