Effects of environment and metacommunity delineation on multiple dimensions of stream fish beta diversity

Introduction Beta diversity represents changes in community composition among locations across a landscape. While the effects of human activities on beta diversity are becoming clearer, few studies have considered human effects on the three dimensions of beta diversity: taxonomic, functional, and phylogenetic. Including anthropogenic factors and multiple dimensions of biodiversity may explain additional variation in stream fish beta diversity, providing new insight into how metacommunities are structured within different spatial delineations. Methods In this study, we used a 350 site stream fish abundance dataset from South Carolina, United States to quantify beta diversity explainable by spatial, natural environmental, and anthropogenic variables. We investigated three spatial delineations: (1) a single whole-state metacommunity delineated by political boundaries, (2) two metacommunities delineated by a natural geomorphic break separating uplands from lowlands, and (3) four metacommunities delineated by natural watershed boundaries. Within each metacommunity we calculated taxonomic, functional, and phylogenetic beta diversity and used variation partitioning to quantify spatial, natural environmental, and anthropogenic contributions to variations in beta diversity. Results We explained 25–81% of the variation in stream fish beta diversity. The importance of these three factors in structuring metacommunities differed among the diversity dimensions, providing complementary perspectives on the processes shaping beta diversity in fish communities. The effect of spatial, natural environmental, and anthropogenic factors varied among the spatial delineations, which indicate conclusions drawn from variation partitioning may depend on the spatial delineation chosen by researchers. Discussion Our study highlights the importance of considering human effects on metacommunity structure, quantifying multiple dimensions of beta diversity, and careful consideration of user-defined metacommunity boundaries in beta diversity analyses.


Introduction
Freshwater systems are threatened by many anthropogenic effects including habitat fragmentation, land cover change, and disruption of the natural flow regime (Vörösmarty et al., 2010;Heino, 2013). Human effects on taxonomic and functional composition of freshwater communities are well documented (Poff, 1997;Allan, 2004) to the point that they are foundational to the paradigm of contemporary environmental assessment (Karr, 1981;Bailey et al., 2004). Landscape level studies have assessed the dependence of measured shifts in fish communities upon current and past land use changes, with some evidence of generalization in thresholds for urbanization and agriculture effects on compositional and richness changes (Chen and Olden, 2020). While land use variables provide the first evidence that urbanization and agriculture change fish communities, pinpointing the specific variables associated with these land uses can give additional insight into what environmental monitors should be measured when completing bioassessments (Waite et al., 2021;Iacarella, 2022). For example, Stoczynski et al. (2021) found that dam density affected how species replaced one another across the landscape, Ortega et al. (2021) identified impervious surfaces as a proxy for urbanization, and Waite et al. (2021) highlighted how contaminants in water and sediment were just one of multiple urbanization and agricultural stressors. Thus, disentangling anthropogenic effects from 'natural' environmental variables in metacommunity analyses may increase our ability to explain factors structuring beta diversity within metacommunities (Duarte et al., 2018;Gianuca et al., 2018). However, few studies have examined stream fish metacommunities from the lens of anthropogenic effects (Erős et al., 2012;Borges et al., 2020;Stoczynski et al., 2021).
Biodiversity is represented by three dimensions: taxonomic, phylogenetic, and functional diversity; patterns of beta diversity can be unique for each dimension. Taxonomic diversity indexes counts and proportions of taxa (typically species or genera). Alternatively, functional diversity represents the range of traits present in a location, regardless of species (Tilman, 2001), while phylogenetic diversity refers to the degree of variation in evolutionary history in a community (Webb et al., 2002). Most analyses of beta diversity have focused solely on taxonomic diversity and have not considered that co-occurring species may be from different phylogenetic positions and/or have different functional roles within the community (Jarzyna and Jetz, 2016). All species have traits which correspond to their roles and habitat preferences within a community. Functional beta diversity describes differences in these traits between communities (Wootton, 2012;Leibold and Chase, 2017). Phylogenetic beta diversity allows us to investigate evolutionary differences among communities (Weinstein et al., 2014;Nakamura et al., 2017Nakamura et al., , 2021Li F. et al., 2021;Li Z. et al., 2021). Thus, metacommunity approaches can also use functional and phylogenetic beta diversity to complement taxonomic beta diversity to provide additional insights and a more complete diversity picture (Anderson et al., 2011;Leibold and Chase, 2017;Cai et al., 2019;Jia et al., 2020).
Riverine fishes provide a unique study system for testing metacommunity concepts. The dendritic structure of river networks restricts movement of fishes to the stream network (Brown and Swan, 2010;Tonkin et al., 2018). This framework creates a dual habitatconnectivity gradient in which (a) some headwater streams contain more isolated communities and thus are more prone to the effects of species sorting and stochastic effects, and (b) other headwaters with similar habitats are well connected to mainstem rivers resulting in a greater influence from dispersal processes (Brown and Swan, 2010;Hitt and Angermeier, 2011;Heino et al., 2015a;Henriques-Silva et al., 2019). Moreover, watersheds are arranged hierarchically; small watersheds are nested within larger watersheds separated by discrete spatial boundaries caused by mountainous and marine areas, creating discrete spatial units for delineating metacommunities as statistical units in large-scale analyses (Jelinski and Wu, 1996;Patrick and Yuan, 2019). While watershed boundaries create distinct dispersal barriers for fishes that naturally delineate metacommunities, data availability or constraints often lead researchers to delineate metacommunities based on other features such as political boundaries or natural geomorphic breaks across the landscape such as ecoregions (Cai et al., 2019;Henriques-Silva et al., 2019;Lech and Willig, 2020). However, few studies have investigated how the choice of spatial delineation may affect the outcome and interpretation of results. Conducting separate analyses on the same dataset with different delineation scenarios representative of existing metacommunity studies will provide important methodological context for conducting metacommunity analyses.
In this study, we used an extensive stream fish abundance dataset covering the state of South Carolina, USA, to quantify variation in the three dimensions of beta diversity attributable to spatial, natural environment, and anthropogenic effects. Our first objective was to quantify variation in stream fish taxonomic, functional, and phylogenetic beta diversity. We expected to see more explained variation in functional and phylogenetic beta diversity because of the association between species traits and common resource utilization. Additionally, relatedness between those traits may be masked by taxonomic diversity in functionally redundant communities (Leibold and Chase, 2017;Gianuca et al., 2018;Cai et al., 2019). Our second objective was to measure the relative variation explained by spatial, natural environmental, and anthropogenic variables. We expected that the use of anthropogenic variables more specific than land use would explain appreciable amounts of variation in all three beta diversity dimensions because of the extensive alterations to stream systems in South Carolina in recent centuries (Allan, 2004;Walters et al., 2005;Wenger et al., 2008;Marion et al., 2015;Langerhans and Kern, 2020). We conducted analyses for these objectives separately for three different spatial delineations: one based on the entire state, one based on a natural geomorphic ecoregion break, and one based on discrete watershed boundaries. In doing so, we highlight the importance of a priori methodological decisions when conducting beta diversity analyses.

