{'fraction_non_alphanumeric': 0.02481026890283894, 'fraction_numerical': 0.015159748899091165, 'mean_word_length': 5.357636406957351, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Analysis of Multiple Stressors on Resident Benthic Communities in a California Agricultural Stream

This 3-year study (2015-2017) was designed to characterize benthic communities (macroinvertebrates) and physical habitat in an agriculturally dominated waterbody in the Central Coast area of California (Santa Maria River). Benthic communities as represented by various metrics that represent richness, composition, tolerance/intolerance and trophic measures were used as response variables for the various stressors described below. Concurrent water quality evaluations, physical sediment parameters (grain size and total organic carbon [TOC]), pyrethroids, bulk metals—including simultaneously extracted metals (SEM) and acid volatile sulfides (AVS) ratios—and nutrients were measured. The relationship of various benthic metrics to physical habitat metrics, pyrethroids, metals, nutrients and sediment characteristics was evaluated for the 3-year data set. Total physical habitat scores in this watershed were considered to be poor. Samples collected for various sediment chemistry measurements were from depositional areas (fine grain areas primarily silt and clay) where hydrophobic chemicals such as pyrethroids could be found if sources exist. Dominant benthic taxa collected were generally considered to be tolerant to moderately tolerant of environmental stressors and rated as impaired based on a benthic index. Potentially toxic sediment concentrations of arsenic, cadmium and nickel were reported at various sites based on a comparison with existing threshold effect levels. Pyrethroid concentrations interpreted by using a highly protective toxics units approach with a laboratory sensitive taxon (Hyalella) suggested potential toxicity at various sites. Nutrient concentrations could not be interpreted within the context of potential impairment because the State of California has not developed nutrient criteria. The results of the stepwise linear regression models comparing benthic metrics with all environmental variables showed that TOC was the most important variable shaping the benthic communities. In contrast, pyrethroids, metals and physical habitat were not shown to be significant factors shaping benthic communities. The summary multivariate canonical correlation analysis indicated that less stressed, more diverse benthic communities tended to be associated more with TOC-rich finer sediments and lower concentrations of phosphorous-based nutrients, and more stressed, less diverse benthic communities tended to be associated with less organically rich, somewhat less fine sediments and higher phosphorous concentrations.


Introduction
In the early 1990s, it was estimated that approximately 30% to 50% of the earth's land surface was impacted by non-point source pollutants. 1 Of the various activities associated with non-point source pollution, agricultural activities have been identified as one of the major contributors to impairment to surface waters of the United States. 2 The United States has more than 330 million acres of row crop agricultural land that produces an abundant supply of food and other products. 3 Multiple stressors that can impact resident biological communities as a result of agricultural activities are associated with soil erosion, feeding operations, grazing, plowing, animal wastes, application of pesticides, irrigation water, and fertilizers. 4 Therefore, decision-makers need to know the impacts of stressors associated with agricultural activities on resident biological communities in water bodies that are in close proximity to agriculture.
The common approaches used to assess impairments in agricultural water bodies due to multiple stressors are chemical monitoring, toxicity testing and biological assessments (bioassessment). Bioassessment, formally defined as a quantitative survey of physical habitat and biological communities of a water body, is a well-established approach for determining the ecological condition of stream and river systems. [5][6][7][8] Bioassessments provide a useful approach for integrating effects from physical, chemical, and biological stressors on aquatic organisms. The underpinnings of bioassessments are that the structure and function of an aquatic biological community can provide critical information about the quality of the surface water and sediment. Bioassessments are extremely valuable for determining the status of aquatic biological communities across large spatial scales and land use types (agricultural and urban). Information on the status of resident biological communities is particularly useful for determining impaired water bodies, developing Total Maximum Daily Loads (TMDLs), and measuring success of voluntary or regulatory actions. Bioassessments serve monitoring needs through three primary functions: (1) screening or initial assessment of conditions, (2) characterization of impairment and diagnosis, and (3) trend monitoring to evaluate improvements 2 Air, Soil and Water Research from mitigation practices or further degradation. In addition, bioassessments also provide a direct means of measuring compliance with the goal of biotic integrity stipulated under the Clean Water Act because assemblages of aquatic organisms (ie, macroinvertebrates) comprise taxa that are differentially responsive to different environmental stressors.
Synthetic pyrethroids are a class of insecticides that have been used for over 35 years to control pests in over 120 economically important agricultural crops in the United States and globally. 9 Ecological risk to sensitive aquatic biota (arthropods) from pyrethroid exposure has been predicted using highly protective approaches such as risk quotients (highest measured concentration/lowest toxicity value) and estimated environmental concentrations using models. 9 However, more realistic approaches using currently available toxicity and monitoring data suggest that aquatic life may not be impacted. 9 One water body in California's Central Coast that has been implicated as a high risk area for pyrethroids based on toxicity testing with a sensitive laboratory taxon (Hyalella) is the Santa Maria River watershed. 10 The California Central Coast Water Board adopted the Santa Maria TMDL for toxicity and pesticides on January 30, 2014. This TMDL includes various pyrethroids. Eight of the top 12 pesticides used for agricultural crops in the Santa Maria River watershed in 2012 were pyrethroids. Since pyrethroids are highly hydrophobic, concentrations of these insecticides in sediments are a potential ecological risk issue.
Since 2006, bioassessment multiple stressor case studies have been conducted in four primarily urban wadeable California streams using benthic macroinvertebrates with concurrent measurements of habitat metrics, metals and pyrethroids to determine which stressors are most important in influencing the condition of benthic communities. 11 The summary results from these four case studies showed that physical habitat metrics were the most important factors influencing benthic community condition while metal concentrations were the second most important factor influencing benthic communities in these streams. Pyrethroids were the least important constituents influencing benthic community conditions in these urban streams and for two of these streams pyrethroids were not a significant stressor to benthic communities when considered in a multiple stressor analysis. The bioassessment multiple stressor approach is ecologically relevant and provides a direct measure of stream health in concert with diagnostic analysis for identifying stressors, within a multiple stressor context, responsible for impairment. However, this approach has not been used for California water bodies located in agricultural areas where other investigators have documented pesticide occurrence and toxicity based on single species laboratory toxicity tests. 12 Therefore, this bioassessment multiple stressor approach was used in an intensive agricultural area of the Santa Maria River watershed in the Central Coast area of California to determine which stressors were most important in shaping benthic communities in this waterbody.
This 3-year study was designed to address the following specific objectives: (1) characterize benthic communities and physical habitat in Santa Maria River watershed during the spring of 2015, 2016, and 2017; (2) measure basic water quality parameters and nutrients in the water column and eight specific pyrethroids, total organic carbon (TOC), grain size, eight bulk metals and simultaneously extracted metals (SEM), and acid volatile sulfides (AVS) in sediment at various sites in concert with the bioassessments; and (3) use univariate and stepwise multiple regressions along with canonical correlation analysis to determine the relationship between various benthic metrics (ie, species richness, abundance) and physical habitat metrics (ie, embeddedness and channel flow status), pyrethroids, metals, nutrients and sediment characteristics from data sets based on sampling annually from 2015 to 2017.

Site selection
The Santa Maria River watershed is approximately 1.2 million acres composed of the large hydrologic areas containing the Cuyama Valley and the Sisque and Guadalupe watersheds. 13 The Guadalupe hydrologic area is located in the lower watershed and is transected by the Santa Maria River which flows east to west from the confluence of the Cuyama and Sisquoc Rivers to the Pacific Ocean. The lower Santa Maria River watershed contains year-round intensively cultivated lands that supports a $3.5 billion/year agricultural industry producing much of the nation's lettuce, artichokes, crucifer crops and strawberries. 14 Our 12 study sites were located in the lower Santa Maria watershed ( Figure 1). These 12 sites were sampled during the first week of March in 2015, 2016, and 2017. Water levels were generally higher in 2016 and 2017 than 2015 (a severe drought year). The lower watershed was adequately represented spatially as follows. Two sites were sampled in the mainstem Santa Maria River. Four sites were sampled in Solomon/Orcutt Creek. One site was sampled in both a channelized ditch at

Water and sediment measurements
The following water parameters were measured at each stream site using procedures described in Kazyak 15 : temperature, pH, salinity, specific conductivity, dissolved oxygen, and turbidity ( Table 1). All instruments were calibrated before and after sampling.
Grain size 16 and TOC 17 were measured on sediment samples collected from each site. Depositional areas-fine grain areas most likely to contain hydrophobic pesticides such as pyrethroids-were specifically sampled at each site when possible and three to five sediment samples from depositional areas were composited for the final sample. A stainless steel spoon (similar to a scoop) was used to collect the top 2-3 cm of sediment from each site. Approximately 1 L of sediment was collected from each site for grain size and TOC determinations as well as pyrethroids and metals. All sampling equipment was cleaned between sites using nitric acid, ethanol and distilled water. Sediment samples were stored in a cooler on ice in the field and later transferred to a refrigerator before shipment to Alpha Analytical Inc. for grain size and TOC analysis.

Physical habitat assessments
Physical habitat was evaluated at each site concurrently with benthic collections, water quality evaluations, nutrient analysis, sediment parameters, pyrethroids, and metals. The physical habitat evaluation methods followed protocols described in Harrington and Born. 18 The physical habitat metrics used for this study were based on nationally standardized protocols described in Barbour et al. 19 A total of 10 continuous metrics scored on a 0 to 20 scale were evaluated (Table 2). Other noncontinuous metrics including percent canopy, percent gradient, and substrate composition were also measured as reported in Hall et al. [20][21][22]

Bulk metals and SEM/AVS analysis
The following bulk metals with existing threshold effects levels (TELs), conservative protective benchmarks, as described by Buchman 23 were measured on composited sediment samples for each site as previously described using EPA method 6020 m: arsenic (As), cadmium (Cd), chromium (Cr), copper (Cu), lead (Pb), nickel (Ni), and zinc (Zn). Mercury (Hg) was also measured on all sediment samples using EPA method 245.7m. The reporting limit (RL) for these metals except Zn ranged from 0.02 to 4.5 μg/g dry weight. The RL for zinc ranged from 14.4 to 21.0 μg/g dry weight.
SEM analysis was conducted for Ni, Cu, Zn, Cd, Pb, and Hg using EPA method 200.8 m. The reporting limits (μmol/dry g) for these SEMs were as follows: Ni (0.02-0.04), Cu (0.001-0.002), Zn (0.02-0.04), Cd (0.001-0.002), Pb (0.01-0.02), and Hg (0.0001). AVS were evaluated on sediment samples from each site using procedures described by Plumb. 24 SEM/AVS ratios were then developed for each site to provide insight on the bioavailability of these metals in sediment. The principle of SEM/AVS is based on the observation that there are some components in sediment that bind certain metals such that they are no longer available and therefore not toxic to benthic organisms. 25,26 Sulfides in sediments have the ability to bind with divalent metals such as cadmium, copper, lead, mercury, nickel  Hall et al 5 and zinc and may render these metals unavailable to the extent sulfides are available. Sediments from the study sites were therefore analyzed for the amount of SEM and for the amount of freely available divalent metals as SEM. Assuming that sulfides would bind with metals on a 1:1 molar basis, dividing SEM by the amount of AVS would suggest that these metals are available when the SEM/AVS ratio is greater than 1.

Nutrients analysis
The following nutrient analysis was conducted on water column samples at each site: ammonium, nitrite, nitrate, dissolved nitrogen, total nitrogen, soluble reactive phosphate, and total phosphorus. [27][28][29] Pyrethroid analysis The pyrethroids bifenthrin, cypermethrin, cyfluthrin, deltamethrin, esfenvalerate, fenpropathrin, lambda-cyhalothrin and permethrin residues were extracted from sediment by shaking with methanol/water mixture and hexane for 1 hour. The sample was centrifuged and an aliquot of the upper hexane layer evaporated to dryness and re-dissolved in a small volume of hexane. The hexane sample was then subjected to a silica solid phase extraction (SPE) procedure prior to residue determination by gas chromatography with mass selective detection using negative ion chemical ionization (GC-MS/NICI). The limit of quantitation (LOQ) of the method was 0.14 to 0.29 ng/g dry weight for all pyrethroids. 30

Benthic macroinvertebrate sampling and identifications
Benthic macroinvertebrates were collected during the first week of March in 2015, 2016, and 2017 from three replicate samples at all 12 sample sites. The sampling procedures were conducted in accordance with methods described in Harrington and Born. 18 Within each of these sample reaches, a riffle was located (if possible) for the collection of benthic macroinvertebrates. A tape measure was placed along the riffle and potential sampling transects were located at each meter interval of the tape. Using a random numbers table, three transects were randomly selected for sampling from among those available within the riffle. Benthic samples were taken using a standard D-net with 0.5 mm mesh starting with the most downstream portion of the riffle. A 30.5 cm x 61 cm section of the riffle immediately upstream of the net was disturbed to a depth of 10.2 to 15.2 cm to dislodge benthic macroinvertebrates for collection. Large rocks and woody debris were scrubbed and leaves were examined to dislodge organisms clinging to these substrates. Within each of the randomly chosen transects, three replicate samples were collected to reflect the structure and complexity of the habitat within the transect. If habitat complexity was lacking, samples were taken near the side margins and thalweg (deepest path) of the transect and the procedures described above were followed. All samples were preserved in 95% ethanol. Due to the physical nature of this water body, it was often difficult to locate a substantial number of riffles to sample. Therefore, all sites were sampled using the non-riffle method. 18 This involved sampling the best available 30.5 cm x 61 cm sections of habitat throughout the reach using the same procedures described above. Best available habitat included general structure such as root wads, tires or other debris. Nine 30.5 cm x 61 cm sections were randomly selected for sampling (ie, stratified random sampling). Groups of three 30.5 cm x 61 cm sections were composited for each replicate for a total of three replicates per site.
All benthic samples were identified to the species level if possible. For taxa such as oligochaetes and chironomids, family and genus level, respectively, were often the lowest level of identification possible. Benthic macroinvertebrate subsampling (resulting in a maximum of 300 individuals) and identifications were conducted by California's Department of Fish and Wildlife (CDFW) in Rancho Cordova, California. The benthic macroinvertebrate samples were subsampled and sorted by personnel at the CDFW Laboratory located at Chico State University. Level 3 identifications (species level identifications) followed protocols outlined in Harrington and Born. 18 CDFW taxonomists conducted the taxonomic identifications. Slide preparations and mounting for species such as midges and oligochaetes followed protocols from the United States Geological Survey National Quality Control Laboratory described in Moulton et al. 31

Statistical analysis
A total of 899 statistical tests were conducted for the combined 2015, 2016, and 2017 data sets. 22 In preparation for statistical analyses, the data for the 10 key benthic metrics were averaged across the three transects sampled for each site in the Santa Maria watershed. These data were merged with data sets of sediment concentrations of pyrethroids, sediment concentrations of metals, simultaneously extractable metals to AVS ratios, sediment characteristics (% TOC, % sand, and % fine sediments), nutrients in the water, and habitat metrics for each site. The sediment concentration data for pyrethroids were converted to toxicity units (TUs) by standardizing them to 1% TOC and dividing by Hyalella LC50 values that were also standardized to 1% TOC. 32 Metal concentrations in sediment were also standardized to their relative toxicities by dividing the dry weight concentrations of each metal by their respective TEL values. 23 The statistical approach used was similar to that used for previous bioassessment/multiple stressor studies in California's urban streams. 11,20,21,33 The potential associations between the benthic metrics and pyrethroids, metals, sediment characteristics, nutrients, and habitat metrics were explored by a series of regression techniques. Prior to this analysis, all data were unit 6 Air, Soil and Water Research deviate standardized to place all dependent and independent variables on the same relative scales, as well as to produce more normalized distributions. Univariate general linear model regressions 34 were conducted to determine whether there were indications of significant relationships (α = 0.01) between benthic metrics and concentrations of pyrethroids (expressed as TUs), metals in sediments (expressed as metals to TELs), SEM to AVS ratios, nutrients, sediment characteristics, and habitat metrics. Then, a series of stepwise multiple regressions 34 were conducted to determine potential relationships between the benthic metrics and combinations of environmental variables. Stepwise regressions were conducted for the benthic metrics versus the toxicants (pyrethroids and metals); the nutrients, sediment characteristics, and habitat metrics; and all variables combined into the same model. A second series of stepwise regressions were conducted for the benthic metrics versus principal components (PCs) of the environmental data that were produced by principal components analyses (PCA) with an orthogonal rotation (Proc Factor, principal components method with a "varimax" rotation 34 ). The confirmation models that have been utilized in other studies 11 could not be utilized in this study because there were fewer samples than the total number of environmental variables. However, separate PCAs were conducted on the toxicant data set (pyrethroids and the metals), the habitat variables (habitat metrics and sediment characteristics), and the nutrients.
The major PCs for each of these analyses (ie, the PCs that explained at least 5% of the variance in the data and having eigenvectors greater than 1.0) and the variables that had the greatest loadings on them were identified. The stepwise regression analyses were conducted on the benthic metrics versus the combined sets of PCs from all of these PCAs. Since multicollinearity could exist between the three sets of PCs, the two confirmation models described previously could not be employed.
A supplementary series of multivariate analyses were conducted to provide additional insight into the relationships between sets of benthic metrics and sets of environmental variables. Stepwise regression analyses were conducted on PCs of benthic metrics versus the entire set of PCs of the toxicants, PCs of nutrients, and PCs of the habitat variables.
The final series of multivariate analyses involved canonical correlation analyses of the PCs of the benthic data versus the PCs of the environmental data. 34 The canonical correlation analyses were conducted on the PCs, because the number of samples collected was relatively small compared to the total number of benthic and environmental variables to be analyzed. Thus, there were insufficient degrees of freedom to permit canonical correlation analysis of the full suite of benthic and environmental variables.

Water measurements
The 3-year mean ranges of water measurements from the 12 Santa Maria sites in Table 1

Physical habitat metrics
The 3-year mean physical habitat metrics (maximum score of 20 per metric) and total score (maximum of 200) for the Santa Maria sites are presented in Table 2. The various metrics as well as the total score were highly variable among the 12 sites. For example, the mean total scores ranged from 16 to 88. The highest total habitat scores were reported at the downstream sites (SM1, SM2, and SM4; see Figure 1). The narrative scoring for total physical scores is as follows: optimal (150-200), suboptimal (100-150), marginal (50-100), and poor (0-49). 18 The mean total physical habitat score across all sites was 42 (ranked as poor). Nine of the 12 sites also had mean total habitat scores less than 49 reflecting poor habitat. The habitat scores reported in the Santa Maria River watershed were the lowest habitat scores we have reported from 11 different water bodies sampled in California since 2000. 22 Observed habitat conditions in wadeable waterbodies such as the Santa Maria River watershed are usually the result of the complex interplay between landscape alterations (ie, agricultural activities) and hydrogeomorphological factors. The physical habitat assessments conducted during this study were used to determine the suitability of the physical environment for aquatic biota such as benthic macroinvertebrates. Impaired physical habitat (including sediment loading) has been identified as a major stressor to aquatic life in California streams. 36,37 Altered physical habitat structure is also considered one of the major stressors of aquatic systems throughout the United States resulting in extinctions, local expirations and population reductions of aquatic fauna. 38 Identifying degraded physical habitat in streams such as the Santa Maria watershed is particularly critical for bioassessments as failure to do so can sometimes hinder investigations on the effects of toxic chemicals, as measured in this study, or other water quality related stressors as discussed later in this paper. Rankin 38 has also reported that there is a small but still significant risk of reporting a water quality related impact when one does not exist (false positive) when habitat assessments are insufficient or absent. Therefore, evaluation of physical habitat in agricultural streams in California's Central Coast is particularly important for multiple stressor analysis due to the intensive agriculture in the Santa Maria watershed.

TOC and grain size
The 3-year mean percent TOC ranged from 1.23% to 2.21% for the 12 Santa Maria sites as presented in Table 3. Percent sand was highly variable among sites ranging from 7.3 to 34.4. All sites were predominately fine grain areas, with silt ranging 7 from 34.5% to 57.8% and clay ranging from 15.2% to 55.6%. Depositional areas-where hydrophobic chemicals such as pyrethroids could accumulate-were specifically targeted for sampling to increase the probability of measuring these insecticides if they were present in the water body.
Sediment mapping studies designed to determine the percent depositional area in this waterbody have not been conducted, but sediment mapping studies have been conducted in another Central Coast stream complex (Salinas streams) that are somewhat similar to the Santa Maria watershed. The results from the sediment mapping studies in the Salinas streams showed that depositional areas were not dominant as only 24% of the bed sediment was determined to be depositional area. 39 Based on field crew observations, it is highly likely that depositional area (silt and clay areas) would not be dominant in the Santa Maria river watershed.

Bulk metals and SEM/AVS
The 3-year mean concentrations for 8 metals measured in the Santa Maria watershed sites are presented in Table 4. TEL exceedances were reported at 8 sites for arsenic and at all 12 sites for cadmium and nickel. Natural sources for these metals in local soil cannot be ruled out as a contributor to the potentially toxic concentrations reported. For example, the cadmium TEL exceedances at all sites are likely related to the presence of shale-based soil in the Santa Maria watershed that have historically high concentrations of cadmium. 40 Serpentine soils in the Santa Maria watershed area, a known source of nickel as reported by other investigators, may be a natural source for the high nickel concentrations reported in this watershed. 41,42 Possible anthropogenic sources of cadmium in the Santa Maria watershed are phosphorous-based fertilizers or micronutrient applications. 43 Anthropogenic sources of arsenic in this watershed may also be related to the use of phosphorous-based fertilizers and micronutrients and possibly diffuse sources such as atmospheric fallout. 43 The SEM/AVS data in Table 5 showed that 10 sites had ratios greater than 1 with at least one metal exceeding at TEL at all sites thus suggesting predicted potential metals toxicity to resident benthic communities at these sites. The upstream site (SM-12) in Osco Flaco Creek had the highest SEM/AVS ratio (15.9), while the two downstream sites in the mainstem of the Santa Maria River had the lowest SEM/AVS ratios (both ratios less than 1).  Table 6). Highest total nitrogen concentrations were reported at upstream site SM11 in Osco Flaco Creek (see Figure 1), while the highest total phosphorus was reported at upstream site SM7. Currently, the State of California is in the process of developing nutrient water quality criteria for the state so at this time there are no California nutrient criteria that can be compared with the nutrient measurements reported in the Santa Maria watershed.

