Assessing the climate suitability and potential economic impacts of Oak wilt in Canada

Species distribution modelling

A total of 1548 occurrence locations were obtained for B. fagacearum from the United States Forest Service (provided to us by Erin Bullas-Appleton of the Canadian Food Inspection Agency in July 2018) and the Global Biodiversity Information Facility¹⁷. These records were filtered at a 300 arcsecond (approximately 10-km) resolution to remove duplicates and reduce spatial clustering, leaving 1401 unique location records (Fig. 1a). Occurrence locations for two key insect vectors of oak wilt, C. truncatus (Fig. 2a) and C. sayi (Fig. 3a), were obtained from GBIF, publications^10,18, and from specimens in the following collections: Atlantic Forestry Centre, Fredericton, NB, Canada; Canadian National Collection of Insects, Arachnids, and Nematodes, Agriculture and Agri-Food Canada Research Centre, Ottawa, ON, Canada; Ontario Forest Research Institute, Sault Ste. Marie, ON, Canada; Gareth S. Powell Collection, Nephi, UT, USA; Reginald Webster Collection, Charters Settlement, New Brunswick, Canada; and the Florida State Collection of Arthropods, Gainesville, FL, USA . After filtering at a 10-km resolution, there were 82 and 58 unique occurrence records for C. truncatus and C. sayi respectively.

Climate estimates were obtained at each occurrence location by interrogating North American climate models (described in McKenney et al.¹⁹) of the 1981–2010 normal period for the following four variables: (1) fall (i.e., September–December) precipitation (FALLPCP); (2) average spring (i.e., March–June) temperature (SPRINGTMP); (3) annual climate moisture index (CMI; a climate-based moisture balance variable; see Hogg²⁰ for details); and (4) Average Minimum Temperature of the Coldest Month (MINTCM). These climate variables were selected based on their reported influences on B. fagacearum (FALLPCP and CMI^21,22) C. sayi and C. truncatus (SPRINGTMP¹⁰), and insect distributions in general (MINTCM²³). Further, none of the selected variables were highly correlated (i.e., r < 0.7 in all pairwise comparisons), thus alleviating concerns regarding the impact of collinearity among environmental variables on species distribution models²⁴.

Future climate habitat maps for the 2011–2040 and 2041–2070 periods were generated using projections of the four climate variables described above from a composite (i.e., average) of four Earth System Models (ESMs): CanESM2, CESM1CAM5, HadGEM2-ES, and MIROC-ESM (see Price et al.²⁵ for details on these models and the downscaling approach used). All projections employed a moderate greenhouse gas emissions scenario (i.e., RCP4.5; van Vuuren et al.²⁶), which incorporates expected reductions in future greenhouse gas emissions and best describes the path of recent emissions²⁷.

Spatial predictions of the potential distribution of B. fagacearum and its two main insect vectors in North America were generated using Maxent²⁸—a machine learning method that estimates the distribution of a species by finding the distribution of maximum entropy subject to a set of spatial constraints defined by the environmental conditions at the occurrence locations. Maxent employs a regularization parameter (set to the default value of 1 for the current work), which determines the smoothness of the resulting models, and a variety of response functions (i.e., linear, product, quadratic, hinge, threshold, and categorical) to model potentially complex occurrence-environment relationships.

The selection of background points is an important component of a Maxent analysis as this provides a null distribution against which the occurrence locations are compared. Phillips²⁹ recommended that background data points be selected from the same general area as the occurrence observations. However, in our experience, selecting background points that are too close to occurrence locations can also produce flawed results as the machine learning algorithm struggles to distinguish between suitable and unsuitable environmental conditions. Given that we are modelling three species whose geographic limits are unknown (and for which occurrence data is likely incomplete), we felt it was appropriate to consider a somewhat wider environmental domain. Thus, we selected 10,000 random background points from within the treed ecosystems of North America by masking out the Arctic Cordillera, Tundra, and North American Desert ecoregions (based on Level 1 ecoregion definitions by the Commission for Environmental Cooperation³⁰) from the domain used for background point selection.

