Modelling the Incursion and Spread of a Forestry Pest: Case Study of Monochamus alternatus Hope (Coleoptera: Cerambycidae) in Victoria

Abstract

Effective and efficient systems for surveillance, eradication, containment and management of biosecurity threats require methods to predict the establishment, population growth and spread of organisms that pose a potential biosecurity risk. To support Victorian forest biosecurity operations, Agriculture Victoria has developed a landscape-scale, spatially explicit, spatio-temporal population growth and dispersal model of a generic pest pine beetle. The model can be used to simulate the incursion of a forestry pest from a nominated location(s), such as an importation business site (approved arrangement, AA), into the surrounding environment. The model provides both illustrative and quantitative data on population dynamics and spread of a forestry pest species. Flexibility built into the model design enables a range of spatial extents to be modelled, from user-defined study areas to the Victoria-wide area. The spatial resolution of the model (size of grid cells) can be altered from 100 m to greater than 1 km. The model allows core parameters to be altered by the user, enabling the spread of a variety of windborne insect species and pathogens to be investigated. We verified the model and its parameters by simulating and comparing the outputs with the 1999/2000 Melbourne incursion, but no establishment of a forestry pest beetle was believed to be Monochamus alternatus Hope (Coleoptera: Cerambycidae). The model accurately predicts the distance and direction of the historic incursion, and the subsequent failure to establish is due to low overall population density of the pest species.

1. Introduction

Australia has a robust, dynamic and evolving plant biosecurity system. Nonetheless, pressure from increasing international trade and travel has identified capacity and capability gaps. An increasing rate of exotic forest pest establishments poses significant challenges to Australia’s forest biosecurity [1]. Quantifying the benefits of pest management actions that reduce spread of invasive species is critical for optimising pest management programs, evaluation of existing programs and to assist in determining future actions. Successful early detection of exotic forest pests significantly contributes to the effectiveness of biosecurity activities and potential outcomes.

Many post-border forestry pest surveillance programs have been piloted in Eastern Australia involving port-environ surveillance, urban tree monitoring and targeted trapping [2,3]. Expansion of these activities into a nationally coordinated program, targeting high-risk sites and priority pests, will maximise the chance of early detection and possible eradication, thereby reducing the likelihood of exotic pests becoming established in Australia’s forests [1,4,5].

Pine beetles (from the families Curculionidae and Cerambycidae) including Dendroctonus and Monochamus species are among the most serious potential pests in forestry in Victoria, Australia [6]. In Australia, 9% of 20,000 crates inspected in Sydney, Brisbane and Melbourne during 2005 have “something of quarantine concern”, such as bark, fungi, live insects, or frass, including Monochamus alternatus [5]. Monochamus is a genus of longhorn beetles (Coleoptera: Cerambycidae) commonly known as sawyer beetles. Their larvae bore into dead or dying trees, especially pines. The nematode Bursaphelenchus xylophilus which causes pine wilt disease is exclusively associated with the genus Monochamus. Pine wilt diseases can cause decline and death in heavily infected trees within one growing season; crowns of infected trees turn from green to reddish brown [7] and symptoms can develop in a matter of weeks. Many Monochamus species have been recognised as Australian national priority pests [8]. The insect and associated nematode have caused extensive damage to pine forests where it has invaded countries with susceptible hosts and with native co-occurring Monochamus beetle species [5]. An Australian pest risk analysis shows that there is a moderate likelihood of pinewood nematode and its primary vector, M. alternatus (Japanese pine sawyer beetle), arriving at, establishing, spreading, and impacting Pinus plantations in Australia [5]. The species of Pinus used in Australian plantations have never encountered this pest and most likely have very little natural defense.

Being able to predict how organisms that pose biosecurity threats will establish and build up population densities and spread across landscapes is essential for designing effective and efficient systems for surveillance, eradication, containment and management of key biosecurity threats.

2. Materials and Methods

We selected the insect pest to be modelled using the following criteria:

• A high-priority pest for the forestry industry.
• An incursion of the pest would require an emergency response.
• Historical data of an incursion in a similar environment to Victoria.
• Sufficient information available on the life history and dispersal of the pest.

