Relationships between crayfish population genetic diversity, species richness, and abundance within impounded and unimpounded streams in Alabama, USA

Parallel processes may structure biodiversity

Community-level taxonomic diversity (e.g., species richness) and population-level genetic diversity are both key components of biological diversity (Allendorf & Luikart, 2007). According to the theory of island biogeography, the balance between colonization and local extinction determines species richness at a given site (MacArthur & Wilson, 1967), whereas the counteracting forces of gene flow and drift determine the standing levels of population genetic variation (Wright, 1940; Kimura, 1983; Nei, 1987). Vellend (2005) proposed that if parallel processes operate at both the community- and population-level, this should result in a positive species-genetic diversity correlation (SGDC). In addition to SGDCs, it has also been proposed that species richness and/or population-level genetic diversity may be positively correlated with the cumulative multispecies abundance of individuals within a community (Lamy et al., 2017). There are several reasons why this can occur. First, under the more individuals hypothesis (MIH; Storch, Bohdalková & Okie, 2018), community-level species richness at a local site may be high if populations of these species are large and stable in size, such that local extinction driven by negative feedback between intrinsic stochastic ecological and genetic processes (e.g., Allee effects, genetic drift, and inbreeding) is negligible. Likewise, direct and indirect species interactions can be beneficial (e.g., mutualism or facilitation), and species-rich communities may also be characterized by greater redundancy among key role players that are critical to ecosystem functioning, making these communities more resilient to environmental perturbations (Waide et al., 1999; Finke & Snyder, 2008; Lamy et al., 2017). Second, according to the abundance-genetic diversity correlation (AGDC) hypothesis (Johansson et al., 2005; Overcast, Emerson & Hickerson, 2019), if effective population size (N_e) and census population size (N_c) scale with one another (which is supported by a meta-analyses (Frankham, 1996; McCusker & Bentzen, 2010)), then environmental conditions that promote high cumulative multispecies abundance of individuals should also translate into relatively weak effects of drift, preventing the loss of allelic diversity (Storch, Bohdalková & Okie, 2018). Large N_e also provides the basis for more efficient natural selection, and balancing selection may play an important role in retaining genetic variation within populations (Chesson, 2000).

Numerous studies have assessed evidence for positive SGDCs, MIH, and AGDCs in diverse groups of organisms (Hurlbert, 2004; Dudgeon & Ovenden, 2015; Xie et al., 2021; Bucholz et al., 2023), but negative and no correlation have also been detected when factors such as habitat fragmentation and disturbance affect species richness, population genetic diversity, and multispecies abundance differently (Scribner et al., 2001; Johansson et al., 2005; Wei & Jiang, 2012; Šímová, Li & Storch, 2013; Seymour et al., 2016; Watanabe & Monaghan, 2017; Reisch & Schmid, 2019). Recent meta-analyses (Xie et al., 2021) and reviews of published studies that investigated evidence for SGDCs (Lamy et al., 2017) and MIH (Storch, Bohdalková & Okie, 2018) found that most studies (80%-SGDC and 72%-MIH) reported positive correlations. Nonetheless, no relationship between species richness, genetic diversity, and abundance was detected in studies conducted in highly disturbed areas or along environmental gradients (Carnicer & Díaz-Delgado, 2008; Wei & Jiang, 2012; Šímová, Li & Storch, 2013; Fan et al., 2021; Petersen et al., 2022). Additionally, negative SGDCs have also been documented when there is an opposite influence of environmental drivers (e.g., altitude, temperature) and/or competition on species diversity, genetic diversity, and abundance (Scribner et al., 2001; Currie et al., 2004; Seymour et al., 2016; Lamy et al., 2017; Marchesini et al., 2018; Ishii et al., 2022). Taken together, these inconsistent results indicate that community- and population-level processes may not operate in parallel ways in some ecological contexts and groups of organisms (Lamy et al., 2013; Storch, Bohdalková & Okie, 2018).

