Sea lamprey wounding in Canadian waters of Lake Huron from 2000 to 2009: Temporal changes differ among regions

Journal of Great Lakes Research 37 (2011) 601–608

Contents lists available at SciVerse ScienceDirect

Journal of Great Lakes Research journal homepage: www.elsevier.com/locate/jglr

Sea lamprey wounding in Canadian waters of Lake Huron from 2000 to 2009: Temporal changes differ among regions David V. McLeod a,⁎, R. Adam Cottrill b, Yolanda E. Morbey c a b c

Department of Biology, University of Western Ontario, 1151 Richmond Street N. London, ON, Canada N6A 5B7 R. Adam Cottrill, Ontario Ministry of Natural Resources, Upper Great Lakes Management Unit, 1450 7th Ave. East, Owen Sound, ON, Canada N4K 2Z1 Yolanda E. Morbey, Department of Biology, University of Western Ontario, 1151 Richmond Street N. London, ON, Canada N6A 5B7

a r t i c l e

i n f o

Article history: Received 30 June 2010 Accepted 7 August 2011 Available online 25 September 2011 Communicated by John Janssen Keywords: Sea lamprey Lake Huron Lake whitefish Lake trout Generalized linear model

a b s t r a c t Lake Huron has undergone a number of substantial changes in recent years, including changes to management of the parasitic sea lamprey, Petromyzon marinus. While control strategies of lamprey involving lampricides have had some success, lamprey spawning in St. Marys River has been a major and persistent problem and has led to intensified treatment beginning in 1998. The objective of our study was to broadly examine lamprey spatial wounding dynamics of lake whitefish (Coregonus clupeaformis) and lake trout (Salvelinus namaycush) within the Canadian waters of Lake Huron from 2000 to 2009. Temporal trends were evident and these differed among regions (North Channel, northern Main Basin, southern Main Basin, northern Georgian Bay, and southern Georgian Bay). There was a monotonic annual increase in probability of wounding for both lake trout and lake whitefish in three of the five regions, with high increases seen in both northern and southern Georgian Bay. The increases in three of the five regions are unexpected given the ongoing treatment of St. Marys River. © 2011 International Association for Great Lakes Research. Published by Elsevier B.V. All rights reserved.

Introduction Lake Huron has undergone a number of recent major ecosystem changes, including the precipitous decline in abundance of the benthic macroinvertebrate Diporeia (Nalepa et al., 2007), the collapse of the offshore demersal fish community in 2006 (Riley et al., 2008), and the significant alteration in habitat use by offshore demersal fish (Riley and Adams, 2010). The direct and indirect effects of these changes on large-bodied piscivores have been substantial. Diporeia was a staple food source for lake whitefish (Coregonus clupeaformis) and its decline caused a significant decrease in whitefish condition and growth (Madenjian et al., 2006; McNickle et al., 2006; Rennie et al., 2009). For lake trout (Salvelinus namaycush), who were virtually extirpated from the lake in the 1960s due to sea lamprey (Petromyzon marinus) parasitism (Morse et al., 2003), a rehabilitation effort which has been in place since the 1970s (Eshenroder et al., 1995) appears to be succeeding (Morbey et al., 2008; Riley et al., 2007). However, despite increases in natural reproduction, lake trout underwent a 21% decline in energy density from 1995 to 2004 (Paterson et al., 2009) and showed a decline in growth and condition (He et al., 2008; He and Bence, 2007). Meanwhile, Chinook salmon (Oncorhynchus tshawytscha), a Pacific salmonid which has been stocked in Lake Huron since 1968, has shifted from a fishery that was once dependent entirely on stocking ⁎ Corresponding author. E-mail addresses: [email protected] (D.V. McLeod), [email protected] (R.A. Cottrill), [email protected] (Y.E. Morbey).

to one that is largely self-regulating (Johnson et al., 2010). In spite of these reproductive gains, there has been a decrease in growth and condition of Chinook salmon within Lake Huron in recent years (He et al., 2008). In addition to these ecological changes, revisions have been made to the management of sea lamprey. Lampreys have exacted an enormous toll on lake whitefish and lake trout populations in Lake Huron since their initial invasion (Applegate, 1950). Conventional lamprey management has focused on chemical treatment of the larval stage of lampreys through the application of the lampricide TFM (3-trifluoromethyl-4nitrophenol) to spawning tributaries. The success of TFM in Lake Huron has been limited however by the amount of spawning habitat in St. Marys River, whose discharge rate precludes cost-effective TFM treatment (Shen et al., 2003). By the late 1990s, St. Marys River was estimated to contribute 88% of adult parasitic phase sea lamprey to the lake (Schleen et al., 2003). Development in the mid-1990s of a granular lampricide Bayluscide® (2′,5-dichloro-4′-nitro-salicylanilide) made treatment of St. Marys River feasible, and application of Bayluscide® to the region has been on-going since 1998. Initial post-treatment assessment indicated a reduction in the larval population by 45% (Fodale et al., 2003). Within 5 years, treatment of St. Marys River was expected to reduce the abundance of parasitic phase lamprey in Lake Huron by 60% (Adams et al., 2003). The effect of the St. Marys River treatment on lamprey-induced wounding of lake trout has been assessed for the Drummond Island Refuge. Declining lake trout wound occurrence indicated that the treatment was achieving success (Madenjian et al., 2008a). To our knowledge

0380-1330/$ – see front matter © 2011 International Association for Great Lakes Research. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.jglr.2011.08.003

602

D.V. McLeod et al. / Journal of Great Lakes Research 37 (2011) 601–608