Nutrients analysis
However, the relationship of nutrients to macroinvertebrate communities has been addressed in various studies as described below. Ashton et al 44 reported that taxonomic and functional metrics of an index of biotic integrity (IBI) were significantly associated with ammonia-N and nitrite-N in Maryland streams. Wang et al 45 reported that 18 of 26 macroinvertebrate measures were significantly correlated with at least one nutrient measure in Wisconsin streams. These investigators specifically reported that EPT individuals and taxa and mean tolerance value metrics had the strongest correlations with most nutrient measures. Maul et al 46 reported that nutrients such as total phosphorus and ammonia were important variables for structuring benthic communities in northwest Mississippi streams. Other investigators also reported that nutrients tended to increase plant-associated macroinvertebrates at low levels with a flattening response at intermediate levels and a decline at high levels. 47 Therefore, based on the above studies, there is certainly background information to suggest that nutrients could have significant relationships to benthic communities in agricultural waterbodies such as the Santa Maria watershed.

Pyrethroids
Ranges of 3-year mean concentration of 8 pyrethroids normalized to 1% TOC presented in Table 7 were as follows: bifenthrin (0.23-31.9 ng/g), cyfluthrin (0.06-3.5 ng/g), cypermethrin    (Table 8). 32 TU concentrations exceeding 1.0 were predicted to be potentially toxic. The sum of pyrethroid TUs using Hyalella indicated that 9 of the 12 sites were predicted to be toxic due to pyrethroids. The highest sum of pyrethroid TU (11.7) was reported at site SM11 in Oso Flaco creek (see Figure 1). The 3-year mean sum of pyrethroid TUs for all 12 Santa Maria site means was 2.7. This value is approximately 50% lower than the 3-year mean of another group of California Central Coast streams (three streams in the Salinas River watershed) where the mean of 13 Salinas stream mean site means of 6.1 was calculated. 48 Therefore, it appears that pyrethroid sediment concentrations are somewhat lower in the Santa Maria watershed when compared with data from other California Central Coast streams. The five most dominant benthic taxa collected during the 3-year study comprising 51% of the total taxa and the reported % of the total were as follows: Tubificidae unidentified immature (oligochaetes)-13.4%; Cricotopus (chiromonid)-10.7%; Physa (snails)-6.8%; Cricotopus bicinctus group (chironomid)-6.5%; and Enchytraeidae (oligochaetes)-6.4% (Table 9). These five taxa are generally considered tolerant to moderately tolerant of environmental stressors. 18 The 3-year mean values for the various benthic metrics were highly variable among the 12 Santa Maria sites. For example, taxa richness ranged from 12 to 23, percent tolerant taxa ranged from 41 to 71, percent collectorsgatherers ranged from 31 to 72, percent predators ranged from 3.3 to 18, and total abundance ranged from 232 to 6506 (Table  10).