SGDC, MIH, and AGDC are not mutually exclusive, but there are several reasons why only a subset (or perhaps just one) of these hypotheses may receive support in a given study. First, extrinsic factors may impact species differently (Kahilainen, Puurtinen & Kotiaho, 2014; Bucholz et al., 2023). This may be due to divergent ecological functions, life histories, ecological optima, phenotypic plasticity and/or capacity to respond to dynamic environmental conditions. Second, contrasting carrying capacities of habitats that support local communities may impose constraints on diversity and abundance of some species but not others (Loreau, 2000). Additionally, rare and specialized species’ genetic diversity often have little correlation to carrying capacities of habitats, while in most communities species diversity and abundance are strongly correlated with carrying capacities (Vellend, 2005).

While the assessment of SGDC, MIH, and AGDC requires a field- and labor-intensive multi-level approach aimed at investigating different levels of biodiversity, as well as their evolutionary ecological drivers within a community, such studies can provide important insights into the “scaling” of processes shaping biodiversity, which in turn have practical implications for conservation (Kahilainen, Puurtinen & Kotiaho, 2014). For example, if all three hypotheses are supported in a given study, then natural resource managers could collect any one of the three types of diversity data (i.e., species richness, genetic diversity, or cumulative multispecies abundance) and make reasonable predictions about the other two (Kahilainen, Puurtinen & Kotiaho, 2014; Overcast, Emerson & Hickerson, 2019). Likewise, if two of the three hypotheses were supported, this information could be used to identify which one of the three data types should be prioritized in biodiversity assessments (e.g., if both SGDC and MIH are true, then species richness data alone could be used to predict genetic diversity and abundance; Bucholz et al., 2023). Although the robustness of such extrapolations should be empirically verified at several sampling sites, in such situations there is potential for this predictive framework to improve the cost-effectiveness of biodiversity monitoring and conservation strategies. Furthermore, making assumptions without assessing biodiversity correlations can lead to sub-optimal or detrimental conservation strategies.

Streams as study systems for assessing parallel processes

When assessing the correlations between biodiversity metrics, having a system that is explicitly linked through dispersal takes into account the spatial effects of migration and dispersal which are key processes involved in diversity dynamics (Vellend & Geber, 2005; Altermatt, 2013; Seymour et al., 2016). As such, a system, such as streams, with a dendritic network, provide a suitable and tractable study system for empirically testing whether parallel processes, operating at hierarchically nested levels, structure biodiversity in similar ways. This is because stream boundaries are clearly demarcated and environmental characteristics (e.g., habitat area and substrate size), and distributions of stream-dependent biota, are usually predictably structured along an upstream-downstream gradient (Vannote et al., 1980; Rimalova, Douda & Stambergova, 2014; Barnett et al., 2022). Furthermore, because most stream-dependent biota are regionally constrained by the network spatial arrangement (Grant, Lowe & Fagan, 2007) and locally restricted to the stream channel, it is possible to temporarily isolate sections so they can be exhaustively sampled (e.g., via block net multi-pass electrofishing). This provides an otherwise rare opportunity for quantitative assessments of stream community species richness and abundance (Hornbach & Deneka, 1996; Ode, Rehn & May, 2005; Hanks, Kanno & Rash, 2018; Barnett et al., 2020).

Potential impacts of stream fragmentation