there has been no further assessment of wounding incidence on a broader geographical scale in Lake Huron, especially important in light of the suite of changes in the lake during the past decade. The objective of our study therefore was to examine changes in lamprey wounding rates during the period of 2000–2009 in Canadian waters of Lake Huron. We focused our analyses on lake trout, as done in similar studies (Rutter and Bence, 2003; Sitar et al., 1999), and on lake whitefish, a species of commercial importance that also is commonly parasitized by lamprey. Methods Data selection Our study focused on the Canadian waters of Lake Huron for the period 2000–2009. Wounding data for lake trout and lake whitefish was extracted from commercial catch sampling data provided by the Ontario Ministry of Natural Resources (OMNR). This program is comprised of OMNR contract biologists who randomly sample a proportion of the commercial fishery catch throughout the year and record a variety of information including fish species, fork length, area of capture, and lamprey wounding details (e.g. Milne, 2003). Recording of lamprey-induced wounds on fish by the OMNR follows specific guidelines. Prior to 2003, the number of scars per fish, wounds per fish, wounds ≥ 25 mm, and wounds b 25 mm was recorded. From 2003 onward, the classification system of lamprey wounding was revised based on King and Edsall (1979). Under this system, lamprey wounds qualify as either a definite opening through skin with muscle visible (A-type), or a B-type where a lamprey attack occurred without penetration to muscle. The key uses a number system to indicate the level of healing (I–IV with I being fresh and IV being fully healed; King and Edsall, 1979). This system also provides estimates of wound size. In our study, wounding data were condensed into binary occurrence data. This approach diverges from other studies that used the number of wounds per fish (Rutter and Bence, 2003; Sitar et al., 1999). However, the probability of wounding and mean wound number would both be substantially less than one, and therefore very similar in Lake Huron where fish having more than one wound are rare. For 2000–2002 data, all wounds ≥ 25 mm were assigned a value of “1”; as were all AI-AIII wounds ≥ 25 mm from the 2003–2009 data. These two classifications are considered equivalent as they both reflect the occurrence of recent lamprey attacks. Other lamprey marks (e.g. AIV and BI-BIV) and those records where no wounds were observed were assigned a value of “0”. Observations where lamprey wounding data were not recorded were excluded from analyses. As most of the OMNR commercial fisheries sampling occurs during the July–October period, data were pooled into this standardized time block to ensure sufficient sample sizes and statistical power to detect any changes that may have occurred over the study period. Additionally, using only AI–AIII wounds≥ 25 captured the bulk of wounds suffered during the previous year. While pooling data in this way mixes wounds from the feeding season occurring at the time of capture and from those of the previous year, and prevents attributing wounds to specific year classes of lamprey (see Schneider et al., 1996), it is appropriate for assessing broad interannual trends in sea lamprey wounding levels. Spatial autocorrelation The finest spatial resolution available was the 5"× 5" (approximately 9.25 km× 6.62 km) grid system used by the OMNR for assessment purposes. The number of fish sampled varied widely among grids and years due to fluctuations in the intensity of fishing and monitoring. In order to obtain adequate sample sizes for spatial and temporal comparison, data from twelve assessment areas were used: two areas in each of the North Channel, northern Main Basin and southern Main Basin, and the remaining six in Georgian Bay. Spatial autocorrelation of wounding

was used to determine how best to aggregate the twelve assessment areas. This was done separately for each species based on values of wounding frequency (the number of wounds divided by the total number of fish sampled) for each assessment area. A spatial lattice arrangement was used to test for spatial autocorrelation. This allows irregularly distributed observations to be grouped into larger, arbitrary areas (Zuur et al., 2007). Distances between pairs of assessment areas were calculated as the shortest distances via water between the grids within each assessment area containing the greatest number of observations for lake trout and lake whitefish (referred to hereafter as observation centers; see Fig. 1). The matrix of pair-wise distances was converted to contiguity matrices (W) using an indicator function, l(x):  IðxÞ ¼

1 if x is true 0 if x is not true

ð1Þ

where the matrix W is of dimension 12-by-12 and its ijth element wij is given by  
ðkÞ wij ¼ wij dk ; dkþ1 ¼ I dk b dij ≤ dkþ1

ð2Þ

In this equation, dk are distances {0, 50, 100, 150, 200, 250, 300, 350 km} and indices i and j refer to the observation centers. An incremental distance of 50 km was chosen because it contained sufficient data to test for spatial autocorrelation. Data on wound frequency were tested for spatial autocorrelation among assessment areas using Moran's I Coefficient (Zuur et al., 2007). Values generated by Moran's I Coefficient were converted to z-scores by taking the difference between the value given by Moran's and the value of the expected score divided by the standard deviation of the value given by Moran's. The expected value (E) is given by. E = − 1 / (n − 1), where n is the number of observations (Zuur et al., 2007). Significant spatial autocorrelation of wounding frequency have z-scores N +1.96 (similarity) or b −1.96 (dissimilarity). Models of probability of wounding Probability of wounding was modeled using a generalized linear model (GLM) with a binomial error distribution and logit link function (i.e. logistic regression). We used two categorical variables (species and aggregated areas) and two covariates (fork length and year [adjusted so that 2000 = 0 and 2009 = 9]). Inclusion of fork length and species in the model was done because lampreys exhibit host-selectivity for larger fish over smaller fish (Schneider et al., 1996; Swink, 2003) and show host specificity (Christie and Kolenosky, 1980; Morse et al., 2003). The fully saturated model considered all explanatory variables, all two-way interactions, all four three-way interactions, and one four-way interaction. The form of the logistic regression function applied was log it ðzÞ ¼ ð1 þ exp½−zÞ−1 z ¼ β0 þ β1 x1 þ β2 x2 þ β3 x1 x2 þ ⋯

ð3Þ

where the βi represent the logistic coefficients while the xi are the variables (or higher order transformation of variables) of interest. The effect of the interaction terms (i.e. x1.x2) varied based upon the classification of the variable in question (categorical or continuous). For interaction terms between categorical variables (species and area), the term effectively meant that the intercept of the relationship varied for every combination of the two categorical variables. For interactions between categorical and continuous variables, the slope parameter for the continuous variable is allowed to vary among categories (i.e. different rates of change in probability of wounding as a function of increasing fork length for lake whitefish versus lake trout). For an interaction

