Land-use and water-quality threats to current and historical Cryptobranchus alleganiensis streams across multiple ecoregions

Freshwater systems have experienced greater biodiversity loss than any other systems on Earth. The Hellbender Cryptobranchus alleganiensis (Daudin, 1803), a fully aquatic salamander species endemic to streams and rivers across the Midwestern and Eastern United States, has experienced drastic declines due to habitat degradation, accelerated sedimentation, aquatic contaminants, and infectious diseases. Although declining water quality is often suggested as a major contributing factor to Hellbender population declines, few studies have evaluated the effects of land use on water quality and Hellbender population status across multiple watersheds and ecoregions. In this study, we sought to understand the relationships between Hellbender population declines with both anthropogenic land use and water quality. We sampled 30 stream sites across a gradient of land use and Hellbender population status in Tennessee and North Carolina, USA, to evaluate relationships between Hellbender populations, water quality, and concentrations of heavy metals and the herbicides glyphosate and atrazine, the latter which we measured via a novel passive sampling technique. We found ecoregion-level differences in atrazine, Cd, pH, and conductivity, with the Interior Plateau ecoregion having the greatest variance in water-quality measurements. Although % watershed development was not different among ecoregions given the sites evaluated, Blue Ridge watersheds were overall less developed than watersheds in the Interior Plateau. We also found relationships between declining Hellbender population status and both increased watershed development and declining estimates of water quality, suggesting that increased anthropogenic watershed disturbance can lead to water-quality declines that can negatively affect the species. Because this research focused on assessing the prevalence of commonly encountered aquatic contaminants, our results and study design are applicable to lotic habitats within the evaluated ecoregions and similar habitats throughout the geographic range of the Hellbender. Our study highlights parameters that should be considered when investigating the effect of anthropogenic disturbance and declining water quality on Hellbender population status and emphasizes the importance of effective land management strategies that reduce anthropogenic effects to freshwater systems.

Freshwater ecosystems have experienced greater biodiversity losses than any other ecosystem on Earth (Dudgeon et al. 2006, Tickner et al. 2020, Vári et al. 2022).Anthropogenic disturbances have degraded the integrity of lotic ecosystems, and these affects tend to be greatest in landscapes with high levels of urban and agricultural development (Martinuzzi et al. 2014).Characteristics of degraded streams include substantial alterations in channel morphology and hydrological patterns, declines in water and habitat quality, and reductions in aquatic species diversity.Biodiversity declines within freshwater ecosystems have been associated with reduction of riparian buffers, fine sediment loading, elevated conductivity, and point and nonpoint source pollution (King et al. 2011).
Nonpoint source pollution includes fertilizers, chemicals, sediment, and bacteria that are transported into waterways through surface runoff, erosion, stormwater discharge, atmospheric deposition, drainage, or seepage.Agriculturerelated nonpoint source pollution (i.e., sediment, manure, and herbicides) represents one of the greatest threats to aquatic biodiversity in North America and globally (Hascic and Wu 2006, Vörösmarty et al. 2010, Rasmussen et al. 2015).Goolsby et al. (1991) reported high annual levels of herbicides in streams and aquifers across 10 Midwestern states, often exceeding recommended levels for drinking water.Atrazine, an herbicide detected by Goolsby et al. (1991) and the 2 nd -most used herbicide in the agricultural sector since 2001 (Atwood and Paisley-Jones 2017), can cause negative morphological effects in amphibians (Rohr andMcCoy 2010, Van Der Kraak et al. 2014).Glyphosate, the most-used herbicide in agriculture since 2001 (Atwood and Paisley-Jones 2017), has also been reported in freshwater systems (Coupe et al. 2012, Annett et al. 2014, Battaglin et al. 2014, Herek et al. 2020) and can be toxic to several freshwater taxa at different concentrations (Annett et al. 2014, Wagner et al. 2016).
Heavy metals can also harm wildlife through acute and chronic exposure and are particularly impactful to aquatic environments and associated biodiversity including fish, amphibians, mussels, and macroinvertebrates (Schuler andRelyea 2018, Ali andKhan 2019).Although naturally occurring, elevated heavy metal concentrations in freshwater systems have been associated with mining and refinement operations, landfills and transportation, and agricultural fertilizers and pesticides (Singh et al. 2011).Toxic levels of heavy metals associated with industrial practices can reach freshwater systems through precipitation, runoff, erosion, or direct discharge (Singh et al. 2011).Most water bodies in the Northern Hemisphere are contaminated by Cd, Hg, Pb, and Cu (Naimo 1995).Land-altering practices and associated pollution of aquatic systems are likely to continue; therefore, greater research is needed to understand landscapescale disturbances and how the resulting pollutants and other environmental impacts are distributed throughout ecosystems.Such data are essential for the creation of effective management plans that can potentially mitigate the negative effects of pollution on aquatic biodiversity.
The Hellbender Cryptobranchus alleganiensis (Daudin, 1803) is a fully aquatic, large-bodied salamander species that has drastically declined across its historical range since the 1980s (Freake andDePerno 2017, Pitt et al. 2017).Impoundments and dams, accelerated sedimentation, emerging pathogens, and compromised water quality are suggested as contributors to these declines (Pitt et al. 2017).Hellbenders are generally restricted to streams and rivers with cool, welloxygenated flowing water with a heterogeneous substrate composed of unembedded rocks of diverse shapes and sizes (Keitzer et al. 2013, Da Silva Neto et al. 2019).Individuals of all age classes use the interstitial spaces between substrate rocks as shelter and breeding sites (Hecht et al. 2019, Unger et al. 2020a).Because Hellbenders are fully aquatic, benthic, and long-lived; have cutaneous respiration; and are sensitive to degradation of aquatic habitats, they can be reliable indicators of stream health where they occur (Keitzer et al. 2013, Pugh et al. 2016).
Despite extensive research on Hellbenders' habitat requirements and demographics of several populations across many states within the USA (Taber et al. 1975, Nickerson et al. 2003, Humphries and Pauley 2005, Burgmeier et al. 2011, Hecht-Kardasz et al. 2012), little is known about the effects of land use and associated chronic environmental stressors on the species across multiple watersheds and ecoregions (Solis et al. 2007, Keitzer et al. 2013, Pugh et al. 2016).In this study, we sought to evaluate the following research objectives across the Eastern Hellbender (Cryptobranchus alleganiensis alleganiensis) range in Tennessee and North Carolina, USA: 1) investigate how ecoregionand watershed-level land use are associated with water quality in streams currently or historically occupied by Hellbenders, and 2) assess the relationships among ecoregion and watershed land use, herbicides, heavy metals, and waterquality parameters and Hellbender population status.Specifically, we expected that locations with greater levels of agricultural land use would have greater herbicide and heavy metal concentrations and more impaired water-quality parameters compared with locations with lower levels of agricultural land use.Further, we predicted that Hellbender population status would mirror patterns of aquatic impairment, with stable and recruiting Hellbender populations in streams with the lowest amounts of anthropogenic disturbance and more robust water-quality measurements.Collectively, assessment of these impacts at the watershed and landscape scales will help researchers identify factors that contribute to Hellbender population declines and ultimately protect additional aquatic biodiversity.