Notwithstanding the aforementioned advantages of streams as study systems, loss of connectivity due to human-mediated habitat fragmentation (e.g., dams and impoundments) is common, and this can impact biodiversity of resident biota by altering stream habitat, community composition, and population genetic structure (Barnett et al., 2020, 2022, 2023). Dams and impoundments are among the most prevalent and extreme alteration on fluvial systems (The Heinz Center, 2002; Liermann et al., 2012; Grill et al., 2015), with river fragmentation and flow regulation being one of the largest biological effects of dams (Stanford & Ward, 2001; Grill et al., 2015; Barnett et al., 2021). Unlike unimpounded streams where organisms experience the natural flow variability and can freely move throughout the stream system, organisms in impounded streams are often isolated to one stream section and natural flow variability is greatly reduced causing changes to habitat composition and accessibility. These changes can impact species richness, abundance, and dispersal throughout the stream system. Furthermore, habitat fragmentation is expected to reduce within-patch species richness and abundance due to decreases in habitat complexity which can reduce suitable habitat for habitat-specialists (Barnett et al., 2022), as well as reduce intraspecific genetic diversity due to restricted dispersal and gene flow leading to isolation and drift (Vellend & Geber, 2005; Hartfield, 2010; Barnett et al., 2020). However, a decoupling of the effects on these two diversity metrics may occur when habitat fragmentation impacts the dispersal ability of some species differently than others (Lamy et al., 2017). For example, crayfish that prefer smaller streams and naturally disperse upstream, up steep slopes and against fast water velocities, may be capable of bidirectional gene flow within fragmented streams, while those preferring larger sized streams may have unidirectional downstream or no gene flow between fragmented sections (Hartfield, 2010; Barnett et al., 2020). Thus, whether changes to species richness, abundance, and population genetic diversity in fragmented systems mimic those in connected systems is an important question for conservation biologists.

Crayfish as focal group of organisms

Nearly 70% of the world’s freshwater crayfish species are found in the United States (US) (Crandall & Buhay, 2008; Richman et al., 2015), with the southeastern US being the major center of diversity (Hobbs, 1989; Richman et al., 2015). These organisms play fundamental roles in stream ecosystem trophic processes (e.g., processing detritus, altering the composition of macrophyte and substrate, transferring energy to higher level organisms), and they are often considered ecosystem engineers due to their modification of the physical habitat (e.g., creating habitat for other organisms through burrow creation, bedform alterations in streams, etc.) (Momot, 1995; Usio, 2000; Statzner, Peltret & Tomanova, 2003; Usio & Townsend, 2004; Reynolds, Souty-Grosset & Richardson, 2013; Krupa, Hopper & Nguyen, 2021). Alarmingly, 48% of North American crayfish species are threatened (Taylor et al., 2007), with extinction rates likely to increase by more than an order of magnitude over the next several decades (Ricciardi & Rasmussen, 1999; Cowie, Bouchet & Fontaine, 2022; Finn, Grattarola & Pincheira-Donoso, 2023). Hence, there is an immediate need for effective crayfish conservation strategies (e.g., Taylor et al., 2019). However, conservation planning at the community-level has been emphasized much more strongly than population-level genetic diversity, as has preservation of specific “units” or phenotypes over the evolutionary processes that generate this diversity (e.g., large self-sustaining populations living in heterogeneous landscapes; Moritz, 2002). Thus, to effectively conserve crayfish diversity, understanding the extent to which similar processes structure biodiversity across different levels of biological organization is key.

In the present article, we assessed evidence for parallel processes, operating at hierarchically nested levels, within connected and fragmented (i.e., unimpounded versus impounded) streams in the southern Appalachian Mountains, Alabama. SGDC and AGDC were each tested using Faxonius validus and F. erichsonianus. These were chosen as our focal species because they share many ecological traits (e.g., life span, mating season, burrowing capabilities, preferred habitat) but differ in their preferred stream size and geographic range (Bouchard, 1972; Williams & Bivens, 2001; Hopper, Huryn & Schuster, 2012; Barnett et al., 2020). Faxonius validus occurs in small intermittent to medium-sized perennial streams in northern Alabama and southern Tennessee (Cooper & Hobbs, 1980; Hobbs, 1989), while F. erichsonianus occurs in medium to large streams in six southeastern states, from Mississippi to Virginia (Hobbs, 1981). For these species, population genetic diversity was measured using several alternative metrics based on mitochondrial DNA sequences (mitochondrial DNA cytochrome oxidase subunit I (mtCOI)), as well as complementary nuclear genetic maker (inter-simple sequence repeat (ISSR)) data. Testing of all three hypotheses (i.e., including MIH) incorporated data from community-level surveys of crayfish species richness and cumulative multispecies abundance in five streams spanning two drainages. To clarify the extent to which habitat fragmentation and environmental characteristics may contribute to such correlations, we also investigated whether stream habitat characteristics (size, connectivity, and habitat complexity) were correlated with species richness, abundance, and population genetic diversity (Table 1). Taken together, outcomes from these analyses should inform strategies for crayfish conservation in a geographic region of high endemism.