D.V. McLeod et al. / Journal of Great Lakes Research 37 (2011) 601–608

603

Fig. 1. Map of Lake Huron showing observation centers used in spatial autocorrelation and the division of areas. Also shown are the locations of TFM-treated tributaries and the number of treatments (denoted by number and shade of circle) during 1998–2008 from Great Lakes Fishery Commission annual reports (e.g. Adair and Young, 2009).

term between two continuous variables (year and fork length) the coefficient represents the multiplicative factor by which the rate of change of probability of wounding with respect to one variable, transformed or otherwise, depended upon the other. Higher-order interactions were removed using likelihood ratio tests in a sequential fashion based upon the highest p-value (N0.05) until only significant interactions and main effects remained (McCullagh and Nelder, 1983). Because a significant three-way interaction involving area (fork length × species× area) remained (p b 0.001) in the initial analysis, separate GLM models were analyzed for each area. In four area-specific models, year was treated as a continuous variable following verification of linearity between wounding rate and year using generalized additive modeling (Hilbe, 2009). For southern Main Basin (see Fig. 1), year did not show linearity with wounding rate and therefore was treated as a factor. Fully-saturated area-specific models were constructed that included all main effects, all two-way interactions (fork length × species, fork length × year, species × year), and one three-way interaction. Non-significant interactions between main effects were sequentially removed starting at the highest pvalue using likelihood ratio tests until only significant interactions and main effects remained (Hilbe, 2009). Although we selected among models by likelihood ratio testing, we note that similar conclusions would be reached by selecting among models based upon AIC (Table 1). A previous study (Rutter and Bence, 2003) observed that marking rates approached an asymptote far below one within the range of observed lengths. This could not occur for a simple logistic regression versus fork length. To allow for nonlinearities, including those that could approximate an asymptotic relationship such as that observed by Rutter and Bence (2003), the method of multivariable fractional polynomials (MFP) was applied (Royston and Altman, 1994; Royston and Sauerbrei, 2008; Sauerbrei and Royston, 1999). Briefly, for each

continuous variable (x) this methodology tests a variety of power (p) transformations (φ(x, p)) to ensure linearity with the link function:  φðx; pÞ ¼

β0 þ β1 xp1 þ β2 xp2 ; p1 ≠p2 β0 þ β1 xp1 þ β2 xp2 lnðxÞ; p1 ¼ p2

ð4Þ

The power values (p) belong to the set p = {−2, −1, −0.5, 0, 0.5, 1, 2, 3}. The βi represent the coefficients of the transformation fitted by maximum likelihood estimation (see Royston and Sauerbrei, 2008). The selection algorithm allows the covariate transformations to increase in complexity by first fitting the most appropriate second degree model using maximum likelihood estimation (by cycling through all appropriate combinations of pi), then sequentially comparing it to the best fitting null, linear, and first degree models. If the pair-wise comparison of models indicates no difference using likelihood ratio testing at a significance level of p = 0.05, the simpler model is chosen as the most appropriate (Royston and Sauerbrei, 2008). As the method makes no underlying assumptions about the relationship between covariates (i.e. year and fork length) and probability of wounding, these transformations will accommodate for an asymptotic or other nonlinear relationship with the link function while avoiding potential biases from pre-specification of a functional form. The method of MFP has been successfully applied in numerous medical studies using regression models (e.g. Kahn et al., 2006; Puhan et al., 2009). To test for interactions, MFP was first applied to a model including main effects only. If transformations of continuous variables were required, then the transformation was divided into its constitutive components (e.g. x p1, x p2) and a full model was refitted using these components and all the interaction terms between these components and the remaining variables in the model (Royston and Sauerbrei, 2008). For example, if a model had only two main effects, both of

604

D.V. McLeod et al. / Journal of Great Lakes Research 37 (2011) 601–608

Table 1 Summary of the stepwise selection of main effects and significant interaction terms for the different regions. AIC and deviance values represent either the selected model with all variables included, or the model with the specified variable eliminated. ΔAIC represents the difference in AIC between focal model and most parsimonious model (i.e. that with the lowest AIC). The likelihood ratio test gives a chi-square value (LR) and associated p-value indicating the contribution of the variable (*** p b 0.001, ** p b 0.01, *p b 0.05). VIF represents the variable inflation factor. Model

Variables

Deviance

AIC

ΔAIC

LR

VIF

North Channel

Final Species Year × Fork length (β1) Year × Fork length (β2) Year × Fork length (β3) Final Fork length Species Year Final Fork length Species Year Final Fork length Species Year (β1) Year (β2) Final Fork length Species Year

4205.9 4264.4 4214.5 4210.9 4212.8 1207.7 1246.4 1222.7 1283.5 913.5 943.2 924.2 934.4 2415.9 2575.6 2431.7 2530.5 2520.6 804.8 824.5 806.0 822.8

4223.9 4280.4 4230.5 4226.9 4228.8 1215.7 1252.4 1228.7 1289.5 921.5 949.2 930.2 940.4 2425.9 2583.6 2439.7 2538.5 2528.6 828.8 846.5 828.0 828.8

– 56.5 6.6 3 4.9 – 36.7 13.0 73.8 – 27.7 8.7 18.9 – 157.7 13.8 112.6 102.7 – 17.7 − 0.8 0.0

– 58.5*** 8.5** 5.0* 6.8** – 38.7*** 15.0*** 75.8*** – 29.7*** 10.7*** 20.9*** – 159.6*** 15.7*** 114.6*** 104.6*** – 19.7*** 1.2 18.0*

– 1.0 – 1.0 1.0 – 2.2 2.4 1.2 – 1.3 1.3 1.1 – 1.1 1.1 1.1 – – 1.5 1.5 1.0