Study area
Located in the southeastern United States, the state of South Carolina contains five level III ecoregions: The northwestern part of the state (thereon referred to as the upstate) comprises the Blue Ridge to the far northwest corner and the Piedmont to the east. The southeastern part of the state (thereon referred to as the lowlands) comprises the Southeastern Plains, Middle Atlantic Coastal Plains, and Southern Coastal Plain (Figure 1; Omernik, 1987). The division Frontiers in Ecology and Evolution 03 frontiersin.org between the Piedmont and Southeastern Plains, referred to as the fall line, divides streams into two distinct types (Paller, 1994): upstate and lowland streams. Upstate streams (above the fall line) have well defined channels with stable flow, which contain heterogenous habitat with riffles, runs, glides, and pools containing a mixture of sand and coarse substrates (Wallace et al., 1992;Denison et al., 2021a,b). A broad transitional zone containing relatively steep elevational change separates the upstate and lowland streams. Lowland streams (below the fall line) tend to be poorly defined with more unstable flow regimes and are more homogeneous with habitats characterized by pools, runs, and fine substrates (Paller, 1994;Rohde et al., 2009;Marion et al., 2015). Because of these geomorphic differences, and the fact that the upland rivers have been isolated from each other over longer evolutionary time, fish metacommunities above and below the fall line differ taxonomically. Due to watershed separation, fish assemblages in different upstate watersheds contain more phylogenetically related species causing them to be more taxonomically different but functionally similar. In contrast, lowland streams contain a greater number of shared species due to watershed connection during colonization following over bank flooding events and more recent hydrologic connectivity associated with sea level rise and fall (Rohde et al., 2009 Figure 2A). Detailed methods can be found in Scott et al. (2009), but briefly, sites were selected randomly from streams across the state to proportionally represent all ecoregions and river basins.
All samples were collected using backpack electrofishing. Sites sampled above the fall line used one electrofishing pass equal to 30 times the average wetted width (m) with most habitats sampled in an upstream direction while an 2.5-3 m seine with mesh size 6.4 mm was used to more effectively capture benthic species in riffles habitats. Sites below the fall line were sampled using three electrofishing passes with a reach length equal to 20 times the average stream width. Average total effort with standard deviation was 5361.04 ± 3791.36 s, 1.95 ± 0.95 number of passes, and 121.92 ± 42.73 meters of sampled distance. Fishes were netted, measured, identified to the species level, and then returned to the stream (Scott et al., 2009). This dataset has been used in other studies to identify statewide community patterns (Marion et al., 2015;Denison et al., 2021a,b;Bower et al., 2022). The dataset contains 350 sites representing 101 species from 51 genera, 17 families, and 14 orders.