Study area description

This study focused on lotic sections of two unimpounded and three impounded streams that are distributed across two drainages with diverse aquatic communities and numerous imperiled species (Allen, 2001; McGregor & Garner, 2003; Phillips & Johnston, 2004; Barnett et al., 2022). In the Bear Creek drainage (Tennessee River Basin), we surveyed one unimpounded (Rock Creek) and two impounded (Little Bear and Cedar creeks) streams (mean stream length: 61.5 km). In the Cahaba River drainage (Mobile River Basin), we surveyed one unimpounded (Shades Creek) and one impounded (Little Cahaba River) stream (mean stream length: 41.6 km) (Fig. 1).

Each impounded stream had one earthen storage dam (17–29 m high), creating impoundments that were 425 to 1,700 ha. The two Bear Creek drainage impoundments were installed for flood control, and water was released from outlets more than 19 m below full-pool levels from November until February and during heavy rain events. Conversely, the Cahaba River drainage impoundment stored water for municipal use. Water was released from outlets more than 10 m below full-pool levels when river water levels were insufficient to meet water municipal demands. Dams in both drainages pass normal inflows via spillways throughout the year.

Site selection

In each of the three impounded streams, we selected sampling sites at set intervals (Data S1) up- and downstream of impoundments, and we mimicked the same sampling design in each of the two unimpounded streams. This approach led to selection of two to five local sites up- and downstream of impoundments, and up- and downstream of the midpoint in the unimpounded stream (hereafter, up- and downstream sections) (collection data same as in Barnett et al., 2020). Overall, our study included 32 sites: 24 in the Bear Creek drainage (10 in Little Bear Creek, six in Rock Creek, and eight in Cedar Creek), and eight in the Cahaba River drainage (four in Shades Creek, and four in Little Cahaba River) (Fig. 1).

Crayfish monitoring and sampling

For community-level assessment of diversity, sampling was conducted in spring (May–July) and fall (September–December) of two sampling years, during which all crayfishes were collected. During the first sampling year (spring and fall 2015), we collected crayfishes from all sites in the Bear Creek drainage. During the second sampling year (fall 2016 and spring 2017), we collected crayfishes from all streams in the Cahaba River drainage. At each site, we sampled one linear reach, 30 times the wetted stream width or a minimum or maximum length of 200 to 500 m, respectively (Barnett et al., 2022). Stream reach lengths remained constant across seasons unless the dry season shortened a reach. Each reach was divided into two subreaches of equal length. Each subreach was simultaneously sampled by electrofishing (3–8 s/m; mean 5 s/m ± 1.2 SD) upstream subreaches and kick seining (20 plots/100 m, 2 m long × 1.5 m wide) downstream subreaches (Barnett et al., 2021). Because these methods are ineffective in pools and deeper waters, crayfish were collected only from riffles and runs with depths less than 1 m (≤15% of each reach). Electrofishing duration and total number of kick seines were calculated based on subreach areas. We recorded the amount of area (m²) sampled by each method once the target sampling effort (electrofishing: 250–2,000 sec/subreach; kick seining: 20–50 kicks/subreach) was reached (Barnett et al., 2021). We used expert knowledge to identify crayfish species in the field, counting the number of individuals for each species collected at each site during each sampling round. We preserved voucher specimens (housed at Mississippi Museum of Natural Science) for each species in ≥70% ethanol and confirmed species identifications in the lab (Hobbs, 1981, 1989). All collections were approved by the State Alabama under Alabama Conservation License numbers 2016064289868680 and 2017092711268680.

For population-level assessment of genetic diversity, the two focal species, F. erichsonianus and F. validus, were both collected from the Bear Creek drainage, but only F. erichsonianus was collected from the Cahaba River drainage. On average, 20 individuals per species (SD = 6.6) were collected per stream section, with a total of 143 F. validus and 173 F. erichsonianus collected) (Table 2). Whole specimens were preserved in 95% ethanol.