N. Main Basin

N. Georgian Bay

S. Georgian Bay

S. Main Basin

which required second degree transformations, then there would be four interaction terms. Non-significant interactions were sequentially removed starting at the highest p-value using likelihood ratio tests until only significant interactions and main effects (transformed or otherwise) remained (Hilbe, 2009). Final models were checked for multicollinearity between explanatory variables using the variable inflation factor (VIF), computed as 1 / (1 − R 2), where R 2 is the value obtained from the regression of the other variables upon the variable of interest. Variable inflation factor values exceeding five indicated multicollinearity (Table 1; Davis et al., 1986). Goodness-of-fit of final models was first assessed through visual examination of marginal model plots (Pardoe and Cook, 2002). Somers'Dxy, c-index and the le Cessie-van Houwelingen test were used to further examine model fit (Table 2). Somers' Dxy provides comparisons of number of concordant versus discordant pairs. The null hypothesis of the le Cessie-van Houwelingen goodness-offit is that the model fits the data, therefore p-values N 0.10 are considered indicative of good model fit (Hosmer et al., 1997). Finally, the c-index represents the predictive value of the model, with higher values indicating better model fit (taking maximum value of 1; Hilbe, 2009). Models were also examined for overdispersion, which occurs when the observed variance differs from the expected variance, and is indicated by a ratio greater than 1.05 (Table 2; Hilbe, 2009). All statistical analyses were conducted in R (R Development Core Team, 2009).

Results Spatial autocorrelation Wounding frequency for lake trout and lake whitefish showed significant positive spatial autocorrelation (p b 0.05) for data collected within 0–50 km (Fig. 2). Based upon the spatial similarity among geographically proximal areas, the initial twelve assessment areas were aggregated into five regions: North Channel, northern Georgian Bay,

Table 2 Goodness of fit statistics for area-specific models showing the number of observations per model (n), likelihood ratio (LR), c-index, Somers'Dxy and le Cessie-van Houwelingen goodness-of-fit (GOF) p-value, and model dispersion value. Model

n

LR

c-Index

Dxy

GOF

Dispers.

North Channel N. Main Basin N. Georgian Bay S. Georgian Bay S. Main Basin

9,654 5,912 3,812 11,746 5,012

197.4 124.43 69.17 296.28 73.04

0.67 0.75 0.70 0.77 0.74

0.32 0.50 0.39 0.54 0.49

0.22 0.68 0.61 0.99 0.39

1.00 1.01 1.04 0.96 0.98

Fig. 2. Spatial autocorrelation of wounding frequency. Moran's I(k) (converted to z-score) shown for the seven distance ranges (k = 50, 100, 150, 200, 250, 300, and 350 km) and for each species. Solid lines (z-score= 0) represent expected values indicating no spatial autocorrelation, whereas dashed lines represent the upper and lower 95% critical values. Observations falling outside the 95% critical values are considered to be significantly spatially autocorrelated at α = 0.05.

D.V. McLeod et al. / Journal of Great Lakes Research 37 (2011) 601–608

605

southern Georgian Bay, northern Main Basin, and southern Main Basin (Fig. 1).

fork length and year in the North Channel were required and took the following form:

Models of wound probability

ϕðxÞ ¼ β1 x þ β2 x ⋅ lnðxÞ

Of all the fish captured, 2.9% (832/29111) of the lake whitefish and 5.4% (377/7030) of the lake trout had at least one sea lamprey wound. Model selection using multivariable fractional polynomials indicated linear relationships for fork length and year in three of the regions: northern Georgian Bay, northern Main Basin, and southern Main Basin. Year and fork length required transformations for the North Channel, and year required transformation for southern Georgian Bay (Figs. 3 and 4). To account for peaks in the probability of wounding in 2007 and in fish of moderately-large size, transformations of

with a significant interaction term (Table 1) between fork length and year given by

3

3

3 3

3

ð5Þ

3

3

3

ψðx1 ; x2 Þ ¼ β1 x1 x2 þ β2 x1 lnðx1 Þx2 þ β3 x1 lnðx1 Þx2 lnðx2 Þ

ð6Þ

where xi was scaled for year as x1 = (year − 1999)/10 and for fork length as x2 = (fork length/1000). The coefficients β1, β2, β3, are provided in Table 3. To account for asymptotic yearly wounding over the study period, Southern Georgian Bay required year to be transformed as εðxÞ ¼ β1 x

−2

þ β2 x

−2

⋅ lnðxÞ

ð7Þ

where x = (year − 1999) / 10. The coefficients β1, β2 are provided in Table 3. The estimated probability of wounding for lake trout was highest in fish from the North Channel during the late 2000s (Fig. 3). However, the probability of wounding varied among years, between species, and with fork length, making generalizations difficult (Figs. 3 and 4). Differences in the probability of wounding among areas and between species can also be seen by focusing on single years (e.g. 2000 and 2009) and a single fork length (Fig. 5). From 2000 to 2009, the probability of wounding showed a monotonic increase for lake trout sized 533 mm of 88.8% in northern Main Basin, 76.1% in northern Georgian Bay, and 86.0% in southern Georgian Bay (Table 3). The probability of wounding peaked in North Channel in 2007 and subsequently declined (Fig. 3). Over the entire period from 2000 to 2009, there was a 12.1% increase in probability of wounding for lake trout sized 533 mm and a 36.0% increase for lake whitefish sized 500 mm (Table 3). The significant interaction term between fork length and year for North Channel fish meant that the pattern of variation in the probability of wounding among years was different for different sized fish. However, a rise in the probability of wounding from 2000 to 2007 followed by a decline from 2007 to 2009 was apparent for all fish sizes. No change in the probability of wounding across years was apparent in southern Main Basin. Fork length affected the probability of wounding in all of the models (Table 1; Table 3). In four of the five areas, probability of wounding was low for smaller size classes of fish, increasing sharply (i.e. exponentially) after approximately 500 mm fork length (Fig. 4). The North Channel deviated from this pattern: the probability of wounding peaked in medium length fish before subsequently declining in larger-sized fish. In comparing species, lake trout showed a higher probability of wounding than lake whitefish in northern Georgian Bay (+ 38.8% for the year 2000 and fish of 500 mm), southern Georgian Bay (+24.9%), and the North Channel (+40.3%) whereas lake whitefish had greater probability of wounding than lake trout in the northern Main Basin (+51.8%; Table 3). Species did not differ in their probability of wounding in the southern Main Basin (Table 3). Discussion