Two pests that fitted these criteria were pine beetles from the genus Monochamus (Coleoptera: Cerambycidae) such as the Japanese pine sawyer, M. alternatus, and pine bark beetles such as the mountain pine beetle, Dendroctonus ponderosae (Coleoptera: Curculionidae). It was decided to focus the model on M. alternatus, as it is a vector of the pinewood nematode, Bursaphelenchus xylophilus. As noted in Akbulut and Stamps [9]: “Studies to date indicate a remarkable similarity among beetle species (Monochamus) around the globe for a variety of life-history traits, including lifespan, adult emergence numbers, flight capability, nematode transmission rates and attraction to pine volatiles.” As such, with little other information available, we can assume the model is applicable across the Monochamus genus. However, it was recognised that the model should have the capacity to accommodate other species with only minor changes to either the programming code or parameters within the model.

2.1. Life Cycle of Monochamus Species and Bursaphelenchus Xylophilus

The life cycles of the Monochamus species and B. xylophilus have been extensively researched, reviewed [5,10,11,12,13] and are summarily illustrated in Figure 1. The emergent adults, carrying the nematode, fly to healthy pine trees and feed during daylight hours on succulent young twigs for 7–12 days. During this feeding and maturation period, the nematode enters feeding wounds in the bark of the smaller branches. The nematodes feed on the lining of the tracheids, disrupting water flow, and susceptible trees wilt and die as a result of the nematode damage to the phloem.

After maturation feeding, sexually mature beetles fly to dying or recently dead pines. Females deposit eggs below the bark and nematodes may enter the oviposition wounds. Both mated and unmated females can produce abundant eggs and transmit nematodes through oviposition scars, although unmated females only lay sterile eggs. The nematode will still proceed through its life cycle and kill the tree even if sterile eggs have been laid. Susceptible trees eventually wilt and die, often rapidly (within 2–3 months). Water stress and high temperatures are consistently associated with expression of pine wilt disease, especially where mean temperatures exceed 20 °C for long periods. Beetles are then attracted to these dying trees or other abiotically stressed trees and oviposit eggs under the bark, also transmitting nematodes.

The beetle egg hatches and the larva feed on the inner bark, cambium and outer sapwood, progressing through 3–8 larval instars before pupation. The insect predominantly overwinters as a larva but can survive as an egg or a pupa. The insect pupates in a chamber for several days where the nematode crawls into the body of the insect via the spiracles and trachea. Adults emerge reproductively immature. Upon emergence they disperse, usually within 100 m, however, long-flight dispersal has been observed with flight mill experiments [12] demonstrating a maximum flight duration of 115 min.

2.2. Model Components

The model was developed mechanistically by including all factors that satisfied two criteria:

• They were considered important influences in the establishment and spread of pine beetles.
• The data were available at a scale that was suitable for inclusion in the model.

The model description follows the Overview, Design concepts, Details (ODD) protocol [14,15]. In summary, the model consists of several modules that interact to simulate the development, maturation and spread of a pine beetle across a real-life landscape:

(a)

Pine beetle life cycle module which simulates the growth from eggs, through the larval stages, pupae and then to adult, dispersal of a proportion of adult beetles to suitable nearby trees, mating and egg production.

(b)

Wind spread module that disperses a proportion of adult beetles across the landscape during spring and summer, depending on actual or averaged wind direction(s) and its interaction with the dispersal kernel.

(c)

Host availability and suitability module that accurately describes the distribution, abundance and health status of trees within the Pineacea family (pine, fir and spruce species), based on local, state and federal government datasets available across Victoria (see [16] for a full description of the host availability and suitability dataset development). Most of the local and state government pine survey data were of high quality and spatially accurate. All data sources were curated for data quality and duplicate records from the diverse data sources were removed. Each grid cell contains a count of healthy, healthy (but infected by the nematode), stressed (but uninfected by the nematode), stressed (after being infected by the nematode) and dead pine trees.

2.3. Model Design