Characterization of population genetic diversity in two focal species

Genomic DNA was extracted from leg tissue of each F. erichsonianus and F. validus individual using a DNeasy blood and tissue kit (Qiagen, Valencia CA, USA), following manufacturer’s recommendations. A 618–640 base pair (bp) region of the mitochondrial DNA cytochrome oxidase subunit I (mtCOI) gene was amplified and sequenced as described in Barnett et al. (2020). Additionally, we performed genotyping using inter-simple sequence repeat (ISSR) markers (Ziętkiewicz, Rafalski & Labuda, 1994) for a minimum of three individuals per focal species per site (mean = 4, SD = 1.07). Whereas mtCOI is a maternally inherited haploid marker, ISSRs are generally presumed to be biparentally inherited nuclear autosomal markers, and as such, they have been used for addressing diverse questions in ecology and evolution (e.g., Wolfe, Xiang & Kephart, 1998; Abbot, 2001; Haig, Mace & Mullins, 2003; Dušinský et al., 2006; Sinn et al., 2022). Using both nuclear and mitochondrial data could provide information on two temporal scales, with mtCOI having a smaller N_e than ISSRs, potentially giving it the ability to detect more recent changes to populations (Moore, 1995). Together, these two types of molecular data should provide a broad overview of population genetic diversity (Garrick, Caccone & Sunnucks, 2010).

For initial assessment of polymorphism and scorability of ISSRs, ten primers (Data S2) were screened using a geographically representative of 10 F. validus and 10 F. erichsonianus individuals. Eight ISSR primers were designed by the author (R.C. Garrick), with the other two primers selected from the UBC Primer Set #9 (University of British Columbia, Canada, available at www.github.com/btsinn/ISSRseq). Polymerase chain reaction (PCR) amplifications were carried out in a final volume of 10 μL, containing the following: 1 μL genomic DNA, 2 μL 5 × buffer (Promega, Maddison WI, USA), 0.8 μL MgCl₂ (25 mM, Promega), 1.6 μL dNTPs (1.25 μM, Promega), 0.5 μL bovine serum albumin (10 mg/μL, New England Biolabs, Ipswich, MA, USA), 3 μL dH₂O, 0.1 μL Go-Taq DNA polymerase (5 U/μL, Promega), and 1 μL of primer (10 μM). Thermocycling conditions were: 95 °C for 2 min (one cycle), 95 °C for 30 s, 48 °C for 30 s, 72 °C for 1 min 30 s (35 cycles), and a final extension at 72 °C for 2 min (one cycle). All PCRs contained a negative control to assess evidence for contamination. Amplified products were electrophoresed on 2% agarose gels at 100 volts for 2 min and 45 volts for 16 h and 40 min in a 4 °C cold room, and then viewed under ultraviolet light, and photographed. Sizes of amplified bands were approximated via comparison to a 100-bp DNA ladder. Preliminary results identified four primers (two per species) that produced informative data (Data S2), and these were selected for population-level screening.

A given allele at each ISSR locus is scored as binary presence (1) vs. absence (0) data, and several loci are typically co-amplified with the same primer. Accordingly, an individual’s gel banding profile usually contains multiple bands of different size (i.e., alleles present at different loci) and represents a multi-locus genotype. To standardize scoring of band sizes across gels and to ensure repeatability of banding profiles, PCR products from each individual were run two to three times with strategic re-ordering of individuals across gels so as to provide key side-by-side comparisons and scoring of each profile performed by two people. Only those loci and individuals that yielded reproducible results were included in downstream analyses.

Habitat characterization

During spring and fall sampling 2015–2017, we measured stream channel characteristics. This included channel wetted width at four evenly spaced transects within our sampled reach. Using Wolman pebble counts procedures (Wolman, 1954; Barnett et al., 2020), we analyzed habitat complexity across the bankfull channel width. Ten zig-zag transects from one bank to the other were sampled at each site, with 10 points in each transect (100 sampling points/site). At each sampling point, we measured the intermediate axis of substrate. Between adjacent sampling points, we visually estimated the percentage of streambed covered by vegetation and counted number of pieces of large woody debris (LWD) (Bain & Stevenson, 1999).

