Effects of an invasive forest pathogen on abundance of ticks and their vertebrate hosts in a California Lyme disease focus.

Invasive species, including pathogens, can have important effects on local ecosystems, including indirect consequences on native species. This study focuses on the effects of an invasive plant pathogen on a vertebrate community and Ixodes pacificus, the vector of the Lyme disease pathogen (Borrelia burgdorferi) in California. Phytophthora ramorum, the causative agent of sudden oak death, is a non-native pathogen killing trees in California and Oregon. We conducted a multi-year study using a gradient of SOD-caused disturbance to assess the impact on the dusky-footed woodrat (Neotoma fuscipes) and the deer mouse (Peromyscus maniculatus), two reservoir hosts of B. burgdorferi, as well as the impact on the Columbian black-tailed deer (Odocoileus hemionus columbianus) and the western fence lizard (Sceloporus occidentalis), both of which are important hosts for I. pacificus but are not pathogen reservoirs. Abundances of P. maniculatus and S. occidentalis were positively correlated with greater SOD disturbance, whereas N. fuscipes abundance was negatively correlated. We did not find a change in space use by O. hemionus. Our data show that SOD has a positive impact on the density of nymphal ticks, which is expected to increase the risk of human exposure to Lyme disease all else being equal. A positive correlation between SOD disturbance and the density of nymphal ticks was expected given increased abundances of two important hosts: deer mice and western fence lizards. However, further research is needed to integrate the direct effects of SOD on ticks, for example via altered abiotic conditions with host-mediated indirect effects.

Introduction

Invasive species can have strong impacts on the biodiversity, structure, and function of invaded ecosystems (Levine et al. 2003). Much attention has focused on the direct effects of invasive species through processes such as competition (Callaway and Aschehoug 2000; Melgoza et al. 1990; Sakai et al. 2001), predation (Goldschmidt 1996; Sax and Gaines 2008), and facilitation (Simberloff and Von Holle 1999). Invasive species can also exert indirect effects, such as apparent competition (Orrock et al. 2008; Sessions and Kelly 2002) or trophic cascades (Flecker and Townsend 1994; Simberloff and Von Holle 1999) within native ecosystems. Increasing reports of community and ecosystem-level perturbations from invasive species attest to their far-reaching, multi-trophic effects (Braithwaite et al. 1989; Lovett et al. 2006; Mack et al. 2000; Schmitz et al. 1997). Damage from an invasive plant pathogen can cause habitat disturbance comparable to other perturbations such as fragmentation, which has been linked to higher disease risk (Allan et al. 2003; Suzan et al. 2008b). Here we assess the impact of an invasive plant pathogen on four species of native vertebrates and a tick vector that are involved in the transmission cycle of B. burgdorferi, an important zoonotic pathogen in California oak woodlands.

Exotic plant pathogens can drastically alter habitats and whole landscapes (Anagnostakis 1987; Hansen 2008; Lovett et al. 2006). The invasive oomycete P. ramorum is the pathogen responsible for extensive habitat fragmentation and the deaths of hundreds of thousands of oak (Quercus spp.) and tanoak (Lithocarpus densiflorus) trees in California and Oregon (Rizzo and Garbelotto 2003; Rizzo et al. 2002). The pathogen was identified in 2001 (Rizzo et al. 2002), and the disease was named sudden oak death (SOD) due to the apparent rapidity of tree-crown death (McPherson et al. 2001; Rizzo et al. 2002; Rizzo et al. 2005). SOD has now spread to 15 counties in central-coastal California and southern Oregon (Rizzo and Garbelotto 2003; Rizzo et al. 2005), spanning a linear distance of 700 km. To date, 68 species in 40 plant genera are known to be susceptible to P. ramorum (Rizzo and Garbelotto 2003). The species at greatest risk of mortality are the stem-canker hosts, oaks (Quercus) and tanoaks (Lithocarpus; Davidson et al. 2005; Garbelotto et al. 2003; Rizzo et al. 2002).

