1. Pine wilt disease, caused by pinewood nematode (Bursaphelenchus xylophilus; PWN), is a damaging and globally distributed insect-vectored forest pathogen. Native forest tree mortality associated with PWN is newly reported from the Front Range of Colorado, but there is no regional information on PWN frequency or biology of local insect vectors, limiting management options.

2. A sampling array was established to survey PWN in native pines (Pinus ponderosa) and longhorn beetles (Monochamus clamator & Monochamus scutellatus) over two years and across natural and urban forest landscapes. We developed flight phenology models and evaluated effects of landscape factors on vector abundance and probability of infection.

3. Flight phenology was similar for vectors; Monochamus flight initiated in mid-July and continued into October for both species. We report the first M. clamator–PWN association in the United States. PWN was distributed in the region at rates lower than reported from its putative native range: 3.6 and 4.2% of sampled pines and beetles, respectively, tested positive for PWN. Many host trees were outwardly asymptomatic; infection frequency in tree populations varied considerably and four epicenters of vector infectivity were identified.

4. Epicenters varied in timing of anomalous infective vector frequency—some epicenters had high abundances of infected beetles early in the growing season whereas others had high abundances of infected beetles late in the growing season, though PWN-positive beetles were captured at all sites. Monochamus populations were found primarily in natural forest stands but migrated to urban areas late in the growing season. The only landscape factor positively correlated with abundances of both Monochamus species was distance to previous wildfire.

5. Synthesis and applications: PWN epicenters in the southern Rocky Mountains exhibit specific temporal windows of vector activity that differ from proximal sites. Urban forests, where the disease was initially observed in the region, do not support vector populations. Our results suggest that natural forest landscapes in the region are important reservoirs of PWN and vector populations are especially abundant near burned stands. Collectively, our findings are important for timing disease management activities appropriately and help to distinguish priority areas for mitigation efforts.

Site measurements and sampling of PWN in trees and beetles

A sampling array was established to regionally survey for PWN in host trees and insect vectors. Sites were established in wildland-urban interface (WUI; N=32; Stewart et al. 2007) and urban (N=12) greenspaces (Table S1). Sites were located 50-500 m from roads where ponderosa pine (P. ponderosa) was the dominant canopy species. Urban greenspaces were located in municipalities of Fort Collins, Loveland, Boulder, and Golden (Colorado, USA) where Austrian pine (P. nigra) or Scots pine (P. sylvestris) were dominant canopy species (3 locations in each municipality). Urban greenspaces were selected to be as close to WUI forests (west) as possible while still having >5 trees within the urban greenspace. Sites elevation ranged from 1725-2567 m (GIS, ARCMAP 10.4, ESRI, Inc.).

At WUI sites several forest measurements were conducted on 0.04 hectare fixed-area plots including tree species and diameters-at-breast-height (DBH; 1.3 m), crown-class (suppressed, intermediate, co-dominant, or dominant), presence/absence of visible fungal infection, and fire damage for all trees with DBH > 2 cm (Figure S2-S2). Site aspect and hillslope were recorded. Landscape variables including distance to recently burned stands (km), distance to edge of ponderosa pine canopy (km), percent canopy cover (250 m radius), and distance to nearest urban area (population >2,000;  km) were derived for each study site using a geographic information system (GIS, ARCMAP 10.4, ESRI, Inc.) (USDA Forest Service). Heat-load index (McCune & Keon 2002), a metric of radiative forcing (MJ·cm^-1·year^-1) incorporating slope, aspect, and latitude, was also calculated for each site. These variables were used to develop predictive models of beetle abundance.