Trait and phylogeny data
We collected 16 available traits using the FishTraits database (Frimpong and Angermeier, 2009; Table 1), and calculated shape factor (body depth/body length) and swim factor (caudal peduncle depth/caudal fin depth) using high quality fish images (Supplementary Table S1 and Figure 2A). Traits were chosen to ensure we were capturing as much of a species' niche dimension as possible. The use of traits covering many niche dimensions provides an overview of community structure (Trisos et al., 2014;Fitzgerald et al., 2017). Traits such as shape and swim factor provide information on habitat niche dimension and dispersal capabilities (Lamouroux et al., 2002;Sagnes and Statzner, 2009;Pease et al., 2012). Any gaps in fish traits were filled in through literature searches of the species in question or the nearest relative (Near et al., 2011), and if literature searches failed to provide information (few cases), then we used the 'missForest' package (Stekhoven, 2011) in R version 4.0.2 (R Core Team, 2021) to impute the missing trait values using those trait values from individuals within the same genus. This approach was only used for a handful of traits within a few species. The 'missForest' package uses random forest models to predict and impute missing data based on relationships within that fish family. Imputations relied on 1,000 iterations and 100 trees to fill in the few missing traits.
The phylogenetic tree used for our study was trimmed from an extensive phylogeny of Actinopterygii, the class containing the majority of fish diversity (Rabosky et al., 2018; Figure 2A). Two species were manually added to the tree, Notropis cummingsae (Myers, 1925) and Pteronotropis stonei (Fowler, 1921), from reading relevant phylogenetic research (Hollingsworth et al., 2013;Mayden and Allen, 2015). This tree contains phylogenetic distances along the branch lengths and can thus be used to calculate phylogenetic distance among species. The trimmed tree is rooted at the ancestor between the outgroup of the Neopterygiian Bowfin (Amia calva; Linnaeus, 1766) and Spotted Gar (Lepisosteus osseus; Linnaeus, 1758), which represent the two most ancestral members of Actinopterygii, relative to other species (Supplementary Figure S1).

Calculating taxonomic, functional, and phylogenetic beta diversity
We calculated a response matrix for fish beta diversity in each diversity dimension (taxonomic, functional, and phylogenetic; South Carolina, USA with sampling site locations and watershed color delineations used in analysis. The fall line crosses through the middle of the state and divides streams into two distinct types: with the upstate to the northwest and lowland to the southeast. Frontiers in Ecology and Evolution 04 frontiersin.org Figure 2B). We first Hellinger transformed fish abundance data (Legendre and Gallagher, 2001). We then calculated Bray-Curtis taxonomic beta diversity for all pairs of sites using the beta.pair.abund function in the 'Betapart' package (Baselga and Orme, 2012). We used Bray-Curtis dissimilarity because this beta diversity metric is widely used with abundance data and is the most robust to taxonomic error and geographic undersampling (Schroeder and Jenkins, 2018).
We calculated functional distances (species x species matrix) with the gowdis function in the 'FD' package (Laliberté et al., 2010) for all trait variables using Gower distance because we have a mixture of continuous and categorical traits. We calculated phylogenetic distances (species x species matrix) using the cophenetic.phylo function in the 'ape' package (Paradis and Schliep, 2019). We weighted the phylogenetic and functional distances by the Hellinger Flowchart for statistical analysis completed by collection of stream fish abundance data (A), preparation of response (B) and explanatory matrices (C), preparation for variation partitioning (D), and final variation partitioning run (E) over our spatial delineations and diversity dimensions. ANTH, anthropogenic; ENV, Natural Environmental; PCNM, Principal coordinate neighbor matrices; EV, Environmental Variable.
Frontiers in Ecology and Evolution 05 frontiersin.org transformed abundance data (site x species matrix) to calculate for pairwise functional and phylogenetic beta diversity (site x site matrix) for all sites using the comdist function from the 'picante' package (Kembel et al., 2010). We conducted principal coordinate analysis (PCoA) for each beta diversity matrix (response matrices) using the pcoa function from the 'ape' package to have all data in a similar reduced dimensionality format for variation partitioning using Bray Curtis dissimilarity . Using all axes of the PCoA could lead to bias interpretation of the results, so we limit variation partitioning to the first and second axes. To visualize beta diversity in the diversity dimensions, red, green, blue (RGB) color plots were produced using the metaMDS function from the 'vegan' package to first perform nonmetric multidimensional scaling (NMDS) followed by the recluster.plot.sites.col function from the 'recluster' package (Dapporto et al., 2013), to place sites in the RGB color space. Similar colors indicate sites with similar species composition, as colors become more distinct the beta diversity between those two sites increases.