Benthic communities
The state of California has recently developed a state-wide index that translates complex data about benthic macroinvertebrates into a measure of stream health. 49 This index is called the California Stream Condition Index (CSCI). The CSCI combines a multimetric index that measures ecological structure and function and an observed to expected (O/E) index that measures taxonomic completeness. Categories of biological condition based on 0-1 scoring are as follows: ⩽0.62 = very likely altered; 0.63-0.79 = likely altered; 0.80-0.91 = possibility altered; and ⩾0.92 likely intact. The 3-year mean value for the Santa Maria sites sampled was only 0.33 with values ranging from 0.19 to 0.53-very likely altered category-so the benthic communities in this watershed would be considered impaired.

Relationship of benthic metrics to all stressors
Interpretation of the relationships of benthic metrics to environmental stressors requires an understanding of how different benthic metrics are expected to respond to impairment. 18,50 In the statistical analyses presented below, the definition of the benthic metrics and expected response of the various benthic Air, Soil and Water Research metrics to impairment is as follows: taxonomic richness-total number of individual taxa (decrease); % dominant taxa-percent composition of the single most abundant taxa (increase); Shannon Diversity-an index that accounts for both abundance and evenness of species present (decrease); tolerance valuevalue between 0 and 10 weighted for abundance of individuals designated as pollution tolerant (higher values) and intolerant (lower values) (increase); % tolerant taxa-% of organisms in sample that are highly tolerant to impairment as indicated by a tolerance value of 8, 9, or 10 (increase); % collectors/filterers-% of macrobenthos that collect or filter fine particulate matter (increase); % collectors/gatherers-% of macrobenthos that collect and gather fine particulate matter (increase); % grazers-% of macrobenthos that graze on periphyton (variable); % predators-% predator individuals (decrease); and abundance-total number of individuals (decrease). The results of the stepwise analyses of benthic metrics versus all environmental variables from the combined 2015, 2016, and 2017 data set are shown in Table 11. Sediment TOC was the most frequently observed significant environmental variable that, in combination with other variables, displayed relationships with benthic metrics: displaying direct relationships with taxonomic richness, % collectors/gatherers, % grazers, and % predators; and displaying inverse relationships with tolerance value, and % collectors/filterers. Some of the benthic metrics that displayed direct relationships with TOC also displayed inverse relationships with various metals: taxonomic richness was inversely related to copper; % collectors/gatherers was inversely related to chromium; and % predators was inversely related to nickel. The % predators metric was also inversely related to Total TUs. Shannon's index, a diversity metric, was inversely related to lambda cyhalothrin, total phosphorous, and directly related to fines, while % dominant taxon, the index that is indicative of lack of diversity, displayed the opposite relationships with these environmental variables. The benthic metric % tolerant taxa, generally indicative of stressed communities when high, was inversely related to the habitat metric bend/riffle frequency, dissolved nitrogen, and directly related to ammonium. Most of the relationships between benthic metrics and the various environmental variables were individually quite weak (most R 2 values <0.20), but some combinations of the environmental variables explained more of the variations of some of the benthic metrics. Thus, with the exception of TOC, no single environmental variable or group of variables appeared to emerge from these analyses as being the dominant one in influencing benthic communities.
The results of the stepwise regression analyses of benthic metrics versus PCs of toxicants, nutrients, and habitat variables are presented in Table 11. While there were fewer significant relationships between benthic metrics and PCs than with combinations of individual environmental variables, several benthic metrics were significantly related to HabPC3, the PC positively loaded by TOC and fine sediments: tolerance value and % collectors/filterers were inversely related to HabPC3, and % collectors/gatherers was directly related to this PC. The metrics % collectors/gatherers and % collectors/filterers were also inversely related to ToxPC4, the PC positively loaded by lambda cyhalothrin, copper, and, to a lesser extent, by bifenthrin. The diversity metric Shannon's index was inversely related to NutPC2, the PC that was positively loaded by total phosphorous and dissolved phosphorous.
Most of the relationships shown to be significant in the stepwise regressions involving combinations of individual environmental variables were not supported by the results from the stepwise regressions of the PCs of the groups of environmental variables. Even those significant relationships that were observed between a few of the benthic metrics and environmental PCs did not consistently make sense ecologically for typical benthic communities. Tolerance value, % collectors/filterers, and % collectors/gatherers, three benthic metrics generally believed to be indicative of stressed communities did not display consistent relationships to HabPC3 (the PC indicative of fine, organic rich sediments), since the former two metrics displayed inverse relationships with it, while the latter metric displayed a direct relationship to it. Moreover, % collectors/filterers, and % collectors/gatherers both displayed inverse relationships to the ToxPC4. If the toxicants positively loaded by this PC (lambda cyhalothrin, copper and bifenthrin) had represented a stress to the benthic communities, one would have expected these benthic metrics to display a direct relationship.
In considering these results, a caveat should be noted. The limited number of samples for the combined 3-year data set prevented the typical confirmation analyses that has been used and previously as described in the statistical analysis section of this paper. Thus, multicollinearity may exist not only between the individual environmental variables, but also between the