Assessing the performance of Maxent models has been a somewhat contentious issue. The area under the receiver operator curve (AUC) statistic, which is the default assessment metric provided with the Maxent software, purports to provide a threshold-independent measure of predictive accuracy based on the ranking of locations³¹. This metric, which ranges from 0 to 1, can be interpreted as the probability that a randomly chosen presence location is ranked higher than a randomly chosen background point. However, AUC has been criticized for producing overly optimistic measures of fit when background points are taken from extensive, ecologically unviable locations³². Alternative metrics, such as the True Skill Statistic (TSS³³) require the definition of a suitable threshold to convert Maxent output to a binary (presence/absence) spatial product, and then employ background points to calculate the standard components of a confusion matrix (e.g., true positives, false positives, etc.). This approach has also been criticized, as confusion matrices are more appropriately constructed using true absence data—particularly for species with low sampling effort and/or those lacking a stable geographic distribution³¹. Recently, Wunderlich et al.³⁴ proposed the repurposing of a metric commonly used in meteorology called the Symmetric Extremal Dependence Index (SEDI) to assess the accuracy of species distribution models. This metric, which makes use of logged confusion matrix components and ranges between − 1 and 1, was shown to outperform TSS under a range of modelling outcomes, particularly when background points greatly outnumber occurrence locations as is the case in the current study³⁴. Thus, we present two performance metrics for the current work: AUC (the most commonly reported method for assessing Maxent performance) and SEDI (a promising new performance metric). These metrics were calculated using a random sample of 25% of the occurrence data that was withheld from the model training process. Note that, in order to produce the strongest final distribution maps, final models were run with the full set of occurrence locations; thus the performance metrics based on withheld data should be considered minimum estimates of the predictive accuracy of the model. In order to convert Maxent predictions into binary outcomes for calculating SEDI, we employed the ‘balanced’ threshold, which is provided as a standard output with each Maxent model. Specifically, grid cells with climate suitability less than 0.05, 0.05, and 0.06 were defined as unsuitable for B. fagacearum, C. truncatus, and C. sayi respectively.

As a complimentary modelling approach, climate envelopes were generated for each organism and time period of interest using the ANUCLIM software system^35,36. This system, which represents an early generation, but still useful, tool for species distribution modelling, provides statistical summaries (i.e., min, max, mean, and various percentiles) for each climate variable of interest based on the distribution of values across occurrence locations. Climate envelopes can be defined for the full range of climate conditions at which a species is known to occur (i.e., using min/max values) or for a core range of conditions (e.g., 5th and 95th percentiles). An overall climate envelope is then defined by intersecting the envelopes for each climate variable of interest. For the current work, climate envelopes were defined based on the minimum and maximum values obtained at occurrence locations for each of the four climate variables described above. Final envelopes were overlaid on the Maxent-based climate suitability maps to provide a further assessment of the reliability of these outputs.

To assist in assessing risk to Canadian forests, the geographic range of oak in Canada was delineated by carrying out a geometric union (in ArcGIS v 10.4.1) of the individual ranges of the ten oak species that occur in Canada. For this analysis, we employed digital versions of Little’s (1971) North American tree range maps³⁷. Projections of oak migration under climate change were not incorporated here as we anticipate minimal tree migration by the middle of the current century³⁸.

Potential impacts on street trees

Street tree information was obtained from a survey that has been described previously³⁹. The original survey involves participants walking (or driving) a number of 0.5-km long spatially randomized routes in an urban area. During a survey, each tree within 10 m of the road is identified to genus (or species if possible) and classified according to height (i.e., small ≤ 5 m; medium = 5–10 m; large ≥ 10 m). The number of survey routes in each community was determined such that tree densities could be estimated within reasonable error bounds: ± 5 trees/km for common species and ± 1 tree/km for uncommon species; typically, routes covered 5–10% of total road length in each community³⁹. Canadian Forest Service staff and volunteers employed this approach to collect information on street tree composition for 53 communities in eastern Canada between 2009 and 2015.

A variation on this approach, in which surveys were carried out using Google StreetView imagery, was introduced in 2016. Surveys in 49 communities have since been completed using this approach. Follow-up ground surveys indicated that the StreetView-based surveys were 89% accurate at the genus level and 66% at the species level (unpublished data). These numbers are comparable to other Google-based tree surveys⁴⁰. Note that, for the current work, the genus-level accuracy is most relevant given that B. fagacearum attacks all oak species.

Based on the surveys described above, we collected tree composition and size data for 106 urban centers in eastern Canada. For each of these centers, we calculated oak frequency per kilometer of urban roadway in each size class. Street tree density values were assigned to the remaining urban centers in the study area using an inverse-distance weighted average of the five nearest communities with survey data. Finally, we estimated the total number of host trees in each size class in each community by multiplying the community-level oak density values by the length of urban roadway in each community.

