Observed and projected functional reorganization of riverine fish assemblages from global change

Climate and land‐use/land‐cover change (“global change”) are restructuring biodiversity, globally. Broadly, environmental conditions are expected to become warmer, potentially drier (particularly in arid regions), and more anthropogenically developed in the future, with spatiotemporally complex effects on ecological communities. We used functional traits to inform Chesapeake Bay Watershed fish responses to future climate and land‐use scenarios (2030, 2060, and 2090). We modeled the future habitat suitability of focal species representative of key trait axes (substrate, flow, temperature, reproduction, and trophic) and used functional and phylogenetic metrics to assess variable assemblage responses across physiographic regions and habitat sizes (headwaters through large rivers). Our focal species analysis projected future habitat suitability gains for carnivorous species with preferences for warm water, pool habitats, and fine or vegetated substrates. At the assemblage level, models projected decreasing habitat suitability for cold‐water, rheophilic, and lithophilic individuals but increasing suitability for carnivores in the future across all regions. Projected responses of functional and phylogenetic diversity and redundancy differed among regions. Lowland regions were projected to become less functionally and phylogenetically diverse and more redundant while upland regions (and smaller habitat sizes) were projected to become more diverse and less redundant. Next, we assessed how these model‐projected assemblage changes 2005–2030 related to observed time‐series trends (1999–2016). Halfway through the initial projecting period (2005–2030), we found observed trends broadly followed modeled patterns of increasing proportions of carnivorous and lithophilic individuals in lowland regions but showed opposing patterns for functional and phylogenetic metrics. Leveraging observed and predicted analyses simultaneously helps elucidate the instances and causes of discrepancies between model predictions and ongoing observed changes. Collectively, results highlight the complexity of global change impacts across broad landscapes that likely relate to differences in assemblages' intrinsic sensitivities and external exposure to stressors.

Global change may affect ecological communities differently across broad gradients because of spatially variable degrees of exposure to environmental stressors. Currently, lowland areas are, in many cases, warmer and more anthropogenically modified, relative to upland areas. However, with ongoing global change, upland regions may become warmer and more modified and may, thus more resemble contemporary lowland environments (Osland et al., 2022;Trew & Maclean, 2021). Importantly, many of these upland regions are characterized by rare, specialized, and functionally unique taxa (Cadena et al., 2011;Myers et al., 2000) that may be disproportionately sensitive to global change. Therefore, assessing which regions or habitats may be most vulnerable to global changes will require understanding of varying degrees of stressor exposure and intrinsic sensitivity across broad landscapes.
In this study, we assessed ongoing and future fish biodiversity responses to climate and LULC in the Chesapeake Bay Watershed (CBW). First, we modeled effects of projected future (2030, 2060, and 2090) conditions ( Figure 1a) on riverine fish habitat suitability and assemblage trait composition. We used this information to address the following questions: (Q1) what functional traits are most sensitive to climate and LULC change? and (Q2) relative to a contemporary (2005) baseline, how will assemblage responses to global changes vary among regions and habitats? For Q1 we expected that a suite of traits including preferences for cold water temperatures, faster lotic flows, silt-free spawning substrates, and lower trophic status would predispose species to greater sensitivity to climate and LULC change and, thus, range contraction (i.e., global change "losers"; Dornelas et al., 2019;Rosset & Oertli, 2011). For Q2, we hypothesized that the magnitude of future assemblage responses would be greatest in upland regions and smaller (headwater) streams because these habitats have been less anthropogenically modified (Moritz et al., 2008) and may support distinctive ecological assemblages, relative to lowland and larger rivers (Meyer et al., 2007).
Finally, we compared future projections with observed timeseries data to ask: (Q3) how do observed environmental and fish assemblage trends (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) relate to model-projected changes over time (2005-2030; Figure 1a)? To do so, we used a subset of sites with long-term sampling records to estimate contemporary trends in environmental conditions and stream fish assemblage responses. We then assessed how time-series trends qualitatively related to modeled patterns approximately halfway through our initial modeling period. This analysis provided a novel perspective on whether climate and land-use model forecasts and subsequent projected effects on ecological assemblages are consistent with ongoing, observed trends through time, as well as through space. This analysis can also indicate specific instances where projecting contemporary trait-environment relationships into future conditions may be improved as future habitat conditions change.