METHODS
To answer our research questions, we conducted a field study in the summer of 2017 at 30 sites distributed throughout the geographic range of the Hellbender in Tennessee and North Carolina.We assessed watershed-level development with geospatial modeling, collected water-quality data, and tested stream sediment for heavy metals (i.e., Cd, Hg, and Pb).In addition, we developed a novel sampling approach to quantify concentrations of the herbicides glyphosate and atrazine at each site.Because these pollutants are commonly associated with anthropogenic activities and are reported to cause harmful effects to aquatic ecosystems and fauna, we used them as indicators of habitat degradation where Hellbenders occur currently and historically throughout the study area.We used a combination of inferential and descriptive approaches and ordinal regression modeling to analyze our data.

Study area
The Hellbender geographic range extends from northern Alabama and Georgia, north through Pennsylvania, and into southern New York (Petranka 1998), but portions of Tennessee and North Carolina provide some of the best remaining habitat for Hellbenders in the eastern USA (Freake andDePerno 2017, USFWS 2018).State boundaries are often not biologically relevant to species distribution patterns, so we performed our analysis considering the study area as the joint species range between Tennessee and North Carolina (Fig. 1).The Hellbender range across both states is composed of ~1300 hydrologic unit code (HUC) 12 watersheds (i.e., the smallest hydrological unit identified by a unique 12-digit number).Of these, ~230 watersheds are in North Carolina, 1005 watersheds are in Tennessee, and 96 watersheds are shared between Tennessee and North Carolina.In addition, the joint range encompasses the following 6 United States Environmental Protection Agency level III ecoregions: Southeastern Plains, Ridge and Valley, Southwestern Appalachians, Central Appalachians, Interior Plateau, and Blue Ridge (Omernik 1987).The species range in Tennessee encompasses all 6 ecoregions, whereas the species range in North Carolina only includes the Blue Ridge ecoregion.The range of the Hellbender in Tennessee includes several large metropolitan areas (e.g., Nashville, Knoxville, Chattanooga) compared with the species range in North Carolina, which lacks urban areas of comparable size.

Watershed development analysis and site selection
To aid in site selection, we used geospatial analysis in ArcGIS (version 10.3,Environmental Systems Research Institute,Redlands,California).Using data from the National Land Cover Database (https://doi.org/10.5066/P97S2IID;USGS 2011), we quantified the percentage of landcover in each HUC12 watershed with combined agricultural and urban development.See Appendix S1 for detailed geospatial methods.We then divided watersheds into 3 groups based on the % of watershed surface area dedicated to combined agricultural and urban development: low (0-25% of watershed surface area), moderate (25-50% of watershed surface area), and high watershed development (agricultural and urban development is ≥50% of watershed surface area; Lu and Weng 2006).
Next, we used a stratified-random sampling design to select sampling sites.We identified streams within the Hellbender range in Tennessee and North Carolina with historical species presence and survey data (i.e., physical surveys and environmental DNA; 185 potential streams: 155 in North Carolina, 30 in Tennessee) that had adequate historical and current data for potential inclusion.We stratified site selection across watershed development categories (low, moderate, high; Table S1), and then randomly selected 5 sites/ category for each state for a total of 30 sites.If a randomly selected site was not accessible (i.e., surrounded by private property without a logistically feasible access point), we removed the site from the list and randomly selected another site.The Central and Southwestern Appalachians ecoregions do not have known Hellbender populations, whereas the Southeastern Plains have only 1 river where they are known to occur.Therefore, our 30 sampling sites were located within the Southeastern Plains (n 5 1), Interior Plateau (n 5 9), Ridge and Valley (n 5 2), and Blue Ridge (n 5 18) ecoregions, which span >550 km from east to west and include a considerable portion of the southern geographic range for the Hellbender (Fig. S1).

Hellbender population status
We assigned Hellbender population status for each site (Fig. S2) based on survey data (on-site surveys, environmental DNA sampling) collected by state agencies and universities in Tennessee and North Carolina and on historical datasets such as Biodiversity Information Serving Our Nation (www.gbif.us/)and the Atlas of Amphibians in Tennessee (www.apsubiology.org/tnamphibiansatlas).We used population status and trend categories proposed in the 2018 United States Fish and Wildlife Species Status Assessment (USFWS 2018) as follows: 1) recruiting when there was evidence of recruitment since the year 2000, 2) declining when there was documentation of declines in abundance or recruitment, 3) functionally extirpated when surveys indicated the presence of adults only, and 4) extirpated when surveys indicated that the species was extirpated or the habitat was no longer deemed suitable.For 50% of sites (15/30), we used expert opinion to determine the appropriate population status when >1 category was possible.For example, if recent surveys at a site identified only a few individuals present and previous surveys indicated a decline in abundance or recruiting over time, we used our knowledge of the site to determine if the population was most likely declining or functionally extirpated (Table S2).