Explanatory matrices preparation 2.4.1. Environmental and anthropogenic data
Water quality measurements, collected at each site on the date of fish sampling, included water temperature (°C), dissolved oxygen (mg/L), conductivity (μs/cm), pH, and salinity (ppm), while a turbidimeter measured turbidity (NTU). Water quality measurements were taken from the area of representative flow and from the middle of the water column (Scott et al., 2009). Instream habitat was also measured on the date of fish sampling; details for the protocols are found in Scott et al. (2009). Briefly, researchers measured depth (m), current velocity (m/s), and substrate type and diameter (mm) along 50 points in the sample transect after community assessment and calculated an average site value for each measurement. These habitat measurements were taken from equal portions of instream habitat types (runs, riffles, and pools). Current velocity was measured at 60% of the depth from the surface and 40% from the bottom. Substrate diameter measurements used the Wolman pebble count (Wolman, 1954). Inorganic material was registered as fine clay or sand, bedrock, or measuring particle to nearest mm, while organic material was categorized into one of five categories based on size and substance: fine particulate organic matter, coarse particulate organic matter, fine woody debris, large woody debris, or aquatic vegetation.
The stream-catchment (StreamCat) database maintained by the U.S. Environmental Protection Agency (USEPA) contains environmental and anthropogenic data for stream segments throughout the United States (Hill et al., 2016). StreamCat is linked with the National Hydrology Dataset Plus Version 2 (NHDPlusV2) to obtain natural environmental and anthropogenic variables for the associated catchment or for the full upstream watershed for millions of stream segments in the United States. We linked the StreamCat database with the 350 SCDNR sites to obtain environmental or anthropogenic variables associated with the stream segment each site was located on. Natural environmental and anthropogenic variable definitions, codes, and units can be found in Table 2, and are summarized in Supplementary Table S2. We performed a collinearity check on all anthropogenic and environmental data by (1) performing a variance inflation factor (VIF) check on all variables and (2) finding variable pairs with a correlation of >0.5, followed by one of those variables being removed (Supplementary Table S3). We used forward selection process using the forward.sel function in the 'adespatial' package (Dray et al., 2018), to select the variables that explained the most variation in the PCoA axis for each of our metacommunities using a two-stop criterion: the adj R 2 value for the global model and the significance value of 0.05 (Blanchet et al., 2008; Figure 2C). Forward selection first finds the explanatory variable that explains the greatest variation ( adj R 2 ). Then, the next best variable is found until the combined adj R 2 values for all explanatory variables reaches the global adj R 2 for each explanatory matrix. Selected natural environment and anthropogenic variables for each metacommunity can be found in Table 3.

Generating spatial variables
Vectors from principal coordinate neighbor matrices (PCNMs) define spatial scales, where the first vectors represent the largest spatial structure, and the last vectors represent finer spatial structure (Borcard  and Legendre, 2002;Griffith and Peres-Neto, 2006). PCNMs chosen with small values represent large geographic, long temporal time scales. As PCNM number increases, these values indicate smaller geographic time scales including dispersal within a watershed or even between nearby stream reaches. Because of the dendritic layout of watersheds, using Euclidian distances for spatial variables would not produce the most realistic scenario when comparing distance among sites. We used an origin-distance cost matrix within ArcMap 10.6 (ESRI, 2018) to calculate the fluvial distance among sites. We then calculated PCNM vectors using the pcnm function in the 'vegan' package (Oksanen et al., 2015). When using PCNM vectors, the use of all vectors can lead to possible misunderstandings in result interpretation, so the appropriate vectors were extracted using forward selection with each of the response matrices in each diversity dimension for the corresponding metacommunity similar to anthropogenic and environmental variables ( Figure 2C). Selected spatial variables for each metacommunity can be found in Supplementary Table S4.

Spatial delineations
We conducted analyses at three different spatial delineations representing logical breaks in the landscape based on natural watershed and geomorphic boundaries, as well as at the whole state level. By delineating our dataset into metacommunities of different spatial sizes, we can understand how stream fish communities are structured by space, natural environment, and anthropogenic variables in metacommunities with different underlying characteristics. The largest delineation analyzed was the whole state, which included all 350 sites and represents an artificially delineated metacommunity because the state boundary of South Carolina crosses both watersheds and ecoregions (Figure 1). While this delineation does not represent a truly natural metacommunity, conducting analyses at this delineation captures any patterns and variability from the entire dataset, offering a regional approach comparable to other studies (Wang et al., 2011;Heino et al., 2018;Cai et al., 2019). We also split the dataset into the four natural metacommunities delineated by watershed boundaries: the Savannah (58 sites); the Ashepoo, Combahee, and Edisto (ACE; 57 sites); the PeeDee (83 sites); and the Santee (152 sites). We also had two metacommunities delineated by the geomorphic properties of the upstate (167 sites) and lowlands (183 sites; Figure 1). All metacommunities designated for the analysis have enough sites to allow for robust analysis (Heino et al., 2015b).

Quantifying importance of space, environment, and anthropogenic influence on beta diversity
We used adj R 2 to determine the number of axes to be used in the variation partitioning ( Figure 2D). If the adj R 2 value for the first axis was larger than the first and second axes, just the first axis was sufficient to use in the final variation partitioning analysis (Anderson and Willis, 2003;Gianuca et al., 2018). Using only the first PCoA axis showed the highest R 2 adj value for all RDAs and was thus used for the variation partitioning analysis except for the metacommunity above the fall line, which required the use of the first and second axes to obtain the maximum R 2 adj value (Table 4). We used variation partitioning to determine the spatial, natural environmental, and anthropogenic variables explaining the most variation in the PCoA axis representing taxonomic, functional, and phylogenetic beta diversity separately for each metacommunity ( Figure 2E). Because we used fluvial distance instead of Euclidean distance, we conducted a distance-based redundancy analysis (db-RDA) using the varpart function in the 'vegan' package to determine explained variation in our metacommunities.

