Aggregation of two imperfectly detected imperilled freshwater fishes: understanding community structure and co-occurrence for multispecies conservation

Characterizing critical habitat for endangered species is the first step for effective conservation and relies on knowledge of the abiotic conditions that promote species occurrence. Knowledge of biotic interactions, a component of species habitat, can also provide vital information for the conservation of endangered species by identifying conditions and/or locations suitable for multispecies conservation. We investigated site-level species co-occurrence patterns for 2 federally protected fish species in Canada, eastern sand darter Ammocrypta pellucida and silver shiner Notropis photogenis, to determine if management measures applied to one species could confer benefits to the other species. Using repeat surveys of fish communities in the Grand River, Ontario, we developed 2-species occupancy models to control for imperfect detection while evaluating species co-occurrence patterns and related those patterns to abiotic habitat covariates. We also used a principal component analysis (PCA) to characterize observed fish community structure and fit abiotic variables to PCA axes to explain observed patterns. Results indicated aggregation between the 2 federally protected species in the Grand River; the probability of occupancy for N. photogenis was higher at sites where A. pellucida was present. Results from the PCA and cooccurrence models confirmed the importance of depth, water clarity, and water velocity for both species. The result of site-level aggregation after correcting for imperfect detection suggests that management measures applied to one species could confer benefits to the other; however, the mechanisms underlying co-occurrence remain poorly understood.