| Study area and spatial framework
The CBW comprises five major river drainages: The Susquehanna, Potomac, Rappahannock, York, and James rivers (Figure 1b; Chesapeake Bay Program, 2022). Collectively, the >83,500 creeks, streams, and rivers in the watershed drain 186,600 km 2 of terrestrial land ( Figure 1b) These land-use changes may be compounded by climatic changes, especially warmer temperatures across the watershed. Clearly, the CBW characterizes a broad diversity of environmental features and potential threats to temperate riverine systems likely making our results widely generalizable in the context of global change.
We used the National Hydrography Dataset Version 2 Plus 1:100,000 resolution (NHD; McKay et al., 2012) as a spatial riverine base layer. The finest resolution in this dataset was confluence-toconfluence stream segments (hereafter, "reaches"). "Catchments" are defined as the terrestrial area draining directly into a specific reach and "watersheds" as the accumulated upstream terrestrial area draining into a reach inclusive of its associated catchment (Hill et al., 2016). We removed catchments classified as lakes, ponds, or reservoirs from our analysis. The total sample size for our projected analysis was 74,490 catchments, following data filtering for complete cases across predictors (described below). Each reach was assigned an Ecoregion (hereafter "region") the quantity of fish sampling locations and characteristics of fish assemblages in these regions that facilitated disaggregation. We distinguished these regions into two broader categories: lowland (CPL and PIED) and upland (SAP and NAP) and classified reaches into habitat sizes based on upstream drainage area as: headwaters, streams, and small, medium, and large rivers (see . These region and habitat size designations were subsequently used as spatial units within which reach-scale model predictions were aggregated to examine variable responses of fish assemblages to global change across the CBW. Throughout the data collection and analysis methods, we use "predictions" for modeled fish responses within contemporary periods, "projections" for modeled future fish responses, and "forecasts" for future scenarios of climate and LULC.

| Predictor variables
We used 2005 as a contemporary baseline period based on available climate and LULC datasets (see below). Predictor variables were summarized at both catchment and watershed scales (Table S1). Despite high correlation between catchment and watershed metrics for many variables (e.g., climate and LULC), the two scales provide independent information in many contexts. For example, a reach can have low values of catchment agriculture but be located downstream of intensive agriculture and, thus, may experience stressors from upstream accumulated agricultural runoff (Allan, 2004). Similarly, while accumulated upstream precipitation can indicate streamflow discharge, local catchment precipitation can reflect the amount of water entering reaches through runoff during rainfall events, particularly in urbanized catchments where much of local precipitation enters into streamflow (Machado et al., 2022). Therefore, we retained both catchment and watershed summaries to account for the many ways in which terrestrial metrics affect local conditions in directional, dendritic systems and selected a modeling approach (random forest models, explained in detail in subsequent sections) equipped to handle multicollinearity (Li et al., 2016).
Landscape and physical habitat variables were downloaded from a stream metrics database (StreamCat; Hill et al., 2016) and the NHD. These metrics characterized topography (elevation, slope), soils (permeability, erodibility, organic matter), baseflows, and anthropogenic modifications (road crossings, dam storage; Table S1). We lacked a basis for forecasting quantitative changes in these variables and so held them constant between contemporary and future periods.
We used climate change forecasts from Hay and McCabe (2019) that included seasonal (based on water years; January-March, April-June, July-September, October-December) temperature (°C) and precipitation (mm; Table S1)   We used a dataset of conterminous US LULC (Sohl et al., 2018) that included 12 anthropogenic and natural LULC classes (Table S1) with 2005 as a contemporary baseline (Krause & Maloney, 2019). and approaches to development and conservation (Nakicenovic et al., 2000). Based on these features, we selected two SRES that represented high (A2; high 2100 population, slow-paced advancement, and economic growth focus) versus low (B2; low 2100 population, medium-paced advancement, and environmental sustainability focus) US development ( Figure S1c; Nakicenovic et al., 2000).
We then designed 10 scenarios of future climate and LULC change based on varying combinations of low, medium, and high values of temperature, precipitation, and high versus low development (Table S2). Four scenarios used the 25th percentile of temperature (i.e., lowest warming): two were paired with the 25th percentile of precipitation (i.e., driest) with one for each of the land-use scenarios (A2 and B2), and two were paired with the 75th percentile of precipitation (i.e., wettest) with one for each of the land-use scenarios (A2 and B2; Table S2). We used two scenarios of moderate climate change based on the 50th percentiles of precipitation and temperature, one for each land-use scenario. The final four scenarios used the 75th percentile of temperature (i.e., greatest warming) where two were paired with the 25th percentile of precipitation (i.e., driest) with one for each of the land-use scenarios (A2 and B2), and two were paired with the 75th percentile of precipitation (i.e., wettest) with one for each of the land-use scenarios (A2 and B2; Table S2).

| Fish dataset
We used an extensive CBW fish sampling dataset that contained over 30,000 sampling events spanning 1969-2019, collected by 20 agencies and programs  Table S3). The objective of sampling varied among events from assemblage surveys that characterized the abundances of all species present to targeted surveys for a specific species (e.g., game fishes). To mitigate potential discrepancies caused by differing sampling objectives, we restricted our analysis to assemblage surveys. Each sampling event was assigned to an NHD reach based on geographic coordinates and location name from the data provider and rated on the confidence in NHD reach assignment. For example, coordinates that fell directly on, and shared a name with, an NHD reach were considered highly confident, whereas those that fell elsewhere in the catchment or were named for a different reach were assigned lower confidence (described in . For purposes of our analysis, we filtered our dataset to high confidence assignment to NHD reaches. Taxonomy was harmonized among samples to species-level identifications using a global database of fishes (FishBase; Froese & Pauly, 2022). To temporally align our survey dataset with baseline predictors, we restricted our fish dataset to 2000-2010 samples and selected the sample nearest to the 2005 baseline when multiple samples occurred from one reach (n = 2442 samples; Figure 1c). These samples were converted into a reach × species abundance matrix.
We used functional trait data (hereafter, "traits") from the Fish Traits (Frimpong & Angermeier, 2010) and National Rivers and Streams Assessment (U.S. Environmental Protection Agency, 2020) to evaluate species' relative sensitives to global change stressors.
These traits included streamflow, substrate, and temperature preferences, and reproductive and trophic guilds ( Table 1) and were used to build a species × trait matrix. This matrix was used in two ways: (1) to assign species to functional groups (FGs) which were, subsequently, used to identify suites of traits that may predispose species to increasing or decreasing habitat suitability with ongoing global change ("winners" and "losers," respectively; Dornelas et al., 2019;Rosset & Oertli, 2011) and (2) incorporated with the reach × species abundance matrix to derive abundance-based functional metrics that characterized assemblage trait composition across all species. For both aims, it was necessary to binary dummy-code all traits. We also assigned weights to each dummy-coded column that summed to one within each categorical trait. For example, the trophic guild trait consisted of four categories (carnivore, herbivore, invertivore, and omnivore), so each dummy-coded column was assigned a weight of 0.25.

| Q1: What traits are associated with sensitivity to global change?
We assigned species into FGs and assumed that species within the same FG would share more similar traits and, thus, be more likely to respond similarly to global change stressors. Although there is likely considerable within-FG variation in species responses to global change, the value of this approach is that it can allow some generalization from predictions made from well-sampled species to less sampled species (Hendrick & McGarvey, 2019;McGarvey et al., 2021;Woods & McGarvey, 2022).
We used dendrogram clustering to assign species into FGs (using the "FD" package in R; Laliberté et al., 2014;Laliberté & Legendre, 2010). This method used the binary dummy-coded species × trait matrix and associated weights to build a species × species Gower dissimilarity matrix (Gower, 1971). A Ward's hierarchical clustering algorithm (Murtagh & Legendre, 2014) was applied to the dissimilarity matrix to produce a functional dendrogram. We visually inspected this dendrogram and determined a pruning height which assigned each species to one of five FGs ( Figure S2). FG assignment depends on the clustering method applied. Therefore, we assessed sensitivity of FGs to an alternative method and performed k means clustering on the first three axes from a principal coordinates analysis (PCoA) applied to the Gower dissimilarity matrix. These two methods confirmed that five FGs were appropriate to characterize the range of multivariate trait space and resulted in similar groupings of species to FGs, which indicated our assignments were robust to clustering methodology ( Figure S3).
Next, we selected one species per FG as a focal species for our analysis based on sample size considerations. Specifically, we calculated the number of presences (i.e., reaches with observed non-zero abundances) per species and selected the most widely distributed species per FG as a focal species for analysis, which ensured a minimum sample size of 350 reaches per species (Table S4). Although these species likely characterized one extreme of multivariate trait space within FGs, this effect would be consistent among FGs and, therefore, not hinder our interpretation of variation in global change sensitivities among FGs.
For each of our five focal species, we modeled future reach-scale habitat suitability. We did not consider range limitation in this analysis (e.g., restriction of modeled estimates within a species range) because the objective of this analysis was to infer the suitable habitat for FGs to assess which functional strategies might be more broadly TA B L E 1 Functional traits used to define functional groupings. sensitive or tolerant to future global changes. However, the inclusion of region and upstream accumulated area as covariates in our models (Table S1) should have accounted for FG requirements related to regional species pools and stream size, respectively. were tuned using default settings (Maloney, Buchanan, et al., 2022;Probst et al., 2018).
Using these reach-scale habitat suitability predictions, we summarized focal species responses to global change in two ways. First, we measured the proportion of CBW reaches assigned as low, medium, or high habitat suitability per FG and modeling period (2005, 2030, 2060, and 2090). Second, we used the 2005 and 2090 predictions and classified total change in suitability at the reach scale (e.g., low suitability to high suitability). We assumed that suites of traits associated with focal species that were projected to gain suitable habitat would indicate FGs favorable to future climate and LULC conditions (i.e., global change "winners") relative to those with declining habitat availability (i.e., global change "losers").

| Q2: How will fish assemblage responses vary among regions and habitats?
We also used functional response metrics that described trait composition across all species in an assemblage (i.e., reach-specific sample) calculated by . These metrics were the percentage of individuals in each assemblage assigned to traits representing specific flow preferences (rheophilic; RHEOPIND), reproductive guilds (coarse substrate [lithophilic] spawners; LITHPIND), thermal preferences (cold-water; COLDPIND), and trophic guilds (carnivores; CARNPIND; Table 2).
We calculated abundance-based functional metrics describing the functional diversity (functional richness; FRIC), and functional redundancy (functional dispersion; FDIS) of fish assemblages (using the "FD" package  (Tucker et al., 2017). MPD quantifies the mean phylogenetic difference among all species pairs in the assemblage. Decreasing MPD indicates species are, on average, more closely related and, thus the assemblage is more phylogenetically redundant (Tucker et al., 2017). We constructed a tree describing phylogenetic relationships between species pairs in our dataset (using the "fishtree" (Chang et al., 2019;Rabosky et al., 2018) package), which was used to calculate the PD and MPD of each assemblage (with the "picante" (Kembel et al., 2010) and "ape" (Paradis & Schliep, 2019) R packages).
We used RF regression models to predict assemblage functional and phylogenetic metrics at all reaches under contemporary and future periods with the "tuneRanger" package. We used an 80-20 training-testing (n = 1953 and 489) split to evaluate predictive performance of RFs with root mean residual squared errors and percentage of variation explained (R 2 ) in training and testing datasets (Table S5). We used the same number of trees and methods for tuning hyperparameters described for the focal species analysis above.
We used modeled reach-scale assemblage responses to assess variable global change impacts among regions and habitat sizes. At  For the time-series analysis, we identified a subset of wellsampled sites from our larger dataset. Specifically, we required that sites be sampled at least once within each of seven periods : 1999-2001, 2002-2004, 2005-2006, 2007-2008, 2009-2011, 2012-2013, and 2014-2016. These periods were selected so that the end year aligned with available LULC data in National Land Cover Dataset  (Figures 1d and 7a). SHEN NP sites lacked LU disturbance, whereas MD sites may be subjected to more intensive LU stress.
We used environmental data to assess observed climate and LULC trends at each site. We used annual monthly catchment temperature and precipitation from the PRISM dataset for summer months (July, August, and September) for each year 1999-2016. These annual metrics were summarized within the same periods as the fish dataset (described above) to derive mean catchment summer temperature and precipitation. These catchment values were paired to the fish dataset for the corresponding time periods. We also paired each observed fish sampling event with watershed values of agriculture and developed LU from the corresponding NLCD release (2001, 2004, 2006, 2008, 2011, 2013, and 2016) for each period.
Next, we quantified observed site-level time-series trends in (1) climate and LU predictors and (2) fish assemblage responses. For assemblage responses, we focused on CARNPIND, LITHPIND, FD, and PD as responses because these variables were representative of biodiversity in SAP and PIED regions and had relatively high explanatory power in RF models (Table S5). For each environmental predictor (n = 4) and assemblage response (n = 4), we ran site-level linear regression models (LMs) with year as a predictor variable (n = 8 LMs for each site). We extracted the slope from each LM to describe the site-level trend in each environmental predictor and assemblage response. We used these slopes to perform a qualitative assessment Forecasted developed LU increases were less extensive in NAP and SAP in both A2 and B2 scenarios.

| Q1: What traits are associated with sensitivity to global change?
The five FGs differed by habitat preferences and resource use ( Figure 3). FG1 (represented by focal species: Lepomis macrochirus) characterized carnivorous, nest-guarding taxa with preferences for vegetated substrate, warm-water, and pool habitats ( Figure 3).
FGs differed in sensitivity to changes in habitat suitability ( Figure 4). Generally, FGs 1 and 4 gained highly suitable habitat over  Table S7). These taxa could be predicted to be "winners" under global change scenarios of warmer temperatures, decreased summer precipitation, and increasing development. In contrast, highly suitable habitat for FGs 2, 3 and 5 decreased over time (Figure 4b,c,e; Table S7). Taxa in these FGs could be considered global change "losers."

| Q2a: How will fish assemblage responses vary among regions?
Projected assemblage responses to global changes varied among regions ( Figure 5). COLPIND decreased in Appalachian regions (SAP and NAP; Figure 5a) but showed little response in more lowland  regions (CPL and PIED; Figure 5a; Table S8). This result is not surprising, given the low percentages of cold-water strategists in lowland regions in contemporary periods. Models projected decreased LITHPIND and RHEOPIND in PIED, SAP, and NAP (Figure 5b,c; Table S8). By 2090, projected LITHPIND losses were greatest in PIED, followed by SAP and NAP, whereas RHEOPIND losses were greatest in SAP, followed by PIED and NAP. Contrary to the above traits, models projected increasing CARNPIND over time (Figure 5d; Table S8). CARNPIND increases were greatest in PIED, followed by CPL, SAP, and NAP.

F I G U R E 3 Functional group assignments for
Projected changes in FRIC and FDIS differed between lowland and upland regions (Figure 5e,f). Models projected increasing FRIC (i.e., greater functional diversity) in Appalachian regions (NAP and SAP) but decreasing FRIC in CPL and PIED (Figure 5e; Table S8).
Models projected decreases in FDIS (i.e., increasing functional redundancy) in CPL and PIED but minimal change in FDIS in SAP and NAP ( Figure 5f; Table S8). Like FRIC, models projected increased PD (i.e., greater PD) in upland regions (SAP and NAP) and decreased PD in lowland regions (PIED and CPL) ( Figure 5g; Table S8). Models projected increased MPD (i.e., less phylogenetic redundancy) in PIED, SAP, and NAP but decreased PD (i.e., more phylogenetically redundant) in CPL (Figure 5h; Table S8).

| Q2b: How will fish assemblage responses vary among habitats?
Model-projected assemblage responses also differed among habitats within regions ( Figure 6). COLDPIND decreases in SAP and NAP were greatest in headwaters relative to large river habitats ( Figure 6a). Projected changes in LITHPIND and RHEOPIND by habitat size varied among regions (Figure 6b), increasing in large rivers relative to smaller habitats in CPL, but decreasing in headwaters relative to larger rivers in PIED and SAP (Figure 6b,c; Table S9). Within NAP, medium-sized habitats (stream-medium rivers) showed larger decreases in LITHPIND, relative to headwaters and large rivers ( Figure 6b; Table S9), whereas medium and large habitats showed greater decreases in RHEOPIND relative to smaller habitats (headwaters-small rivers; Figure 6c; Table S9).
Model-projected increases in CARNPIND varied by habitat size, and the relationship differed among lowland and upland regions ( Figure 6d). In lowland regions (CPL and PIED), larger increases in CARNPIND were projected in smaller habitats (headwaterssmall rivers), relative to larger habitats (medium-large rivers) ( Figure 6d; Table S9). In contrast, in upland regions, larger habitats were projected to gain more CARNPIND, relative to smaller habitats ( Figure 6d).  Table S9). In SAP, models generally projected increased FRIC, FDIS, PD, and MPD (i.e., increasing diversity and lower redundancy) in smaller habitats (headwatersstreams), whereas larger habitats (medium-large rivers) were projected to become less diverse and more redundant (Figure 6e-h; Table S9). In NAP, models projected the greatest increases in FRIC, FDIS, PD, and MPD in medium-sized habitats (streams-medium rivers) (Figure 6e-h; Table S9).