Diversity metrics

Species richness was measured following Chao (1984), using the Chao-1 metric, which extrapolates the probability of undetected species within each site from the number of rare species detected (i.e., singletons). Chao-1 species richness was calculated using the “chao1” function of the fossil package in R software (version 4.2.1; R project for Statistical Computing, Vienna Austria) (R Core Team, 2022).

To quantify cumulative multispecies abundance, we counted the number of crayfish collected after electrofishing and kick seining and used total number of individuals and area sampled to calculate the number of crayfishes collected/100 m². The cumulative multispecies abundance was summed across two sampling rounds per site.

Within-population genetic diversity was calculated using several different metrics. For mtCOI sequence data, we used DnaSP v.5.10.01 (Librado & Rozas, 2009) to calculate haplotypic diversity (hd; Nei, 1987). Notably, hd treats mtCOI haplotypes as a multi-state unordered variable. To capture information on how different haplotypes within a population were from one another, we also calculated nucleotide diversity (π; Nei, 1987). One of the limitations of both hd and π is that they are calculated using only those haplotypes present within a given local population, and so these metrics can suffer from issues of small sample sizes. Accordingly, we also explored the utility of Faith’s (1992) Phylogenetic Diversity (PD) as means of jointly considering all of the available mtCOI sequence data for a given focal species when calculating population-specific diversity values. Briefly, PD is the sum of branch lengths on phylogenetic tree uniting all “taxa” (i.e., haplotypes present within a location/population), back to the root of the tree. For both F. erichsonianus and F. validus, rooted maximum-likelihood (ML) phylogenetic trees were estimated in MEGA v.10.2.6 (Kumar et al., 2018). Given that the use of an appropriate outgroup is important for PD, yet we are not aware of a published genus-level phylogeny for Faxonius, a phylogenetic analysis was conducted using mtCOI sequence data from 71 formally recognized species available in NCBI’s nucleotide database (see Data S4). This identified F. spinosus as the sister taxon of F. erichsonianus, whereas a clade including both F. cooperi and F. pagei was sister to F. validus. For each of the two focal species and their associated outgroup(s), the best-fit model of nucleotide evolution was determined via Akaike information criterion (AIC) model selection, and an ML tree was estimated using the following search settings: missing data = partial deletion (cut-off: 95%), maximum parsimony starting tree, nearest-neighbor-interchange branch swapping, and branch swap filter = moderate. Node support was assessed via 500 bootstrap replicates. The resulting tree was exported in Newick format containing tree topology plus estimated branch lengths, and PD was then calculated using the picante package (v1.8.2; Kembel et al., 2010) in R.

For the ISSR data, within-population genetic diversity was calculated as the proportion of polymorphic loci (PPL) (i.e., number of loci that were polymorphic among individuals at local site, divided by total number of loci screened for the focal species). Because there was no correlation between PPL and number of individuals within a population sample (Pearson correlation: F. validus: r = 0.27, P = 0.22; F. erichsonianus: r = 0.29, P = 0.16), subsequent rarefaction correction of PPL was not applied (Table 2).

Species richness, abundance, and genetic diversity correlations

For the focal Faxonius species at all study sites, collectively, we investigated SGDC assessing the relationships between species richness and each of the four genetic diversity metrics separately (i.e., hd, π, and PD for mtCOI sequences, and PPL for ISSR markers). To do this, we calculated the Pearson correlation coefficient and asymptotic confidence intervals based on Fisher’s Z transformation using the “cor.test” function of the stats package in R to determine significance (R Core Team, 2022). At each of the same 32 sampling sites for F. validus and F. erichsonianus (Fig. 1), we tested the MIH by assessing the correlation between species richness and cumulative multispecies abundance of all crayfishes, again using Pearson correlation coefficient and confidence intervals to determine significance, calculated in R. AGDC was assessed in the same way as SGDC, except that species richness was replaced by cumulative multispecies abundance.

