DOI: 10.5061/dryad.sn02v6x32
title: Probability of occurrence and phenology of pine wilt disease transmission by insect vectors in the Rocky Mountains

description: 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.
description: <i>Site measurements and sampling of PWN in trees and
 beetles</i> 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 (<i>P. ponderosa</i>) was the
 dominant canopy species. Urban greenspaces were located in municipalities
 of Fort Collins, Loveland, Boulder, and Golden (Colorado, USA) where
 Austrian pine (<i>P. nigra</i>) or Scots pine (<i>P.
 sylvestris</i>) 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 &gt;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 &gt;
 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 &gt;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<sup>-1</sup>·year<sup>-1</sup>)
 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<i> </i>trees per site (DBH &gt; 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
 <i>Monochamus</i> 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. <i>Data
 analysis</i> 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 <i>M.
 clamator</i> and <i>M. scutellatus</i> 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 (<i>M. clamator</i> N=30 sites, <i>M.
 scutellatus</i> 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
 <i>t</i>-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 (<i>P</i>=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 (&gt;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.
description: 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
 &gt;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
 <em>M. clamator</em> captures and 'mosc'
 <em>M. scutellatus</em> 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
 &gt; 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.