Results
Variation partitioning explained 25-81% of variation in stream fish beta diversity, depending on the metacommunity delineation and dimension of beta diversity (Figure 3 and Table 4). All db-RDA analyses resulted in significant explained variance by the explanatory variables (all p-values <0.001; Table 5). We observed a strong phylogenetic signal in our fish traits for each metacommunity (Supplementary Table S5). We found that overall explained variation in stream fish beta diversity, the relative contribution of the three suites of variables to beta diversity, and the contribution of individual variables, all differed depending on how metacommunities were delineated-whether the entire state, above and below the fall line, or based on watershed boundaries. Breaking up our study region into the upstate and lowland metacommunities resulted in the least explained variance; however, the explained variance for anthropogenic portions alone was much higher than when metacommunities were distinguished based on major watersheds or the whole state ( Figure 3 and Table 4).
Frontiers in Ecology and Evolution 09 frontiersin.org
Frontiers in Ecology and Evolution 10 frontiersin.org Red, green, blue (RGB) plots of fish beta diversity over taxonomic (A), functional (B), and phylogenetic (C) diversity dimensions. Similar colors reflect communities with similar species compositions, while differing colors reflect increased differences in species composition between sites.
(61%), followed by the functional (56%), and phylogenetic dimensions (46%). The whole state was the only spatial delineation in which we observed no influence of solely anthropogenic variables on stream fish beta diversity for any of the dimensions. All three diversity dimensions were structured by stream order (Table 3). Functional and taxonomic beta diversity were additionally structured by base flow index, while taxonomic beta diversity was also structured by stream temperature and percent of mixed forest within the watershed (Table 3). Spatial vectors explaining beta diversity across each diversity dimension ranged between vectors explaining large and small spatial scales (Supplementary Table S4).

Fall line delineation
RGB plots showed two distinct regions of taxonomic, functional, and phylogenetic beta diversity roughly separated by the fall line ( Figure 4). However, the differences between the upstate and lowland regions of South Carolina in phylogenetic beta diversity were smaller than between the functional or taxonomic dimensions.
The explained variation for the upstate and lowlands was less than the other spatial delineations (Figure 3 and Table 4). However, we observed the largest significant explained variance for just anthropogenic variables on the functional (Adj R 2 = 0.22; F df = 13.75 4 , p = 0.001) and phylogenetic (Adj R 2 = 0.19; F df = 21.87 2 , p = 0.001) dimensions in the upstate compared to other metacommunities ( Figure 3B). Significant fractions of explained variance were observed across each diversity dimension for anthropogenic variables (Tax upstate : Adj R 2 = 0.06; F df = 4.87 5 , p = 0.002; Tax lowland : F df = 2.94 4 , p = 0.016; Fun lowland : F df = 5.00 5 , p = 0.001; Phy lowland : F df = 8.06 2 , p = 0.001). In the upstate, beta diversity for all three diversity dimensions was explained by the amount of manure application to agricultural land (Table 3). In total, taxonomic beta diversity was explained by five variables, functional beta diversity by four variables, and phylogenetic beta diversity by two variables (Table 3). In the lowlands, beta diversity in all three dimensions was explained by dam density and road crossings. Functional beta diversity was explained by five variables, taxonomic by four, and phylogenetic by two variables (Table 3).
Only variation in beta diversity in the upstate showed any influence from natural environmental variables. Taxonomic beta diversity (Adj R 2 = 0.02; F df = 4.28 2 , p = 0.013) was structured by stream temperature and stream order, and functional beta diversity (Adj R 2 = 0.03; F df = 7.54 1 , p = 0.009) was structured by runoff (Table 3). While two environmental variables were picked for the lowlands, neither of these variables had a significant influence on beta diversity in any dimension.
Spatial variables showed greater importance within the lowlands compared to the upstate, though there was a large spatial fraction for the taxonomic dimension in the upstate (Adj R 2 = 0.26; F df = 4.60 2 , p = 0.01). Variance explained by spatial factors was greatest for the Frontiers in Ecology and Evolution 11 frontiersin.org taxonomic dimension (Adj R 2 = 0.34; F df = 6.30 23 , p = 0.001), followed by the functional (Adj R 2 = 0.29; F df = 7.08 17 , p = 0.001) and phylogenetic (Adj R 2 = 0.24; F df = 8.29 11 , p = 0.001) dimensions in the lowlands ( Figure 4C). The lowlands had the most PCNM vectors selected to explain beta diversity variation compared to all other metacommunities analyzed. The upstate only had four vectors selected between the three diversity dimensions with PCNM1 and PCNM8 explaining taxonomic beta diversity, PCNM3 selected for but not significant in explaining functional beta diversity, and PCNM16 explaining phylogenetic beta diversity (Supplementary Table S4).