| DISCUSS ION
Global change effects on biodiversity are complex. Species vary in sensitivity to stressors based on functional traits that confer the environmental conditions required for species persistence and, spatially, habitats differ in their degree of exposure to stressors. We integrated these species-and habitat-specific aspects of vulnerability and projected global change impacts on CBW riverine fishes.
We found that suites of traits including preferences for warm water, pool habitats, and fine or vegetated substrates, nest-guarding reproductive strategies, and carnivores may gain suitable habitat in the future. We also found that assemblage responses to global change varied among regions and habitats because of differing stressor exposure. We also assessed how well model-projected biodiversity responses agreed with ongoing time-series trends halfway through the initial future modeling period. Although seldom applied in global change studies, this approach allowed us to better understand circumstances where observed trends differed from modeled patterns and potential causes of these discrepancies.

| How do fish functional traits relate to climate and LULC change?
Commonly, analyses project global change impacts from speciesspecific perspectives (e.g., species distribution models). Although useful, taxonomic approaches overlook functional aspects that can cause species to differ in their sensitivity to global change stressors. Such information is useful to generalize inferences of global change impacts to less-sampled taxa or regions. We addressed this challenge and used functional trait information to compare differences in projected habitat suitability, which were then compared among species to elucidate suites of traits that may be particularly sensitive to global change stressors. Generally, carnivorous species that preferred warmer water, finer substrates, and pool habitats and exhibited higher parental care (guarder reproductive strategists) were projected to gain suitable habitat with ongoing global change. Species that exhibited this "winning" functional strategy included many species of the Centrarchidae and Esocidae families (Table S6). Increased habitat suitability for these generalist predators could alter biotic interactions and thus, ecosystem functioning. For example, expansion of native or nonnative centrarchids (e.g., basses) may locally extirpate native prey (e.g., cyprinids) and predators (e.g., walleye) due to increased predation or competition, respectively (Alofs & Jackson, 2015;Hansen et al., 2017;Hughes & Herlihy, 2012;Jackson & Mandrak, 2002).
In contrast, taxa that may be disproportionately sensitive to global change preferred cooler temperatures, coarse, clean spawning substrates, and faster flow velocities. These findings are well supported by literature as traits that predispose freshwater taxa to negative global change impacts in many instances (Comte et al., 2013;Heino et al., 2009;Lynch et al., 2022;Rosset & Oertli, 2011). These traits also characterize many functionally unique montane fishes (Colvin et al., 2019;Freeman et al., 2007), prominently, salmonids that provide many ecosystem services (Watz et al., 2022) and are particularly at risk of global change-induced range contractions (Chu et al., 2005;Lynch et al., 2022;Sharma et al., 2007).
Although these "winning" and "losing" trait assignments are broadly confirmatory with prior traits-based assessments of global change sensitivity, our analysis approach provides benefits over many prior species-specific assessments. Integration of functional traits allowed us to form and assess hypotheses of species sensitivities while accounting for multidimensional trait composition. This approach is advantageous compared to analyses that fail to incorporate functional traits or rely on post-hoc descriptions of species sensitivities with respect to a single trait (e.g., thermal sensitivity).
Furthermore, a priori assignment of species to FGs (Table S6) can provide managers with information on potential sensitives of species that lack sufficient data for modeling efforts.