Relationships between habitat characteristics, community- level diversity, and population-level diversity

To clarify the extent to which stream fragmentation and environmental characteristics (size, connectivity, and habitat complexity) may impact SGDC, MIH and/or AGDC, we examined whether species richness, abundance, and genetic diversity metrics showed any relationships with stream channel characteristics (see Habitat characterization, above) using linear models. If stream channel characteristics affect the two levels of diversity in a different way (e.g., positive relationship between stream size and genetic diversity vs. negative relationship between stream size and species richness), we would expect this to lead to no or negative correlation between different diversity metrics, thus resulting in no support for SGDC, MIH and/or AGDC.

To assess impacts of stream channel characteristics on diversity correlations, we calculated the median wetted width, percent vegetation, substrate size (i.e., D50), and LWD from spring and fall sampling for each site. All sites within unimpounded streams were characterized as connected, and all sites within impounded streams were characterized as fragmented. Because dams and associated impoundments can have far-reaching effects up- and downstream of impoundments (Falke & Gido, 2006; Johnson, Olden & Zanden, 2008), all sites sampled within impounded streams have the potential to be fragmented (e.g., isolation of upstream sites, alterations of seasonal flow patterns downstream and habitat modification both up- and downstream) (Yeager, 1993). We fit linear models with least squares estimates using the ‘lm’ function in with stats package in R. In these models, stream characteristics were treated as the independent variables (potential predictors), whereas the different metrics for species richness, cumulative multispecies abundance, or genetic diversity were treated as the dependent variable (response). We included 2-way interactions of stream characteristics and binary stream type classification (i.e., connected vs. fragmented). We used the MuMIn R package (Barton & Anderson, 2002) to analyze all possible models. Model selection was based on corrected Akaike information criterion (AICc) because sample sizes were small relative to the number of estimated parameters (Burnham & Anderson, 2004). We compared alternative models by weighting their level of data support (Hurvich & Tsai, 1989), with delta AICc values ≤2 representing the best-supported models. We calculated relative variable importance (RVI; number of models predictor variable appears in/number of total models) scores for each predictor variable, based on variables appearance in the AICc-best models. Predictors with RVI >0.5 were considered most important. If there were significant stream characteristics by stream type (connected or fragmented) interactions, pairwise comparisons between each stream type and stream characteristic were done. Tukey HSD P-value adjustment approach (Sokal & Rohlf, 1981) was used to correct for the effect of multiple comparison on the family-wise error rate.

Crayfish collections, and characterization of genetic diversity

Across all sites, we collected 12 crayfish species, with six and eight species collected in the Bear Creek and Cahaba River drainages, respectively (Table 3). Additionally, nine crayfish species were collected in both impounded and unimpounded streams. Cumulative multispecies abundance of individuals varied greatly between sites (N crayfish/100 m² = 0.001–0.282), with an average of 0.032 crayfish collected per 100 m² (Datas S5 and S6). The highest densities of crayfishes were collected in Rock Creek (0.723) and lowest in Shades Creek (0.016) (Table 3).

For F. validus, we successfully sequenced 143 individuals, obtaining a 618-bp mtCOI alignment, with 25 polymorphic sites and 28 unique haplotypes (Table 2; Data S3) (data from Barnett et al., 2020 assessed; GenBank accession numbers MN053979–MN054006, Data S3). For F. erichsonianus, we obtained a 640-bp mtCOI alignment, with 68 polymorphic sites and 42 haplotypes from 173 individuals (data from Barnett et al., 2020 assessed; GenBank accession numbers MN054007–MN054048, Data S3). We obtained ISSR data from 109 F. validus and 95 F. erichsonianus individuals. Faxonius validus and F. erichsonianus included in ISSR analyses were a subset of those used in mtCOI assessments. We assessed a minimum of three individuals per site (mean four individuals/site) and only used individuals that yielded reproducible results. Nelson & Anderson (2013) showed that genetic diversity estimates were similar when using five compared to 10 individuals per site, suggesting that our sample sizes are likely reasonable. ISSR primers yielded 24 and 34 polymorphic loci for F. validus and F. erichsonianus, respectively (Datas S5 and S6). While this is a relatively low number of polymorphic loci (Nelson & Anderson, 2013), studies have shown that using 20–30 dominant markers can yield acceptable results for population genetic assessments (Vandergast et al., 2009; Guasmi et al., 2012; Nelson & Anderson, 2013).