Fig. 3. Conditional effects of year upon wound probability by area. Line-of-best fit calculated from model where year was treated as a continuous variable (scaled so that 2000 = 0 and 2009 = 9) for lake trout of fork length of 533 mm and individual points represent same model format with year treated as a categorical variable (e.g. allowing it to vary freely without any supposition of a underlying pattern). Φ, ψ, and ε functions are given by Eqs. (5), (6), and (7) respectively using year scaled as (year − 1999) / 10 (i.e. (2005–1999)/10) while fork length (FLEN) was scaled as FLEN/1000.

Sea lamprey numbers vary at both temporal and spatial scales (Young et al., 2003), and large regional variations in probability of wounding are to be expected (Rutter and Bence, 2003). Accordingly, we observed significant spatial autocorrelation of wound frequency of lake whitefish and lake trout at distances of 0–50 km. Despite the large scale treatment of St. Marys River, probability of wounding monotonically increased during 2000–2009 in three of five areas of the Canadian waters of Lake Huron. Furthermore, there was spatial

606

D.V. McLeod et al. / Journal of Great Lakes Research 37 (2011) 601–608

Fig. 4. Conditional effects of fork length upon wound probability by area. Points represent the observed mean wound probability binned by size class for specified year while plotted line is calculated from parameter estimates from Table 4. Mean wounding points were not calculated for size classes with less than 10 fish. Φ, ψ, and ε functions are given by Eqs. (5), (6), and (7) respectively using year scaled as (year − 1999) / 10 (i.e. (2005–1999)/10) while fork length (FLEN) was scaled as FLEN/1000.

variation in how year, fork length, and species affected the probability of wounding. Temporal effects There are a number of likely factors that contributed to the annual increase in probability of wounding in both species. First, it is plausible that there has been an increase in sea lamprey numbers within Lake Huron. Unfortunately, parasitic lamprey population estimation is a difficult task (Bergstedt et al., 2003; Young et al., 2003) and depending upon the magnitude, a population increase or decrease may be impossible to detect. Furthermore, given the large year-toyear variation in estimated numbers (Young et al., 2003) detecting potential trends would be difficult. Second, the recent deterioration in condition of both lake whitefish and lake trout (He et al., 2008; Paterson et al., 2009; Rennie et al., 2009) could have decreased their ability to evade lampreys, leading to a higher incidence of attacks. Furthermore, lampreys spend less time feeding on smaller, poor condition hosts (Bence et al., 2003), potentially resulting in each lamprey engaging in more attacks. Third,

there was a substantial constriction of the depth range frequented by offshore demersal fish, including both lake whitefish and lake trout, during the study period (Riley and Adams, 2010). The potential increase in habitat overlap and decreases in the ranges of these host species may have increased host encounter rates by lamprey. Among regions, some of the greatest annual increases in probability of wounding occurred within Georgian Bay (Table 3). Georgian Bay has historically received relatively low amounts of lampricide treatment (Fig. 1) because tributaries of Georgian Bay have not historically supported a significant spawning population relative to other areas in Lake Huron. However one of the recent major changes in Lake Huron ecology has been the naturalization of the Chinook salmon, the majority of which appear to spawn in southern Georgian Bay tributaries (Johnson et al., 2010). While Chinook salmon continue to be stocked in Lake Huron, survivorship of stocked fish is low and most fish are feral. It is possible that the Chinook salmon are acting as a vector for parasitic phase sea lamprey as the salmon return to spawn in the rivers in the southern Georgian Bay region. If this is the case, lamprey population dynamics within the lake could be undergoing a substantial change, as the presence of parasitic lampreys

D.V. McLeod et al. / Journal of Great Lakes Research 37 (2011) 601–608 Table 3 Parameter estimates for the most parsimonious models of wounding probability (see Table 3). For significant factors (species or year in southern Main Basin), one level was chosen as a baseline (lake trout for species and 2000 for year). Model North Channel

Variable

Intercept Fork length (β1) Fork length (β2) Year (β1) Year (β2) Species(Lake whitefish) Year × Fork length (β1) Year × Fork length (β2) Year × Fork length (β3) N. Main Basin Intercept Fork length Year Species(Lake whitefish) N. Georgian Bay Intercept Fork length Year Species(Lake whitefish) S. Georgian Bay Intercept Fork length Year (β1) Year (β2) Species(Lake whitefish) S. Main Basin Intercept Fork length Species(Lake whitefish) Year (2001) Year (2002) Year (2003) Year (2004) Year (2005) Year (2006) Year (2007) Year (2008) Year (2009)

Estimate

SE

z-value

pvalue

− 13.23 − 57.52 − 210.76 2.02 − 68.05 − 0.92

1.97 14.57 43.48 0.54 21.50 0.11

− 6.72 − 3.95 − 4.85 3.77 − 3.17 − 8.09

b 0.001 b 0.001 b 0.001 b 0.001 0.001 b 0.001

− 12.32 − 370.11 − 1265.50 − 11.37 0.01 0.25 1.14

4.35 164.22 484.04 1.17 0.002 0.03 0.31

− 2.83 − 2.25 − 2.61 − 9.74 6.00 8.33 3.63