The model operates across a landscape populated with real data, and can accommodate any spatially defined local, regional or statewide extent. The model creates an initial master grid based on a user-defined spatial resolution (>100 m × 100 m) and extent. All subsequent spatial layers are then re-projected and resampled to align to this master layer. The model is initialised with a user-defined pest infestation site, such as pre-defined assigned importation-approved arrangement (AA) site, a point selected interactively within the model interface, or user-defined co-ordinates.

The model operates on a daily time-step. The time frame of the model run is defined by the user. For any point in time (day) that the model runs, each cell contains the following information:

• The number and health status of pine trees,
• Average meteorological wind data, relevant to the time frame of the model run,
• The presence, sex, and life stages of any beetles,
• The nematode infection status of any beetles, and
• Degree day information relevant for beetle development.

Based on meteorological conditions and the distribution of pine trees, the model calculates the dispersal of beetles over time, the infection of trees by the pathogen, and the beetle’s population growth (or decline).

The sub-models, which describe how each module works (life cycle, host availability and suitability, and dispersal), are summarised below and in more detail in our technical report [16]. The report also details the parameter values used in each module, and the sources and rationale for the data used. We utilised MatLab^® (Mathworks, Inc., Natick, MA, USA) as a modelling platform as it supports both live and static climatic data feeds.

2.4. Life Cycle

The life cycle of the pinewood nematode and its vectors have been extensively researched (as reviewed by [5,9,10,11,12,13]). As the insect is not known to have established in Australia, we utilised the European and American research to estimate the thermal degree day thresholds for the key stages in the life cycle (i.e., from eggs to adult emergence, to sexual maturity, and to death). Daily climate data were sourced at a 5 km gridded resolution across Victoria [17], enabling statewide estimates of the likely date of occurrence of key phenological stages. Within the model construct itself, each beetle was identified as an agent with assigned properties such as; “sex”, “mated”, “location”, “degree day counter” and “phenological stage” (

Table 1. Agent properties for the Monochamus beetle model.

	Property	Description
	X	X-coordinate (meters) of agent
	Y	Y-coordinate (meters) of agent
	PhenoStage	Current phenological stage of agent [1,2,3,4,5]
	DegDayCount	Current cumulative degree count for agent
	AdultDayCount	Current number of days since emergence
	Male	Male/Female (true/false)
	Mated	Whether agent has mated (true/false)
	Infected	Is agent infected with nematode (true/false)

). Agent behaviour was associated with the underlying gridded base-layers of the model, which included degree day data for phenological stages, wind data for long-distance dispersal function and tree-health information.

2.5. Feeding, Mating and Ovipositioning

During spring and summer, Monochamus beetles emerge from trees infected and killed the previous year, carrying pinewood nematodes in their tracheal system. This component of the model also updates the number and health status of any pines. The adult beetles fly to healthy pine trees and feed on young shoots and twig bark (Figure 2).

In general, mating occurs around 10 days after the adults have emerged, fed, and possibly infected healthy trees. Females then lay an average of 60 to 109 eggs in a stressed tree [10,11]. However, mating can only occur if both males and females are present both spatially and temporally (i.e., within a localised region when the female is ready to mate). The model accounts for this, by stipulating mating occurring when there are both male and female beetles within the same cell after dispersing from a feeding tree. At low densities, mate limitation reduces the reproduction potential and decreases the rate of population growth, which may lead to inbreeding depression and lead to failure of an incursion to establish [18]. In addition to mated females, unmated females can also lay eggs and subsequently infect host trees with the pathogen, which is captured in the model. Associated with each phenological stage is a beetle mortality rate. During oviposition, eggs are laid by mated females at a predefined daily rate. All beetles progress through their annual life cycle until they die and are then removed from the model.

2.6. Phenological Development

Degree days were calculated as a function of minimum and maximum temperatures and a base threshold, Tbase (Equation (1)).

$$\max \left(\frac{Tmax+\max \left(Tmin,Tbase\right)}{2}-Tbase,0\right),\mathrm{where}Tbase=12.5°\mathrm{C}$$

(1)

Phenological thresholds were applied to the cumulative degree data. The relationships between minimum and maximum temperatures and degree day are illustrated in Figure 3.

