Almost half of the mussel species in North America are imperiled, and eight species found in the eastern Gulf Coastal Plain drainages were recently federally listed. Information regarding the status of known populations of these species is either limited or outdated. Three sites in the Choctawhatchee River watershed (southeast Alabama), where federally threatened mussel species were known to occur, were sampled for mussels eight times each over 4 mo. Three federally threatened species, Fusconaia burkei, Hamiota australis, and Pleurobema strodeanum, and one common species, Elliptio pullata, were individually tagged and released using a robust mark-recapture sampling design. Each species-site combination having sufficient sample sizes was analyzed using a set of six candidate mark-recapture models chosen a priori, and estimates of apparent survival, detectability, and density were derived using the computer program MARK to average models. A total of 820 mussels, 427 of which are listed as federally threatened or endangered, were tagged over eight sampling occasions at three sites. Apparent survival of E. pullata varied among sampling occasions (0.96–0.99), while threatened species tended to have nearly constant survival. Detectability increased with mussel length for E. pullata at all sites (0.07–0.82), but with the exception of P. strodeanum at 8M1, length did not affect detectability of threatened species (0.11–0.52). Densities of threatened species (0.05–1.0 individuals/m2) were typically lower than those of E. pullata (0.15–1.78 individuals/m2) at each site. These data offer insights into the current status of known populations of threatened species at three sites in the Choctawhatchee watershed and will serve as a baseline against which the future status of these populations can be measured. These data also demonstrate the potential viability of using these methods for long-term monitoring of these populations.
Unionid mussels are among the fastest declining groups of freshwater organisms in North America (Vaughn and Taylor 1999). About 30 of the approximately 300 species now recognized are believed extinct. Of those remaining, 65% are either endangered, threatened, or vulnerable (Haag and Williams 2014). In November 2012, eight mussel species endemic to the Gulf drainages of southeast Alabama were listed under the Endangered Species Act (U.S. Fish and Wildlife Service 2012). Margaritifera marrianae Johnson, 1983; Fusconaia rotulata (Wright, 1899); Ptychobranchus jonesi (van der Schalie, 1934); and Obovaria choctawensis (Athearn, 1964) were listed as endangered, and Fusconaia burkei (Walker in Ortmann and Walker, 1922); Fusconaia escambia Clench and Turner, 1956; Hamiota australis (Simpson, 1900); and Pleurobema strodeanum (Wright, 1898) were listed as threatened. Information regarding the local population status of many of these protected species indicates areas where conservation efforts should be focused to prevent further decline and extirpation. Such information is also useful when comparing current population levels to those historically reported and serves as a baseline against which future assessments of these populations can be compared. Further, known populations of threatened mussels may serve as future sources of gravid females for captive propagation of these species, making knowledge of their status of vital importance. The population parameters of interest in this study were apparent survival, detectability, and population size, all of which can be estimated using mark-recapture models (Hart et al. 2001).
Detectability is the probability of detecting (or finding) an individual or species at a site if present (MacKenzie et al. 2006). Imperfect species detectability is a major source of bias in estimating parameters such as population size and survival (Mazerolle et al. 2007). Mark-recapture models are able to account for detectability when estimating population parameters of interest, producing relatively unbiased estimates compared to other methods. Mark-recapture models generally fall under one of two categories: (1) closed population models or (2) open population models. Sampling in a closed population model occurs over a short enough period that changes in population due to births, deaths, or migration are assumed to be negligible (Kendall 1999). Closed population models are primarily used to estimate population size and generally have only one parameter: the probability p that an individual is detected given that it is available for capture (i.e., detectability). In open population models, sampling occurs over a period during which the population is vulnerable to change. Open models have an additional parameter (Φ), which is the probability that a marked individual survives between sampling periods (Amstrup et al. 2010). Open models, such as Cormack-Jolly-Seber, are commonly used to estimate apparent survival and recruitment (Meador et al. 2011). Apparent survival is the probability of a marked individual both surviving in the interval between primary sampling periods and not emigrating from the sample area.
Pollock's robust design is a combination of closed and open population models (Pollock 1982). The robust design consists of a number of secondary sampling occasions (closed), each taking place over a short period (usually consecutive days). These secondary sampling occasions are nested within primary sampling periods (open) with much longer intervals. Survival, population size, and capture probabilities (detectability) can be simultaneously estimated by combining open and closed models.
To estimate parameters such as survival, a model with a given set of assumptions (e.g., constant recapture probability) is chosen. This model is applied to the data, and estimates of the desired parameters, such as population size or survival, are calculated. Estimates will vary depending on the model chosen. To ensure that an appropriate model is used, a set of candidate models are chosen a priori and tested to determine which best explains the underlying process given the data collected. Akaike's Information Criterion is a method of choosing the model that best fits the data using the fewest number of parameters. The most likely model is one that best balances bias (fit of data) and variance (number of parameters) (Burnham and Anderson 2002). Akaike weights can then be used to determine the relative level of support for each model in a given set of candidate models. When no single model is well supported, parameter estimates can be derived from the entire candidate model set using multimodel inference (Burnham and Anderson 2002). A mean of each individual parameter can be calculated using Akaike weights of all models in the candidate set.
The objective of this study was to use a robust mark-recapture design to sample three federally threatened and one common mussel species at three sites in the southeastern Coastal Plains (Choctawhatchee River watershed). We used an initial set of six candidate models to calculate estimates of detectability, apparent survival, and density for the federally threatened species Fusconaia burkei, Hamiota australis, and Pleurobema strodeanum, and a common species, Elliptio pullata. Results were then compared among species and among sites. Few studies have used these methods on mussels in general (Villella et al. 2004; Meador et al. 2011), and an indepth search of the literature failed to reveal any studies that used the robust design on federally threatened mussel species in particular.
METHODS AND MATERIALS
The Choctawhatchee River watershed lies in the Southeastern Plains level III ecoregion and covers about 12,297 km2 (Heath et al. 2010). Three sites from the Choctawhatchee River watershed were selected based on previous knowledge of species composition (Pilarczyk et al. 2006; Reátegui-Zirena et al. 2013) and were sampled from June to October 2012 (Fig. 1). The first site (BS, 31°39′49.6″N, 85°30′18.8″W) was located on the West Fork of the Choctawhatchee River and is a fourth-order stream near Blue Springs State Park, in Barbour County, Alabama. The remaining two sites (8M1, 30°58′50.3″N. 86°10′45.5″W; 8M2 , 30°58′46.7″N, 86°10′45.4″W) were located at Eightmile Creek, a second-order stream in Walton County, Florida.
We captured and tagged mussels between June and October 2012 using Pollock’s robust design (Pollock 1982). We sampled each site during four primary periods, each of which was approximately 1 mo apart. Each primary period consisted of two secondary sampling days as close to each other as possible (usually consecutive). The length of intervals between primary sampling periods varied slightly both within and among sites. A combination of tactile and visual searching was used to sample the stream, generally following Georgia/ Florida qualitative/semiquantitative protocols (Carlson et al. 2008). We did not excavate quadrats because of its time-intensive nature, rarity of the species of interest, and greater disturbance to the habitat. We sampled the same reach for approximately five man-hours at each site per sampling event (Metcalfe-Smith et al. 2000). Mussels were identified to species on site using Williams et al. (2008). All listed mussels encountered in the study sites were tagged during the first sampling event, and all un tagged (not previously collected) individuals encountered on subsequent sampling events were also tagged. Because of their abundance, the first 30 E. pullata encountered within each site were tagged during the first sampling event. Subsequent sampling for this species was limited to the segment of each site where these individuals were found. Untagged E. pullata found within that segment were tagged in subsequent sampling events. Elliptio pullata found outside this segment of each site were counted but were not tagged and were not included in modeling. Listed mussels were tagged throughout the entire site reach. Tag numbers of any previously captured mussels were recorded for each sampling event.
Target mussels were tagged using Hallprint glue-on shellfish tags (available from www.hallprint.com). No juveniles of federally threatened species and few juveniles of E. pullata were captured during this study, which were not tagged due to their small size. Consequently, estimates of apparent survival, detectability, and population size calculated for these populations apply only to adult mussels that exceeded a minimum size (26.2 mm in length). We removed a small section of the periostracum with a file and cleaned and dried the area with a cotton swab and isopropyl alcohol. Tags were attached using forceps and cyanoacrylate glue (superglue). Mussels were kept out of the water for at least 2 min to allow the glue to dry (Lemarie et al. 2000). We returned mussels to the general area from which they were collected by an individual not involved in the collection to avoid bias during subsequent sampling events. None of the current federally threatened species we sampled were listed at the time of sampling, before November 2012 (U.S. Fish and Wildlife Service 2012).
A set of six candidate models were developed for each site and species combination. These models contained the following parameters:
Si = apparent survival during primary period i
γ′ = probability of not being available for capture during primary period given that an individual was not available for capture during primary period i–1 (i.e., the probability of not immigrating back into the study area)
γ″ = probability of not being available for capture during primary period given that an individual was available for capture during period i–1 (i.e., the probability of temporarily emigrating)
Pij = probability of being captured during secondary sampling occasion j of primary period i
Cij = probability of being recaptured during secondary sampling occasion j of primary period i.
Reference table for models used to estimate apparent survival, detectability, and population size for mussel species at three sites in the Choctawhatchee River watershed. Capture probabilities are denoted by p and recapture probabilities are denoted by c, (t) indicates a given parameter varies with time, (·) indicates that a parameter is constant across all sampling periods, and (length) indicates that that a parameter varied with the length of an individual mussel. Models where p and c are equal indicate no behavioral response to initial capture; models where p and c are not equal indicate a behavioral response.
All models assumed that capture probability was constant within a primary period but varied among primary periods (i.e., [p11 = [p12] ≠ [p21 = p22). Temporary emigration was assumed to be constant and random (i.e., γ′[·] = γ″[·]). Finally, all models assumed that apparent survival varied with time (S[t]).
Model 1 was the most general model, allowing both initial capture (p) and recapture (c) probabilities to vary with time between primary periods (interval between primary sampling period) and not be equal to each other between secondary sampling occasions within each primary period (i.e., a behavior response to being captured initially) (Table 1). Model 2 still allowed capture and recapture probabilities to vary with time (interval between primary sampling periods), but they were equal between secondary sampling occasions within a primary period (i.e., no behavior response). Capture and recapture were constant between primary sampling periods in models 3 and 4, but model 3 had no behavior response and model 4 had a behavior response. Finally, models 5 and 6 allowed capture and recapture to vary based on length of individual mussels. Model 5 had no behavior response, and model 6 had a behavior response.
To ensure that models adequately fit the data, some form of goodness-of-fit testing is necessary. In cases where the most general model in a candidate model set does not fit the data, a correction factor, ĉ, can be estimated and applied to the model set (Burnham and Anderson 2002). We tested goodness of fit for each site-species combination by collapsing encounter histories to primary sampling periods and using the live encounter Cormack-Jolly-Seber model in the computer program MARK. One thousand bootstrap simulations were run using a fully time-dependent model. We estimated ĉ using two methods: (1) by dividing the deviance of the time-dependent model by the mean deviance of the simulations and (2) dividing the ĉ of the model by the mean ĉ of the simulations. We used the higher ĉ as a correction factor when analyzing our candidate model set (Table 1). In cases where the estimated ĉ was less than one, we used a ĉ of one for analysis.
We analyzed our candidate model sets using MARK to determine the model with the highest likelihood (Villella et al. 2004; Meador et al. 2011). The same set of candidate models was used for all site-species combinations. Likelihood estimates were based on Akaike's Information Criterion (AIC) (Akaike 1973) modified for small sample sizes (AICc) (Sugiura 1978):
where L() is the likelihood of the parameter estimates, given the data, K is the number of parameters, and n is the sample size. In cases where we used a ĉ correction, likelihood estimates were based on QAICc by dividing the likelihood by the ĉ correction (Burnham and Anderson 2002):
where c is the calculated correction factor (see above). Akaike weights (w) are used to normalize the relative likelihoods of models in a candidate set. For a given candidate set of models, Akaike weights sum to 1 and can be thought of as the weight of evidence for each model.
After analysis of the initial candidate set, we added subsets of the top models where apparent survival was constrained to be constant to see if there was support for survival being constant across our primary sampling periods. For example, if model 3 of a given site-species combination was one of the top models, we ran the same model with constant apparent survival. If this model had a substantially lower AIC, then there was evidence that apparent survival was indeed constant rather than varying among primary sampling periods. This also helped improve our estimates and the precision of parameter estimates in cases where the data were sparse. If these post hoc models were not improvements, we removed them and used the original model set for subsequent analyses.
There is often substantial uncertainty when selecting the best model from a candidate set (Burnham and Anderson 2002). Ignoring this uncertainty can overestimate the precision of parameter estimates. Consequently, we estimated parameters using the model averaging function of MARK to obtain estimates of apparent survival, probability of capture and recapture, and population size during each primary period for each species. To facilitate comparisons among sites, population sizes were converted to density using the length of stream sampled and the mean width of the stream.
Number of tagged mussels by species and site for mussel species at three sites in the Choctawhatchee River watershed.
A total of 820 mussels, 427 of which were either federally threatened or endangered, were tagged over eight sampling occasions at three sites (Table 2). Relative abundance was highest at 8M2 with 462 total tagged mussels and 296 threatened mussels tagged, and lowest at BS (172 total and 38 threatened tagged). Based on poor goodness-of-fit results (possibly due to low recapture rates) we did not analyze H. australis at any site or F. burkei at 8M1. Nor did we analyze F. burkei or O. choctawensis at BS due to low sample sizes (n = 1) (Table 2). Finally, estimates of parameters of P. strodeanum at BS were highly suspect (e.g., high standard errors of population size, poor convergence on estimates of S), even for the top models (i.e., those with the lowest AIC values).
The top models (i.e., those with lowest AIC values) for E. pullata at 8M1 were models 5 and 6. The combined weight of these models (ω ≈ 0.99) strongly suggests that length affected capture probability of this species at this site. Although the top model (5) showed no effect on capture probability due to sampling (i.e., behavior response) (ω ≈ 0.60), the second best model (6) suggested otherwise (ω ≈ 0.40). The high ranking of this second model may be due more to the importance of length as a variable, although there was evidence for a small behavioral response. Constrained versions of models 5 and 6 with constant apparent survival had less support (AIC) than the original models. These two post hoc models were therefore discarded before averaging models for parameter estimation.
The top models for P. strodeanum at 8M1 were 5, 3, and 6. Initial post hoc analysis of models 5, 3, and 6 suggested strong evidence for apparent survival being constant. Therefore, we included a total of six post hoc models corresponding to the original model set but with apparent survival constrained as constant. This resulted in the final model set used for analysis having 12 models. Based on this model set, the strongest evidence was for no behavior effect (ω ≈ 0.77), followed by constant survival (ω ≈ 0.67), and then capture rate varying by length of individual (ω ≈ 0.59).
We found models 5, 3, and 6 to be the top models for E. pullata at 8M2 (ω ≈ 0.84). A post hoc analysis of models 5 and 3 with survival constrained to be constant showed these models had lower QAICc values. However, the estimates of the temporary emigration parameters for these models were suspect (zero with zero standard error), so we discarded these models and used the original model set for subsequent analysis. Using these models, there was evidence for no behavior effect (ω ≈ 0.71). Length had the most effect on capture probabilities (ω ≈ 0.57), and there was very little evidence for a time effect on capture (ω ≈ 0.07).
Model 2 was the only model with any support for both F. burkei (ω ≈ 0.93) and P. strodeanum (ω ≈ 0.99) at 8M2, suggesting that capture and recapture probabilities varied with time and that there was no behavior effect. Models with length as a covariate of capture had the least support. Constraining apparent survival in model 2 to be constant resulted in a lower QAICc in both cases, suggesting that survival was constant over the course of our study. This model was included in the final model set for parameter estimation of both species.
Models 5 and 6 had almost all the support for E. pullata at BS (ω ≈ 0.99), which suggests that length as a covariate affecting capture was the most important variable (ω ≈ 0.99) followed by no behavior effect (ω ≈ 0.56). We ran a post hoc analysis on models 5 and 6, constraining survival to be constant. Both models had substantially higher AICc values and were thus discarded. We used the original model set for parameter estimation.
The lowest apparent survival estimates were 0.96 (95% CI [0.94, 0.97]) for E. pullata at BS during the second sampling interval and 0.98 (95% CI [0.95, 0.98]) for E. pullata at 8M1 during the third sampling interval (Figs. 2 and 3). Otherwise, apparent survival rates for all target species were in the 0.98–0.99 range with generally more precise estimates for threatened species at 8M2 (Fig. 4).
The lowest densities were found at BS. The estimated densities of E. pullata were less than 1.0/m2 on all sampling occasions, and no other species at this site occurred in sufficient numbers to estimate density (Fig. 5). Estimated densities were higher at 8M2 than 8M1 for both E. pullata and P. strodeanum, although estimates at 8M2 had wider confidence intervals (Figs. 6 and 7). Presumably 8M2 had a higher density of F. burkei as well, since this species did not occur in sufficient numbers to estimate density at 8M1. Overall, densities of E. pullata (0.15–1.78) were higher at all sites than any threatened species (0.1–1.0).
Detectability increased with mussel length for E. pullata at all three sites. This relationship was strong at 8M1 and BS but weaker at 8M2 (Figs. 8, 9, and 10). This is due to the lower support for length as a covariate at 8M2 (ω ≈ 0.57) as compared to 8M1 and BS (ω ≈ 0.99). Detectability increased with mussel length for P. strodeanum at 8M1 as well (Fig. 11), although the relationship was also not as strong as E. pullata at 8M1 and BS because of lower support for length as a covariate (ω ≈ 0.59). Overall, detectability ranged from 0.07–0.82 for E. pullata and from 0.11 to 0.52 for threatened species.
Length had no effect on detectability of either P. strodeanum or F. burkei at 8M2. There was strong support that detectability varied with time in these cases (ω ≈ 0.99). Fusconaia burkei appeared to have slightly higher detectability on a given occasion than P. strodeanum (Fig. 12). Detectability was highest on the first and last sampling occasions with a decrease during the two intervening occasions.
Apparent survival is the probability of both surviving between primary sampling periods and not permanently emigrating (i.e., remaining in the superpopulation). Thus, apparent survival estimates will typically underestimate actual survival. As permanent emigration increases, apparent survival estimates become more negatively biased compared to true survival estimates (Gilroy et al. 2012). However, as our estimates were all relatively high (all but two were at least 0.98), permanent emigration was likely not very high at our sites. Given the relatively sedentary nature of mussels, this should not be surprising. Villela et al. (2004) found evidence of very little movement of mussels in their study. Another study found that Elliptio complanata moved a mean of 27 cm downstream over the course of 1 yr (Balfour and Smock 1995). Additionally, Reátegui-Zirena et al. (2013) were able to recover 17% of tagged mussels released 7 yr before from a previous study located at our Eightmile sites.
We found evidence that apparent survival of E. pullata varied with time at all sites, and apparent survival was constant for P. strodeanum at 8M1 and 8M2 and F. burkei at 8M2. Villella et al. (2004) found that apparent survival varied with time in three species of mussels (including two Elliptio species) in the Cacapon River, West Virginia, although their study took place over 3 yr, with sampling occurring throughout the year, while our study took place over only 4 mo. To extrapolate these estimates to annual survival, we would have to assume that apparent survival is the same during the rest of the year in which we did not sample. Given that Villela et al. (2004) found much lower survival estimates during the fall and winter months, this seems unwarranted. Most of the species in our study were in the Tribe Pleurobemini, which often live at least 20 yr (Haag and Rypel 2011). Because of this, apparent survival estimates are likely to be high for all of our species over such a short period, regardless of actual differences in survival. Indeed, the lowest survival estimate we found was 0.95 (95% CI [0.94, 0.97], while most other estimates were at least 0.98. Despite these shortcomings, we were able to successfully estimate apparent survival over the time frame of our study. This suggests that using these methods for long-term monitoring of these populations, with a larger interval between primary periods, could potentially allow more robust estimates of annual apparent survival, allow trends in apparent survival to be tracked over time, and allow differences in apparent survival among species to be found.
Detectability is seldom 100% and may vary over time and space due to choice of sampling method, habitat, weather conditions, experience of collectors, and other variables (Bailey et al. 2004; Wisniewski et al. 2014). Reproductive condition can also affect detectability, as a larger proportion of the adult population is at the surface during the breeding season (Balfour and Smock 1995; Villella et al. 2004). We found that detectability increased with mussel length for E. pullata at all three sites, although this relationship was relatively weak at 8M2. Length seemed less important for P. strodeanum at 8M1 and not important at all for P. strodeanum and F. burkei at 8M2. Meador et al. (2011) found that length had a positive effect on detectability of mussels in all habitat types. Length's lack of influence on detectability at 8M2 maybe due to environmental factors being far more diverse and important at this site, as well as 8M2 being generally more difficult to sample due to depth. Because of rain, water levels were also higher during the second and third primary sampling periods at 8M2, resulting in reduced detectability. Wisniewski et al. (2013) found that detectability of E. nigella decreased with increasing depth, and Schwalb and Pusch (2007) found that the surface density of three species of mussels decreased as river discharge increased due to vertical mussel migration. The time effect on detectability observed at 8M2 may have been a result of these environmental variables (either by causing vertical migration or by decreasing searcher efficiency) and could have masked any length effect (if present) on P. strodeanum and F. burkei and made length less important for E. pullata compared to 8M1 and BS. It is also possible that length is less important for P. strodeanum and F. burkei in general because of their smaller maximum lengths. Elliptio pullata had a maximum length of 81.7 mm in this study and thus a larger range of lengths over which detectability could vary. In contrast, the maximum lengths of P. strodeanum and F. burkei were 57.1 mm and 61.7 mm, respectively. In addition to finding that detectability increased with length, Meador et al. (2011) found that habitat type affected capture. Specifically, detectability was higher in slackwater habitats (zones with low velocities and deposition) than pool and swiftwater habitats. This is consistent with the lower detectability at 8M2, which was deeper with greater habitat heterogeneity, compared to 8M1.
Differences in collector experience may have affected detectability during this study. For example, sites were sampled with four people during the first primary sampling period and with only two people on all subsequent sampling events. Further, although the same two people sampled on all occasions, their efficiency likely increased over time. However, even if detectability varied due to these issues, the final parameter estimates should not have been affected. Detectability in our models was explicitly estimated based on actual data collected. If detectability varied from one primary period to another, parameter estimates such as density and apparent survival accounted for this, even if the causes (collector experience, environmental conditions, etc.) of differences in detectability were not explicitly modeled.
In conclusion, the apparent survivability and detectability calculated from our models serve as repeatable baseline data. These models provide robust data that can be used in longterm monitoring studies to evaluate the status of these species over time. These sites were previously qualitatively sampled. For example, Pilarcyzk et al. (2006), in a single sample, found 13 individuals of Hamiota australis at BS, but the most we found in a single sampling event was five (although they sampled 150 m instead of our 100 m). Pilarcyzk et al. (2006) also found 47 P. strodeanum and 18 O. choctawensis at BS, but we found only 11 P. strodeanum and one O. choctawensis during a single sampling event. This superficially suggests that Eightmile Creek is a stronghold and species are declining at other sites. However, previous data were qualitatively collected and could be biased by detectability (sampling conditions, reproductive condition, collector experience) and are not directly comparable. Using the mark-recapture methods suggested in this study in the future could provide a robust means of determining population trends of these species over time. However, the 4 mo sampling period used in this study is fairly short compared to the long life span of these species. For example, P. strodeanum can live at least 70 yr (Reátegui-Zirnea et al. 2013). Testing these methods over several years is needed to evaluate their efficacy.
We thank Evelyn Reátegui-Zirena and Miluska Olivera Hyde for their assistance in the field. Financial support for this project was provided by the ALFA Fellowship.