Watershed delineation
Significant portions of beta diversity were explained by anthropogenic variables alone in each watershed metacommunity (Figure 3). No anthropogenic variables were significant in explaining variation of phylogenetic beta diversity within any of the watershed metacommunities. Functional beta diversity showed the largest portion of explained variation from just anthropogenic variables for the Savannah (Adj R 2 = 0.07; F df = 5.14 3 , p = 0.006) and ACE (Adj R 2 = 0.05; F df = 5.41 2 , p = 0.005) watersheds ( Figures 3E,G). Stream biological condition, percent cropland and road crossings contributed toward explained variation in beta diversity of the Savannah watershed, while maximum dam storage and road crossings explained the beta diversity variation in the ACE watershed (Table 3). Taxonomic beta diversity showed the largest portion of explained variance from just anthropogenic variables for the ACE watershed (Adj R 2 = 0.05; F df = 4.26 2 , p = 0.02).
Stream fish beta diversity explained only by natural environmental variables was seen in each watershed metacommunity. We observed the largest fractions of variation explained by only environmental variables in all three diversity dimensions for the Santee (Tax: Adj R 2 = 0.26; F df = 18.27 10 , p = 0.001; Fun: Adj R 2 = 0.17; F df = 14.30 8 , p = 0.001; Phy: Adj R 2 = 0.20; F df = 16.18 6 , p = 0.001) and ACE watersheds (Tax: Adj R 2 = 0.36; F df = 6.99 7 , p = 0.001; Fun: Adj R 2 = 0.19; F df = 10.01 5 , p = 0.001; Phy: Adj R 2 = 0.38; F df = 15.68 3 , p = 0.001). The Savannah showed significant environmental signal for the functional (Adj R 2 = 0.15; F df = 5.67 6 , p = 0.001) and taxonomic (Adj R 2 = 0.13; F df = 5.81 7 , p = 0.001) dimensions, while only the taxonomic dimension (Adj R 2 = 0.05; F df = 4.81 5 , p = 0.001) showed an environmental signal for the PeeDee watershed ( Figure 3 and Table 4). In the Santee watershed many variables were similar in explaining beta diversity, but the taxonomic and functional beta diversity dimensions were additionally explained by percent aquatic vegetation and base flow index, while taxonomic dimension also included percent clay substrate and stream order as influential variables (Table 3). Environmental variables differentially explaining beta diversity within the ACE watershed include stream depth, percent total wetlands, and percent clay substrate for taxonomic beta diversity, fine woody debris for functional beta diversity, and percent other substrates for phylogenetic beta diversity (Table 3).
Spatial variables contributed toward the explained variation for all watershed metacommunities. The highest explained variance from just the spatial fraction was found in the taxonomic dimension of the PeeDee watershed (Adj R 2 = 0.11; F df = 4.23 9 , p = 0.001), followed by the phylogenetic dimension of the ACE watershed (Adj R 2 = 0.05; F df = 3.40 3 , p = 0.031). PCNM2, which represents a large spatial scale was selected for all watershed metacommunities except the ACE. The ACE watershed was unique because no spatial variables explained taxonomic beta diversity. Additionally, the PCNMs selected for the functional and phylogenetic dimension in the ACE basin mostly related to fine spatial scales (Table 4). The PeeDee watershed had the most PCNMs selected to explain each of its diversity dimensions, however the PeeDee was one watershed where variation explained by just space, environment, or anthropogenic factors was only seen in the taxonomic dimension. Instead, distinguishing between how space, natural environment, and anthropogenic variables explained beta diversity in the PeeDee was unsuccessful due to the large fraction of explained variation from the overlap of all three matrices ( Figure 3D). Large portions of explained variation due to overlap between two or more matrices were also observed for all three dimensions in the Savannah and Santee watersheds.

