Working Title: Understanding landscape influences on aquatic fauna across the central and southern Appalachians

: Understanding influences of multiple stressors across the landscape on aquatic biota is important for conservation, as it allows for an understanding of spatial patterns and informs stakeholders of significant conservation value. Data exists for land use/landcover (LULC) and other physicochemical components of the landscape throughout the Appalachian region yet biological data is sparse. This dearth of biological data relative to LULC and physicochemical data creates difficulties in making informed management and conservation decisions across large landscapes. At the HUC12 watershed scale we sought to create a single score for both abiotic and biotic values throughout the central and southern Appalachian region. We used boosted regression trees (BRT) to model biological responses (fish and aquatic macroinvertebrate variables) to abiotic variables. Variance explained by BRT models ranged from 62-94%. We categorized both predictor and response variables into themes and targets respectively to better understand large scale patterns on the landscape that influence biological condition of streams. We combined predicted values for a suite of response variables from BRT models to create a single watershed score for aquatic macroinvertebrates and fish. Regional models were developed for fish but we were unable to develop regional models for aquatic macroinvertebrates due to the low number of sample sites. There was strong correlation between regional and global watershed scores for fish models but not between fish and aquatic macroinvertebrate models. Use of such multimetric scores can inform managers, NGOs, and private land owners regarding land use practices; thereby contributing to largescale landscape scale conservation efforts.

Abstract: Understanding influences of multiple stressors across the landscape on aquatic biota is important for conservation, as it allows for an understanding of spatial patterns and informs stakeholders of significant conservation value. Data exists for land use/landcover (LULC) and other physicochemical components of the landscape throughout the Appalachian region yet biological data is sparse. This dearth of biological data relative to LULC and physicochemical data creates difficulties in making informed management and conservation decisions across large landscapes. At the HUC12 watershed scale we sought to create a single score for both abiotic and biotic values throughout the central and southern Appalachian region. We used boosted regression trees (BRT) to model biological responses (fish and aquatic macroinvertebrate variables) to abiotic variables. Variance explained by BRT models ranged from 62-94%. We categorized both predictor and response variables into themes and targets respectively to better understand large scale patterns on the landscape that influence biological condition of streams.
We combined predicted values for a suite of response variables from BRT models to create a single watershed score for aquatic macroinvertebrates and fish. Regional models were developed for fish but we were unable to develop regional models for aquatic macroinvertebrates due to the low number of sample sites. There was strong correlation between regional and global watershed scores for fish models but not between fish and aquatic macroinvertebrate models.
Use of such multimetric scores can inform managers, NGOs, and private land owners regarding land use practices; thereby contributing to largescale landscape scale conservation efforts.