INTRODUCTION
Characterizing species−habitat associations is fundamentally important when developing recovery strategies for endangered fishes because knowledge of abiotic conditions that support fish populations allows managers to identify and protect areas of critical habitat (Rice 2005, Barletta et al. 2010. In Canada, critical habitat is defined under the Species at Risk Act (SARA; Government of Canada 2002) as 'habitat that is necessary for the survival or recovery of a listed wildlife species and that is identified as the species' critical habitat in the recovery strategy or in an action plan for the species', where habitat for aquatic species is defined as 'spawning grounds and nursery, rearing, food supply, migration and any other areas on which aquatic species depend directly or indirectly in order to carry out their life processes, or areas where aquatic species formerly occurred and have the potential to be reintroduced' (Government of Canada 2002, p. 4). Therefore, fish habitat consists of all resources and conditions -abiotic and biotic -that support a given organism (Hall et al. 1997). While abiotic associations have been the subject of most habitat-related research for endangered fishes and represent an important science need to inform recovery strategies, understanding biotic interactions can also provide vital information for conservation managers (Angermeier & Winston 1999, Wenger et al. 2011.
Biotic interactions are characterized by a direct relationship between 2 organisms at a given space and time. The relationship can occur between 2 individuals of the same species or between different species, and can be mutualistic (+,+), predatory (−,+), competitive (−,−), commensalistic (+, 0), or amensalistic (−, 0). Despite the importance of identifying biotic interactions for aquatic restoration and species conservation (Halpern et al. 2007, Preston et al. 2008, biotic interactions have been relatively understudied compared to abiotic associations because they are often more difficult to measure in the wild, especially for endangered fishes. For example, it is difficult to observe nest associations or feeding behavior of small-bodied fishes that are found in turbid environments with non-lethal methods. These challenges preclude direct measurement of the mechanisms underlying biotic interactions for many aquatic organisms. Consequently, researchers often use species co-occurrence patterns as a proxy for biotic interactions by identifying whether species co-occur more than, less than, or as frequently as expected based on a hypothesis of independence (Peres-Neto 2004, Araújo et al. 2011. Although it is difficult to derive mechanisms of interactions from species co-occurrence patterns, analysis of co-occurrence provides insight into the direction of relationships between species, thereby informing future experimental work or providing evidence of shared or disparate abiotic preferences.
In addition to providing a better understanding of the local ecology of endangered fishes, testing for site-level patterns of positive or negative species cooccurrence, hereafter referred to as aggregation and segregation, respectively, has implications for the management of both species. For example, evidence of aggregation between 2 species would indicate that factors influencing the habitat patches of one species -whether restorative (i.e. habitat restoration) or destructive (i.e. habitat threats or other stressors)would preferentially influence the other species (Halpern et al. 2007). Alternatively, evidence of segregation or independence would indicate that the 2 species require different habitat characteristics or that a negative interaction exists whereby the mechanism is unknown. For example, evidence of segregation could suggest a predator−prey relationship between 2 species. Nevertheless, understanding patterns of aggregation, segregation, and/or independence among co-occurring species, particularly en dan gered species, has direct implications for multi species conservation planning (Halpern et al. 2007, Preston et al. 2008, DeCesare et al. 2010. Southwestern Ontario has the highest diversity of freshwater fishes in Canada (Chu et al. 2008), but also contains the greatest number of SARA-listed freshwater fish species (n = 20 species; Staton & Mandrak 2006, Hutchings & Festa-Bianchet 2009). Most of these SARA-listed fishes are small bodied (< 20 cm total length), are at the northern limit of the species' range (Glass et al. 2017), and co-occur or cooccurred with other imperilled species in the region. For ex ample, eastern sand darter Ammocrypta pellucida is a benthic percid listed under SARA as Threatened with populations distributed discontinuously at their northernmost extent within Ontario and Québec (COSEWIC 2009, Fisheries and Oceans Canada [DFO] 2011, 2012. A major threat to A. pellucida in Ontario is increased sedimentation leading to the loss of habitats containing sand and fine gravel substrates (DFO 2012), the favored habitat of A. pellucida (Daniels 1993, Drake et al. 2008, Dextrase et al. 2014. One of the few drainages supporting A. pellucida populations in Canada is the Grand River, Ontario, which also supports a population of silver shiner Notropis photo genis, a species listed as Threatened under SARA (DFO 2013, 2019 1 ). N. photo genis is a surface and midwater column feeding species known to be associated with moderately deep riverine habitats containing relatively high water clarity and moderate to fast water velocity (McKee & Parker 1982, Glass et al. 2016. Although both species are extant in the Grand River drainage, it is unknown if A. pellucida and N. photogenis cooccur within small habitat patches (i.e. sites) more frequently than expected based on chance, and/or if these SARA-listed species display directional associations. For example, the burrowing and bottomdwelling behavior of A. pellucida could create a bottom disturbance, suspending particles and insects into the water column (e.g. Teresa & Carvalho 2008), 1 www.dfo-mpo.gc.ca/species-especes/profiles-profils/silvershiner-mene-miroir-eng.html benefiting the midwater column feeding N. photogenis. Alternatively, A. pellucida may compete for food resources with N. photogenis, particularly the juvenile stage, which does not display the same affinity for surface feeding.
The goals of this study were to evaluate co-occurrence patterns of 2 SARA-listed freshwater fish species in Canada (A. pellucida and N. photogenis) and relate patterns of co-occurrence to environmental covariates to assess the potential for multispecies conservation planning, including habitat restoration and the susceptibility to threats. We addressed 2 primary objectives related to abiotic associations and potential biotic interactions for A. pellucida and N. photogenis while accounting for imperfect detection to determine (1) the direction and magnitude of the association between A. pellucida and N. photogenis occupancy in the Grand River watershed and (2) the relative importance of sand and fine gravel substrates, depth, and water clarity for A. pellucida and N. photo genis occupancy. Overall, our study provides novel insights into the joint abiotic and biotic associations for 2 federally listed fishes in Canada.

Data collection
Fish and habitat surveys were performed at 151 sites along the main stem of the Grand River, On ta rio, Canada . Point measurements of depth (cm), substrate size (mm), and water clarity (cm) were taken at each site, along with water velocity at 60% depth (m s −1 ; Table 1). Substrate composition was assessed using a modified pebble count (Wolman 1954) and included the mean substrate size at the site as well as the proportion of pebble counts in the substrate categories: clay/silt (< 0.063 mm), sand (0.063−2.0 mm), fine and medium sand (0.063−0.05 mm), fine gravel (2.1− 8.0 mm), sand and fine gravel combined (0.063− 8.0 mm), and coarse material (> 8.0 mm; Dextrase 2013).
Sampling of fishes was performed by conducting 3 consecutive 10 m downstream hauls at each site with a 9.2 m bag seine (hence 92 m 2 sites). After each haul, captured fishes were placed into separate holding aquaria, identified, and counted before release. For the purposes of occupancy modelling, we converted count data into detection and non-detection data ( Fig. S1 in the Supplement at www. int-res. com/ articles/ suppl/ n040 p123 _supp.pdf provides a plot of A. pellucida and Notropis photogenis naïve relative abundance at sites where either species was captured; n = 49 sites). The 3 repeat seine hauls were used as repeated surveys at each site for modelling occupancy and detection, whereas the 3 hauls were combined into a single observation of abundance for multivariate analyses.

Multivariate analysis
We performed a principal component analysis (PCA) on the observed fish community abundance data and fit environmental vectors onto the PCA ordination plot to visualize community and speciesspecific habitat associations. The fish community data were Hellinger transformed prior to performing the PCA, and nontrivial axes were retained through permutations (Legendre & Gallagher 2001, Lamothe et al. 2018. We permuted each column of the transformed community data matrix and conducted a PCA 9999 times; a component was retained if the proportion of variance explained in the observed data was greater than 95% of the permuted PCAs for that component (Peres-Neto et al. 2003, 2005. Pearson correlations were calculated to ensure a lack of correlation between abiotic variables. Mean substrate size, depth, and water clarity measures were scaled and centered (i.e. z-scores) prior to fitting them to the significant PCA axes. Squared correlation coefficients (r 2 ) were calculated for each environmental variable and the significant PCA axes, and their significance was assessed using 9999 permutations.

Occupancy models
We developed maximum-likelihood based, singleseason 2-species occupancy models for A. pellucida and N. photogenis with notation de scribed by Richmond et al. (2010). With a 2-species occupancy modelling approach, sampling sites could be characterized as one of 4 mutually exclusive states: (1) A. pellucida and N. photogenis are both present, (2) A. pellucida is present and N. photogenis is ab sent, (3) A. pellucida is absent and N. photogenis is present, and (4) A. pellucida and N. photogenis are both absent.
We compared 9 possible detection models that in -cluded 1 or 2 detection covariates and used the best detection model to examine patterns of occupancy. Variables considered in the detection models included a species-specific detection parameter, a variable characterizing species removal after each seine haul, and scaled and centered site-specific mean substrate size (Table 1). The removal variable consisted of a row vector of survey-specific species removal history, whereby a 0 was assigned for a haul when the species was not detected in a previous haul, and a 1 was assigned when the species was captured in a previous seine haul ( We developed 17 occupancy models to consider combinations of site-specific habitat covariates (scaled and centered water clarity, depth, water velo city, and the proportion of sand and fine gravel), speciesspecific effects, and interactions between spe cies whereby the presence of A. pellucida has an effect on the presence of N. photogenis. Three models failed to converge and were removed from the analysis. We compared candidate models using Akaike's information criterion (AIC), where the best models had the lowest AIC scores and a ΔAIC of less than 2 units.
A species interaction factor (ϕ) can be derived from occupancy probabilities to identify positive, negative, or independent associations between species using: where ψ A = probability of occupancy for A. pellucida, ψ BA = probability of occupancy for N. photogenis with A. pellucida present, and ψ Ba = probability of occupancy for N. photogenis with A. pellucida absent. When species occur independently, ϕ = 1. When ϕ < 1, N. photogenis is less likely to co-occur with A. pellucida than expected under a hypothesis of independence (i.e. segregation). In contrast, when ϕ > 1, N. photogenis is more likely to co-occur with A. pellucida than expected under a hypothesis of independence (i.e. aggregation; Richmond et al. 2010). Note that significant aggregation or segregation patterns do not imply direct species interactions but instead suggest that when imperfect detection is accounted for, A. pellucida and N. photogenis cooccur more or less frequently than expected under a hypothesis of independence.

Fish community description
In total, 49 species representing 34 genera were observed in the Grand River (Table S1 in the Supplement). Spotfin shiner Cyprinella spiloptera was the most prevalent species across sites (90.1% of sites), whereas many species were only observed at a single site (brindled madtom Noturus miurus, brook stickleback Culaea inconstans, creek chub Semotilus atromaculatus, mooneye Hiodon tergisus, northern pike Esox lucius, northern sunfish Lepomis peltastes, river chub Nocomis micropogon, and white perch Morone americana). Eastern sand darter Ammocrypta pellucida was observed at 34 sites (22.5% of sampled sites), silver shiner Notropis photogenis was observed at 29 sites (19.2% of sampled sites), and the 2 species were observed co-occurring at 12 sites (7.9% of sampled sites).
Four axes were extracted from the community PCA explaining 25.68% of the total variation in observed community composition (Fig. 2). Based on random permutations, water clarity, water velocity, and site depth were significant (p < 0.05) predictors of community composition (Fig. 2), with r 2 values ranging from 0.17 (depth) to 0.47 (water velocity). Notably, the proximity of A. pellucida and N. photogenis in ordination space indicates that patterns of relative abundance in the Grand River are similar and that both species relate to areas with higher water clarity, depth, and water velocity (Fig. 2).

Occupancy models
The model-averaged probability of detection across sites was 0.605 ± 0.029 SE for A. pellucida and 0.480 ± 0.030 SE for N. photogenis. The best detection model included mean substrate size, where the effect of mean substrate size differed between species (Table 2, Fig. 3). As substrate size increased, detection probability for both species decreased (Fig. 3). The model-averaged probability of occupancy in the Grand River was 0.327 ± 0.021 SE for A. pellucida and 0.455 ± 0.008 SE for N. photogenis. In  Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ n040 p123 _ supp. pdf. WC: water clarity (cm); WV: water velocity (m s −1 ). Circled are eastern sand darter Ammocrypta pellucida (ESD) and silver shiner Notropis photogenis (SS) the absence of A. pellucida, modelaveraged N. photo genis occupancy was 0.330 ± 0.008 SE, whereas in the presence of A. pellucida, model-averaged N. photogenis occupancy was 0.715 ± 0.010 SE. The best occupancy model inclu ded the proportion of sand and fine gravel, where the effect of the proportion of sand and fine gravel on occupancy differed between species, as well as an interaction term (INT; Table 2). The interaction term indicates that the presence of A. pellucida had an effect on the probability of occupancy for N. photogenis, but the effect was constant with respect to the proportion of sand and fine gravel at a site. The species interaction factor (ϕ) for the top occupancy model indicated a positive association between A. pellucida and N. photogenis across sites (ϕ = 2.055 ± 0.104 SE; Table 2); similarly, the remaining models with interaction terms all estimated a positive association (ϕ > 1; Table 2).
We used the best models that included either the proportion of sand and fine gravel, depth, water clarity, or water velocity to investigate patterns of occupancy across habitat covariates (Fig. 4). The proba-bility of occupancy for A. pellucida increased with higher proportions of sand and fine gravel, depth, water clarity, and water velocity (Fig. 4). Similarly, the probability of occupancy for N. photogenis increased with depth, water clarity, and water velocity but showed a unimodal response to the proportion of sand and fine gravel (Fig. 4). Across all proportions of sand and fine gravel, depths, and water clarity values, the probability of N. photogenis occupancy was always higher when A. pellucida was present than when A. pellucida was absent (Fig. 4) -consistent with expectations from ϕ calculations (Table 2) and the results of multivariate analysis (Fig. 2) Table 2. Summary of occupancy models selected using Akaike's information criterion (AIC). LL: log likelihood; npar: number of parameters in the model; ψ: occupancy probability; R: removal; SP: species-level effect; Substrate: mean substrate size (mm); (.) intercept model; INT: occurrence-level interaction between species; sfg: proportion of sand and fine gravel; p: detection probability; WV: water velocity (m s -1 ); WC: water clarity (cm); Depth: mean depth (m)

DISCUSSION
Aggregation patterns between 2 federally listed species, whereby the 2 species show similar patterns of occupancy across abiotic covariates, have important conservation implications given the challenges associated with species protection. Based on the positive association between Ammocrypta pellucida and Notropis photogenis in the Grand River and the species' similar responses to abiotic features (i.e. water clarity, depth, and water velocity), management to improve water clarity would likely return twice the expected benefits from a single restoration measure. A long-standing debate exists about whether imperilled species should be protected through speciesspecific management or management at the ecosystem or landscape scale (e.g. Simberloff 1998, Franklin 1993, Poos et al. 2008). An intermediate approach that blends species-and ecosystem-based approaches is the protection of an umbrella species, or a species for which species-specific pro tection would confer protection to a large number of co-occurring species (Launer & Murphy 1994). Based on the ordination of naïve community observations, which showed consistent patterns with our oc cupancy models, N. photogenis shared habitat as sociations with A. pellucida and sev eral minnow species, including mi mic shiner N. volucellus, spot fin shiner Cypri nella spiloptera, and rosyface shiner N. rubellus. Perhaps N. photogenis represents an umbrella species, whereby management to im prove N. photogenis ha bitat (i.e. water clarity and flow regime) would lead to benefits for several other species. Alternatively, the presence of these common species (i.e. C. spiloptera, N. rubellus, and N. volucellus) may provide an indication of potentially suitable habitat to support N. photogenis or habitat that could be selected for restoration.
Although species co-occurrence was used here as a proxy for biotic interactions, the mechanism behind the positive association is still unknown. It is possible that the scale of sampling was larger than actual microhabitats selected by A. pellucida and N. photogenis, and therefore the positive association was an artefact of sampling. Based on an allometric relationship (Minns 1995), the home range for an average-sized adult A. pellucida (Ontario average: 60 mm; Holm et al. 2009) is approximately 47 m 2 , whereas the home range of an average-sized adult N. photogenis (100 mm; Holm et al. 2009) is approximately 109 m 2 . The sampling approach employed in this study focused on erosional and depositional bars of river reaches, the anticipated habitat for A. pellucida, and was performed using a bag seine over 92 m 2 sites (Dextrase 2013). Seining is an approach that covers the entirety of the water column and therefore does not allow differentiation between how species are vertically distributed. However, A. pellucida is a benthic species that feeds on small bottom-dwelling invertebrates (Fisheries and Oceans Canada 2012), whereas N. photogenis is an opportunistic feeder that capitalizes on aquatic invertebrates throughout the water column, including jumping out of the water to capture airborne insects (COSEWIC 2011). As such, positive co- occurrence patterns could indicate correlated areas of foraging, whereby the smaller A. pellucida is residing in the benthic zone and the larger N. photogenis is residing in mid-column areas.
A finding of segregation between federally listed species, contrary to results presented here, would have more challenging conservation implications and would likely require protection of unique habitats. For example, we recently demonstrated negative co-occurrence patterns between A. pellucida and johnny darter Etheostoma nigrum and rainbow darter E. caeruleum in Ontario, 2 species widely present across the southern Ontario landscape . Based on our present multivariate analysis, E. nigrum and E. caeru leum were associated with areas of lower flows and shallower depths (Fig. 2). If the general protection of darters was a key management objective, then it is possible that the scale of management would change by focusing on multiple habitat types, with different flow characteristics, rather than focusing on water clarity alone.
Our study also demonstrates the importance of repeat sampling to account for imperfect detection when characterizing patterns of species distributions, co-occurrence, and community assembly. Imperfect detection commonly impacts monitoring efforts for imperilled species (MacKenzie et al. 2002), as their rarity often requires unique sampling approaches. In this study, N. photogenis was detected infrequently in the Grand River (19.2% of sites) but had an average detection probability of approximately 0.5. As a result, the estimated mean occupancy probability was more than twice that from the naïve sampling data. Nevertheless, the results of our community analysis performed on the observed data showed similar habitat association results to those from our species co-occurrence models that incorporated the probability of detection. For example, N. photogenis was found near the tip of the water velocity vector within the community ordination (Fig. 2), indicating an association with habitats having higher water velocity, which is consistent with the literature (Glass et al. 2016) and the occupancy model results (Fig. 4D).
Understanding and quantifying patterns of species distributions, co-occurrence, and community assembly is fundamental for protecting endangered species. While traditional species distribution models based on abiotic relationships provide important insight into habitat components critical for sustaining populations (Anderson & Martínez-Meyer 2004, Williams et al. 2009), knowledge of imperilled species co-occurrence and community assembly patterns can help identify sites or biotic assemblages where the species is, or is not, likely to occur (Bailey et al. 2009, Cove et al. 2018. This is particularly important in diverse systems such as the Grand River, one of the most speciose rivers in Canada, and also helps to manage against impending human-mediated changes (the Ontario population is projected to grow to ~18 million people by 2041; Ontario Ministry of Finance 2018). However, the Grand River is also subject to a suite of environmental stressors threatening species persistence (e.g. Balasingham et al. 2018, Dean et al. 2018, Raab et al. 2018, Hanief & Laursen 2019. The results of our study indicate that improving water clarity in this watershed will benefit fish species that are most in need of protection, supporting recommendations in federal recovery documents (DFO 2012) and the longstanding goals of many aquatic restoration programs (e.g. Allison & Meselhe 2010). (2008)