Herbicide exposure
Passive sampling protocol development To assess the potential exposure of Hellbenders to herbicides at our sam-pling sites, we developed a novel sampling approach for quantifying the concentration of atrazine and glyphosate in the water.Our goals were 1) to develop a biologically relevant sampling technique that would expose the sampling units to the same environmental conditions that an individual Hellbender would experience in natural conditions and 2) to replicate the diffusion of analytes across the skin surface.To do so, we used a passive sampling technique that relied on diffusion to collect neutral and lipophilic contaminants diluted in water (Fig. S3).The principle of passive sampling is based on the ability of analyte molecules to flow freely (i.e., nonassisted) from the sampled medium into a receiving phase in a sampling device because of the difference between analyte concentrations of the 2 media (Seethapathy et al. 2008, Salim andGórecki 2019).Analytes flow into the receiving phase (i.e., deionized water inside the passive sampler) until the analyte concentration reaches equilibrium between the receiving phase and the sampled medium (i.e., stream water; Seethapathy et al. 2008).We designed the passive samplers to mimic the process in which active compounds in herbicides are absorbed through plant and animal cell membranes via simple diffusion (Sterling 1994).Amphibians are especially susceptible to environmental pollutants because of their permeable skin, which they use for gas exchange and water uptake (Brühl et al. 2011, Larsen 2021).Therefore, we assumed that the concentration of analytes in the receiving phase of the passive sampler represented the concentration that an individual Hellbender would be exposed to in situ.We acknowledge that absorption of herbicides through animal skin is a complex process that may not be fully replicated using passive samplers.See Appendix S2 for detailed methodology for protocol development, laboratory assessment and results, calibration, and statistical analysis.
Herbicide sampling design To measure herbicide contamination, we deployed 2 passive samplers at each of the 30 stream sites.We placed each passive sampler in a protective plastic mesh tube (ALS Global, Brisbane, Australia), attached an 85.05-g weight at each end of the mesh tubing, and closed both ends with zip ties.To mimic the benthic behavior of Hellbenders, we deployed passive samplers under or next to potential cover rocks and as close as possible to the substrate.We secured each protective tube to a 0.91-m piece of rebar sunk 0.61 m into the substrate with 31.8-kg(70-lb) braided fishing line.
To account for seasonal variation in herbicide concentration, we deployed the 2 sampling units/site each month during June, July, and August.We deployed samplers every 30 d and usually took 2 d to deploy samplers at all sites because of travel time among sites.Each group of sampling units was deployed for a 48-h exposure period to allow the concentration of analytes between the sampled medium (i.e., river water and substrate) and the receiving phase (i.e., deionized water inside passive sampler) to reach equilibrium.
Sampler deployment was targeted to avoid effects from environmental variables such as heavy rain, which can dilute analytes and alter analyte exchange dynamics (Seethapathy et al. 2008).We recovered sampling units in the order they were deployed in the respective streams.Immediately after recovering the sampling units, we made an incision at the end of the passive sampler and stored the contents of the receiving phase in 50-mL centrifuge tubes (Thermo Fisher Scientific, Waltham, Massachusetts).The tubes were stored on ice and later frozen at 2207C until subsequent processing at the Biological Small Molecule Mass Spectrometry Core at the University of Tennessee, Knoxville, Tennessee.

Herbicide detection
We analyzed herbicide contamination using reverse-phase ultra-performance liquid chromatography coupled to a high-resolution mass spectrometry (Q Exactive ™ Plus Hybrid Quadrupole-Orbitrap ™ Mass Spectrometer; Thermo Fisher Scientific).We injected a total of 10 lL from each sample onto an Accucore ™ Phenyl-X HPLC column (Thermo Fisher Scientific), which measured 100 Â 2.1 mm and contained 2.6-lm particles.We eluted analytes with a 2-solvent system (A 5 0.1% formic acid in water, B 5 methanol) with the following gradient: 95% A and 5% B at 0 min, 50% A and 50% B at 2.5 min, 10% A and 90% B at 3 min, 10% A and 90% B at 7 min, 95% A and 5% B at 10 min, and 95% A and 5% B at 12 min.To ionize samples, we used heated electrospray with a sheath gas flow of 35, auxiliary gas flow of 10, spray voltage of 3.6 kV, capillary temp of 3007C, S-lens radio frequency level of 50, and auxiliary gas heat of 3057C.We performed mass spectrometry with a polarity-switching method, which acquired in negative mode from 0 to 3.6 min and in positive mode from 3.6 to 12 min.All other parameters were constant for the duration of the run (resolution 5 140,000, Automatic Gain Control target 5 1e 6 , max inject time 5 250 ms).Glyphosate was detected as the [M -H] -with a retention time of 1.5 min, and atrazine was detected as the [M 1 H] 1 with a retention time of 5 min.After collecting these data, we converted .RAW files to .mzmlfiles and used Metabolomics Analysis and Visualization EngiNe to integrate the peaks (Melamud et al. 2010, Clasquin et al. 2012).

