Mussel monitoring data are abundant, but methods for analyzing long-term trends in these data are often uninformative or have low power to detect changes. We used a dynamic occurrence model, which accounted for imperfect species detection in surveys, to assess changes in species occurrence in a long-term data set (1986–2011) for the Tar River basin of North Carolina, USA. Occurrence of all species decreased steadily over the time period studied. Occurrence in 1986 ranged from 0.19 for Utterbackia imbecillis to 0.60 for Fusconaia masoni. Occurrence in 2010–2011 ranged from 0.10 for Lampsilis radiata to 0.40 for F. masoni. The maximum difference between occurrence in 1986 and 2011 was a decline of 0.30 for Alasmidonta undulata. Mean persistence for all species was high (0.97, 95% CI = 0.95–0.99); however, mean colonization probability was very low (<0.01, 95% CI = <0.01–0.01). These results indicate that mussels persisted at sites already occupied but that they have not colonized sites where they had not occurred previously. Our findings highlight the importance of modeling approaches that incorporate imperfect detection in estimating species occurrence and revealing temporal trends to inform conservation planning.
INTRODUCTION
Mussel survey and monitoring data often are not collected or analyzed in a manner that allows strong inference about population trends over time (Downing and Downing 1992; Strayer and Smith 2003). A particular weakness of traditional approaches is an inability to account for imperfect species detection, which is inherent in all survey methods (MacKenzie et al. 2003; Royle and Kery 2007; Dorazio et al. 2010). Species detection has a large random component, and nondetection does not necessarily indicate species absence; failure to account for the probability of detection can lead to faulty conclusions about long-term assemblage changes (MacKenzie et al. 2002, 2003; Royle and Kery 2007; Wisniewski et al. 2013).
More robust analytical approaches that explicitly estimate detection probability, such as occupancy modeling, provide more useful inference from survey data (Meador et al. 2011; Shea et al. 2013; Wisniewski et al. 2013). A related approach, dynamic occurrence modeling, can provide information on the processes that underlie population trends (Betts et al. 2008; Walls et al. 2011; Frey et al. 2012). In addition to accounting for imperfect detection, these models provide estimates of local colonization and extinction probabilities, which are often the focus of long-term monitoring projects (MacKenzie et al. 2003; Royle and Kery 2007).
We used recent and existing survey data to assess changes in mussel occurrence in the Tar River basin, North Carolina, USA, over a 26-yr period with a dynamic occurrence model that estimated and incorporated detection probability. We then estimated persistence and colonization rates to examine processes driving changes in mussel occurrence.
METHODS
Study Sites
We conducted mussel surveys at 20 sites in the Tar River basin and analyzed previous survey data from those sites (Fig. 1). The Tar-Pamlico River basin is the fourth largest basin in North Carolina, with a 14,090-km2 watershed and about 3,790 km of streams. Land use in the basin is primarily forest and wetland with areas of agriculture and urban development (NCDENR 2008). This river basin is among the most species rich of North Carolina, supporting a diverse mussel community of 24 species, 13 of which are imperiled (Bogan 2002). Sites were located among three subbasins (upper Tar River, Swift Creek, and Fishing Creek) and were selected to span a range of environmental conditions and include known occurrences of two U.S. federally endangered species, Alasmidonta heterodon and Elliptio steinstansana.
Mussel Surveys
We conducted timed snorkel and tactile search mussel surveys at all 20 sites during summer 2010 and at eight of those sites during summer 2011. Surveyed stream reaches were accessible at bridge crossings and coincided with known mussel beds or apparently suitable mussel habitat. Reaches extended for 200–500 m, approximately 20 times the average stream width of each site. A minimum of 6 person-hours of effort was expended at each site.
We compiled additional freshwater mussel survey data from the 26-yr period spanning 1986–2011 from the North Carolina Wildlile Resources Commission (NCWRC) database. We included data from surveys that occurred within 500 m upstream or downstream of our 2010–2011 survey sites. We included only survey data that were collected or vetted by NCWRC staff to ensure appropriate survey techniques and accurate species identification. Sampling effort for these surveys was not recorded in the database. When database records included count data, we converted these records to detection/nondetection for each species.
Table 1.
Number of sites surveyed and number of sites with mussel detections for 14 freshwater mussel species in the Tar River basin, North Carolina, USA, from 1986 to 2011.
Our survey data combined with NCWRC database records provided data from 127 surveys among our 20 sites (Table 1). We included all species detection/nondetection data from these surveys with the following exception. We excluded data for three Elliptio species—E. complanata, E. icterina, and E. congaraea—because these species are difficult to identify, and we were uncertain about the consistency of identifications in previous surveys. We did not pool these records (as Elliptio spp.), because E. complanata is extremely common in the Tar River basin and consistently present at all sites, and the lack of absence data would interfere with the ability of the occurrence model to accurately estimate parameters. We included all data for four other Elliptio that are more easily identified: E. fisheriana, E. lanceolata, E. roanokensis, and E. steinstansana.
Dynamic Occurrence Model
We developed a dynamic occurrence model using the detection/nondetection data from all 127 surveys to evaluate changes in species occurrence over the 26-yr period. We adopted a state-space representation of the model wherein we described two component processes: a submodel for the observations conditional on the unobserved state process and a submodel for the unobserved or partially observed state process. We followed an approach developed by Dorazio et al. (2010) and Walls et al. (2011) wherein the single-species model of Royle and Kery (2007) was extended to account for variation in model parameters among ecologically similar species. We modeled the entire species assemblage where each species' individual estimates influence the parameter estimates of every other species in the assemblage and inferences about one particular species are borrowed across all species. Essentially, the parameter estimates for one species are a compromise between the individual species estimates and the mean estimate of those parameters for the assemblage. This is referred to as “shrinkage” in the statistical literature (Gelman et al. 2003) because each species-specific estimate is shrunk in the direction of the estimated mean parameter value. The amount of shrinkage depends on the amount of information for each species and how closely it resembles the overall mean effect for a particular parameter. A major benefit of shrinkage is the ability to estimate parameters for species that are rarely detected. Such species may be critically imperiled species and are an important component of assemblage dynamics, but if analyzed individually there would be too few data to make relevant inference. Instead, these species can be included in the analysis, and species-specific estimates can be obtained. Similarly, parameters may be estimated for years in which data are limited. For example, in years with no survey data, we obtained occurrence estimates by borrowing information from all other years and species through the estimated global mean and the prior distribution that we assigned to that year. In this case, we had noninformative priors indicating that the occurrence probabilities could take any value between 0 and 1, with most of the information drawn from other species and years.
Because sites were not surveyed multiple times per year, we were unable to estimate detection probability directly. However, a related study designed to explicitly estimate detection probability in the same streams found little evidence of species-specific detection rates (Pandolfo et al. 2016). Therefore, we focused on accounting for temporal (i.e., among-year) variability in detection rates. For each year, we randomly drew a single value for detection probability from a normal distribution with a mean of 0.42 (SD = 0.03), estimated previously from 14 Tar River basin mussel species (Pandolfo et al. 2016). This equated to using an informative prior on detection rates and allowed us to accommodate annual variability in detection. Although variable sampling effort across the 26-yr study period would be expected to contribute increased variability in occurrence estimates, our model accounted for this by allowing detection probabilities to vary from year to year. Therefore, unless there has been a systematic trend or bias through time that we have not accounted for explicitly, our estimates reflect the variable sampling effort across the study period.
Following Walls et al. (2011), we then specified a model using the randomly generated detection rates and conditional on the binary occurrence state (detected or not detected). We defined detection state as yikt for each combination of site (k), year (t), and species (i), where each binary observation indicates whether the species was detected (yikt = 1) or not detected (yikt = 0). We defined the occurrence state as zikt for species i, site k, and year t, such that zikt = 1 indicated species presence and zikt = 0 indicated absence of the species. It is noteworthy that if we observed no detections, there is ambiguity in defining the occurrence state because the site could be occupied and we failed to detect the species or the site could be unoccupied. Therefore, we defined the model for each element of the data as follows: yikt | zikt,pt ∼ Bernoulli(ziktpt), where pt denotes the probability of detecting a species in year t given that the species is present. This implies that if the kth site is unoccupied by species i in year t, then yikt = 0 with probability 1 and otherwise the species is detected with probability pt.
We modeled changes in occurrence state for each species by using a first-order Markov process (Royle and Kery 2007). We assumed the initial occurrence state for the ith species at site k is modeled as zik1 | ψi1 ∼ Bernoulli(ψil), where ψi1 denotes the probability of occurrence for species i in year 1. Using a recursive relationship wherein occurrence states in subsequent years (t+1, t+2, ..., T) depend on the occurrence states 1 yr earlier, occurrence in subsequent years can be written as folows: zik,1+1 | ∼ Bernoulli(zik,t,ϕit+γit (1 –zikt)), where γit = Pr(zik,t+1 = 1 | zikt = 0) denotes the probability of local colonization (i.e., a site unoccupied at time t will become occupied at time t + 1), and (ϕit = PR(zik,t+1 = 1 | zikt = 1) denotes the probability of local persistence (i.e., the probability of an occupied site at time t staying occupied at time t + 1). We defined the probability of local extinction, εit, as the probability of an occupied site at time t becoming unoccupied at time t + 1 and defined this as the complement of local persistence probability: εit ≡ 1 – ϕit. Colonization and persistence were fixed between years, with one colonization probability and one persistence probability modeled for each species and applied yearly throughout the 26-yr period.
We used a multivariate normal prior distribution to model species-specific deviations from the mean group-level parameter values (Dorazio et al. 2006; Kery and Royle 2008). We estimated parameters using a Bayesian approach with Markov chain Monte Carlo (MCMC) implemented using R statistical software (with the R2WinBUGS package; Sturtz et al. 2005) and WinBUGS (Lunn et al. 2000) using flat priors for each of the group-level parameters. The MCMC approach allowed us to explicitly measure uncertainty in parameter values by examining a posterior distribution for each parameter. We ran three chains of each model for 20,000 iterations, thinned by 5, after a burn-in of 10,000 iterations (resulting in 12,000 posterior samples for each parameter), and we assessed model convergence by examining trace plots and Gelman–Rubin statistics by using package CODA in R (Gelman et al. 2003). We estimated occurrence, persistence (1 – extinction), and colonization probabilities for the entire mussel assemblage.
RESULTS
Detection probabilities ranged from 0.39 to 0.45 among all species and years. The modeled overall occurrence for all 14 mussel species over 26 yr was 0.35 (95% CI = 0.20–0.51). Initial occurrence rates (1986) ranged from 0.19 for Utterbackia imbecillis to 0.60 for Fusconaia masoni (Table 2). Every species exhibited a decline in occurrence from 1986 to 2011, regardless of initial occurrence (Fig. 2). In 2011, occurrence ranged from 0.10 for Lampsilis radiata to 0.40 for F. masoni. The maximum difference between occurrence rates in 1986 and 2011 was a decline of 0.30 for Alasmidonta undulata. The mean persistence for all species was high (0.97, 95% CI = 0.95–0.99) and ranged from 0.93 for L. radiata to 0.98 for A. heterodon, E. fisheriana, E. lanceolata, E. roanokensis, F. masoni, and V. constricta. However, the mean colonization probability was very low (<0.01, 95% CI = <0.01–0.01). The modeled colonization probability for all 14 species was <0.01.
DISCUSSION
The dynamic modeling approach showed that the occurrence of all 14 mussel species in the study area declined steadily over the 26-yr period. This finding is consistent with qualitative assessments of the region's fauna, which portray steep declines in mussel abundance and species richness (Alderman 1997). Although persistence probability was high among all species for every year in the study, it never reached a value of 1.0. This indicates that every year, there was at least one site where a mussel species was extirpated. Because persistence probability was accounted for annually, the effects of less than total persistence were compounded over the 26-yr period. This, combined with extremely low colonization probabilities, resulted in decreases in mussel occurrence probabilities over time.
Table 2.
Parameter estimates and SDs of occurrence, persistence, and colonization probabilities for 14 freshwater mussel species in the Tar River basin, North Carolina, USA, from 1986 to 2011.
Figure 2.
Estimated occurrence (Ψ) for 14 freshwater mussel species in the Tar River basin, North Carolina USA, from 1986 to 2011.
A major advantage of our modeling approach was that it incorporated imperfect detection of mussels, and therefore we avoided basing inferences on biased data (MacKenzie et al. 2003; Kery and Schmidt 2008). For example, a cursory inspection of survey data (Table 1) might suggest that the federally endangered A. heterodon occurred at survey sites more frequently in later time periods than in earlier periods. However, dynamic modeling indicated that occurrence of A. heterodon declined steadily, and this finding is in agreement with regional mussel experts who infer that populations of A. heterodon are declining in the Tar River basin (R.B. Nichols, North Carolina Wildlife Resources Commission, personal communication). The apparent increase in the number of surveys in which A. heterodon was detected is most likely related to sampling effort and intensity or species detection issues. The species was federally listed as endangered in 1990, which led to expanded research interest in this species (Strayer et al. 1996). It is likely that survey efforts became more frequent and targeted for this species and surveyors became more adept at locating and identifying this small mussel (Strayer and Smith 2003; Meador et al. 2011). The apparent increase in species detections at survey sites could be misinterpreted as an increase in occurrence and perhaps misinform conservation planning.
Despite our finding of a steady decrease in mussel occurrence, the high persistence probabilities in the Tar River basin indicate that the majority of sites where mussels previously occurred have remained occupied. However, our analysis did not include a measure of abundance, and we have no information about changes in population size. Similarly, our analysis did not include individual size and provides no information about other changes in population status such as age structure or strength of recent recruitment.
The extremely low colonization probabilities that we modeled for mussels in the Tar River basin cannot be conclusively attributed to a particular cause. Other studies found that colonization rate depends on mussel density and distribution and larval dispersal traits (McClain and Ross 2005; Vaughn 2012). Habitat requirements of newly settled and established juvenile mussels may also influence colonization success. For example, juvenile mussels may not be able to settle during high velocity or shear stress conditions (Payne and Miller 2000). We have no information about changes in mussel density or habitat conditions that may influence colonization in the Tar River basin. However, it is likely that, to a large extent, unoccupied sites are not being colonized simply because mussels have reached a steady state in which all species suited for the habitat and resources at a particular site have already colonized that site.
Our analyses were constrained by the data available in the existing database. Most records in the database did not report sampling effort, which limited our analyses to presence/absence. Including sampling effort in survey records can greatly enhance their utility by allowing assessment of long-term trends in abundance. Also, our data were not collected specifically for application in an occurrence model, so results should be interpreted with a degree of caution. For example, methods of estimating detection probability are data intensive, and we were unable to model it using survey data from the database. Instead, we relied on detection probability estimates derived from a complementary study in the same river basin (Pandolfo et al. 2016). These detection probabilities were estimated for the same species at the same sites examined in this study, and our model included annual variation around the mean detection probability from the complementary study. However, we were unable to empirically estimate changes in detection probability over time. Therefore, our modeled occurrence probabilities over the study period may not reflect actual changes in detection probability that may have occurred during that time. For example, our model cannot conclusively address the presumption that detection of A. heterodon has increased during the study period due to increased focus on this species (see above). In addition, because our modeling approach uses shrinkage (i.e., it borrows data across all species), the parameter estimates are influenced by the mean estimates for the entire assemblage (Gelman et al. 2003). Thus, if one species is more data rich than others in the assemblage, it may influence the parameter estimates of the other species.
Despite the limitations of our data set, our dynamic occurrence modeling approach incorporated imperfect detection to generate parameter estimates for an entire mussel assemblage, including rare species that are more data limited. This approach enabled us to document gradual declines in occurrence for all species in this region since 1986. The specific causes of these declines are unknown, but species life history traits, agricultural land use, and stream power influence occurrence of mussels in this region (Pandolfo et al. 2016). In addition, this region is experiencing intensive climate and land use changes, rendering the aquatic fauna vulnerable (Ingram et al. 2013). A wealth of mussel survey data exist among individuals, agencies, and universities, and more effective analytical approaches can increase the value of these data for assessing long-term trends in mussel populations and the causes of mussel declines (Burgman et al. 1995; Reichman et al. 2011). Our ability to assess long-term trends can be enhanced further by recording sampling effort and population measures such as individual size in future surveys.
ACKNOWLEDGMENTS
We thank T. Black, B. Jones, E. Buttermore, D. Weaver, J. Archambault, A. White, B. Cope, and M. Fisk for mussel field survey assistance and support and A. Drew for a constructive manuscript review. This research was funded by the U.S. Geological Survey, National Climate Change and Wildlife Science Center, through Research Work Order 171. The North Carolina Cooperative Fish and Wildlife Research Unit is jointly supported by North Carolina State University, North Carolina Wildlife Resources Commission, U.S. Geological Survey, U.S. Fish and Wildlife Service, and Wildlife Management Institute. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