0.004 0.02 0.009 b 0.001 b 0.001 b 0.001 b 0.001

− 8.70 0.01 0.17 − 0.83

1.10 0.00 0.04 0.10

− 7.90 5.34 4.46 5.34

b 0.001 b 0.001 b 0.001 b 0.001

− 9.03 0.01 − 0.38 − 0.15 − 0.51

0.60 0.00 0.04 0.02 0.13

− 15.12 12.58 − 9.65 − 9.27 − 4.03

b 0.001 b 0.001 b 0.001 b 0.001 b 0.001

− 8.22 0.01 − 0.38

1.26 0.00 0.34

− 6.54 4.14 − 1.13

b 0.001 b 0.001 0.26

− 0.16 0.23 − 0.96 0.68 − 0.75 − 0.49 0.34 0.42 0.63

0.42 0.60 0.78 0.38 0.67 0.68 0.46 0.48 0.46

− 0.37 0.38 − 1.24 1.78 − 1.11 − 0.73 0.76 0.87 1.38

0.71 0.71 0.22 0.07 0.27 0.46 0.45 0.38 0.17

could increase the incidence of spawning lampreys in Georgian Bay tributaries. The North Channel had the highest probability of wounding in comparison to all other regions. Additionally, lake trout and lake whitefish in this area experienced an annual increase in the odds of wounding from 2000 to 2007 (Fig. 3, Table 3). This increase was

Fig. 5. Estimated probability of wounding of lake trout (fork length = 533 mm) and lake whitefish (fork length = 500) for 2000 and 2009.

607

apparent for fish of all sizes. The observed increase in the probability of wounding in the North Channel was unexpected given the scale and duration of the Bayluscide® treatment of St. Marys River. However, in light of the greater annual increases seen in three of the other regions of Lake Huron, the situation in the North Channel may have been much worse if the St. Marys River was not treated. Effects of fish size Unlike in previous studies (Bence et al., 2003; Rutter and Bence, 2003), there was no asymptote in wounding rates over biologically conceivable fork lengths (e.g., up to 900 mm). However, recently Madenjian et al. (2008a) reported an absence of an asymptote in lamprey wounding for lake trout post-treatment of St. Marys River (2000–2005) for the Drummond Island Refuge. A potential contributing factor to the lack of asymptote in our study could be the lower sample sizes of larger fish. The one exception to the monotonically increasing relationship between fork length and observed wounding rate was the North Channel, where the occurrence of one or more lamprey wounds started at noticeably smaller sizes than in other areas, peaked at approximately 500 mm fork length, and then dropped rapidly with increasing size (Fig. 4). It is unknown why this parabolic relationship existed. The number of individuals in the larger size classes also was lower in the North Channel than in the other areas. Additionally mortality from other sources (i.e. fishing) undoubtedly influenced the lake trout and lake whitefish in the North Channel as substantial commercial and recreational fisheries exist in the North Channel, both of which are likely to influence the demographics of the lake trout and lake whitefish populations there. The small number of large lake trout in the North Channel with relatively few lamprey wounds may be at least partially explained through immigration into the North Channel from other parts of lake. Coded wire tags (CWT) recovered during the study period support this assertion. The majority (50/58) of CWTs recovered from large lake trout (≥ 600) were stocked outside of the North Channel, while the majority (140/223) of lake trout less than 600 mm total length were originally stocked at one of two locations in the North Channel (OMNR, unpublished data). While similar CWT data do not exist for lake whitefish, a lake-wide tagging study of lake whitefish in Lake Huron indicated that there were similar, frequent movements of adult lake whitefish into the North Channel (Ebener et al., 2010). Comparison of lake trout and lake whitefish Species was not a significant predictor of probability of wounding in southern Main Basin (Table 4). In three of the areas and consistent with previous studies (Christie and Kolenosky, 1980; Morse et al., 2003), wounds were more common on lake trout than on lake whitefish. Notably, lake whitefish had higher probability of wounding than lake trout in the northern Main Basin (Table 4). The reason for the discrepancy in the wounding rates between species in this part of the lake is not clear. One contributing factor may by the prevalence of Seneca strain lake trout in northern Lake Huron. Seneca strain lake trout are thought to be less vulnerable to lamprey predation than other strains commonly stocked into the Great Lakes (Eshenroder et al., 1995; Madenjian et al., 2004; Schneider et al., 1996) and as such, large numbers of Seneca strain lake trout are stocked by the US Fish and Wildlife Service in northern Lake Huron (Madenjian et al., 2008a). The prevalence of Seneca strain lake trout is reflected in the proportion of CWTs recovered from the commercial fishery in the northern part of the lake. Between 2000 and 2009, 73% (199/271) of the recovered CWTs from the northern main basin were from Seneca strain lake trout compared to 24% to 32% (n=1342 recovered CWTs) in other region considered in this study (OMNR, unpublished data). Consequently, management agencies should consider using measures that incorporate more than one species, or even community level

608

D.V. McLeod et al. / Journal of Great Lakes Research 37 (2011) 601–608

indicators to assess the effectiveness of lamprey control measures. To date, the vast majority of reports examining lamprey wounding rates in the wild have focused on almost exclusively on lake trout (including backcross and splake) without considering alternative prey species (e.g. Madenjian et al., 2008a; Madenjian et al., 2008b; Sitar et al., 1999).

Conclusions Sea lamprey wounding of lake trout and lake whitefish differed substantially among regions and between species in Canadian waters of Lake Huron for the period of 2000–2009 and a single lake-wide trend does not adequately convey the dynamics of lamprey and their host populations through time. The largest annual increases in probability of wounding occurred in Georgian Bay, a region where wounding has historically been the lowest. As the patterns in probability of wounding differed between species among regions, multiple species or community level indicators should be considered when wounding data is used as a proxy for lamprey abundance.