An overview of phenological stages for agents in the Monochamus model is presented in Figure 4. Degree day thresholds for the phenological development stages were obtained from the United States Department of Agriculture Animal and Plant Health Inspection Service (USDA, APHIS) model [19].

We determined the length and completion dates of the phenological stages across Victoria and Melbourne Region. As expected, there is spatial and temporal variability in the average adult emergence date (utilising 2010–2018 meteorological data) across Victoria, with most beetles emerging between early December and mid-March. Median results are presented in

Table 2. Degree day thresholds (DD 12.5 is the cumulative days above 12.5 °C) [19], typical length and completion dates of phenological stage across (1) Victoria and (2) the Melbourne Region.

Days Required to Complete Each Phenological Stage at DD 12.5		Victoria		Melbourne
		Completion of Stage (median)	Stage Length (days)	Completion of Stage (median)	Stage Length (days)
Eggs	80	12 September	73	25 September	86
Larvae	705	2 January	112	18 January	115
Pupae	865	19 January	17	8 February	21
Pre-oviposition	965	30 January	11	21 February	13
Oviposition	1333	20 March	49	7 May	75

.

2.7. Host Availability and Suitability (Pine Model)

Pine tree locations were based on two types of spatial data; polygons of softwood plantations areas and point data representing individual or groups of trees surveyed. Several government, council and departmental datasets were used to create this layer (see [16] for full details). Data for all susceptible Pineacea species were included, representing pines, fir and spruce. The most common species were Pinus. Original data sources were used to develop a statewide pine tree density layer, calculated using a 100 m grid cell. This was used as an initial layer in the model.

The host availability and suitability sub-model iteratively updates the abundance and health status of all pine trees in each of the pixels. Each grid cell contains a count of the number of pine trees within the following categories:

• Healthy: tree is growing robustly (this category is included in the initial input dataset and is based on observations recorded in source data).
• Healthy and infected by the nematode: tree is not showing signs of stress but has been infected by the nematode after beetle feeding (this category is updated following beetle dispersion, feeding and egg laying events within the model).
• Stressed but uninfected by the nematode: tree is showing signs of stress but has not been infected by the nematode (this category is included in the initial input dataset and is based on observations recorded in source data).
• Stressed and infected by the nematode: tree is showing signs of stress and has been infected by the nematode (this category is updated following beetle dispersion, feeding and egg laying events within the model).
• Dead: tree is dead but still present (this category is included in the initial input dataset and is based on observations recorded in source data).

Even with large numbers of beetles entering a stand of healthy trees, not all trees will be infected. Under the assumption that each beetle randomly selects a healthy tree, the expectation that a single tree, $T_i$, will be infected if B beetles enter the stand (Equation (2) and Figure 5a) can be described as:

$$E\left(T_i=infected\right)=1-{\left(\frac{T-1}{T}\right)}^B$$

(2)

It follows that the total number of trees infected is the expectation of a single tree being infected multiplied by the total number of trees in the stand, T (Equation (3) and Figure 5b).

$$Infectedtrees=T\left(1-{\left(\frac{T-1}{T}\right)}^B\right)$$

(3)

2.8. Dispersal

The distance and direction of adult beetle dispersal is determined by the status of pine trees and the wind dispersal sub-model (Figure 6). Where suitable trees are present within a cell or adjacent cell, most beetles will stay within an immediate area. A small proportion will disperse further based on wind direction and strength (wind dispersal sub-model). If there are no suitable trees within the immediate area, beetle dispersion will also be determined by the wind dispersal sub-model which is a function of the number of trees in the initial cell, wind speed and maximum beetle flight time.

To inform wind dispersion, we utilised wind direction and speed data available from the Australian Bureau of Meteorology (BOM). Averaged hourly data for separate weather stations were calculated in eight cardinal directions: N, NE, E, SE, S, SW, W, NW. For each weather station/month/hour/direction combination, the 5th, 25th, 50th, 75th and 95th wind speed percentiles were calculated. These percentiles were then averaged by an hourly dispersion pattern, derived from [20], who found that Monochamus species flight activity was highest between 9 p.m. and 12 a.m.