Air, Soil and Water Research
PCs calculated for each group of variables separately. Thus, the significant, albeit rather weak, patterns that were observed may have been confounded by this multicollinearity (i.e. other correlated environmental variables, PCs or combinations may have been responsible for the relationships). The definitive confirmation analyses could not be conducted to determine whether significant patterns between benthic metrics and habitat variables persisted when the potential effects of toxicants were taken into account and vice versa. Thus, the results of the stepwise statistical analyses should not be over-interpreted.

Multivariate analysis of the PCs of benthic metrics versus PCs of environmental variables
The results of the stepwise regressions of PCs for the benthic metrics versus the PCs for all of the environmental variables for the combined 2015, 2016, and 2017 data sets are presented in Table 12. BenPC1 (the PC that was positively loaded by the benthic metrics % collectors/gatherers, taxonomic richness, % predators, and abundance, and negatively loaded by % collectors/filterers) was directly related to HabPC3 (the PC that was positively loaded by the sediment characteristics TOC and Fines) and was inversely related to ToxPC2 (the PC that was positively loaded by the metals nickel, cadmium, and chromium). Thus, the stepwise regressions for the PCs tended to confirm some of the patterns observed for the benthic metrics associated with BenPC1 in the analyses of individual benthic metrics with combinations of environmental variables shown in Table 11.
BenPC2 (the PC that was positively loaded by the diversity metric Shannon's index and was negatively loaded by % Dominant taxon) was inversely related to NutPC2 (the nutrient PC that was positively loaded by total phosphorous and dissolved phosphorous) and to ToxPC4 (the toxicant PC that was positively loaded by lambda cyhalothrin, copper, and, to a lesser extent, by bifenthrin). These relationships are similar to those observed for the relationships between the Shannon's index and % dominant taxon, and the environmental variables lambda cyhalothrin and phosphorus (see Table 11).
For the combined 2015, 2016, and 2017 data set, a canonical correlation was conducted on the benthic metric PCs versus the PCs for the environmental variables (ie, PCs for toxicants, nutrients, and habitat metrics). The results from the canonical correlation analysis are shown in Figure 2. There was a significant and moderately strong (R 2 = 0.65) relationship between the Canonical Variate for benthic PCs (CVBen; plotted on the y-axis) and the Canonical Variate I for PCs for environmental data (CVEnv; plotted on the x-axis). The CVBen for benthic metrics was positively correlated to BenPC1 and BenPC2, and displayed a somewhat weaker negative correlation to BenPC3. The benthic PCs that were positively loaded by benthic metrics associated with more diverse and healthier benthic communities tended to increase going up the y-axis, while those loaded on metrics associated with less diverse, more stressed communities tended to increase going down the y-axis. The CVEnv was directly correlated to HabPC3, the PC that was positively loaded by the sediment characteristics TOC and Table 11. Results of stepwise multiple linear regression models for the Santa Maria watershed in the combined 2015, 2016, and 2017 data set: (a) benthic metrics versus toxicity units for pyrethroids, habitat metrics, nutrients, sediment characteristics, and metals to TEL ratios; and (b) benthic metrics versus principal components of the environmental data; only variables that were significant at α = 0.01 were included in the models (NS = not significant). The direction of the relationship for each significant variable is indicated (+ = direct; -= inverse), as is the contributed R 2 value. Fines, and displayed a somewhat weaker negative correlation to NutPC2, the PC that was negatively loaded by total phosphorous and dissolved phosphorous. Samples higher on the x-axis tended to be from sites characterized by more organic rich, fine sediments and lower phosphorous concentrations.
In looking at the trend line relationship in Figure 2, data points found in the upper right hand quadrant of the figure tended to be from samples that were characterized by somewhat healthier, more diverse benthic communities, finer and more organic rich sediments, and lower phosphorous concentrations at these sites. Conversely, data points found in the lower left quadrant tended to be from more stress impaired, less diverse benthic communities with less organic rich sediments and greater phosphorous concentrations.
The direct relationship we have reported between TOC and diversity measures dominated by tolerant benthic taxa (oligochaetes) in the Santa Maria watershed has been addressed in other studies. 51,52 Grumiaux et al. 51 reported that polluoresistence benthic taxa abundance increased with an increase in organic content in rivers and canals in northern France. Since our sites in the Santa Maria watershed are dominated by tolerant benthic taxa, such as oligochaetes, the direct relationship we have reported with TOC and diversity measures for benthic assemblages dominated by tolerant taxa is supported by these authors. In another study conducted in the San Francisco estuary, Thompson and Lowe 52 reported that tolerant taxa such as oligochaetes increased with elevated concentrations of TOC. The results from the Thompson and Lowe 52 study also support our findings showing that TOC is an important factor shaping benthic communities dominated by tolerant taxa such as oligochaetes.
Overall, several observations may be made from the results of the canonical correlation analysis comparing the overall patterns of combinations of benthic metrics to combinations of environmental variables. First of all, the ToxPCs related to toxicants did not emerge as appearing to be significant to the BenPCs in the 3-year Santa Maria River data. The nitrogenbased nutrients also did not appear to be significant in the analyses of the 3-year data set. Rather, higher concentrations of phosphorous-based nutrients (NutPC2) appears to be somewhat important in characterizing the environmental conditions associated with more stressed, less diverse benthic communities (BenPC2 and BenPC3). The habitat metrics that characterize physical stream quality conditions did not appear to be significant in relationship to the benthic communities in the "big picture" canonical correlation analysis. This observation is quite different from the results of previous studies in California urban streams 11 in which some of the habitat metrics represented key environmental factors that were shown to be associated with the relative health of benthic communities as indicated by the benthic metrics. Table 12. Results of stepwise regression analyses of principal components for benthic metrics (BenPC a ) versus the principal components for environmental data (see footnote for Table 11 for loadings of the ToxPCs, HabPCs, and NutPCs) for the Santa Maria watershed in the combined 2015, 2016, and 2017 data set. Only variables that were significant at α = 0.01 were included in the models (NS = not significant). The direction of the relationship for each significant variable is indicated (+ = direct; -= inverse), as is the contributed R 2 values. BENTHIC Tables 11 and 12 for loadings of variables for PCs). The PCs for benthic biological metrics (y-axis) and the PCs for the environmental (x-axis) that were most highly correlated with the canonical variates are shown (r ⩾ |0.5|; moderate correlations at r ~ |0.5| are shown in parentheses), along with the direction of their relationships (shown by arrows). The regression trend line and 95% confidence limits are indicated by solid and dashed lines, respectively.