Acknowledgments We are indebted to the Ontario Ministry of Natural Resources and the Upper Great Lakes Management Unit for provision of the data set. We thank N. Keyghobadi and G. McDonald for useful input into the design of the study and interpretation of the results. We also thank three anonymous reviewers for their helpful comments. Salary support for D.V.M. was provided by an NSERC Discovery Grant to Y.E.M.

References Adair, R.A., Young, R.J., 2009. Integrated management of sea lampreys in the Great Lakes 2008. Annual Report to the Great Lakes Fishery Commission, Ann Arbor, MI. Adams, J.V., Bergstedt, R.A., Christie, G.C., Cuddy, D.W., Fodale, M.F., Heinrich, J.W., Jones, M.L., McDonald, R.B., Mullett, K.M., Young, R.J., 2003. Assessing assessment: can the expected effects of the St. Marys River sea lamprey control strategy be detected? J. Great Lakes Res. 29 (Suppl. 1), 717–727. Applegate, V.C., 1950. Natural history of the sea lamprey (Petromyzon marinus) in Michigan. US Fish and Wildlife Service Special Scientific Report: Fisheries, p. 55. Bence, J.R., Bergstedt, R.A., Christie, G.C., Cochran, P.A., Ebener, M.P., Koonce, J.F., Rutter, M.A., Swink, W.D., 2003. Sea lamprey (Petromyzon marinus) parasite–host interactions in the Great Lakes. J. Great Lakes Res. 29 (Suppl. 1), 253–282. Bergstedt, R.A., McDonald, R.B., Mullett, K.M., Wright, G.M., Swink, W.D., Burnham, K.P., 2003. Mark-recapture population estimates of parasitic sea lampreys (Petromyzon marinus) in Lake Huron. J. Great Lakes Res. 29 (Suppl. 1), 226–239. Christie, W.J., Kolenosky, D.P., 1980. Parasitic phase of the sea lamprey (Petromyzon marinus) in Lake Ontario. Can. J. Fish. Aquat. Sci. 37, 2021–2038. Davis, C.E., Hyde, J.E., Bangdiwala, S.I., Nelson, J.J., 1986. An example of dependencies among variables in a conditional logistic regression. In: Moolgavkar, S.H., Prentice, R.L. (Eds.), Modern Statistical Methods in Chronic Disease Epidemiology. Wiley, New York, pp. 140–147. Ebener, M.P., Brenden, T.O., Wright, G.M., Jones, M.L., Faisal, M., 2010. Spatial and temporal distributions of lake whitefish spawning stocks in Northern lakes Michigan and Huron, 2003–2008. J. Great Lakes Res. 36, 38–51. Eshenroder, R.L., Payne, N.R., Johnson, J.E., Bowen, C., Ebener, M.P., 1995. Lake trout rehabilitation in Lake Huron. J. Great Lakes Res. 21 (Suppl. 1), 108–127. Fodale, M.F., Bergstedt, R.A., Cuddy, D.W., Adams, J.V., Stolyarenko, D.A., 2003. Planning and executing a lampricide treatment of the St. Marys River using georeferenced data. J. Great Lakes Res. 29 (Suppl. 1), 706–716. He, J.X., Bence, J.R., 2007. Modeling annual growth variation using a hierarchical Bayesian approach and the von Bertalanffy growth function with application to lake trout in southern Lake Huron. Trans. Am. Fish. Soc. 136, 318–330. He, J.X., Bence, J.R., Johnson, J.E., Clapp, D.F., Ebener, M.P., 2008. Modeling variation in mass-length relations and condition indices of lake trout and Chinook salmon in Lake Huron: a hierarchical Bayesian approach. Trans. Am. Fish. Soc. 137, 801–817. Hilbe, J.M., 2009. Logistic Regression Models. Chapman & Hall/CRC Press, Boca Raton, FL. Hosmer, D.W., Hosmer, T., le Cessie, S., Lemeshow, S., 1997. A comparison of goodnessof-fit tests for the logistic regression model. Stat. Med. 16, 965–980. Johnson, J.E., DeWitt, S.P., Gonder, D.J.A., 2010. Mass-marking reveals emerging self regulation of the Chinook salmon population in Lake Huron. N. Am. J. Fish. Manage. 30, 518–529. Kahn, J.M., Goss, C.H., Heagerty, P.J., Kramer, A.A., O'Brien, C.R., Rubenfeld, G.D., 2006. Hospital volume and the outcomes of mechanical ventilation. N. Engl. J. Med. 355, 41–50.