Within the model, wind dispersion occurred at emergence or the end of the pupae stage of development. The dispersion vectors were randomly sampled from a probability density function (PDF) that was based on the density of healthy trees in a cell, the wind speed and wind direction and the maximum flight time of the beetle. A maximum potential travel distance was used to determine a gamma function for distribution of distances travelled. This distribution was randomly sampled to determine the direction and distance of beetle wind dispersion. Cumulative rainfall was also generated, to restrict dispersal during rain events. A full description of this dataset development is given in our technical report [16].

The number of beetles dispersed to a particular cell was not limited. In cells where beetles were present, the status of a percentage of healthy trees was changed progressively from “healthy” to “healthy but infected” to “stressed and infected” to “dead”. After the pre-oviposition period, the beetles would disperse to either “stressed but not infected” or “stressed and infected” trees.

Numerous studies have reported or modelled spread rates for pinewood nematodes following invasion of a new country [14], which range from 2 to 15 km per year. The rate of spread is influenced by forest density and vector behaviour with the majority of vector flight being shorter distances in contiguous stands and longer distances where hosts are more dispersed [5].

2.9. Study Area

Two study area extents were defined: an area of approximately 160 km × 100 km encompassing the Greater Melbourne area and the State of Victoria (an area of approximately 600 km × 300 km). The Greater Melbourne study area was delineated based on the locations of AAs, as an area would initially be subject to pest incursion response. The datasets were developed statewide to allow the potential spread of a pest species to be modelled across a large region as well as smaller landscapes. This capability was tested to assess model performance using large datasets.

Publicly available Class 1 AA addresses were downloaded from http://agriculture.gov.au/import/arrival/arrangements/sites and converted to a point data layer for use within the model as an incursion source point. To verify the location of these listed businesses, their addresses were cross-referenced to online data sources such as telephone directories and google earth (see [16] for full details of this dataset development).

2.10. Verification of Model

Historical Incursion

To test the model, we simulated a historical incursion using the suspected origin (Docklands) of the nematode incursion in 1999 [21]. We compared the dispersal pattern derived from the model to the 1999/2004 data based on the distance and direction of observed and modelled dispersion.

In late November 1999, dying pine trees were observed near the docks in Melbourne. The cause was eventually identified as a little-known nematode Bursaphelenchus humanensis, which for the sake of this case-study is assumed to have a similar pathogenicity to the more common pinewood nematode Bursaphelenchus xylophilus, vectored by what was assumed to be the Monochamus (Coleoptera: Cerambycidae) beetle. In January 2000, the death of a mature pine tree at the botanic gardens near Williamstown was reported with symptoms similar to those of the pine wilt disease [22]. An initial sampling scheme focused on all trees within a 1 km radius of the Williamstown site, with no nematodes detected. In May 2000, another single dying Pinus spp. was reported, 10 km from the original site this time. The sampling radius was extended. Through aerial surveys and public reporting, a total of 285 trees were sampled of which 31 nematode-infected trees were detected, the furthest 60 km from the original tree (Figure 7). All infected trees were removed and destroyed. The source was presumed to be at the docks. Surveillance was maintained for a further 18 months; however, no further trees were found to have the B. humanensis nematode. Although 31 infected trees were initially detected, no beetles were ever collected, and despite extensive surveys, no further evidence of beetle establishment or infected trees was found.

3. Results

The figures below (Figure 8 and Figure 9) illustrate the outputs of the beetle dispersal model based on the number of adults initially released during September 1999. We ran 2 scenarios using;

(1)

10,000 simulated beetles for comparative analysis against the observed data;

(2)

100 simulated beetles as an indication to what was thought to happen in 1999/2000.

The probability dispersion function is presented in Figure 8, along with examples of dispersal distributions of 100 and 10,000 simulated beetles during October 1999. The Monochamus beetle dispersion probability density function is based on the density of healthy trees in a cell, the average wind speed and wind direction, and the maximum flight time of the beetle. Locations where pine tree death was positively linked with nematode presence based on observed data in Hodda et al. [21] are shown in Figure 9.