Tree removal and replacement costs (and ranges) were obtained through consultation with tree care professionals and municipal foresters (Table 1). We developed a spreadsheet model to calculate costs for each urban center and for the study area as a whole. Uncertainty in parameter values was explored using @Risk, a spreadsheet add-on that enables detailed Monte Carlo simulations⁴¹. Specifically, variation in tree removal and replacement costs was explored using a triangular distribution with parameters shown in Table 1, while variation in tree replacement was explored using rates of 0%, 50%, and 100%.

Timber-related losses

One approach for placing an economic value on potential forest losses is to employ standing timber (or stumpage) values, which are the fees paid by forest companies (typically to a provincial government in Canada) for the rights to harvest trees on a given land base. To carry out such a valuation, we used national forest attribute maps¹³ to derive gross merchantable oak volumes for four age classes (0–20, 20–40, 40–60, and > 60 years) for each province in our study area. Focusing on current/near-term oak timber stocks (since we are not estimating B. fagacearum spread), we multiplied merchantable volume over 40 years old—roughly the age at which oak becomes harvestable in Ontario⁴²—by average provincial stumpage values.

Estimating stumpage values was somewhat challenging due to inter-provincial variation in stumpage systems and reporting of stumpage fees. For the province of Québec, we obtained oak-specific stumpage fees for 191 harvest zones for the period April 1, 2019 to March 31, 2020 (Bureau de mise en marché⁴³). Since stumpage fees in Québec vary by wood quality class (i.e., A, B, and C), we further obtained information on the proportion of wood harvested in each class over the same period (unpublished dataset, Bureau de mise en marché des Bois). We then calculated the average stumpage fee for the province by averaging across harvest zones and quality classes, while weighting by the proportion of wood in each quality class. For Ontario, we obtained stumpage fees for two quality classes of hardwoods (i.e., Class 1 and Class 2) and four oak-related product types (Veneer, Sawlogs, Composite, and Firewood) for January 1 to December 31, 2019 (Ontario Ministry of Natural Resources and Forestry⁴⁴). Given that oak is typically considered a higher value hardwood, and in lieu of information on how oak is partitioned across product types in Ontario, we calculated oak stumpage fees as an average across product types for the Class 1 hardwood category. Note that the stumpage rates employed here include the Renewal and Forest Futures fees that are part of the Ontario stumpage system. Finally, stumpage values for the province of Nova Scotia were obtained for a single hardwood quality class for the period April 1, 2017 to March 31, 2018 (Province of Nova Scotia⁴⁵). As in Ontario, stumpage fees were averaged across product types. The stumpage values for Nova Scotia were applied to the relatively small amount of oak in the neighbouring Maritime Provinces of New Brunswick and Prince Edward Island (PEI).

Alternatively, gross domestic product (GDP) can provide an estimate of the total economic activity associated with a given industry. Annual GDP estimates for broad categories (e.g., forestry and logging industry, and wood product manufacturing industry) are available for each province at Natural Resources Canada’s forestry statistics website (https://cfs.nrcan.gc.ca/statsprofile/overview/ca). In order to estimate GDP specifically for oak-related timber products, we first multiplied these provincial broad category GDP values by the proportion of the total provincial harvest that was composed of hardwoods (multipliers obtained from published provincial data sources as detailed in the Results section below). This value was then further refined by multiplying by the proportion of hardwoods in the province that was composed of oak species. These estimates were obtained from forest attribute grids¹³, by summing merchantable volume of (1) oak and (2) all broadleaf species within the industrial forestry limits of each province. Spatial summaries were carried out using the raster and rgdal packages in r. Though admittedly coarse, we felt this approach was the best available given the dearth of readily available economic data for individual tree species/genera; similar approaches have been used previously to estimate economic impacts of invasive species^46,47.

These two approaches (i.e., stumpage-based and GDP-based) provide different perspectives on oak-related timber values at risk. The stumpage approach attaches a basic price to standing timber resources, but does not consider downstream economic activities associated with harvest, such as wages, equipment purchases, and capital expenditures. This approach implicitly assumes that substitution possibilities (e.g., other tree species) can fully replace oak-related contributions to the economy with minimal adjustment costs and, as such, is a conservative estimate of potential timber value losses. Alternatively, the GDP approach attempts to include all downstream economic contributions and assumes little or no opportunity for substitution, such that oak timber losses would be accompanied by a proportional reduction in economic activity. We present both estimates here to provide policy-makers with a range of possible impacts. The value of costs through time is generally arrived at using economic discounting; however, here we have no estimates of the timeline associated with oak wilt spread and hence have chosen to report gross, undiscounted values. See Aukema et al.⁴⁸ for further discussion.

Source: Ecology - nature.com