Measurement of heavy metals
We estimated heavy metal (Cd, Hg, and Pb) concentrations because these metals are known to have negative impacts on stream biota (Bergeron et al. 2010, Ali and Khan 2019Friberg et al. 2019, Eagles-Smith et al. 2020).To do so, we used a modified sampling method described by the United States Geological Survey (Nelson et al. 2014) to collect 1 stream sediment sample/site (n 5 30).At each site, we selected areas with accumulated fine sediment in slow-moving water and collected 50-mL of sediment into 50-mL centrifuge tubes (Thermo Fisher Scientific).If sediment was com-posed of mud or silt, we collected samples by pressing a 50-mL centrifuge tube 2 cm into the sediment surface.If sediment was composed mainly of sand or loose material, we used a small plastic scoop to collect the top 2 cm of sediment and placed the substrate into a sterile 50-mL centrifuge tube.In areas with larger substrate structure (i.e., pebble, cobble, and boulders), we moved those larger substrate materials until we reached fine sediment (>2 cm in depth).We stored each sample on ice during transport and then froze and shipped samples overnight for heavy metal analysis at the Center of Environmental Science and Engineering at the University of Connecticut, Storrs, Connecticut.
To evaluate Cd and Pb contamination in streambed sediment, we brought samples to room temperature and placed 0.35 g of each sample into a separate hot block tube.We added 5 mL of concentrated trace-metal-grade nitric acid to each tube, placed each sample on a hot block, and refluxed each sample for 4 h at 957C.We cooled each sample and added 2 mL of deionized water and 3 mL of tracemetal-grade hydrogen peroxide.We then heated the sample in the hot block until effervescence subsided.We cooled samples and filled them to a final volume of 25 mL with deionized water.We analyzed samples on an Inductively Coupled Plasma Mass Spectrometer (DRCe; PerkinElmer ® , Waltham, Massachusetts) with standard protocols.We also analyzed interference check solutions (A and A 1 B; High-Purity Standards, Charleston, South Carolina) with all sample runs to compensate for any matrix effects that would interfere with sample analysis.
To evaluate the concentration of Hg in streambed sediment samples, we measured inorganic Hg (total Hg) instead of methylmercury (MeHg).Although MeHg may be more relevant for wildlife toxicity, inorganic Hg is found at higher concentrations in aquatic systems (Eagles-Smith et al. 2018, Obrist et al. 2018).MeHg is created when inorganic Hg is methylated by aquatic microorganisms (Eagles-Smith et al. 2018, Hsu-Kim et al. 2018, Obrist et al. 2018).In aquatic systems, the conversion process is complex and is influenced by several factors, such as chemical composition of the sediment, pH and dissolved O 2 (DO) of water, and microbial community structure (Hsu-Kim et al. 2018, Obrist et al. 2018).We assumed that we would be more likely to detect inorganic Hg at the study sites, which would indicate potential exposure to MeHg if methylation conditions occurred.
In the laboratory, we brought samples to room temperature and placed ~0.25 g of each sample in a separate hot block tube.We added 2.5 mL of deionized water and aqua regia (3:1 HCl:HNO 3 ) to each tube and placed each sample in the hot block at 957C for 5 min.We removed samples from the hot block, cooled them to room temperature, and added 5 mL of deionized water and 7.5 mL of potassium permanganate solution.We heated each sample to 957C for 30 min.We added 1.5 mL of hydroxylamine hydrochloride to each sample and filled samples with deionized water to a final volume of 50 mL.We analyzed samples on a Perkin Elmer FIMS Cold Vapor Atomic Absorption Spectrometer using standard protocols.
We used standard quality-assurance procedures for all samples, which included analysis of method blanks, post-digestion spiked samples, and laboratory control samples.Sample mass was too small to duplicate analyses.We used a calibration verification standard and blank to evaluate instrument response initially, every 10 samples, and at the end of an analytical run.

Water-quality variables
During deployment and recovery (June, July, and August) of passive samplers, we recorded pH, specific conductivity (lS/cm), temperature (7C), and DO (%) at each site with a YSI ProDSS multiparameter water-quality unit (Yellow Springs Instruments, Yellow Springs, Ohio).