Habitat characterizations

Crayfish were collected in medium (median wetted width = 11.81 m) sized streams with mostly pebble substrate (median substrate size = 25 mm) and a relatively wide range of LWD (2–25 pieces of LWD; median = 10.0) and percent vegetation (6–32%; median = 15.6) (Table 4; Datas S5 and S6).

Species richness and genetic diversity estimates

Regarding community-level diversity, Chao-1 species richness ranged from one to seven species, with an average of four species per site (Datas S5 and S6). Regarding population-level diversity, F. validus mtCOI haplotypic diversity and nucleotide diversity were typically lower than that of F. erichsonianus (mean = 0.57 [SD = 0.20] vs. 0.63 [0.17]; mean = 0.002 [0.001] vs. 0.003 [0.002], respectively Table 2). Likewise, phylogenetic diversity (PD) values, measured in branch length units of substitutions per site, were generally lower for F. validus than F. erichsonianus (mean = 0.006 [0.004] vs. 0.015 [0.014], respectively; Table 2). Conversely, the ISSR-based proportion of polymorphic loci (PPL) showed close equivalence between the two focal species (mean PPL = 0.87 [0.11] vs. 0.85 [0.07] for F. validus and F. erichsonianus, respectively; Table 2).

Species richness, abundance, and genetic diversity correlations

Pearson correlation tests showed no significant correlations between species richness and genetic diversity for any of the mtCOI- or ISSR-based diversity metrics for F. validus or for F. erichsonianus (Tables 5A, 5B). Additionally, Pearson correlation tests showed no significant correlations between cumulative multispecies crayfish abundance and species richness (P = 0.46; Table 5C). As such, neither SGDC nor MIH were supported.

An AGDC was evident in one focal species (Table 5A). For F. validus, Pearson correlation tests showed a negative relationship between cumulative multispecies abundance and genetic diversity for two of the three mtCOI-based metrics (hd: r = −0.45, P = 0.03; π: r = −0.36, P = 0.09), although only the abundance-hd correlation was significant at the 0.05-level. However, the Pearson correlation test showed no significant correlation between abundance and ISSR-based PPL (Table 5A). For F. erichsonianus, Pearson correlation tests showed no significant relationship between cumulative multispecies abundance and genetic diversity for any mtCOI (hd, π, PD) or ISSR (PPL) diversity measure (r = −0.10–0.27, P > 0.05; r = 0.07, P = 0.75, respectively) (Table 5B).

Association between environmental characteristics, community-level diversity, and population-level diversity

Species richness was significantly correlated with stream width, explaining 36% of species richness variation (Table 6). However, this correlation varied depending on stream type (fragmented vs. connected) (Fig. 2A). There was a positive relationship between species richness and stream width in unimpounded streams, but a negative relationship in impounded streams.

Cumulative multispecies abundance was significantly correlated with stream width, LWD, and substrate size, which explained 97% of crayfish abundance variation (Table 6). Cumulative multispecies abundance increased with decreasing stream width, and its relationship with LWD and substrate size varied depending on stream type (Figs. 2B, 2C). There was a significant negative relationship between cumulative multispecies abundance and amount of LWD in unimpounded streams, but no relationship with LWD in impounded streams (Fig. 2B). Additionally, there was a significant positive relationship between cumulative multispecies abundance and substrate size in unimpounded streams, but no relationship in impounded streams (Fig. 2C).

Stream width was significantly positively correlated with F. validus nucleotide and haplotype diversity, explaining 19% and 22% of their variation, respectively (Table 6). There was no significant relationship between ISSR-based PPL and stream characteristics. Additionally, there were no significant relationships between any F. erichsonianus genetic diversity metrics and stream characteristics (Table 6).