To estimate PWN infection frequency in host trees at sample sites, branch and sawdust samples were taken from a subset of 6-10 randomly selected ponderosa pine trees per site (DBH > 10 cm). A 20-cm section proximal to the bole from each of 2 branches was taken from each selected tree using a pole pruner. Sawdust was collected from 2 holes drilled on opposing N and S aspects at 1.3 m height on the bole with an auger-style drill bit (15 mm) to a depth of 6 cm. Tissues samples from each tree were homogenized into a composite sample and nematodes were extracted using the Baermann funnel (Viglierchio and Schmitt 1983); extracted nematodes were stored at -20°C until molecular testing. To sample insect vectors, black crossvane traps were centrally placed in each plot and (described in Morewood et al. 2002) supplied a diffuse pesticide (No. Pest 2 Strips; Dichlorvos; 18.6% 2,2-dichlorovinyl dimethyl phosphate; Hot Shot Corp., St. Louis MO) to kill captured insects. Traps were baited with lures containing host tree volatiles, ethanol, monochamol and ipsdienol (Monochamus lite combo lure - lot #546371; Synergy Semiochemicals, Victoria BC). In 2018 and 2019, traps were visited weekly after commencement of the beetle flight season (July) until flight termination (October) (2018: N=13 weeks; 2019: N=15 weeks). Urban sites were sampled only during 2019. Collected specimens were stored at -20°C until molecular testing to preserve nematode DNA.

All captured Monochamus beetles and Baermann extracts of wood tissues (WUI N=289,  Urban N = 42 samples) were subsequently analyzed for the presence of PWN using a molecular assay. Beetles were bisected longitudinally and homogenized using a sterile micropestle prior to analysis. Baermann extracts of tree tissues and homogenized beetle tissues were tested for PWN using a loop-mediated isothermal amplification (LAMP) assay (Bx Detection Kit, Lot #’s 29000H-L, Nippon Gene Co., Tokyo Japan) according to methods of Kikuchi et al. (2009). Samples were resolved to the individual level (i.e., all sampled trees and captured beetles were tested). This molecular assay is commonly employed in the invasive range of PWN and is 1,000 times more sensitive than traditional PCR approaches (Kikuchi et al., 2009). Presence/absence data was recorded via this assay as opposed to attempting to quantify nematode load/beetle for two reasons; 1) the captured vectors were far too numerous, and 2) the aim of the study was to identify areas where PWN was present rather than evaluate vector competency.

Data analysis

All analyses were performed in the R statistical programming environment and unless otherwise stated use a Type I error rate of α=0.05 for assigning statistical significance (R Core Team 2019).

Flight phenology for M. clamator and M. scutellatus was modeled using a 2-parameter logistic regression (function ‘nplr’) with ordinal day as the independent variable and cumulative proportion of captures as the response variable. Initiation, peak, and termination of flight were approximated using 10%, 50%, and 90% cumulative capture for each species and site × year combination to evaluate differences in phenology between vectors and years (Figure S7, S8). Flight synchrony was estimated by solving growth rate of logistic curves at 50% capture—greater flight synchrony is consistent with more rapid logistic growth (Dell & Davis 2019). Only sites with 10 or more captures recorded for each species were considered reliable for informing species-level flight phenology models (M. clamator N=30 sites, M. scutellatus N=18 sites). Phenology thresholds (mean dates of flight initiation, peak, and termination) and flight synchrony were compared between beetle species using a 2-sample Student’s t-test.

To evaluate effects of landscape factors on vector abundances, beetle capture abundances were modeled for each species using a multiple-regression model selection  using distance to burned stands (burned since year 2000), elevation, heat-load index, distance to ponderosa pine cover boundary, distance to city edge, and canopy cover as predictors. Trap-capture data were root-transformed where necessary in order to meet assumptions of normality and heteroscedasticity. Sample year was included as a random effect, and models were selected via minimization of AIC (Akaike 1974) with a ΔAIC threshold of 2 (function ‘dredge’).

Vector beetle species and sex ratios were compared using Chi-squared tests. The probability of vector association with PWN was evaluated using a log-likelihood mixed-effects modeling approach. Factors considered included both vector species (N=2 factor levels) and sex (N=2 factor levels), as well as day-of-year of capture (continuous effect), and  capture year (N=2 factor levels). Site (N=44) was included as a random effect and evaluated using a likelihood ratio test (P=1).

The recency of first reports of PWN in Colorado indicate that the disease may not be established uniformly in the region, with infections radiating from central locations or (i.e., epicenters). This may be observable via differences in patterns of disease incidence throughout the growing season. Here, disease epicenters were defined as sites with an occurrence of spatiotemporal outliers in the frequency of infection in vector captures during either year (Kitron 1998). Identification of epicenters was made using a scanning statistic to identify sites or aggregate zones where the rate of infection is dissimilar to proximal areas. Epicenters were identified using the function ‘scan_eb_poisson’ with 999 Monte-Carlo iterations (package ‘scanstatistics’). This function computes an expectation-based Poisson scan statistic useful for identifying anomalous spatiotemporal clusters of disease incidence and is commonly employed in human epidemiological studies for a similar purpose (Kulldorff et al. 2005). The method compares all possible temporal windows for each group in each zone list to test a null hypothesis of spatiotemporal randomness using a likelihood ratio statistic. This analysis allows for the possibility of sampling from within a population of vectors twice by considering similarities in infection frequencies between nearest-neighbor groups and across the flight season. The site list used all possible levels of nearest-neighbor combinations for sites grouped within an area while excluding combinations that would include a nearest-neighbor from a geographically discrete (>5 km distance) area based on reported vector flight capacity (Akbulut & Linit 1999; Togashi & Shigesada, 2006). To validate findings, a second model using sites with an identified spatiotemporal anomaly with an interactive site-by-date term was used to test the hypothesis that the likelihood of capturing infected vectors is higher early in the vector flight season. This pattern is consistent with results reported from the southeastern United States where PWN is long-established (Pimentel et al. 2014). Observing significant interactive effects with date would serve as further evidence that previously identified epicenters reflect patterns observed in established systems.

The datasets uploaded here correspond to the following methods section from "Probability of occurrence and phenology of pine wilt disease transmission by insect vectors in the Rocky Mountains (Atkins, Davis, Stewart, 2021). A brief description of the analyses is appended following the methods.

All analyses were performed in the R statistical programming environment and unless otherwise stated use a Type I error rate of α=0.05 for assigning statistical significance (R Core Team 2019).

Flight phenology for M. clamator and M. scutellatus was modeled using a 2-parameter logistic regression (function ‘nplr’) with ordinal day as the independent variable and cumulative proportion of captures as the response variable. Initiation, peak, and termination of flight were approximated using 10%, 50%, and 90% cumulative capture for each species and site × year combination to evaluate differences in phenology between vectors and years (Figure S7, S8). Flight synchrony was estimated by solving growth rate of logistic curves at 50% capture—greater flight synchrony is consistent with more rapid logistic growth (Dell & Davis 2019). Only sites with 10 or more captures recorded for each species were considered reliable for informing species-level flight phenology models (M. clamator N=30 sites, M. scutellatus N=18 sites). Phenology thresholds (mean dates of flight initiation, peak, and termination) and flight synchrony were compared between beetle species using a 2-sample Student’s t-test.

To evaluate effects of landscape factors on vector abundances, beetle capture abundances were modeled for each species using a multiple-regression model selection  using distance to burned stands (burned since year 2000), elevation, heat-load index, distance to ponderosa pine cover boundary, distance to city edge, and canopy cover as predictors. Trap-capture data were root-transformed where necessary in order to meet assumptions of normality and heteroscedasticity. Sample year was included as a random effect, and models were selected via minimization of AIC (Akaike 1974) with a ΔAIC threshold of 2 (function ‘dredge’).

Vector beetle species and sex ratios were compared using Chi-squared tests. The probability of vector association with PWN was evaluated using a log-likelihood mixed-effects modeling approach. Factors considered included both vector species (N=2 factor levels) and sex (N=2 factor levels), as well as day-of-year of capture (continuous effect), and  capture year (N=2 factor levels). Site (N=44) was included as a random effect and evaluated using a likelihood ratio test.

Flight phenology data contained within the file 'All_Flight_10_or_more_per_site.csv'. This file contains cumulative flight capture data for all sites that recorded >10 captures during the 2018 and 2019 field season. Column 'site' indicates the site for which the captures occurred (can be compared to supplemental table for coordinates), the column 'doy' references the day of the year (1-365), and 'cp' is the cumulative proportion of flight captures recorded (0-1). These data were analyzed using the 'nplr'package as described above.

Landscape modeling data is contained within the file 'PWN Landscape data.csv'. the column 'site' again corresponds to the site as listed in the supplemental table. 'total' is the total number of beetles captured across both years, while 'mocl' refers to M. clamator captures and 'mosc' M. scutellatus captures. 'pipo' is the distance to the nearest eastern edge of ponderosa pine cover as estimated using the above described GIS data. 'f09' and 'f16' were distances to fire (km) that had occured from 2000-2009 and 2010-2019, respectively, while 'fire' is distance to burned area of any age. 'city' indicates distance to neared population center (km) derived as described above. 'canavg' is the average canopy cover (%) within 250m of the site as estimated using the above described GIS data. 'x' and 'y' are UTM coordinates, elevation is the square-root of elevation (m), 'year' is the capture year, 'hli' is the heat-load index derived as described in McCune et al 2002, 'pos' is the total abundace of beetles that tested positive for PWN.

All tree-health data are recorded in 'PWN Site-Lvel Data.csv'. The full stem inventory (all growth > 1in DBH) is similarly recorded in the file 'stem inventory.csv'. This includes the survey data containing 'site', as described above, 'spec', the 4-letter FIA abbreviation for each tree specis, 'DBH' diameter at breast height (inches), 'CRN CL' crown class (dominant, co-dominant, intermediate, supressed) CRN DB% crown dieback observed (%), flagging is the presence/absence (1/0) of any dead branch tip sections, 'CANKER' is the presence/absence (1/0) of any fungal cankers on the tree bole, DW. MIST is the presence/absence (1/0) of dwarf mistletoe infection, MECH DMG. is the presence/absence (1/0) of mechanical damage, including large broken limbs or bole scarring, fire is the FIA burn severity rating (0-3), Slope is the slope of the site in degrees, aspect is the aspect of the site in degrees, ASPECT W/C indicates whether a site is a Warm (91-270) or Cold (271-90) aspect. 'PWN Site-Lvel Data.csv' contains only the samples that were screened for PWN, as noted in the PWN (1/0) column.

Beetle screening data are contained in 'beetle.csv' which also contains date and site of capture, beetle species (4-letter code) and whether the beetle was positive/negative (1/0) for PWN.

'map_data_pub_3panel.csv' contains all the input data necessary to create the map from the publication, besides the above referenced GIS layers which are publicly available.