In California, SOD has caused the greatest damage in mixed oak woodlands composed of California bay (Umbellaria californica) and Quercus spp., although U. californica is not killed by the pathogen (Davidson et al. 2005). Over the course of 14 years, SOD-related mortality rates in California were 4–55% for coast live oak (Quercus agrifolia) and 20–70% for tanoak (L. densiflorus; (Brown 2007; Swiecki and Bernhardt 2002). Relative to control plots, mortality on SOD-impacted plots has doubled for Q. agrifolia and tripled for L. densiflorus (Swiecki and Bernhardt 2002). Although the effects of SOD on plant communities and abiotic factors are being characterized (Brown 2007; Garbelotto et al. 2003; Rizzo and Garbelotto 2003; Rizzo et al. 2002), the impacts of this exotic plant disease on vertebrate communities have not been well studied (but see Apigian 2005; Monahan and Koenig 2006).

Oak woodlands of California and Oregon host a diverse assemblage of vertebrates that are important in food web structure and ecosystem functioning. In California, some of the most abundant and conspicuous vertebrates in oak woodlands are involved in the maintenance of either the tick vector (I. pacificus) or the bacterial agent (Borrelia burgdorferi) of Lyme disease. Columbian black-tailed deer (Odocoileus hemionus columbianus) are important blood-meal hosts for adult I. pacificus. Small mammals such as P. maniculatus and N. fuscipes are reservoirs of B. burgdorferi (Brown and Lane 1992; Peavey and Lane 1995), whereas lizards such as S. occidentalis are important hosts for juvenile stages of I. pacificus. In contrast to the rodents, however, these lizards actively suppress the prevalence of B. burgdorferi by cleansing infected ticks of their infection (Kuo et al. 2000; Lane and Quistad 1998; Wright et al. 1998). Changes in the abundances of these species from SOD could have important implications for Lyme disease risk and public health.

Reports of human diseases emerging directly as a result of habitat disturbance are increasing [e.g., Lyme disease (Allan et al. 2003; Ostfeld et al. 2006), malaria (Diuk-Wasser et al. 2007; Vanwambeke et al. 2007), hantavirus pulmonary syndrome (Lehmer et al. 2008; Suzan et al. 2008a)]. We expected that SOD-related mortality could act similarly to other forest disturbances such as logging or fire, which have been demonstrated to have important impacts on vertebrate populations (Bolger et al. 1997; Lee and Tietje 2005; Suzan et al. 2008a). We assessed whether SOD changed the abundances of key vertebrate hosts and the vectors of Lyme disease with potential consequences for Lyme disease transmission.

Methods

Field sites

Field work was conducted in Marin County, California, north of San Francisco. This area is characterized by a Mediterranean climate with warm, dry summers and mild, wet winters. The first reported case of SOD occurred in Marin County (Rizzo et al. 2002), and it is still one of the areas most severely impacted by SOD (Kelly et al. 2008). Plots were established in two sites, China Camp State Park (CCSP) (elevation 159 m; 38°0′9.50″N; 122°28′2.53″W) in San Rafael and Marin Municipal Water District (MMWD) Sky Oaks headquarters (elevation 430 m; 37°58′5.39″N; 122°36′15.20″W). At each of the two sites, CCSP and MMWD, seven 1 hectare (ha) plots were established. All 14 plots were established in mixed-evergreen continuous canopy forests composed of coast live oak (Q. agrifolia), California bay laurel (U. californica), Pacific madrone (Arbutus menziesii), black oak (Q. kelloggii), Douglas fir (Pseudotsuga menziesii), and occasionally coastal redwood (Sequoia sempervirens).

The presence of P. ramorum infection has been confirmed at both sites (Brown and Allen-Diaz 2009; Kelly et al. 2008; McPherson et al. 2005). Despite the extensiveness of SOD in the far-western US, the local distribution is patchy (Kelly and Meentemeyer 2002; Swiecki and Bernhardt 2002). Plots were established to capture a range of SOD impacts. Several plots were chosen from those established at the two sites by other researchers to monitor SOD (Apigian 2005; McPherson et al. 2005; Brown and Allen-Diaz 2009), and the remainder were selected based on Google Earth (Google, Mountain View, CA, USA) satellite imagery and field ground-truthing. Sampling stations were arrayed in a 10 × 10 grid on each plot, spaced 11 m apart, for vegetation, vertebrate, and tick surveys.

Quantifying SOD impact

A quantitative assessment of the impact of SOD on vegetation was conducted using several metrics commonly used to characterize SOD-impacted forests, including stem composition, stem mortality estimation using stand reconstruction methods, and coarse woody debris (CWD) volume (Waddell 2002; Brown and Allen-Diaz 2009). Stem mortality and CWD were identified (when possible) to genus or species and classified based on decomposition stage. Early decomposition stages were assumed to be recent cases of mortality due to SOD because background mortality levels in oak woodlands are minimal (Brown 2007; Brown and Allen-Diaz 2009; Swiecki and Bernhardt 2002). Within each 1 ha plot, five circular subplots of radius 15 m (0.07 ha) were randomly selected from the 100 sampling stations and assigned for SOD assessment.

Within each subplot, all stems with diameter at breast height (DBH) ≥ 5 cm were identified to species and measured with a diameter tape. These measurements were then converted to basal area. In addition, all dead standing and dead fallen stems that could be identified to genus or species were included in vegetation surveys as dead stems. These data were used to estimate the impact of recent mortality on stem basal area.

Coarse woody debris (CWD), defined as downed and dead woody matter greater than 12.5 cm in diameter, was estimated using five center-radiating transects per subplot, following the protocol detailed in Waddell (2002). Estimates for the entire plot used volumetric estimates of all encountered downed wood.

Canopy cover was estimated with digital photographs taken with a Nikon Coolpix 995 and FC–E8 Fisheye converter. All photos were taken at sunrise or sunset to ensure diffuse light conditions. Thirteen photos were taken at regular intervals throughout each plot (ranging from 23.3 to 30 m apart) to capture the plot canopy. The area covered by each photo differed depending on the canopy height, but all photos were visually checked to assure the absence of overlap. Photos were analyzed with Gap Light Analyzer v.2.0 to calculate canopy openness values (Frazer et al. 1999).

Estimating vertebrate abundances

Mark-recapture methods were employed to estimate the local abundances of lizards and small mammals and the survival rate of N. fuscipes in each plot. Small mammals were sampled once each year for three consecutive nights between March and May from 2006 to 2008 for a total of 9 sampling occasions per plot and 25,200 trap-nights from 2006 to 2008. Extra-large Sherman live traps (7.6 × 9.5 × 30.5 cm; H.B. Sherman Traps, Tallahassee, FL, USA) were pre-baited for three consecutive nights prior to trapping with rolled oats and peanut butter. At each of the plots’ 100 sampling stations, two Sherman traps were set each night. All captured animals were given an individually numbered fingerling eartag (National Band and Tag Co., Newport, KY, USA), identified to species, sex, reproductive status (i.e., lactating, pregnant, nonlactating, testes descended or undescended), and weighed before being released at the point of capture.

Columbian black-tailed deer (Odocoileus hemionus columbianus) use of the plots was estimated using standardized methods of pellet-group counts (>4 pellets in a cluster) (Rowland et al. 1984; White and Eberhardt 1980). Six pellet-sampling subplots were randomly located per plot (4 × 22 m). In the fall of 2007, all subplots were surveyed and cleared of pellets. In the spring of 2008, cleared subplots were re-sampled and new pellet groups were quantified. Pellet-group data were then summed across all subplots to obtain a relative measure of deer use per day on the plots.

Western fence lizard population surveys were conducted along every other plot gridline, for a total of five transects, each spaced 22 m apart. Each transect line was visually surveyed. Sighted lizards were sprayed on their dorsum with a diluted latex paint mixture using an Idico hand tree-marking gun (Idico Products Co., Miami, FL, USA; Eisen et al. 2004). Lizard surveys took place over three consecutive days, with a different paint color used each day to determine a lizard’s encounter history. Each marked lizard was recorded as adult or juvenile; only adults were used in population abundance estimates. Each plot was surveyed once for three consecutive days from April to June of 2007 and 2008, with the exception of one plot (MW4) that was sampled for two days in 2008.

Nymphal tick sampling and infection testing

Lyme disease is vectored by the nymphal stage of I. pacificus, so disease risk focused on assessing the density of nymphal ticks (DON) and the density of infected nymphs (DIN). Questing nymphal ticks were sampled by dragging a 1 m2 flannel cloth along marked transects of woodland understory. Sampling was stopped every 15 m to check and remove ticks from the drag cloth. Sampling was conducted between 10 am and 4 pm during the peak questing activity of juvenile I. pacificus. On each plot, five randomly selected 100 m transects were sampled for ticks. Each plot was sampled twice per year from 2006 to 2008.

All nymphal tick samples were extracted and tested for infection with B. burgdorferi sensu stricto. Due to the high volume of I. pacificus samples, ticks were first screened with a real-time PCR protocol that targets the outer surface protein A (OspA) gene of B. burgdorferi sensu lato. We used an ospA PCR protocol developed for California-specific strains of B. burgdorferi sensu lato that excludes the relapsing fever group of Borrelia (B. coriaceae and B. miyamotoi; Ullmann et al. 2005). To amplify the OspA gene, the forward primer FORospA (5′-AGC AAA ATG TTA GCA GCC TTG AC-3′) and the reverse primer REVospA (5′-CTT TCA TTT CAC CAG GCA AAT CT-3′) were used. An Applied Biosystems (Foster City, CA, USA) Taqman MGB probe, OspA MGB (6FAM-AGA AAA ACA GCG TTT CAG), was developed to maximize the specificity of this assay. Each reaction contained 300 nM FORospA primer, 900 nM REVospA primer, 200 nM Osp-A MGB probe, and 12.5 ul of Taqman Universal master mix (Applied Biosystems, Foster City, CA, USA) for a total reaction volume of 25 ul. Real-time PCR standards were serially diluted from type strain B31 obtained from American Type Culture Company (ATCC, Manassas, VA, USA). Standards ranged from 1 to 1000 genomic copies per 25 ul reaction and were run in triplicate to generate a standard curve on a 7300 real time PCR machine (Applied Biosystems, Foster City, CA, USA). All positive samples were sequenced following the protocol provided in Lane et al. (2005).

Statistical analyses

Local abundances of small mammals and lizards were estimated using Program Mark (White and Burnham 1999). A Huggins closed capture model was used. Each species was analyzed separately by plot and year. Four a priori models were used: (1) capture probability (p) and recapture probability (c) varying with time, (2) p varying with time, c not varying, (3) p not varying with time, c varying with time, and both p and c not varying with time. Model selection was based on Akaike weights, derived from AICc values.

Survival (S) was estimated for N. fuscipes only. Other small mammals (e.g., P. maniculatus) were not recaptured between years, and S. occidentalis lizards were not individually marked and therefore interannual data were not available for them. Survival was estimated using the Huggins closed robust design model in Program Mark (White and Burnham 1999). Six a priori models were used to estimate survival (S): (1–2) Markovian movement with and without variation in S between sampling occasions, (3–4) random movement with and without time variation in S, and (5–6) no movement (i.e., γ′ and γ″ equal to zero) with and without time variation in S. Model selection was based on Akaike weights, derived from AICc.

SOD has varying impacts on the forest. Because animals may be responding uniquely to these impacts, we assessed SOD by quantifying the following vegetation measurements: Quercus spp. basal area (cm3), Quercus spp. mortality (proportion), canopy openness (proportion) and coarse woody debris (CWD) volume (cm3). We then summarized these variables using principal components analysis (PCA) to account for autocorrelation between the variables. Variables were transformed to achieve normality. Quercus spp. basal area was log transformed, and Quercus spp. mortality was arcsine square root transformed, canopy openness was square root transformed, and CWD was log transformed. Transformed variables were then z-score transformed to scale all measurements to a mean of 0 and σ = 1. PCA was performed using a correlational matrix. Components were considered significantly loaded by a variable if the loading values were greater than the absolute value of 0.5 (Hair et al. 1987). All four principal components were retained for generalized linear mixed-effect model analysis (Gauch 1982).

PCA scores were used as independent variables in generalized linear mixed-effect models (GLMM) to assess the impact on P. maniculatus, N. fuscipes, O. hemionus columbianus, and S. occidentalis populations. Vertebrate abundance was regressed against the four principal components generated in the PCA analysis with “site” (CCSP vs. MMWD) included as a fixed effect. Year and plot were included as random effects to account for repeated measures of the same plots over the 3-year sampling period of this study. The analyses for P. maniculatus, N. fuscipes, and S. occidentalis abundances implemented a GLMM with a Poisson error distribution for count data. O. hemionus columbianus plot use and N. fuscipes survival rate utilized a Gaussian error distribution.

GLMM models were also used to assess the response of the density of nymphal ticks (DON) and the density of infected nymphs (DIN). All statistical analyses, unless otherwise noted, were performed in R (Team 2008).

Results

Vegetational responses to SOD

Quercus spp. mortality due to SOD ranged from 6.2 to 75.3%, with a mean prevalence of 28.3% across all plots. Stand structure at CCSP and MMWD was dominated by Q. agrifolia, U. californica, A. menziesii, Q. kelloggii and P. menziesii, each comprising at least 5% of the total basal area at one or both sampling sites. Coarse woody debris (CWD) volume varied tenfold across our plots from 173.75 m3 to 1756.49 m3, with a mean of 781.91 m3 across all plots. Mean canopy openness across all plots and both sites was 13.5%. In general, canopy openness was higher at CCSP than MMWD (Fig. 1).

Fig. 1

Summaries of vegetation data at each site. a Total living basal area (cm2 ha−1) of species by site are shown in quantile box plots. Only species that comprised a minimum of 5% of the total living basal area at each site are included in the summaries. b Boxplots showing quantiles of CWD volume (m3 ha−1) at both sites. c Boxplots and quantiles of canopy openness (%) at the two sites, CCSP and MMWD

Mammal abundances and survival

For small mammals, plots at CCSP were numerically dominated by N. fuscipes (n = 701), whereas plots at the MMWD were dominated by P. maniculatus (n = 1032). Other small mammals captured included California voles (Microtus californicus; n = 45), shrews (Sorex spp.; n = 21), western harvest mice (Reithrodontomys megalotis; n = 9), black rats (Rattus rattus; n = 8), pinyon mice (P. truei; n = 7), and broad-footed moles (Scapanus latimanus; n = 2). The greatest abundance of the two most dominant rodent species occurred in 2006, whereas both species declined substantially in 2007. Abundances rebounded strongly in 2008, particularly for P. maniculatus [Electronic Supplementary Material (ESM) Table 1]. Within CCSP, the plots showed relatively stable numbers of N. fuscipes averaged over the three years of the study versus slightly more variation in P. maniculatus numbers (ESM Table 1). In MMWD, however, average P. maniculatus abundances varied but were consistently high across all plots. N. fuscipes populations were consistently low and occasionally absent (ESM Table 1).

All plots surveyed showed signs of deer use. Pellet group totals on the subplots ranged from 22 to 103, with a mean of 55.4. Deer-use days, a quantitative measurement of average deer presence on each plot, was derived from pellet group counts. The number of deer-use days was calculated based on deer pellet group counts, and these ranged from a low of 33.0 deer use days per year at CC1 to 198.3 deer use days per year at MW7. Mean number of deer-use days per year was 94.3.

The mean annual survivorship rate of N. fuscipes varied little in CCSP (Fig. 2b). In MMWD, survivorship rate was greater in 2006–2007 than in 2007–2008. This was largely due to 0% survival in 2007–2008 in plots MW3 and MW4 (Fig. 2c). However, overall survivorship was near 30% throughout the three years of this study, with a high of 50% at both sites (Fig. 2).

Fig. 2

Survival estimates of N. fuscipes. Survival estimates are shown for two winters at CCSP, summarized by a plot and b year. Survival estimates of N. fuscipes at MMWD summarized by c plot and d year. Survival rates are presented as the proportions surviving from 2006 to 2007 and from 2007 to 2008. Error bars represent ±1 SE

Lizard abundance

S. occidentalis abundances were generally greater at CCSP than MMWD, with mean abundances per plot over the years 2007–2008 of 54.31 ± 8.14 (SE) and 18.13 ± 4.48 (SE), respectively (ESM Table 1). There were no clear demographic trends between sampling years. We observed low year-to-year variation in the mean abundance of lizards per plot: 35.82 ± 7.53 (SE) in 2007 and 36.52 ± 8.93 (SE) in 2008.

SOD impact on vegetation

Principal components analysis (PCA) produced 4 principal components with more than 90% of the variance in the first 2 components (Table 1). A principal component was considered significantly loaded by any input variable if its loading value was greater than |0.5| (Hair et al. 1987). PC1 represented a positive measure of SOD intensity. PC2 was a positive measure of canopy closure. PC3 was a measure of standing Quercus spp. basal area and coarse woody debris (CWD), which is primarily Quercus spp. biomass, and thus can be considered a measure of pre-SOD Quercus spp. presence. Finally, PC4 was a positive measure of Quercus mortality.

Table 1

Loading values for principal components (PCs) 1–4 summarizing vegetation characteristics related to sudden oak death impact

	PCA loading values
	PC1	PC2	PC3	PC4
Proportion of variance ϕ	61.52%	28.60%	6.41%	3.46%
CWD	0.556	0.301	0.660	−0.406
Quercus spp. basal area	−0.543	0.368	0.584	0.478
Quercus spp. mortality	0.569	−0.326	0.126	0.744
Canopy cover	−0.267	0.817	0.456	−0.231
Component summary	SOD	Canopy cover	Pre-SOD oak biomass	Oak mortality

Values in bold indicate variables that dominate a particular principal component (loading > |0.5|). The proportion of the variance for each principal component is shown in the top row, and a summary of each principal component for interpretation is provided

Vertebrate mixed-effects model

The generalized linear mixed-effects model revealed a significant negative effect of Quercus spp. mortality (PC4) on N. fuscipes abundance (Table 2). N. fuscipes survival rate was not correlated with any of our SOD vegetation parameters (Table 2). In contrast to the negative effects of SOD on woodrat abundance, SOD-related factors were significantly correlated with an increase in P. maniculatus abundance (Table 2). SOD intensity (PC1), canopy openness (PC2), and oak mortality (PC4) were all positively correlated with increased mouse abundance (Table 2).

Table 2

Generalized linear mixed-effect model (GLMM) results of principal components 1–4 and the abundances of N. fuscipes, P. maniculatus, and S. occidentalis

	N. fuscipes estimate (P)	P. maniculatus estimate (P)	S. occidentalis estimate (P)	O. hemionus estimate (P)	N. fuscipes survival estimate (P)
Intercept	4.06 (P < 0.001)***	2.05 (P < 0.001)***	3.44 (P < 0.001)***	4.03 (P < 0.001)***	0.27 (P < 0.001)***
PC 1 (SOD)	−0.04 (P = 0.58)	0.10 (P = 0.005)**	0.12 (P = 0.36)	0.08 (P = 0.33)	0.29 (P = 0.18)
PC 2 (canopy cover)	−0.05 (P = 0.70)	−0.17 (P = 0.02)*	−0.61 (P = 0.03)*	0.14 (P = 0.43)	−0.04 (P = 0.23)
PC 3 (pre-SOD Quercus spp. biomass)	−0.09 (P = 0.63)	−0.09 (P = 0.15)	−0.01 (P = 0.98)	0.01 (P = 0.95)	0.03 (P = 0.63)
PC 4 (Quercus spp. mortality)	−0.64 (P < 0.008)**	0.33 (P = 0.02)*	0.06 (P = 0.91)	0.004 (P = 0.99)	−0.15 (P = 0.12)
Site (MMWD)	−3.11 (P < 0.001)***	1.42 (P < 0.001)***	−0.82 (P = 0.14)	−0.18 (P = 0.62)	NA

GLMM results of SOD and O. hemionus plot use, as measured by deer-use days, and N. fuscipes overwinter survival rates are also shown. Results are shown as model estimates and P value. Levels of significance for significant variables are: * ≤ 0.05, ** ≤ 0.01, and *** ≤ 0.001

S. occidentalis abundance estimates were significantly negatively correlated with canopy closure, PC2 (Table 2). Western fence lizard abundances were not significantly affected by other measures of SOD impacts.

Deer-use days showed no correlation with SOD vegetation parameters or their principal components (Table 2). There also did not appear to be any differences in the distributional abundance of pellet groups between sites (Table 2).

The model-predicted abundances were plotted against a single principal component to convey a sense of the response. These responses are displayed in ESM Figs. 1 and 2.

When only considering the range of oak mortality measured on our study plots, GLMM prediction of N. fuscipes abundance was 65% lower on the plots with highest oak mortality compared to the lowest mortality (ESM Fig. 1). The abundance of P. maniculatus is predicted to increase by 83% when comparing the plot with the lowest amount of SOD damage relative to the plot with the highest damage (ESM Fig. 1).

Tick mixed-effect model

A total of 4316 nymphal I. pacificus were collected and tested for infection in this study. From 2006 to 2008, tick infection prevalence with B. burgdorferi sensu stricto was 8.18%. We found a significant negative correlation between the density of nymphal ticks (DON) and canopy cover (PC2), so, conversely, nymphal tick density is positively correlated with canopy openness. When just considering canopy openness, the GLMM model predicted that DON would increase 140% in the plots with the most open canopy relative to plots with the least open canopy (ESM Fig. 2a). DON was also negatively correlated with pre-SOD Quercus spp. biomass (PC3; Table 3). Taken together, these results indicate that the density of nymphal I. pacificus increases with SOD. Similarly, we found a significant negative correlation between the density of infected nymphal ticks (DIN) and pre-SOD Quercus spp. biomass (PC3; Table 3). This negative correlation of pre-SOD oak biomass and DIN suggests an increase in the density of infected ticks as oak biomass decreases, although we did not find a direct correlation of SOD-related forest impact with DIN (Table 3). When regressing the change in DON against the range of canopy openness (PC2) in our data, DON more than doubles from plots with the least to the most canopy openness, while DIN is predicted to nearly triple (ESM Fig. 2b).

Table 3

Generalized linear mixed-effect model results for principal components 1–4 and the density of nymphal (DON) I. pacificus as well as the density of infected nymphal ticks (DIN)

	Density of nymphal ticks (DON) estimate (P)	Density of infected nymphal ticks (DIN) estimate (P)
Intercept	4.48 (<0.001)***	1.51 (<0.001)***
PC 1 (SOD)	−0.01 (0.78)	−0.05 (0.56)
PC 2 (canopy cover)	−0.24 (0.003)**	−0.23 (0.14)
PC 3 (pre-SOD Quercus spp. biomass)	−0.20 (0.02)*	−0.41(0.007)**
PC 4 (Quercus spp. mortality)	−0.03 (<0.84)	−0.25 (0.44)
Site (MMWD)	0.30 (<0.08)	0.61 (<0.054)

Results are shown as model estimates and P values. Levels of significance for significant variables are: * ≤ 0.05, ** ≤ 0.01, and *** ≤ 0.001

Discussion

SOD has affected the stand composition, structure, and resource base in oak woodlands in northern California. The most obvious impact of SOD on Californian woodland habitat is the reduction of up to 75% of Quercus spp. biomass and the loss of oaks as the dominant woodland species. Using 14 1 ha plots that span a gradient of SOD intensity, our study found that SOD has led to an increase in the local abundances of deer mice (P. maniculatus) and western fence lizards (S. occidentalis) and a decrease in dusky-footed woodrats (N. fuscipes). No effect on space use by Columbian black-tailed deer was detected. These species are involved in the ecology of Lyme disease in California: N. fuscipes and P. maniculatus as hosts for immature I. pacificus and as reservoirs of B. burgdorferi (the pathogen causing Lyme disease), and S. occidentalis as an important host of I. pacificus immatures. Changes in the abundances of these species may have a critical impact on the risk of Lyme disease in SOD-impacted forests. In fact, our study found that SOD is positively correlated with the density of nymphs (DON). We also found that the density of infected nymphs (DIN) was negatively correlated with pre-SOD Quercus spp. biomass.

Vertebrate responses

The positive effect of SOD on deer mice abundance is consistent with other studies that found high resiliency of generalist species such as Peromyscus spp. to habitat degradation and disturbance (Allan et al. 2003; Tallmon et al. 2003; Tevis 1956). Mice may benefit from a disturbance that is negatively affecting a resource competitor, in this case, N. fuscipes (Aloiau et al. 1997). Our study found a negative impact of Quercus spp. mortality on N. fuscipes abundance. This negative impact was consistent with studies that find N. fuscipes are generally quite sensitive to habitat disturbance (Castleberry et al. 2001; Tevis 1956) and rely heavily on oaks for food and shelter (Linsdale and Tevis 1951; Tevis 1956).

The increase in S. occidentalis abundance with SOD damage is probably a consequence of changes in abiotic conditions, due to the creation of canopy gaps and an increase in basking sites. Sites with more canopy gaps are characterized by more light penetration, higher temperatures, and lower humidity (Killilea and Swei, unpublished data). In dense-canopied woodlands, basking sites may be an important limiting resource (Adolph 1990). Forest fragmentation from SOD increases the frequency and size of forest gaps, providing more basking sites for lizards. This finding is consistent with prior studies that reported the association of a Sceloporus lizard with timber-cleared sites (Bateman et al. 2008).

We did not find a correlation between plot use of O. hemionus columbianus and any of the vegetation variables. Deer have home ranges larger than the patch size of SOD at our sites (Innes et al. 2009) (Livezey 1991). Due to high vagility, deer probably experience habitat impacts over larger spatial scales than considered in our study, and the effects of SOD on their population densities will require larger plot sizes. However, our data indicate that they are neither concentrating their activity within nor avoiding areas within oak woodlands that are affected by SOD.

Tick response to SOD

We found that DON was positively correlated with SOD. As SOD intensifies, it generates new canopy gaps that provide more forest edges, potentially a suitable microhabitat for ticks (Talleklint-Eisen and Lane 2000). There is evidence from studies on I. scapularis in the northeastern United States that nymphal tick density itself is a reliable indicator of Lyme disease risk (Falco et al. 1999; Stafford et al. 1998), suggesting that Lyme disease risk may rise as SOD progresses. Increased tick density in our plots is likely generated by the increase in lizard abundance that we found in our study. Lizards are the most important host for juvenile I. pacificus, hosting up to 90% (Casher et al. 2002; Eisen et al. 2004) of larval ticks. An increase in the local abundance in lizards, more so than any other vertebrate in this community, should be a strong driver of elevated nymphal tick abundance.

In assessing the impact on Lyme disease risk, we found that DIN was negatively correlated with pre-SOD oak biomass. However, this is not evidence of a direct relationship between SOD-related forest impacts and DIN. Surprisingly, DIN was not affected negatively by canopy cover, as was the case with DON. This discrepancy between how the total density of nymphs (DON) and the density of infected nymphs (DIN) responded to canopy cover suggests that increases in tick density with canopy openness were offset by decreases in tick meals taken from reservoir hosts (rodents), or increases in meals taken from non-reservoirs (lizards). Nevertheless, because we found a negative correlation between DIN and pre-SOD oak biomass, we might expect the continuing SOD epidemic to result in overall increases in Lyme disease risk. Given the important impacts of SOD on the community of reservoir and non-reservoir hosts, it is clear that a much more detailed, mechanistic understanding of the ecological roles played by these vertebrates in supporting ticks and pathogens is needed. Furthermore, data on the direct impact of SOD on ticks and infected tick density via abiotic pathways is needed, and would contribute to our understanding of the impact of SOD-related disturbance on the Lyme disease vector.

The movement of pathogens globally is accelerating, so it will be increasingly important for ecologists to study the potential interactions between diseases. This study attempts to address the issue by looking at the interaction between a devastating invasive pathogen and a human zoonotic disease. We found compelling evidence for disease interaction, and that SOD is affecting the abundances of vertebrates and invertebrates involved in the transmission of B. burgdorferi, with a net positive effect on the tick vector in this oak woodland ecosystem.

Below is the link to the electronic supplementary material.

Supplementary material 1 (DOC 542 kb)