Statistical analysis
Ecoregion-level watershed characteristic differences To evaluate how watershed size, development, contaminant levels, and water quality varied among ecoregions, we tested for differences in ecoregion means for each variable.Because there were only 2 sites in the Ridge and Valley and 1 site in the Southeastern Plains ecoregion, we removed those 3 sites from the ecoregion-level analysis and assessed these variables only for the Blue Ridge and Interior Plateau ecoregions (n 5 27).Prior to analysis, we assessed whether response variables were normally distributed via frequency distributions, quantile-quantile plots, and the Shapiro-Wilk test for normality of residuals.When response data were not distributed normally, we used log10, natural log, and square-root transformations to approximate normality.When normality was verified, we used independent sample t-tests to evaluate ecoregion-level differences between Blue Ridge and Interior Plateau ecoregions in watershed size (ha), relative % of anthropogenic watershed development, herbicide, heavy metal concentrations (atrazine, Cd, Pb), and water-quality variables (pH, DO, specific conductivity).In the instance when response variables (atrazine and Hg concentration) were not normally distributed (even after transformation), we used the Mann-Whitney U test to evaluate differences in mean ranks.We used the stats package in R (version 4.0.2;R Project for Statistical Computing, Vienna, Austria) for these analyses.
Land-use and water-quality relationships with Hellbender population status To address our questions regarding relationships among land-use and water-quality variables and their associations with Hellbender population status at watershed and ecoregion scales, we conducted 2 analyses.First, we used principal components analysis (PCA) to explore multivariate relationships among 9 landscape and waterquality variables at the watershed scale, which included watershed size (ha), % watershed development, herbicide and heavy metal concentrations (atrazine, Cd, Hg, Pb), and water-quality variables (pH, DO, specific conductance).As described above, we screened predictor data for deviations from statistical normality and used appropriate statistical transformations to approximate normality.We conducted PCA via the FactoMineR package (version 2.3; Lê et al. 2008) in R and evaluated scree plots of resulting eigenvalues and variable relationships based on component loadings for each principal component.We then plotted principal component scores from the first 2 principal components by Hellbender population status to evaluate variable correlations within individual principal components to understand their relative importance for the ordination of population status.
Next, we used a multiple hypothesis approach to evaluate the associations of watershed size (ha), % watershed development, heavy metal/herbicide concentrations (Cd, Hg, Pb, atrazine), and water-quality variables (pH, DO, specific conductance) with Hellbender population status.Prior to analysis, we developed a series of 26 candidate models (Table S3) to evaluate the singular and additive effects of landscape condition and water-quality variables on Hellbender population status.We used cumulative link mixed models (CLMMs) via the ordinal package (version 2019.1.10;Christensen 2015) in R to evaluate the effects of the 9 predictor variables (as fixed effects) on Hellbender population status (response variable) while controlling for nonindependence of sample stream selection within the ecoregions (Blue Ridge, Ridge and Valley, Interior Plateau, and Southeastern Plains; random effect).Briefly, CLMMs are used when response data are ordinal (i.e., data have natural ordered categories, but the distance among categories is not known), and there is a lack of independence in observations because of clustering, blocking, and repeated measurements that necessitates the use of a mixed structure.Prior to fitting models, we evaluated the data for multicollinearity using Pearson's correlation tests among fixed effects.We found no correlations >0.46 among predictor variables, so we did not remove any variables from model fitting.
We used Akaike's Information Criterion adjusted for small sample sizes (AICc) to evaluate relative support of individual models.We considered models competitive when DAICc estimates were ≤2.0 and used model averaging, as recommended by Burnham and Anderson (2002), when model terms were contained in multiple top models.We report model-averaged beta coefficients, SE, and 95% CI for all fixed effects that occurred in multiple competitive models as conducted in Sutton et al. (2014).

Agricultural and urban development within HUC12 watersheds
We found that HUC12 watershed size (total ha) of sampled sites was different among ecoregions (t 14.23 5 22.18, p 5 0.04), with greater watershed size in the Interior Plateau (mean 5 10,691 ha ± 3355 SD) compared with the Blue Ridge ecoregion (mean size 5 7756 ± 2526; Table 1, Fig. 2A-F).Although the mean HUC12 % watershed development of sampled sites was lower in the Blue Ridge ecoregion (mean 5 0.32 ± 0.20) compared with the Interior Plateau (mean 5 0.43 ± 0.21; Table 1), evidence for this difference was weak (t 15.98 5 21.29, p 5 0.21).

Hellbender population status along development and water-quality gradients
Of the 9 principal components generated via PCA, 4 components reflected ecological gradients that we considered to be plausible drivers of Hellbender status and that collectively explained ~83% of the variance in the dataset.Principal component axis 1 (PC1; 39% of variance) described a gradient of water quality and watershed development, where conductivity, heavy metal concentration (Cd, Hg, and Pb), pH, and atrazine were greater in streams with in-creased watershed development (Fig. 3).PC2 (18% of variance) described a gradient where DO decreased as watershed development increased (Fig. 3A, B).PC3 (14% of variance) described a gradient where stream pH was positively correlated with watershed development.PC4 (12% of variance) described a site gradient based on HUC 12 watershed size.Hellbender population status was ordinated primarily along PC1, with recruiting populations in watersheds with lower conductivity, heavy metal concentration, pH, and watershed disturbance.In contrast, functionally extirpated and extirpated populations were in watersheds with progressively increasing quantities of water-quality and watershed disturbance variables.In addition, sites where Hellbenders were classified as extirpated and functionally extirpated ordinated along the portion of PC2 at sites with relatively lower DO and greater watershed disturbance.
Overall, we found support for a relationship between the decline in Hellbender population status (i.e., recruiting vs declining, functionally extirpated, and extirpated) when % watershed development in a HUC12 watershed exceeded ~25% and when specific conductivity at the study site exceeded an estimated ~50 lS/cm (Fig. 4A, B).

DISCUSSION
In this study, we investigated the relationships between agricultural and urban development and water quality in streams currently and historically occupied by Hellbenders in Tennessee and North Carolina.We also asked if ecoregion, % development of watersheds, herbicide and heavy metals, and water-quality parameters were related to Hellbender population status at our sites.Although our study cannot establish direct causality, our results indicate there are relationships between declining Hellbender population status and declining environmental conditions throughout the study region.Specifically, we found that watersheds with declining Hellbender populations were associated with increased watershed development and impaired water quality.In contrast, recruiting Hellbender populations were associated most strongly with a lower percentage of watershed development and lower conductivity, as well as lower concentrations of heavy metals.Aquatic species population declines    1.
are often directly associated with urban and agricultural development through physical habitat alterations and, indirectly, through water-quality impairment and increased pollution.A better understanding of land-use effects and water quality on Hellbender populations may help identify potential causes of population declines, reveal further research needs, and guide strategic watershed-level habitat management.