| How do stream fish assemblage responses to climate and LULC vary among regions and habitats?
Although species sensitivities to global change stressors can provide valuable information on global change impacts, a less frequently examined aspect is how ecological assemblage vulnerability may vary across space owing to differing degrees of exposure to environmental stressors. We addressed this aim and used trait-environment relationships to project future suitability for key assemblage functional and phylogenetic metrics following exposure to climate and land-use change.
Our models projected broad declines in habitat suitability for cold-water, rheophilic, and lithophilic species and increases for carnivorous species throughout the CBW. However, the direction and magnitude of these assemblage responses varied among regions and habitats. The Southern Appalachian region and smaller habitats (headwaters) within this region were projected to experience the largest decreases in suitable habitat for cold-water and rheophilic individuals. Two, non-exclusive mechanisms may underlie these projections. First, future stressor exposure from climate change (i.e., warming temperatures) may be greatest in upland regions and smaller habitat sizes (Moritz et al., 2008). Second, taxa inhabiting upland regions and smaller habitats characterized by these traits may be disproportionately sensitive to global changes owing to intrinsic characteristics (e.g., low critical thermal maxima ;Frishkoff et al., 2015;Nowakowski et al., 2017;Sunday et al., 2014). In contrast, projected habitat suitability increases for carnivores and decreases for lithophiles were greatest in the Piedmont region and smaller habitat sizes within these regions. One habitat feature that may drive this pattern is forecasted intensive urban development in lowland regions ( Figure 2) that can facilitate the upslope expansion of generalist predators (Riley et al., 2005;Santana Marques et al., 2020) and reduce lithophilic habitat via increased sedimentation (Sutherland et al., 2002). In addition, expansion of agriculture (particularly irrigated agriculture) in lowland regions can compound urbanization effects, further modify stream physical habitat and degrade ecological condition (Herlihy et al., 2020;Kaufmann, Hughes, Paulsen, Peck, Seeliger, et al., 2022a, 2022b. When combined with the above information on species and trait sensitivities, functional and phylogenetic metrics can be used to understand how trait-based environmental filtering may differentially affect assemblage composition between upland and lowland regions. Lowland regions (and medium to large habitat sizes within these regions) were projected to become less diverse and increasingly redundant (i.e., more homogenized). In contrast, upland regions (and smaller habitat sizes within these regions) were projected to become more diverse and less redundant, despite model-projected declines in cold-water, rheophilic, and lithophilic species. These results suggest different processes may underlie ecological assemblages' responses to global changes among regions. Lowland fish assemblages may become more homogenized but changing environments in upland regions may facilitate the expansion of lowland generalists, thus increasing assemblage diversity (Buisson et al., 2008;Scott, 2006;Scott & Helfman, 2001).
However, caution is warranted in our interpretation of lowland homogenization owing to relatively low model performance for redundancy metrics (Table S5), differing patterns between observed trends and projected changes (see Section 4.3), and limited support for taxonomic homogenization despite generalist expansion facilitated by land-use change in Southern Appalachian fish assemblages (Petersen et al., 2021).
Our results highlight the variable vulnerabilities of ecological assemblages to global change among regions and habitats. For example, upland regions may be disproportionately sensitive with respect to some assemblage responses (e.g., losses of cold-water rheophiles).
These results may be related to differences in forecasted stressor exposures of greatest warming in upland regions compared to lowland regions which may experience greater increases in developed landuse ( Figure 2). Several studies support the idea that upland regions are disproportionately vulnerable to global change impacts (Moritz et al., 2008;Osland et al., 2022), but our results support a more nu-

| How do observed trends relate to modeled patterns?
Our time-series analysis indicated that CBW riverine environments and fish assemblages are, indeed, changing. Both Southern Appalachian and Piedmont sites showed trends toward warmer and drier summer conditions over a 16-year period. However, Piedmont sites additionally experienced increasingly developed land-use and, subsequently, fish assemblages at these sites showed more extensive changes compared to Southern Appalachia. These results may indicate that ecological assemblages exposed to multiple stressors are more affected compared to those exposed to fewer stressors (Brodie, 2016;Comte et al., 2022;Outhwaite et al., 2022). Alternatively, these results could indicate that landuse change exerts greater pressures on fish assemblages relative to climate changes over the period examined (Brodie, 2016;Leprieur et al., 2008;Outhwaite et al., 2022;Tedesco et al., 2013).

| CON CLUS ION
Our CBW analysis has implications for stream fish biodiversity globally. For example, Comte et al. (2022) analyzed stream fish abundance time series and found that assemblages shifted functional composition toward greater representation of warm-and slowwater individuals through time in response to climate and land-use changes, globally. Our findings from a more limited spatial extent (but greater sampling density) are consistent with these findings and additionally project that these trajectories may continue in the future. Furthermore, our work highlights the importance of additional niche axes including substrate preferences, reproduction, and trophic level in understanding stream fish vulnerability to global change.
However, there are several caveats to our findings. Although we used a broad suite of climate, LULC, and landscape features as predictors, we were unable to include instream variables including streamflow and water temperature. Our coarser air temperature metrics are likely suitable surrogates for water temperature because they provide similar information in large-scale assessments (McGarvey et al., 2018). But, given the potentially strong role of streamflow in structuring fish assemblage composition (Freeman et al., 2022), future studies could make use of mechanistic hydrologic models to forecast climate and LULC impacts on streamflow and, subsequently stream biota (e.g., Andres et al., 2019;Cheng et al., 2022;. Although our study assessed potential changes in average climate conditions, global change may also increase the occurrence of extreme weather events (e.g., floods, heat waves; Alifu et al., 2022;Tassone et al., 2022), which may require forecasts that account for changing frequencies of anomalous discharge and water temperature events. Likewise, potential upland and upstream expansion of generalist predators may require the presence of instream habitat feature such as deep pools (Harvey & Stewart, 1991) and hydrologic connectivity for dispersal (Radinger & Wolter, 2015). Therefore, more detailed information on stream channel morphology or river fragmentation could indicate our projections of predatory individual expansion are unrealistic. Habitat fragmentation can also limit the native species dispersal and models that incorporate detailed information on instream barriers to fish migration could yield different insights, particularly with respect to the accessibility of potential upland refugia (Ebersole et al., 2020). Also, there are caveats based on model performance and our qualitative trend assessment has limits. For example, the percentage of variation explained by our RF models was moderate (Table S5) and coupled with uncertainty in climate and LULC forecasts, indicates caution in applying model predictions. With increasing quantity and quality of biological and environmental data and advancements in powerful, data-intensive statistical methods, we expect improvement in projecting climate and LULC impacts on natural systems.
Our work identified lines of inquiry for future studies in the CBW and beyond. One key direction for global change research is to investigate the implications of functional reorganization for biotic interactions. Predator-prey interactions are strong structural mechanisms in ecological systems and the extinction or introduction of top predators can initiate trophic cascades (Natsukawa & Sergio, 2022). Increasing habitat suitability for carnivores (and, potentially, associated predation pressure) can propagate within stream food webs via linkages to prey and competitor organisms (Jackson et al., 2020) and to terrestrial systems by altered aquatic-terrestrial linkages (Baxter et al., 2005). Our assemblagescale functional assessment provides an additional perspective on global change impacts compared to species-specific approaches.
However, understanding of the scale and scope of anthropogenic alterations on biological assemblages will likely require knowledge of how changes in one ecosystem component affect the associated web of interconnected components among regions and globally (Romero et al., 2018).

ACK N OWLED G M ENTS
We thank members of U.S. Geological Survey's CBW Regional Assessment team, Than Hitt, and James McKenna Jr. for insightful comments during the project. We also thank all the data providers whose data were the backbone of this study. Support was provided by the Fisheries Program of the U.S. Geological Survey's Ecosystems Mission Area and Priority Ecosystems Science (Chesapeake Bay study unit) programs. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare that there is no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The fish assemblage metrics , fish sampling data , and climate and land-use/ land-cover summaries (Krause & Maloney, 2019)