Air, Soil and Water Research
A few final caveats should be made regarding the multivariate analyses. While the stepwise analyses involving the PCs for various groups of environmental PCs (toxicants, nutrients, habitat metrics, and sediment characteristics) did address issues associated with multicollinearity between variables within the groups, there may still be multicollinearity between the PCs of the different groups of environmental variables. Thus, one PC could have been selected as being significant during the stepwise process, others that were correlated with it would have been discarded, even though they may be ecologically significant. The canonical correlation analysis tends to overcome this issue, since important environmental PCs identified as being correlated to the CVEnv by the model could be mutually related to the CVBen and the BenPCs to which it was correlated. Therefore, the relatively modest number of samples for even the 3-year data set prevents taking the ideal statistical approach of using PCs for the entire set of environmental variables (in which case, the PCs would not be correlated to each other) or even the use of the entire individual environmental variables in the canonical analyses. 11 Furthermore, the relatively small data set may affect the power of the multivariate models to detect more subtle relationships.

Results overview of the statistical analysis
There were relatively few significant relationships between benthic metrics and toxicants in the sediments in the stepwise analyses of the data from the Santa Maria watershed. Most of the statistically significant relationships that were observed between benthic metrics and pyrethroids were largely influenced by an outlier sample (Site SM11) observed in 2017 that had very high concentrations of bifenthrin and to a somewhat lower lambda-cyhalothrin TU, although still greater than 1. However, these significant relationships disappeared with these two pyrethroids in the stepwise analyses for combined 2015, 2016, and 2017 data set analyses, suggesting that the impact of the outlier data point was no longer a factor in the analyses of a somewhat larger data set.
There were also few relationships between the benthic metrics and metals. For the analyses of the combined 2015, 2016, and 2017 data set, the benthic metric % predators was inversely related to nickel. Nickel was quite high in the sediments with the TEL value being exceeded at all sites (Table 4). While the relationships between this benthic metric and nickel made sense ecologically, the strength of this relationship only explained a low to moderate proportion of the variance for this metric (R 2 = 0.21). However, none of the benthic metrics were significantly related to SEM to AVS ratios for any of the metals, including nickel (ie, the benthic metrics were not significantly related to the forms of the metals generally believed to be most biologically available), so it is difficult to say that nickel in sediments represent a major factor influencing the benthic communities of Santa Maria watershed.
There were more apparent relationships between benthic metrics and nutrients. The benthic metric % tolerant taxa displayed a direct relationship with ammonium in the combined 2015-2017 data set. This relationship persisted in the stepwise analyses of all environmental data, as well as those involving PC loaded by nitrogen nutrients. In contrast, the phosphorousbased nutrients did not display significant relationships in the stepwise analyses involving only the nutrients, but they were significant in relationship to diversity-related metrics in combination with other environmental variables in the stepwise analyses for all variables and in the multivariate analyses: Shannon's index was inversely related to total phosphorous, while % dominant taxon was directly related to it. These relationships persisted in the multivariate analyses involving benthic and environmental PCs.
In contrast to the results of our previous bioassessment multiple stressor studies, the benthic metrics did not have many significant relationships with traditional habitat metrics in the Santa Maria data set. 11 Tolerance taxa was inversely related to Bend/Riffle frequency in the stepwise analyses of the combined data set, but it did not display a significant role in the big picture multivariate analyses (ie, the regressions involving the PC of the habitat metrics or the canonical correlation analysis).
TOC in sediment displayed a significant relationship with benthic metrics in the stepwise regression analyses and the multivariate analyses. The metrics tolerance value and % collectors/filterers were inversely related to TOC in sediments, while % grazers was directly related to it. These benthic metrics displayed similar relationships to HabPC3, the PC that was positively loaded by TOC in the sediment and fine sediments.
The summary multivariate canonical correlation analysis indicated that less stressed, more diverse benthic communities tended to be associated more with TOC-rich finer sediments and lower concentrations of phosphorous-based nutrients, and more stressed, less diverse benthic communities tended to be associated with less organically rich, somewhat less fine sediments and higher phosphorous concentrations. Neither toxicants nor habitat metrics were shown to be important factors in this analysis. The extremely poor habitat metric scores observed consistently at all sites throughout the Santa Maria watershed may be responsible for this group of variables not being a major factor that would shape spatiotemporal patterns of the benthic communities throughout the waterbody.