Watershed development and water-quality relationships with Hellbender population status
Our finding that % watershed development was the most important predictor of Hellbender population status at the watershed scale was unsurprising because watershed development often degrades stream systems and poses a threat to aquatic wildlife (Martinuzzi et al. 2014).Channelization, dredging, loss of riparian buffer, flow alternation, and excess pollution through surface runoff and concentrated deposition are often associated with urban and agricultural development (Wang et al. 2001, King et al. 2011, Martinuzzi et al. 2014).These impacts to aquatic systems often result in alteration or loss of physical and chemical stream functions (i.e., hydrologic, hydraulic, geomorphologic, physiochemical) that negatively affect aquatic fauna (Revenga et al. 2000, Wang et al. 2001).Other studies have indicated that loss of forest at the catchment level contributes to waterquality declines, which negatively affects Hellbender populations (Keitzer et al. 2013, Pitt et al. 2017, Bodinof Jachowski and Hopkins 2018, Wineland et al. 2019).Likewise, our results suggest that recruiting Hellbender populations are associated with lower levels of watershed development (<25%), lower concentrations of heavy metals and atrazine, and lower conductivity.Notably, Da Silva Neto et al. ( 2020) suggested that the greatest area of suitable Hellbender habitat available in Tennessee was located within the Blue Ridge ecoregion, which during our study had lower conductivity levels as well as lower heavy metal and herbicide concentrations than the Interior Plateau and Ridge and Valley ecoregions.Similarly, Bodinof Jachowski et al. ( 2016) predicted that sites with lower percentages of agricultural lands within catchments were more likely to be occupied by Hellbenders within the Blue Ridge and Ridge and Valley ecoregions, although the effects of development on occupancy remained uncertain.
Although we did not find a difference in % watershed development among ecoregions, conductivity, pH, atrazine, Cd, and Hg were all greater in Interior Plateau ecoregion streams compared with Blue Ridge ecoregion streams.Conductivity is often used as a general measure of water quality, though it varies naturally across the landscape because of differences in stream characteristics, climate, geology, and land-use practices among ecoregions.Griffith (2014) reported variable reference values for specific conductivity among the Blue Ridge (4.6-151.5 lS/cm), Interior Plateau (103.0-658lS/cm), Ridge and Valley (10-824.2lS/cm), and Southeastern Plains (10.7-2310.0lS/cm) ecoregions.However, these differences do not account for the influence of anthropogenic disturbance on variation within a particular ecoregion.Previous studies indicate that in some areas, high levels of conductivity (>140 lS/cm) are often associated with anthropogenically caused deforestation, loss of riparian buffers, and elevated erosion and runoff, all with overall detrimental impacts to water quality (Dow and Zampella 2000, Jackson et al. 2017, Bodinof Jachowski and Hopkins 2018, Hutton et al. 2020).
Notably, elevated levels of conductivity have been reported to negatively affect several freshwater taxa including fish, macroinvertebrates, and amphibians (Kimmel and Argent 2010, Chambers 2011, Johnson et al. 2015, Clements and Kotalik 2016, Hutton et al. 2020).In the current study, mean conductivity in sample sites with stable and recruiting Hellbender populations was 44.5 lS/cm, whereas mean conductivity at sampled sites with declining, functionally extirpated, and extirpated Hellbender populations was 97.8 lS/cm, 155.0 lS/cm, and 192.3 lS/cm, respectively (Fig. 4).These results are congruent with previous results by Keitzer et al. (2013), where mean conductivity across 31 streams occupied by Hellbenders in West Virginia was lower (36 lS/cm) than sites where Hellbenders were absent (230 lS/cm).Similarly, Wineland et al. (2019) reported lower mean conductivity (42.33 lS/cm, range 5 0-216 lS/cm) at Hellbender-occupied sites than where the species was extirpated (162.51 lS/cm, range 5 70-546 lS/cm) in the Ohio River drainage in West Virginia.Pitt et al. (2017) also reported a negative association between high levels of conductivity (278 lS/cm) and Hellbender presence in the Susquehanna River drainage in Pennsylvania.Finally, Bodinof Jachowski and Hopkins (2018) reported decreased abundance of subadult/adult Hellbenders and lower proportions of young adults at stream reaches with low catchment-wide riparian forest cover and associated increased conductivity across six stream reaches in Virginia.
Although atrazine concentration was not included in the best-fit models of Hellbender population status, this variable was included in PC1, which was negatively related to Hellbender population status.We found atrazine at 7 sites (6 in Tennessee and 1 in North Carolina) with Hellbender populations assigned as declining, functionally extirpated, or extirpated (Table S4).The mean concentration of atrazine (0.05 ppm) across all 7 sites was 500 Â greater than concentrations that have been shown to cause developmental abnormalities in frogs (0.0001 ppm; Hayes et al. 2002), though 90 Â lower than the 96-h LC50 (4.5 ppm) for Rainbow Trout (Helena Chemical Company 2016).Atrazine is highly mobile through soil, is stable in water, and will precipitate after 24 h of entering an aquatic environment (Helena Chemical Company 2016).The long lifespan of Hellbenders, in combination with these characteristics, may increase the potential for chronic exposure.In addition, atrazine may indirectly affect Hellbender populations and other aquatic organisms through a bottom-up effect.Atrazine is highly toxic to aquatic plant communities that are a food source to invertebrates, which are, in turn, food sources for other aquatic fauna including larval Hellbenders (Helena Chemical Company 2016, Unger et al. 2020b).For example, experimental atrazine additions caused a reduction of algae available to tadpoles exposed during a closed-system, cattle-tank-pond experiment (Boone and James 2003).In addition, Forson and Storfer (2006), Kerby andStorfer (2009), andSifkarovski et al. (2014) suggested that chronic exposure to ecologically relevant concentrations of atrazine can increase susceptibility of some amphibians to pathogens.We did not investigate the direct impacts of atrazine on Hellbender populations, but it is likely that atrazine at the concentrations we observed has the potential to affect Hellbender population health through chronic impacts to individuals and community trophic levels.In addition, elevated levels of atrazine, along with elevated conductivity, in Interior Plateau streams suggests that this ecoregion is experiencing high levels of stream impairment, which is likely because of its relatively high amount of agricultural land use.
Like atrazine, elevated heavy metal concentrations can also negatively affect aquatic organisms, and we detected Cd, Hg, and Pb in stream sediment samples at most of the sample sites.The presence of Cd is often associated with mining and industrial waste, and it has been shown to cause negative effects to amphibians that are chronically exposed to higher concentrations (>0.1 ppm; James andLittle 2003, Friberg et al. 2019).Gross et al. (2007) exposed Northern Leopard Frog (Rana pipiens) tadpoles to multiple concentrations of Cd for 100 d and reported altered larval growth and development at lower Cd concentrations (i.e., 0.005 and 0.02 ppm) than the mean Cd concentration we detected in the Interior Plateau sites (0.52 ± 0.06 ppm).Little is known about the toxicity of Hg to amphibians, although a few studies have reported abnormal embryo development and tadpole mortality associated with exposure to different concentrations of HgCL (Wolfe and Sulaiman 1998).It is important to note that we measured the concentration of inorganic Hg (total Hg), and not MeHg, which is the relevant form of Hg for wildlife exposure.Although our results demonstrate that inorganic Hg is indeed present in several sampled sites and could potentially be methylated, our results do not permit a discussion of true Hg exposure potential.Finally, Pb is a common heavy metal that can be found in concentrations of, on average, 0.004 ppm in water and 0.0018 ppm in the sediment of lakes and rivers across the United States (Berzins and Bundy 2002).The mean Pb concentration we detected in the Interior Plateau is ~1000 Â the mean concentration of Pb found in water and 200 Â the mean concentration of Pb found in sediments of lakes and rivers of the United States.Xenopus laevis tadpoles exposed to Pb concentrations of 0.04 ppm, which is 150Â lower than the concentration we detected on average in the Interior Plateau, experienced considerable decreases in body mass (Berzins and Bundy 2002).Further, Pb can be acutely and chronically toxic to humans and several wildlife species (Pain et al. 2019, Eagles-Smith et al. 2020).The half-life of Pb in the environment is ~20 y; therefore, Pb may pose a threat to Hellbender populations through chronic exposure (Berzins and Bundy 2002).Although our results do not directly link heavy metal concentrations with Hellbender population declines across the study area, we provide baseline data regarding the presence of Cd, Hg, and Pb at several historical and current Hellbender sites.The potential negative effects of these pollutants in combination with other environmental stressors on Hellbender populations cannot be discounted.