<A>Introduction
Assessing biological condition of rivers and streams has been performed at various spatial scales and grains. As data availability and computing power both increase, the ability to perform spatial analysis across large geographic extents and at small grain size is increasing. Examples of large scale spatial analyses across the United States include; Esselman et al. (2013), who used fish as indicators, the US Environmental Protection Agency's Wadeabe Stream Assessment (2006) which included benthic macroinvertebrates and fish as indicators, and Hill et al. (2017) (Harris and Silveira, 1999;Hering et al., 2004;Pont et al., 2006;Pont et al., 2009). Multimetric indices use taxonomic and functional metrics that are known to be sensitive to gradients of anthropogenic disturbance as a means of assessing the health of aquatic systems.
Landscape-scale analyses are important as they cross jurisdictional boundaries and make efforts to more fully understand anthropogenic stressors on biological and ecological needs of species, populations, and communities across broad spatial extents which are often not represented in more regional analyses. We define landscape scale as the spatio-temporal extent that allows for multigenerational population processes to occur. Such large scale analyses require extensive data that is representative across taxonomic groups through both space and time and physicochemical conditions must be represented that range from natural to extensive anthropogenic disturbance (Stoddard et al. 20008). Difficulties exist to ensure data sets are consistent across regions, representative taxa (response variables), and physicochemical factors (predictor variables). It is often difficult to ensure such consistency due to biological and environmental datasets being expensive to obtain, varying goals and objectives from entities collecting data, and varying monetary support for such collections. This often results in either poor spatial coverage and spatially clumped data or both. To overcome these difficulties data sets are often combined from various collection programs. However, this can result in data existing at varying taxonomic levels, problems with standardization of sampling protocols, and biases associated with site selection. For example, site densities may be more dense in some areas than others resulting in biases towards trends that are prevalent in heavily sampled areas. Angermeier and Smogor (1995) found that sampling effort influenced accuracy and precision of estimates of biological community attributes.
Anthropogenic induced stressors often correlate strongly with biological conditions of aquatic systems and thus landscape-scale analyses typically aim to use land cover and land use practices over large geographies as surrogates for local environmental conditions to provide insight into biological conditions of streams and rivers. Use of such analyses allow broad geographic scale understanding of the existing gradients of poor-to-good biological condition and the environmental factors that affect them. This aids in identifying landscape-scale stressors that further the ability of conservation practitioners, managers, and regulators to take actions which positively impact aquatic fauna through a greater understanding of large-scale stressors.
However, over such great spatial extents the biological responses to anthropogenic stressors can be regionally variable and the result of local biogeographical interactions with natural and anthropogenic stressors over expansive time periods (Esselman et al. 2013).
Indices for evaluation of biological condition across large geographies have developed standard procedures enabling transferability across regions through the use of statistically robust methods (Bramblett et al. 2005;Esselman et al. 2013;McCormick et al. 2001;Stoddard et al. 2008). Biological condition indices should sufficiently cover the range of environmental gradients, exhibit low temporal variability to repeated visits to the same site(s), and demonstrate high responsiveness to environmental degradation. Indicators that do not meet these requirements should be excluded from analyses (Stoddard et al. 2008;Whittier et al., 2007). Use of the screening methods provided by Whittier et al. (2007)  Cooperative (App LCC) region (see methods for full description) using methods developed by Stoddard et al. (2008) andWhittier et al. (2007). Additionally, we aimed to identify large scale themes that structure biological communities and thus biological condition across this large geography. The App LCC encompasses biologically important aquatic ecosystems (e.g., Tennessee and Mobile Rivers) and is predicted to undergo continued land conversion, resulting particularly in increased urbanization (Terando et al. 2014). Using broad scale land use information to develop models to predict biological condition of streams and watersheds provides the opportunity to predictively map spatially explicit biological condition across large geographic regions (Hill et al., 2017;Carlisle et al. 2009). An understanding of the anthropogenic threats that exist in the Appalachians is an important endeavor, as more than half of the U.S. population lives within a day's drive of some part of the Appalachians (Ford 2015).
Thus, modelling and mapping the current biological condition and developing an understanding of threats to conservation across the region is an important step towards identifying conservation opportunities and needs so that conservation practitioners can act accordingly.

<B>Methods
Our study focused on the rivers, streams, and landscape encompassed by the Appalachian Landscape Conservation Cooperative (App LCC). The App LCC is comprised of 592,129 km 2 , 266,307 HUC 12 watersheds, and intersects 15 states ranging from central New York to the north, and central Georgia, Alabama, and Mississippi to the south. The eastern most extent reaches western New Jersey in the north and the western most portions of North and South Carolina, while reaching as far west as southern Illinois. The topography and habitat of the region are complex over such a large land area which is comprised of six Level II ecoregions (U.S. EPA), 11 Freshwater Ecoregions, and seven major river basins and has been greatly influenced by the lack of glaciation and ancient geology which gave rise to high levels of aquatic biodiversity. The App LCC region encompasses some of the most biologically diverse freshwater fish resources in the United States. This is particularly true in the southwestern portion of the region. As the 19 th and 20 th centuries have progressed, aquatic diversity has been impacted by anthropogenic activities ranging from agricultural practices, land conversion, and extractive industries.