King, E.L., Edsall, T.A., 1979. Illustrated field guide for the classification of sea lamprey attack marks on Great Lakes lake trout. Great Lakes Fisheries Commission Special Publication. Madenjian, C.P., Desorcie, T.J., McClain, J.R., Woldt, A.P., Holuszko, J.D., Bowen, C.A., 2004. Status of lake trout rehabilitation on Six Fathom Bank and Yankee Reef in Lake Huron. N. Am. J. Fish. Manage. 24, 1003–1016. Madenjian, C.P., O'Connor, D.V., Pothoven, S.A., Schneeberger, P.J., Rediske, R.R., O'Keefe, J.P., Bergstedt, R.A., Argyle, R.L., Brandt, S.B., 2006. Evaluation of a lake whitefish bioenergetics model. Trans. Am. Fish. Soc. 135, 61–75. Madenjian, C.P., Ebener, M.P., Desorcie, T.J., 2008a. Lake trout population dynamics at Drummond Island Refuge in Lake Huron: implications for future rehabilitation. N. Am. J. Fish. Manage. 28, 979–992. Madenjian, C.P., Chipman, B.D., Marsden, J.E., 2008b. New estimates of lethality of sea lamprey (Petromyzon marinus) attacks on lake trout (Salvelinus namaycush): implications for fisheries management. Can. J. Fish. Aquat. Sci. 65, 535–542. McCullagh, P., Nelder, J.A., 1983. Generalized Linear Models. Chapman and Hall, New York, NY. McNickle, G.G., Rennie, M.D., Sprules, W.G., 2006. Changes in benthic invertebrate communities of South Bay, Lake Huron following invasion by zebra mussels (Dreissena polymorpha), and potential effects on lake whitefish (Coregonus clupeaformis) diet and growth. J. Great Lakes Res. 32, 180–193. Milne, S.W., 2003. Commercial Catch Sampling Program 2002 Summary Report. Upper Great Lakes Management Unit Lake Huron: Ontario Ministry of Natural Resources Queens Printer Report PS_LHA_CF02_001. Morbey, Y.E., Anderson, D.M., Henderson, B.A., 2008. Progress toward the rehabilitation of lake trout (Salvelinus namaycush) in South Bay, Lake Huron. J. Great Lakes Res. 34, 287–300. Morse, T.J., Ebener, M.P., Koon, E.M., Morkert, S.B., Johnson, D.A., Cuddy, D.W., Weisser, J.W., Mullett, K.M., Genovese, J.H., 2003. A case history of sea lamprey control in Lake Huron: 1979 to 1999. J. Great Lakes Res. 29 (Suppl. 1), 599–614. Nalepa, T.F., Fanslow, D.L., Pothoven, S.A., Foley, A.J., Lang, G.A., 2007. Long-term trends in benthic macroinvertebrate populations in Lake Huron over the past four decades. J. Great Lakes Res. 33, 421–436. Pardoe, I., Cook, R.D., 2002. A graphical method for assessing the fit of a logistic regression model. Am. Stat. 56, 263–272. Paterson, G., Whittle, D.M., Drouillard, K.G., Haffner, G.D., 2009. Declining lake trout (Salvelinus namaycush) energy density: are there too many salmonid predators in the Great Lakes? Can. J. Fish. Aquat. Sci. 66, 919–932. Puhan, M.A., Garcia-Aymerich, J., Frey, M., ter Riet, G., Anto, J.M., Agusti, A.G., Gomez, F.P., Rodriguez-Roisin, R., Moons, K.G.M., Kessels, A.G., Held, U., 2009. Expansion of the prognostic assessment of patients with chronic obstructive pulmonary disease: the updated BODE index and the ADO index. Lancet 374, 704–711. R Development Core Team, 2009. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria http://www.Rproject.org. Rennie, M.D., Sprules, W.G., Johnson, T.B., 2009. Factors affecting the growth and condition of lake whitefish (Coregonus clupeaformis). Can. J. Fish. Aquat. Sci. 66, 2096–2108. Riley, S.C., Adams, J.V., 2010. Long-term trends in habitat use of offshore demersal fishes in western Lake Huron suggest large-scale ecosystem change. Trans. Am. Fish. Soc. 139, 1322–1334. Riley, S.C., He, J.X., Johnson, J.E., O'Brien, T.P., Schaeffer, J.S., 2007. Evidence of widespread natural reproduction by lake trout Salvelinus namaycush in the Michigan waters of Lake Huron. J. Great Lakes Res. 33, 917–921. Riley, S.C., Roseman, E.F., Nichols, S.J., O'Brien, T.P., Kiley, C.S., Schaeffer, J.S., 2008. Deepwater demersal fish community collapse in Lake Huron. Trans. Am. Fish. Soc. 137, 1879–1890. Royston, P., Altman, D.G., 1994. Regression using fractional polynomials of continuous covariates: parsimonious parametric modeling. Appl. Stat. 43, 429–467. Royston, P., Sauerbrei, W., 2008. Multivariable Model-building. Wiley, New York, NY. Rutter, M.A., Bence, J.R., 2003. An improved method to estimate sea lamprey wounding rate on hosts with application to lake trout in Lake Huron. J. Great Lakes Res. 29 (Suppl. 1), 320–331. Sauerbrei, W., Royston, P., 1999. Building multivariable prognostic and diagnostic models: transformations of the predictors by using fractional polynomials. J. R. Stat. Soc. A. 162, 71–94. Schleen, L.P., Christie, G.C., Heinrich, J.W., Bergstedt, R.A., Young, R.J., Morse, T.J., Lavis, D.S., Bills, T.D., Johnson, J.E., Ebener, M.P., 2003. Development and implementation of an integrated program for control of sea lampreys in the St. Marys River. J. Great Lakes Res. 29 (Suppl. 1), 677–693. Schneider, C.P., Owens, R.W., Bergstedt, R.A., O'Gorman, R., 1996. Predation by sea lamprey (Petromyzon marinus) on lake trout (Salvelinus namaycush) in southern Lake Ontario, 1982–1992. Can. J. Fish. Aquat. Sci. 53, 1921–1932. Shen, H.T., Xu, Q., Yapa, P.D., 2003. Modeling lampricide transport in the St. Marys River. J. Great Lakes Res. 29 (Suppl. 1), 694–705. Sitar, S.P., Bence, J.R., Johnson, J.E., Ebener, M.P., Taylor, W.W., 1999. Lake trout mortality and abundance in southern Lake Huron. N. Am. J. Fish. Manage. 19, 881–900. Swink, W.D., 2003. Host selection and lethality of attacks by sea lampreys (Petromyzon marinus) in laboratory studies. J. Great Lakes Res. 29 (Suppl. 1), 307–319. Young, R.J., Jones, M.L., Bence, J.R., McDonald, R.B., Mullett, K.M., Bergstedt, R.A., 2003. Estimating parasitic sea lamprey abundance in Lake Huron from heterogeneous data sources. J. Great Lakes Res. 29 (Suppl. 1), 214–225. Zuur, A.F., Ieno, E.N., Smith, G.M., 2007. Analysing Ecological Data. Springer, New York, NY.