The results from modelling 10,000 beetle dispersal events showed strong similarity in both distance and direction from the suspected 1999 origin (release point) (Figure 10 and Figure 11). A comparison between the modelled and observed beetle dispersal as a distance from the central AA (Figure 12) shows no significant difference between the two datasets (two-sample t-test statistic F = 1.74, 28 degree’s of freedom; p Value = 0.02). Similarly, a 2- tailed Mann-Whitney test indicated that the angle from the assumed release point was not significantly different (U = 124,934, p Value = 0.19, significance > 0.05).

4. Discussion

Our forest pest dispersal model accurately simulated the distance and direction of the historic 1999/2000 incursion. In the 100 simulated beetle scenarios, pines were infected with the nematode; however, there was no successful mating of any of the females, and the population died out. This appears to be what happened with the historical incursion as well (D. Smith pers. comm. June 2018).

Our model also enabled us to determine the number of adults in the initial release that would be required to have a potential mating couple within a cell (area: 100 m × 100 m). For example, for 100 adults dispersing on the median October wind speed, there is less than a 20% chance that at least one mating couple (male and female) would disperse to the same cell. If 400 beetles were released, it would be likely that at least three cells would end up with at least one mating couple (Figure 13).

Despite multiple interceptions of M. alternatus at Australia’s border [6], there is no evidence that the beetle has established [5]. Our model not only models the dispersal of M. alternatus, but also the life cycle, population growth, and status of host plants. Although in the modelled historic incursion no successful matings occurred and there was no future generation, the model has the capacity to capture this scenario if plausible. The model can be adapted relatively easily for other forestry pest insect species (e.g., changing the degree day model, flight times) and can account for multi-voltine species (many generations per year).

Although the model utilises historic wind data to produce an average monthly wind speed and direction, the model has the capacity to incorporate real-time as well as predictive climatic data inputs. The model thus has the capability to be used for real-life incursion responses. The dispersal model has only been verified against one species and at one location and time in Victoria. In addition, although the model has used the extensive international research on this species to develop the life cycle and degree day model, this component has not been fully verified.

Although the model accurately simulated the maximum dispersal distance observed in the 1999/2000 incursion (60 km), this dispersal distance is excessive compared to the observed dispersal in Europe and North America, which range from 2 to 15 km per year [14]. Our model appears to have a conformational bias based on the knowledge of the Victorian incursion. Further refinement could include evaluating the model against European incursion information and modifying the gamma dispersal function as needed. However, this is dependent on the availability of suitable high-quality quantitative data on yearly spread of the target pest.

One component which has not been included in the model was upwind movement of the beetles during host searching [23]. The lack of quantitative data of this movement for Monochamus limited our ability to include it in the dispersal sub-module.

The establishment of any new pest into Australia is likely to add substantial costs to industry, as has been the case for the sirex wood wasp, Sirex noctilio [24]. Thus, early detection and initial incursion management is far more cost-effective than managing them once they are well established. Having a predictive dispersal model in place will enhance Victoria’s ability to respond and manage new forestry pest incursions.

Future work and development of the model could focus on validating the life cycle and degree day sub-model, possibly against overseas data and locations. The dispersal sub-model could also be used to illustrate the spread within a high-pine-density region (such as a plantation), or the “escape” of pre-mated beetles from an AA. In addition, modelling other pest species (including pathogen spread) could be investigated or developed for different regions of Australia. This model could also be used to evaluate the success of any forestry pest trapping system and compare the efficiency of different trapping grids (fixed grid, random, clustering around “hot spots” or high-priority-risk areas).

5. Conclusions

This modelling platform has a strong potential to enhance our ability to respond to and manage new forestry pest incursions, in this case, at a local, regional or state level. With minor modifications to parameters within the model, it can be applied to other windborne insect species and pathogens. It can support managers in planning and evaluating their responses to potential incursions of forestry pests and be used to evaluate the effectiveness and efficiency of different trapping systems and grids. It also has the potential to support real-time incursion management and monitoring.