Novel passive sampling approach to atrazine detection
Herbicides accounted for ~60% of total United States pesticide use in 2012, with ~64 to ~74 million kg of atrazine used (Atwood and Paisley-Jones 2017), which emphasizes the need for pesticide monitoring approaches.Our pilot and field study results demonstrated a successful novel passive sampling approach to quantify the concentration of atrazine at sampled sites.We did not detect a difference in the concentration of atrazine between passive samplers and known standards, which suggests that passive samplers are an effective tool to assess concentration of atrazine in aquatic environments.In addition, there was low variance among replicates and a nonlinear increase in atrazine detected in each dilution level.This suggests that this approach is precise, and the dilution curve can be used to estimate atrazine concentration from field samples.In contrast, we were unable to effectively use our passive sampling approach to detect glyphosate.The pilot study demonstrated that our analysis method did not detect glyphosate dilutions that were <300 lL/L (Fig. S4).Glyphosate has previously been detected at lower concentrations in freshwater streams using other analytical methods; however, our current method was designed to screen for multiple herbicides that ionized in both positive and negative mode and was not optimized for detection of glyphosate alone (Coupe et al. 2012, Battaglin et al. 2014).Fate and concentration of glyphosate in the environment can vary greatly depending on the characteristics of the chemical compound (e.g., half-life and surfactants), soil and sediment composition, microorganism community, rate of precipitation and runoff, presence and quality of riparian buffer, and compliance with herbicide label recommendations (Coupe et al. 2012, Wagner et al. 2013, Annett et al. 2014, Battaglin et al. 2014).Our inability to identify glyphosate in samples within sites does not mean that glyphosate was not present.It does, however, highlight the limit of sensitivity of our analytical methods and the necessity for a modified and improved protocol.
A key component of our study was our ability to estimate herbicide concentrations that wild Hellbender individuals may be exposed to in situ.Studies using both in vitro and in vivo models have indicated that cutaneous absorption can be a major pathway for herbicide absorption in amphibians (Willens et al. 2006, Brühl et al. 2011).Although our passive sampling approach is not a morphologically identical representation of Hellbender skin, our approach is an effective proxy to estimate relative herbicide exposure concentration of Hellbenders located within sampled streams.In addition, our approach permits investigation of potential herbicide exposure of individual Hellbenders without the need for lethal sampling.We recommend that future research that uses this approach to assess herbicide presence in streams occupied by Hellbenders incorporate a greater number of samples and increase the frequency of passive sampler deployment to account for seasonal variation in herbicide application and delivery.