Landscape data used as predictors
Through the aid of a steering committee comprised of 60 regional experts that include academic scientists, state wildlife agency managers, non-governmental organization scientists, and conservation practitioners we identified 52 predictor variables comprising six themes that are known to influence fish and aquatic macroinvertebrate communities (Appendix Table 1A).
Themes included streamflow, geomorphic condition, connectivity, water quality, non-point source pollution, and point source pollution and were established and evaluated to understand large scale thematic processes that influence aquatic biota across the App LCC. Predictor data were compiled from regional and national datasets at the National Hydrography Dataset version 2 (NHD+V2) HUC12 watershed level (Final predictors used in BRT models are shown in Table   1; Comprehensive list of predictor data with source provided in Appendix A). Riparian corridors are known to be influential for developing ecological diversity (Naiman et al., 1993) and in structuring aquatic communities; hence, we used The Nature Conservancy's Active River Area (Smith et al., 2008) as a buffer into which land cover variables were clipped for modeling. The active river area framework provides an ecologically sound modelling mechanism to capture important variables that structure aquatic ecosystems (Smith et al., 2008). We centered and scaled all predictor variables to overcome instances where predictor variables had wide ranges resulting single predictors dominating models.

Fish and macroinvertebrate data
Our fish and aquatic macroinvertebrate community data were assembled from federal databases of national datasets (United States Environmental Protection Agency's National Rivers and Streams Assessment and the National Fish Habitat Partnership). We used these community datasets as predictors to assess the status of biological communities within HUC12 watersheds (Table 2). Sample point locations were spatially joined to HUC12 watersheds by projecting into a conic equal area projection in Quantum GIS (Q GIS Development Team, OSGF). Our fish community dataset was comprised of 2,991 sample sites across the region and was somewhat sparse (2,991 locations within 266,309 HUC12 watersheds = 1.1%) but exhibited good spatial distribution ( Figure 2). Our aquatic macroinvertebrate dataset was comprised of only 194 sampling locations and had poor spatial coverage (194 locations of 266,309 HUC12 watersheds = 0.07%) and distribution across the region (Figure 2).
For fish we calculated taxonomic richness and diversity and several functional group metrics using the fish traits database provided by Frimpong and Angermeir (2009) and the tolerance database of Stoddard et al. (2005). In total we calculated 45 metrics for further analysis of fish responses to our predictor dataset. Our final analysis used 9 of these metrics and these were further aggregated into three groupings (i.e., targets) ( Table 2). For aquatic macroinvertebrates the metrics modeled as response variables were % EPT (i.e., percent Ephemeroptera/Plecoptera/Tricoptera), % of total abundance comprised by the five most dominant taxa, % of total abundance that were tolerant taxa, and % of the total abundance that were intolerant taxa (Table 2).

Boosted regression tree models and multimetric indices
Prior to model development we successively screened all predictor variables (e.g., Stoddard et al. -2008-,Esselman et al. -2013.-. After which, we employed boosted regression trees (BRTs) to predict fish and aquatic macroinvertebrate response variables at the watershed scale (Table 2). Boosted regression trees are robust to missing data, have both explanatory and predictive power, and are capable of handling complex relationships that are both non-linear and interactive without the need for prior data transformation. By iteratively fitting and combining simple regression models BRT models improve model structure and predictive performance.
Tree complexity and learning rate were systematically altered in order to identify the optimal model structure. We reduced model complexity by removing redundant variables (r > 0.80) and variables with minimum variation among sites. Global models were developed with the remaining variables and then further simplified using scree plots of the predictor relative influence. Simplified models were retained if their cross-validation was greater than or equal to that of the global model. BRT models were developed with methods and code provided by Elith et al. (2008). Initial analyses were performed across the entire geography without accounting for regional differences (i.e., ecoregions) but due to a clear need for the inclusion of regional Lakes River Basins), and Tennessee (Tennessee and Mobile Bay River Basins). Regional models were successfully run for fish but due to the lack of spatial coverage, we were unable to model aquatic macroinvertebrates regionally and we therefore retained models for the entire geography for aquatic macroinvertebrates.
Predictions of response variables from BRT models were centered and scaled. As a multimetric approach the predicted responses (biological metrics in table 1) were then averaged within their target group (see table 1) and further condensed via averaging into overall fish and aquatic macroinvertebrate scores. Boosted regression trees provide the relative influence of predictor variables in the modelling of responses and we averaged these relative influences within our thematic framework (see above and table 1) for each fish and aquatic macroinvertebrate response variable.

<C>Results
The screening process reduced the total number of predictor variables to 37 from 52. When predictor variables were correlated, we attempted to retain the variable most responsive to our models; however, this was not formally tested. The themes with the greatest number of predictor variables were non-point source pollution and point source pollution (each with eight variables).
While the theme with the least number of predictor variables was geomorphic condition (two variables) (Table 1).
Below, we present results from the boosted regression tree models where the models were run with all ecoregions together for both fish and aquatic macroinvertebrates. We refer to this as a "Global" model (see below). After which, we present results where separate models were run for each Freshwater Ecoregion (defined in the materials and methods section). We refer to these models as regional models and because there were not enough data points within the regional data sets for aquatic macroinvertebrates, they were developed for fish only.

Global models
Boosted regression tree models were developed for 45 % (28/62) of the response variables; however, we only retained 12 (eight and four response variables for fish and aquatic macroinvertebrates respectively) for our final multimetric index of biotic condition (Table 2).
Response variables were collapsed into targets where for the fish category we retained three targets (Shannon diversity, functional group, and taxa quality) and for aquatic macroinvertebrates we retained macroinvertebrate taxa quality as the only target.
The range of deviance explained by the BRT models was 48 % for percentage of invertivore taxa to 85 % for fish taxa diversity (x̄ = 73 %) for fish and 76 % for percentage of the five dominant taxa to 89 % for percentage of tolerant taxa (x̄ = 81 %) for aquatic macroinvertebrates.
For the fish BRT models, the cross-validated deviance ranged from 32 % for percentage of invertivore taxa to 69 % for fish taxa diversity (x̄ = 52 %), while for aquatic macroinvertebrates it ranged from 24 % for percentage of the five dominant taxa to 76 % for percentage of tolerant taxa (x̄ = 41 %). The lower cross validated deviances for aquatic macroinvertebrates relative to those of fish, are likely due to the lack of spatial data coverage in the aquatic macroinvertebrate data (Figure 1 b).
Based on the results of the BRT models, the predictor variables that had the greatest relative influence were freshwater ecoregion and network percent impervious surface for fish and aquatic macroinvertebrates respectively. Those with the least relative influence were CONF_CL and watershed mine density for fish and aquatic macroinvertebrates respectively. The theme with the greatest cumulative relative influence for fish was geomorphic condition (46.6 %) and for aquatic macroinvertebrates was water quality (41.7 %), while point source pollution was the least important theme for both fish (1.4 %) and aquatic macroinvertebrates (0.8 %) ( Table 4).
The relative thematic influence was variable between response variables; where, for example the flow theme was the most influential for % EPT (27.4 %) but for % intolerant taxa the most influential theme was water quality (87.9 %).

Regional models
Regional differences were evident in the results of the relative influence of predictor variables on responses. For example, the variable with the greatest relative influence for richness in the Atlantic Slope data was July temperature (22.8 %), for the Cumberland River Basin was ground flow recharge (17.6 %), for the Ohio River Basin was average catchment elevation (21.4 %), and for the Tennessee River Basin it was dam storage capacity (8.2) ( Table 6). Additionally, the relative influence varied widely within a given region for the same predictor variable but for different response variables. For example, k-factor (a measure of soil erodibility) had the 24 th most influence on diversity in the Atlantic Slope but was the 4 th most influential variable in the % herbivore taxa model. The standard deviation associated with the relative influence of predictor variables for a specific region was also variable; where, for example, for invertivore diversity the standard deviation was 2.4, 2.6, 5.1, and 4.3 for the Atlantic Slope, the Cumberland River Basin, the Ohio River Basin, and the Tennessee River Basin respectively.
The mean relative influence for themes were also somewhat variable across regions; where, for example, flow was the most influential theme for each region for the BRT model for the % of species that prefer coarse sediment. However, for the % of tolerant fish species the flow theme represented the 4 th most influential theme for the Atlantic Slope (x̄ = 1.7 %), 3 rd most for the Cumberland River Basin (x̄ = 2.5 %), 2 nd most for the Ohio River Basin (x̄ = 2.7 %), and 1 st most for the Tennessee River Basin (x̄ = 3.6 %).
Generally, when averaged across all response variables, geomorphic condition was the most influential theme followed by flow, non-point source pollution, water quality, connectivity, and point source pollution. With little variation, the themes were consistent in terms of their relative influence across regions; for example, non-point source pollution was the 4 th most influential theme for the Atlantic Slope, the Ohio River Basin, and the Tennessee River Basin but the 5 th most influential for the Cumberland River Basin.

Fish models
Model results suggest there are distinct areas of the App LCC region that exhibit characteristics of poor biological condition (i.e., watershed score) (Figure 3). For example, western Tennessee (i.e., the Southeastern Plains) and the norther-most portion of the region (i.e., the Atlantic Highlands). Additionally, large rivers tended to exhibit poor biological condition.
This was especially true in areas with high agricultural land use (i.e., east Tennessee). Regions where the models displayed relatively good biological condition were largely in areas free from agricultural land use. However, even in areas with good biological condition, large rivers exhibited poor biological condition.
Pearson correlations between the global and regional models were variable and ranged from 0.33 (% intolerant taxa) to 0.99 (% invertivore taxa) for pairwise comparisons. The overall indices for watershed scores (i.e., biological condition) were moderately correlated (0.69). The targets were highly correlated for functional group (0.98) but were relatively poorly correlated for taxa quality (0.34) ( Table 8) ( Figure 5).

Aquatic macroinvertebrate model
We were unable to develop regional models for aquatic macroinvertebrates due to the data set being sparse in some regions. The global model for aquatic macroinvertebrates showed poor ability to identify areas with good and bad biological condition. Only very large geographic regions were distinguishable from one another. For example, the western portion of the App LCC region was identified as having relatively higher biological condition when compared to the eastern and northern portions of the region. Where the fish models identified large rivers as having poor biological condition, the aquatic macroinvertebrate models were unable to identify biological condition at this scale ( Figure 4). Additionally, the Pearson correlation coefficient between the aquatic macroinvertebrate watershed score and those of the global (Figure 5) and regional fish watershed scores was 0.01 for both.

<D>Discussion
Freshwater ecosystems are highly threatened worldwide. In North America the Appalachian LCC region is a particularly valuable aquatic ecosystem due in large part to its high aquatic biodiversity. It is also a region of great conservation concern due to current and projected rapid land cover change (i.e., urbanization), biodiversity loss, habitat fragmentation, and climate change (Terando et al. 2014, Leonard et al 2017. Models such as the ones present here and others (e.g., Esselman et al., 2013 andHill et al., 2017), are valuable resources for identifying areas of greatest concern for restoration activities and those areas in greatest need of increased and/or continued protection.
Multimetric indices may be beneficial in identifying biological condition across large geographies. Our results suggest that multimetric scores (i.e., watershed scores) corroborate one another at the global and regional scales (see table 8), indicating that the predictor variables we used were responsive in a consistent manner across scales. This is an important finding as the ability to rely on a consistent set of landscape variables enables cross-regional and cross-scale comparisons; which can be particularly important when data sets are assembled from multiple sources (Esselman et al. 2013).
The relative influence of individual predictors to our boosted regression tree models was variable between regions in our regional models. This indicates the high variability between regions regarding what environmental factors influence community structure. These findings are not surprising as they are congruent with other large-scale biological condition assessments (Esselman et al. 2013;Hill et al. 2017). While individual predictors were relatively variable between regions, our models suggested themes (i.e., combined predictors) were more consistent yet still exhibited some discrepancies between regions. The themes that exhibited the greatest influence on fish watershed scores across the entire region were geomorphic condition and flow regime. Geomorphic condition in our analysis was a measure of soil stability and erosivity. This finding is unsurprising as aquatic taxa are known to be sensitive to sedimentation. Surprisingly, point source pollution and connectivity were least influential. This is likely due to these two factors being measured somewhat sparsely (i.e., few watersheds have NPDES sites) on the landscape. For example, NPDES density, a measure of point source pollution, is likely to exhibit strong influence on aquatic biota in the immediate area but these sites are scattered across the landscape and thus across a large region such as the App LCC region, they were less influential in our models than predictors that are found in nearly every watershed (e.g., soil erosivity).
High correlation between the global and regional fish multimetric indices suggests the predictor variables used were responsive across both scales. While there was high correlation between global and regional fish multimetric indices (i.e., watershed score) there was very low correlation between multimetric indices from either of the fish models and the global aquatic macroinvertebrate model. The low correlation between the fish watershed scores and aquatic macroinvertebrate watershed scores is most likely due to the low number of sites (N=334) resulting in poor spatial coverage for our aquatic macroinvertebrate data set, as compared to the good spatial distribution and relatively high number of sites (N=8,355) in our fish data set.
Correlations between fish and aquatic macroinvertebrate indices are likely to be higher when the same sites are sampled for both. For example, at 18 sites in the River Raisin watershed in Michigan, correlations between fish and aquatic macroinvertebrate indices ranged from 0.21 to 0.51 (Lammert and Allen 1999). Additionally, stressors have been found to affect fish and aquatic macroinvertebrates similarly. Free flowing stream sections of the Fox River, IL exhibited higher index scores for both fish and aquatic macroinvertebrates when compared to sections of the river above low-head dams (Santucci et al., 2005).
Esselman et al. (2011) suggested regional stratification for analyses that span large geographies, while also including multiple biological response variables, and using statistical modeling techniques that are non-linear would lead to improved multimetric indices. This is precisely the approach we took modeling biological condition across the 15 state App LCC region. However, one major difference between our study and previous studies is that previous work to developed multimetric indices across broad geographies included regional influences by using Omernick ecoregions (Whittier et al., 2007;Esselman et al. 2013, andHill et al. 2017).
Omernick ecoregions are a spatial framework which includes aquatic and terrestrial ecosystems (Omernick 1987;USEPA 1996); however, we used Freshwater Ecoregions of the World for our analysis because they were developed on a framework that focused exclusively on aquatic ecosystems with fish as the focal taxonomic group (Abell et al. 2008).
Another difference between our work on previous work is that we aggregated both predictor and response variables into themes and targets respectively. By aggregating predictor variables into themes we were able to identify categories of predictors that were influential in structuring fish and aquatic macroinvertebrate communities in aggregate (e.g., watershed score for fish), as well as for individual response variables (e.g., % EPT). Response variable aggregation allows for a greater understanding of how individual and aggregate predictors (i.e., themes) structure fish and/or aquatic macroinvertebrate communities in a more functional manner. For example, two targets (e.g., fish functional group and fish taxa quality) may respond differently to the same stressor. Understanding which higher order stressors (theme) affect our targets (fish functional group) may lend insight into more effective regional management and conservation strategies. However, thematic importance was variable between regions and care should be taken when using large scale modeling efforts such as ours to make management recommendations across large regions that span more than one ecoregion.
While we were able to develop regional models for fish indices, we were, unfortunately, unable to develop regional aquatic macroinvertebrate models due to low numbers of sample sites. As aquatic macroinvertebrates are known to be good indicators of biological condition (Gerritsen et al. 2000;Pond et al. 2008) we believe there needs to be increased systematic sampling efforts that result in an increased number of sites that are spatially dispersed across multiple freshwater ecoregions and that includes known areas of both high and low biological integrity. With an extensive sampling approach that includes both fish and aquatic macroinvertebrates the ability to develop accurate models of large geographies will improve.
Such improved spatial models will lend greater confidence in such indices which will only improve the ability of managers, conservation practitioners, and scientists alike to address and meet regional conservation goals.         non-point source pollution (% agricultural land cover, % impervious surface, and % natural land cover in the watershed, active river area, and catchment), and point source pollution (CERCLA site density, permit compliance system site density, toxic release site desnity, coal mine, wind turbine and natural gas well density, and mine density). Basin, and T=Tennessee River Basin. Table 7. Summary of the regional average relative influence for each theme on response variables modeled with the boosted regression trees for fish. The average thematic % relative influence is the percent relative contribution of each theme to the BRT model. Themes are comprised of flow condition (NID storage, dam density, altered streamflow, agricultural and industrial water withdrawals), geomorphic condition (erosive and resistive forces), connectivity (dam and road crossing density at the watershed and catchment levels), water quality (total nitrogen load, total phosphorous load, and dissolved organic carbon), non-point source pollution (% agricultural land cover, % impervious surface, and % natural land cover in the watershed, active river area, and catchment), and point source pollution (CERCLA site density, permit compliance system site density, toxic release site desnity, coal mine, wind turbine and natural gas well density, and mine density). Table 8. Pearson correlation coefficients for the global and regional models for all response variables. Correlations for the same response variables for each model are highlighted in gray.