Conclusion
Resident benthic communities in the Santa Maria River watershed were reported to be impaired based on a benthic index developed by the State of California and the dominant presence of tolerant taxa. A suite of possible constituents that could be responsible for the benthic community impairment, such as pyrethroids, metals, nutrients, sediment characteristics, and altered physical habitat conditions, were analyzed concurrently in order to identify which stressor or stressors were most influential. Predictive approaches based on laboratory-derived toxicity data for toxicants (pyrethroids and metals) suggested that these toxicants could be toxic to resident benthic communities. However, a thorough statistical analysis of all stressors based on field data that yields an observed response did not support this prediction. The consistently poor habitat conditions at all sites in the Santa Maria waterbody was not identified in the statistical analysis to be a major contributor to the impaired benthic communities possibly due to a lack of variation in habitat conditions among the sites. This finding is in sharp contrast to other studies conducted in urban streams in California using a similar sampling design where habitat conditions were determined to be a key stressor.
The key constituents identified as major contributors in shaping benthic communities were TOC and phosphorusbased nutrients. In summary, the less stressed, more diverse benthic communities highly populated by tolerant taxa such as oligochaetes tended to be associated more with TOC-rich finer sediments and lower concentrations of phosphorousbased nutrients while the more stressed, less diverse benthic communities tended to be associated with less organically rich, somewhat less fine sediments and higher phosphorous concentrations. Field approaches conducted over multiple years in a waterbody designed to determine the relationships of multiple stressors to resident benthic communities provide the most ecologically relevant assessment of which stressors are most important in shaping benthic communities. Relying solely on single species laboratory toxicity data to predict the impact of specific toxicants on resident benthic communities may provide misleading conclusions if confirmation from a bioassessment approach addressing possible multiple stressor impacts is not determined.