Study limitations
There are several limitations to our study that must be acknowledged.First, our assessment of Hellbender population status was not entirely based on abundance data and population counts.Instead, we used multiple lines of evidence including occupancy surveys, verified observations, environmental DNA, and expert opinion.Because of imperfect detection during surveys and the impracticability of surveying all habitat available to the species, we acknowledge that the best available data may not reflect true population status.In addition, we used expert opinion as one of multiple methods to assign Hellbender status to HUC12 watersheds, and there are limitations associated with expert opinion, especially when influenced by assumptions regarding species-habitat associations.Research conducted by Bodinof Jachowski et al. (2016) suggests that reliance on specieshabitat associations to infer Hellbender occupancy may lead to misclassification of habitat suitability and may lower presence estimates.Although we do not fully understand the extent that expert opinion influences population status assignment, the use of multiple lines of evidence (i.e., physical surveys and expert opinion) to make informed decisions is a reasonable approach to assign population status.
It is important to note that assessment of true absorption rates due to herbicide and heavy metal exposure falls outside the scope of this study.Exposure potential and total absorption rates are difficult to assess because they are a product of complex interactions among physiological characteristics (e.g., skin thickness and cell membrane biochemistry), species natural history (e.g., foraging behavior), and environmental factors (e.g., water velocity, precipitation rate, soil texture; Willens et al. 2006), which all vary among taxa and life stages.
We did not detect a difference in the amount of development between ecoregions, but this lack of difference may be related to spatial variation in development that we did not assess.Watersheds in the Interior Plateau ecoregion were larger than watersheds in the Blue Ridge; therefore, differences in contaminant concentrations and water quality may be associated with variation in spatial distribution of development within each watershed across ecoregions (Carle et al. 2005, Homer et al. 2015).However, our analysis considered the combined impacts of agricultural and urban development at the ecoregion level, and we did not assess the impact of spatial variation in development within each watershed.In addition, riverine systems in middle and eastern Tennessee and western North Carolina have been affected directly or indirectly by point and nonpoint source pollution since the early 1900s (Howells 1990, Brooks andSouthworth 2011).It is important to consider historical land alteration and legacy chemicals when assessing current water-quality parameters and the presence and concentration of pollutants in freshwater environments.
Here we investigated how ecoregion and watershed development is associated with water quality, herbicide presence, and heavy metals and how these variables are associated with Hellbender population status at 30 sites across the species range in Tennessee and North Carolina.Most importantly, our findings indicate that agricultural and urban development and higher levels of conductivity are negatively associated with Hellbender population status.Specifically, medium to high watershed development and specific conductivity >~50 lS/cm were negatively correlated with Hellbender population status.Overall, our study highlights numerous stressors potentially affecting Hellbender populations, but with unknown population-level effects that should be evaluated further to better understand relationships among watershed disturbance, declining water quality, and Hellbender population status.We recommend long-term assessment of water quality and herbicide presence at known Hellbender locations to permit a quantitative assessment to document the effect of land-use change on water quality and potential chronic effects of these changes on Hellbender populations.Finally, our study emphasizes the importance of effective land management and implementation of best management practices (Broadmeadow andNisbet 2004, Yates et al. 2007) across the landscape to reduce anthropogenic effects to freshwater systems.

Figure 1 .
Figure 1.Map of sampling sites located across a % development gradient within hydrologic unit code 12 (HUC12) watersheds across the Hellbender (Cryptobranchus alleganiensis alleganiensis) range in Tennessee and North Carolina, USA.

Figure 3 .
Figure 3. Principal component 1 (PC1) site scores vs PC2 site scores by Hellbender (Cryptobranchus alleganieneis alleganiensis) population status (A) and mean (±SE) PC1 scores vs mean (±SE) PC2 scores by Hellbender population status (B) in an ordination diagram of a principal component analysis for 9 variables (pH; dissolved O 2 ; conductivity; concentrations of atrazine, Cd, Pb, and Hg; watershed size; and % watershed development) and 4 Hellbender population status categories (RE 5 recruiting, DE 5 declining, FX 5 functionally extirpated, EX 5 extirpated) at 30 sampled sites within the species range in Tennessee and North Carolina, USA.Component loadings for variables (in parentheses) indicate the correlation between a variable and a principal component.

Figure 2 .
Figure 2. Boxplots displaying differences and trends in concentrations of atrazine (A), Cd (B), and Hg (C); watershed area (D); pH (E); and conductivity (F) between Blue Ridge (BR) and Interior Plateau (IP) level III ecoregions at 27 stream sites in Tennessee and North Carolina, USA.Boxes represent the interquartile range, and the center line represents the median.Whiskers represent maximum and minimum values, excluding outliers.Dots outside of the whiskers represent outlier values.Mean values are indicated with an * within each ecoregion.A breakdown of values can be found in Table1.

Figure 4 .
Figure 4. Boxplots of % watershed development (A) and conductivity (B) vs Hellbender (Cryptobranchus alleganieneis alleganiensis) population status (RE 5 recruiting, DE 5 declining, FX 5 functionally extirpated, EX 5 extirpated) in 30 stream sites within 4 level III ecoregions in Tennessee and North Carolina, USA.These variables were included in a single top model that explained the effects of land use, aquatic contaminants, and water-quality parameters on Hellbender population status.Boxes represent the interquartile range, and the center line represents the median.Whiskers represent maximum and minimum values, excluding outliers.Dots outside of the whiskers represent outlier values.Mean values within population status categories are represented by *.

Table 1 .
Mean (±SD) for landscape, heavy metal and atrazine concentrations, and water-quality variables from 30 stream sites within Tennessee (n