Discussion
We investigated how spatial, natural environmental, and anthropogenic variables influence taxonomic, functional, and phylogenetic beta diversity for stream fish communities across South Carolina, USA. Our measured variation explaining beta diversity ranged from 25-81%, which is high relative to other studies investigating similar questions (Gianuca et al., 2018;Burdon et al., 2019;Benone et al., 2020). Many freshwater studies investigating drivers of beta diversity using variation partitioning have explained around 30% of the variation for a variety of taxa including stream macroinvertebrates, ostracods, fish, and macrophytes (Alahuhta and Heino, 2013;de Campos et al., 2018;Huang et al., 2019;Krynak et al., 2019). Accounting for biogeographical features of metacommunities while investigating taxonomic and functional beta diversity of Amazonian stream fishes resulted in explained community variance above 60% (Benone et al., 2020). Studies that added anthropogenic effects also increased the explained variance, similar to our results, where researchers described over 60% of the variance in freshwater macroinvertebrate community structure in Switzerland (Burdon et al., 2019). Variables including road crossings, impervious surfaces, fertilizer application, and forest loss were selected over generic land use variables. Our study, in conjunction with previous studies, indicates that the unexplained variance of past studies was due partly to information that the inclusion of anthropogenic effects provides. Our research highlights the importance of including anthropogenic variables in beta diversity studies and how these variables can provide better context for understanding human effects on stream fish beta diversity.
Quantifying beta diversity in functional and phylogenetic dimensions more fully captures the factors driving beta diversity by introducing traits and relatedness among taxa that is not captured by using taxonomy alone. We explained more variation using taxonomic beta diversity relative to the other diversity dimensions for all but one metacommunity. Large regional species pools provided a greater difference in species composition among sites to account for the increased explained variance of the taxonomic beta diversity compared to differences in traits or evolutionary relatedness, indicating higher redundancy of traits and evolutionary history for stream fish compared to other taxa. Substantial differences in explained variation between the diversity dimensions has been observed among macroinvertebrates (Gianuca et al., 2018; Frontiers in Ecology and Evolution 12 frontiersin.org Bispo et al., 2021;Li Z. et al., 2021) and macrophytes (Garcia-Giron et al., 2019). Gianuca et al. (2018) observed up to three times more explained variation in functional and phylogenetic beta diversity compared to taxonomic diversity of zooplankton metacommunities over an urbanization gradient. Therefore, the inclusion of more than one diversity dimension can not only provide more explained variation but can offer differential insights on the factors structuring beta diversity. For example, functional diversity can provide insights into how environmental stressors effect community structure in stream fish (Jia et al., 2020;Zhang et al., 2020). The only other study to our knowledge which also investigated all three diversity dimensions in stream fish observed that instream habitat variables had a weak relationship with taxonomic beta diversity, while phylogenetic beta diversity values failed to show deterministic structuring processes (Roa-Fuentes et al., 2019). In contrast, our study showed phylogenetic beta diversity was best explained by varying degrees of spatial, natural environmental, or anthropogenic factors based on the spatial delineations within our study. Thus, the amount of explained variation in beta diversity is context dependent because even similar taxa can show different outcomes. Similar studies conducted on other regions of the world with varying climates, geography, and evolutionary histories would improve our understanding of context dependency. We conducted analyses separately for three spatial delineations that are representative of commonly applied approaches-delineations based on political boundaries, natural geomorphic breaks, and watershed boundaries. Our results indicate that the choice of spatial delineation affects the outcome of and interpretation of variation portioning analyses, resulting in alternative interpretations of overall variance explained in metacommunity structure, as well as relative variable contributions to changes in beta diversity. Analyses based on the whole state delineation (based on political boundaries) emphasized the importance of natural environmental factors and spatial gradients indexing large-scale dispersal. Alternatively, analyses based on the fall line delineation provided different inference in which the importance of anthropogenic and spatial factors far outweighed the importance of environmental factors. This delineation also shed light on the fact that the relative contributions of anthropogenic and spatial dispersal gradients differ appreciably between the two regions, indicating different assembly rules structure stream fish metacommunities in different parts of the state. Finally, analyses based on watershed boundary delineations highlighted the synergistic effects of anthropogenic factors, spatial dispersal gradients, and natural environmental variables for structuring stream fish metacommunities. The latter analyses exposed differences between the three dimensions of stream fish beta diversity that were not evident in analyses based on the other two delineations. Many metacommunity studies are often applied to datasets that were collected without the expressed intent of quantifying beta diversity across the landscape, and are subject to variation according to methodological attributes of the dataset such as sampling density and intensity (Stoczynski et al., 2021). Likewise, our results show that user-defined choices of how metacommunities are defined a priori have important bearing on model results and ecological inference, and should be considered carefully before conducting metacommunity analyses.
Differences in the results of analyses based on different spatial delineations bring to the forefront how past geologic change affects community assembly processes across the state. The fall line represents the farthest inland the coast has been since the Cretaceous period when the sea levels were at their highest. The Appalachian uplands including the majority of the Piedmont ecoregion has been continuously above sea level since that time, which predates the modern ostariophysan fauna (Cavender, 1998). Accordingly, the faunas of river basins above the fall line have had an essentially uninterrupted evolutionary history of isolation and speciation, with geologic events such as stream capture permitting some dispersal across basins (Swift et al., 1986). This uninterrupted evolutionary history is evident in the occurrence of numerous sister species across these basins (e.g., in the families Leuciscidae and Percidae), and is reflected in the fact that factors affecting the three dimensions of beta diversity differed among basins and between the Upstate and Lowlands. In contrast, the lowering of sea level during periods of glaciation including the Pleistocene would have allowed gene flow in coastal plain drainages as rivers merged further out on the Atlantic Slope, with flood events in these low topographic relief areas enhancing interbasin dispersal. Rohde et al. (2009) noted that similarity was higher in the faunas of the Santee and Pee Dee basins, and in the Savannah and ACE basins than between those two basin pairs, indicating these watershed pairs shared more recent connections permitting dispersal. As recently as 20,000 years ago, potential basin connections occurred as sea level retreated to roughly 120 m below current sea level, placing the coast near the edge of the Continental Shelf (Miller et al., 2020). Fluctuations in shoreline location through geologic time resulted in observed differences in taxonomic beta diversity, including evidence that anthropogenic effects are greater in the upland region due to less historic mixing of communities compared to below the fall line on the Atlantic Coastal Plain.
The rise and fall of the coast resulted in different geological properties found above and below the fall line. Providing context to the geographical properties of our study region is important for allowing potential comparison between studies in different regions of the world. The South Carolina lowlands are characterized by a mosaic of highly populated areas and unpopulated more agricultural areas (Homer et al., 2015). Lowland streams are dominated by harsh environmental conditions, frequent exposure to extreme weather, and poorly defined stream channels. These streams are also connected more closely to large mainstem rivers and reservoirs and share more species among watersheds due to historical sea level change (Marion et al., 2015;Miller et al., 2020;Denison et al., 2021b). Moreover, the species inhabiting the lowland streams tend to have more dispersal-associated traits and have populations that respond more to connectivity than local-scale abiotic factors (Midway and Peoples, 2019). Humans also historically connected lowland watersheds during the canal boom of the 1790s-1830s that saw over 3,200 km of cross-watershed canal construction, as well as widespread channelization and woody debris removal (Kapsch, 2010). These disturbances have likely contributed to a landscape legacy of increased connectivity, which we observe in the spatial signals in the lowland metacommunity across each diversity dimension (Robinson et al., 2002;Wenger et al., 2008). When comparing rivers, geomorphic variables including the width and gradient of streams dictated community structure with constrained river valleys having more varied rocky substrate and log jams to provide fish cover and wide river valleys showing increased riparian vegetation, slower water velocities, and less woody debris (Walters et al., 2003;Shields et al., 2021). In contrast to the lowlands, beta diversity in the upstate was affected more by anthropogenic and natural environmental factors. In upland regions of the eastern United States, stream fish community structure is often Frontiers in Ecology and Evolution 13 frontiersin.org dominated by more specialized species in less disturbed watersheds, which transition to communities dominated by generalists as watersheds become more disturbed (Scott and Helfman, 2001;Petersen et al., 2021), which would result in high beta diversity between the two ends of this continuum. Within the upstate, increased disturbance associated with higher gradient streams can result in increased sediment transport, flash flooding, channel incision, and nutrient or pollution pulses detected by the anthropogenic signal (Marion et al., 2015). The disparity in the anthropogenic signal between the upstate and lowlands may be due to the sensitivity of upstate fishes as well as the lower connectivity for recolonization after disturbance events. The geomorphic characteristics describing streams above and below the fall line may account for the lack of a large environmental signal observed in the fall line spatial delineation compared to the statewide and watershed spatial delineations (Warren et al., 2000;Smogor and Angermeier, 2001). While we explained a large portion of variation in stream fish beta diversity, we still had some unexplained variance, as well as anthropogenic variables with poor explanatory power. Unexplained variation could be due to unmeasured variables such as biotic interactions and behavior (D' Ambrosio et al., 2009), as well as natural and anthropogenic variables not included in our study such as local hydrologic variation (Marion et al., 2015). We may have seen a lack of an anthropogenic signal in some watershed metacommunities due to the branching complexity within river networks and the environmental heterogeneity provided by geomorphic variables offering a natural defense against human induced environmental changes (Terui et al., 2021). Unexplained variation could be due to homogenization from fish introductions, resulting in disruptions to 'natural' riverine processes structuring species traits and relatedness (Peoples et al., 2020). However, delineating metacommunities by geomorphic variables may allow for the anthropogenic mechanisms driving community structure to become more apparent.
In conclusion, we found that human modification of the land and riverscape is evident and important for structuring stream fish beta diversity. These effects operate both independently and in concert with the natural ecological variables that are commonly included in metacommunity analyses, as well as spatial gradients representing potential for species dispersal and homogenization of stream fish community structure across the landscape. Future studies should consider these anthropogenic effects alongside natural environmental effects and dispersal to obtain a clearer picture of regional beta diversity. We also found that user-defined a priori metacommunity delineations have important consequences for the outcomes and interpretation of metacommunity analyses. Analyses based on different spatial delineations highlighted the importance of synergistic spatial and environmental effects, regional differences in assembly rules, and basinspecific differences in the factors affecting different dimensions of stream fish beta diversity. Researchers should give careful consideration to this important methodological aspect, and seek spatial delineations that are ecologically relevant, given ecology of their focal taxa.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Ethics statement
Ethical review and approval was not required for the animal study because the data collection requiring IACUC approval was not completed in this study. The authors were given a dataset that on stream fish abundance.

Author contributions
LS: conceptualization, formal analysis, writing-original draft, and visualization. LB: resources and writing-review and editing. MS: writing-review and editing and resources. BP: supervision and writing-original draft. All authors contributed to the article and approved the submitted version.

Funding
This work was supported in part by the NIFA/USDA, under project number SC-1700599, and represents technical contribution number 7136 of the Clemson Experiment Station. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement of the U.S. Government.
Frontiers in Ecology and Evolution 14 frontiersin.org