Blocksom_Johnson_2009_Ecological-Indicators_withErratum_0001.pdf

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/223411069

Development of a Regional Macroinvertebrate Index for Large River

Bioassessment

Article in Ecological Indicators · March 2009

DOI: 10.1016/j.ecolind.2008.05.005

CITATIONS

READS

2 authors:

Some of the authors of this publication are also working on these related projects:

Examining impacts of unconventional oil and gas infrastructure development on stream communities View project

Ohio River Biological Monitoring and Assessment View project

Karen A Blocksom

United States Environmental Protection Agency

53 PUBLICATIONS 1,554 CITATIONS

SEE PROFILE

Brent R Johnson

United States Environmental Protection Agency

32 PUBLICATIONS 1,079 CITATIONS

SEE PROFILE

All content following this page was uploaded by Karen A Blocksom on 17 March 2020.

The user has requested enhancement of the downloaded file.

https://www.researchgate.net/publication/223411069_Development_of_a_Regional_Macroinvertebrate_Index_for_Large_River_Bioassessment?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_2&_esc=publicationCoverPdf https://www.researchgate.net/publication/223411069_Development_of_a_Regional_Macroinvertebrate_Index_for_Large_River_Bioassessment?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_3&_esc=publicationCoverPdf https://www.researchgate.net/project/Examining-impacts-of-unconventional-oil-and-gas-infrastructure-development-on-stream-communities?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_9&_esc=publicationCoverPdf https://www.researchgate.net/project/Ohio-River-Biological-Monitoring-and-Assessment?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_9&_esc=publicationCoverPdf https://www.researchgate.net/?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_1&_esc=publicationCoverPdf https://www.researchgate.net/profile/Karen-Blocksom?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_4&_esc=publicationCoverPdf https://www.researchgate.net/profile/Karen-Blocksom?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_5&_esc=publicationCoverPdf https://www.researchgate.net/institution/United_States_Environmental_Protection_Agency?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_6&_esc=publicationCoverPdf https://www.researchgate.net/profile/Karen-Blocksom?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_7&_esc=publicationCoverPdf https://www.researchgate.net/profile/Brent-Johnson-5?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_4&_esc=publicationCoverPdf https://www.researchgate.net/profile/Brent-Johnson-5?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_5&_esc=publicationCoverPdf https://www.researchgate.net/institution/United_States_Environmental_Protection_Agency?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_6&_esc=publicationCoverPdf https://www.researchgate.net/profile/Brent-Johnson-5?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_7&_esc=publicationCoverPdf https://www.researchgate.net/profile/Karen-Blocksom?enrichId=rgreq-2854abcf4b16bf9d32877023039d56d0-XXX&enrichSource=Y292ZXJQYWdlOzIyMzQxMTA2OTtBUzo4NzAwNzg2Nzc1MzY3NjlAMTU4NDQ1NDMzMDE3MQ%3D%3D&el=1_x_10&_esc=publicationCoverPdf e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8 Development of a regional macroinvertebrate index for large river bioassessment

K.A. Blocksom *, B.R. Johnson

U.S. Environmental Protection Agency, National Exposure Research Laboratory, 26 West Martin Luther King Drive, Cincinnati, OH 45268, United States a r t i c l e i n f o

Article history:

Received 7 December 2007

Received in revised form

6 May 2008

Accepted 8 May 2008

Keywords:

Non-wadeable

Upper Midwest

Least disturbed condition

Biological indicator a b s t r a c t

Large river bioassessment protocols lag far behind those of wadeable streams and often rely on fish assemblages of individual rivers. We developed a regional macroinvertebrate index and assessed relative condition of six large river tributaries to the upper Mississippi and

Ohio rivers, Midwest USA. In 2004 and 2005, benthic macroinvertebrates, water chemistry, and habitat data were collected from randomly selected sites on each of the St. Croix, Wisconsin, Minnesota, Scioto, Wabash, and Illinois rivers. We first identified the human disturbance gradient using principal components analysis (PCA) of abiotic data. From the

PCA, least disturbed sites showed strong separation from stressed sites along a gradient contrasting high water clarity, canopy cover, habitat scores, and plant-based substrates at one end and higher conductivity and nutrient concentrations at the other. Evaluation of 97 benthic metrics identified those with good range, responsiveness, and relative scope of impairment, as well as redundancies with other metrics. The final index was composed of

Diptera taxa richness, EPT taxa richness, Coleoptera taxa richness, percent oligochaete and leech taxa, percent collector-filterer individuals, predator taxa richness, percent burrower taxa, tolerant taxa richness, and percent facultative individuals. Each of the selected metrics was scored using upper and lower thresholds based on all sites, and averaging across the nine metric scores, we obtained the Non-wadeable Macroinvertebrate Assemblage Condi-tion Index (NMACI). The NMACI showed a strong response to disturbance using a validation data set and was highly correlated with non-metric multidimensional scaling (NMDS) ordination axes of benthic taxa. The cumulative distribution function of index scores for each river showed qualitative differences in condition among rivers. NMACI scores were highest for the federally protected St. Croix River and lowest for the Illinois River. Other rivers were intermediate and generally reflected the mixture of land use types within individual basins. Use of regional reference sites, though setting a high level of expectation, provides a valuable frame of reference for the potential of large river benthic communities that will aid management and restoration efforts.

Published by Elsevier Ltd.

avai lab le at www.sc iencedi rec t .com journal homepage: www.e lsev ier .com/ locate /eco l ind

1. Introduction

Biological monitoring is well-developed for wadeable streams of the upper Midwest, and many states have developed

* Corresponding author. Tel.: +1 513 569 7139.

E-mail address: blocksom.karen@epa.gov (K.A. Blocksom).

1470-160X/$ – see front matter. Published by Elsevier Ltd.

doi:10.1016/j.ecolind.2008.05.005 indicators for one or more assemblages (i.e., algae, macro-invertebrates, fish) (USEPA, 2002). In some states, indicators are available for non-wadeable sites as well (e.g., Ohio EPA, 1987; IDEM, 2002; Niemela and Feist, 2002), but many are based mailto:blocksom.karen@epa.gov http://dx.doi.org/10.1016/j.ecolind.2008.05.005 e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8314 only on fish assemblages and are only applicable to individual rivers (e.g., Emery et al., 2003). The use of at least two assemblages has been suggested for more robust biological assessment of condition, as each assemblage may respond differently to potential stressors (Yoder and Rankin, 1995).

However, few states in this region currently have multimetric indices specifically designed for macroinvertebrate assem-blages in large rivers. Even in states where multimetric indices are planned, there are relatively few river miles available for index development. A regional approach to developing a biological condition index allows for representation of a wider range of conditions than might be found within a single state.

Multi-state indicators for fish, algae, or macroinvertebrate assemblage condition have been developed for wadeable streams in some regions of the U.S. (Angermeier et al., 2000;

Hilletal.,2000;McCormicketal.,2001;Klemmetal.,2003;USEPA, 2006; Whittier et al., 2007). These studies on wadeable streams involved randomly sampling a large number of sites across the regionandassumedawiderangeofconditionswererepresented among the sites. To support indicator development, additional sitesweresometimessampledtoaugmentsitesfallingatthetwo extremes of condition (i.e., very high or very low quality).

However, assessment of non-wadeable rivers presents issues that may notapply to wadeable streams.Non-wadeablerivers in the U.S. represent only a small fraction of the total number of miles of flowing water, yet a single river may represent a large proportion of non-wadeable river miles in a given state. Thus, in order to develop a regional indicator for non-wadeable rivers, multiplesitesalongagivenriverwouldlikelyneedtobeincluded in any random sampling design. Because states may want to reportontheconditionofaspecificriver,manysitesalongariver would be required to achieve adequate precision of estimates

(based on proportions of river miles in a particular condition)

(Agresti, 1996, pp. 8–12). Although some types of anthropogenic disturbance may affect biological condition on a limited, local scale, we recognize that exposures to stressors are often cumulative,sothat theoverall levelofdisturbance inaparticular section of river is the sum of local-scale and upstream disturbances. Thus, spatial autocorrelation of both disturbance and condition is expected among sites within a river. However, because sites were randomly selected in each river and were not equidistant longitudinally, we assumed independence among samples from an individual river for statistical purposes. For the purposes of developing a biological indicator, capturing a wide range of conditions in a single river is difficult or impossible because conditions within a river may be relatively homo-geneous. For this study, our goals were two-fold: to both develop a macroinvertebrate indicator and to report on the condition of multiple rivers within a region. We therefore identified rivers of varying predominant land use, and randomly selected a large number of sites (i.e., 25 sites) along each river in an attempt to capture a wide a range of conditions.

The rivers sampled for this study were all tributaries to the upper Mississippi and Ohio rivers, but they do not include all large river tributaries of the states included (Illinois, Indiana, Minnesota, Ohio, and Wisconsin). We selected rivers for this study based on land use/land cover in an attempt to include a gradient of human disturbance, and we sampled enough sites on each river to estimate condition for individual rivers.

Inherent in this site selection process is the assumption that all of the rivers in this region have similar natural character-istics. However, because the level of human disturbance might be affected by natural gradients, we also evaluated the relationships between physical site variables representing natural variation and disturbance.

We followed a straightforward process for evaluating metrics to develop a multimetric index for macroinverte-brates, although we relied on designation of reference (least disturbed) conditions using an ordination of sites based on abiotic variables related to disturbance. After developing the index, we performed ordination using macroinvertebrate data to verify the relationship of the index with the stressor gradient and determine how well it reflected variation in benthic macroinvertebrate data. This process reflects the primary goal of this research, to provide an inclusive starting point for measuring condition in non-wadeable rivers in the upper Mississippi and Ohio river drainages.

2. Materials and methods

2.1. Study area and site selection

We selected six large river tributaries to the upper Mississippi and Ohio rivers in U.S. Environmental Protection Agency

(USEPA) Region 5 based partially on basin level land use data, evaluated using National Land Cover Data (NLCD) (Homer et al., 2004). Rivers were also targeted to cover a large range of basin size and locations within the region. The overall goal in selecting specific rivers was to include in the study as much of the gradient of conditions as possible. We ultimately selected three northern rivers (the St. Croix, Wisconsin, and Minnesota rivers) and three southern rivers (the Illinois, Wabash, and

Scioto rivers) for inclusion in the study. On each river, we randomly selected 25 sites from the entire length of river having a watershed area larger than approximately 2590 km2

(1000 mi2) (Fig. 1). An oversample site list of>200% (60 sites per river) was provided to account for sites in the original list that could not be sampled due to flooding, lack of site access, or safety hazards. Only when a site on the original list of 25 could not be sampled for the reasons listed was the next site on the oversample list sampled. Only eight sites (32% of targeted total) were sampled on the Minnesota River in both 2004 and

2005 because flooding in both years prevented collection of representative benthic samples.

As intended, the six rivers sampled reflect a variety of land uses and physiographic regions (Table 1). The Illinois River includes the Chicago metropolitan area in its upper watershed and passes through land dominated by row crop and pasture before entering the upper Mississippi near Alton, IL. Similarly, the Minnesota River runs through predominantly row crop and pasture land in southwest Minnesota, although it also passes through a few larger municipalities before emptying into the Mississippi River. The Wisconsin River runs south-ward through a mixture of forest and agricultural land, including the Wisconsin Dells, to the Mississippi River. Forest predominates in the upper part of the St. Croix watershed, with some residential and agricultural land interspersed, and agricultural land use increases in the lower part of the basin before the river meets the Mississippi. The St. Croix is

Fig. 1 – Sites sampled in each of the six study river basins.

e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8 315 designated as a National Scenic River. The Scioto River drains agricultural land in its uppermost reaches, passes through

Columbus, Ohio, and then runs through agricultural land or is flanked by row crops amid forest before emptying into the

Ohio River. The land use in the Wabash River is primarily cropland and pasture, with pockets of residential land cover and some forestland along the river corridor before it reaches its confluence with the Ohio River in the south (WRHCC, 2004).

2.2. Field sampling methods

Sampling was conducted from June through mid-October during 2004 and 2005. Crews followed the macroinvertebrate sampling method, as described in Flotemersch et al. (2006), with some modifications. Briefly, they sampled within a 500 m reach in which both banks were sampled over 6 evenly spaced transects, for a total of 12 sampling locations throughout the reach. At each transect, crews collected macroinvertebrates within a 10-m sampling zone (5 m on each side of transect) that extended to the midpoint of the river or until depth exceeded 1 m. Within each sampling zone, they performed six sweeps, each 0.5 m in length, using a 0.3 m wide D-frame net (500 mm mesh). The six sweeps were divided proportionally based on available habitats (e.g., snags, macrophytes, cobble) and were collected by kicking or agitating available substrates by hand. If water at a site was>1 m deep at water’s edge, six sweeps were collected along the bank from theboat. All samples from a reach were composited into a single sample, and the total sampling area was approximately 10.8 m2 (0.15 m2 per sweep � 6 sweeps per sampling zone � 12 sampling zones). The use of six standardized sweeps represents a modification from Flote-mersch et al. (2006), which includes two kicks and a timed multiple habitat sample. The use of only a 500 mm mesh D-frame net rather than a rectangular net for kicks and a D-frame net for multiple habitat sampling, both with 595 mm mesh also reflects a modification. All of these changes were designed to increase standardization in the field and reflect equipment commonly used by state agencies.

To characterize the abiotic condition of each site, crews collected habitat and water chemistry data. The habitat

Table 1 – River basin characteristics for six large river tributaries to the Ohio and upper Mississippi rivers (modified from Delong, 2005 and White et al., 2005)

Parameter River

St. Croix Wisconsin Minnesota Scioto Wabash Illinois

Relief (m) 319 300 85 307 275 45

Basin area (km2) 20,018 30,000 27,030 16,882 85,340 75,136

Mean discharge (m3/s) 131 261 125 189 1001 649

Land use (%)

Forest 65 54 1 21 22 7

Agriculture 27 27 95 69 65 87

Urban 2 3 1 9 5 5

River order 6 8 7 6 7 9

Annual precipitation (cm) 78 85 66 98 96 92

Water temperature (8C) 10.4 11.8 10.5 15 15 10.4

Physiographic provinces Superior Upland, Central Lowland

Central Lowland, Superior Upland

Central

Lowland

Central Lowland, Appalachian

Plateau, Interior

Low Plateau

Central Lowland, Interior Low Plateau

Central

Lowland

Major mainstem dams 2 16 6 1 1 5

Population density

(people/km2)

13 46 6 108 62.7 97 e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8316 protocol was a combination of a modified subset of USEPA’s

Environmental Monitoring and Assessment Program for Sur-face Waters (EMAP-SW) physical habitat protocols for non-wadeable streams (Kaufmann, 2000) and modification of a habitat quality index created for non-wadeable rivers in

Michigan (Wilhelm et al., 2005). Quantitative habitat variables relevant to indicator development included depth and width measurements, distance to riparian vegetation, water velo-city, canopy cover on banks, and substrate composition in the sampling zone. The habitat quality index included qualitative assessments of riparian width, large woody debris quantity

(modified from Wilhelm et al., 2005), aquatic vegetation, bottom deposition, bank stability, thalweg substrate, and quantity of off-channel habitats.

At each site, crews collected physico-chemical data and water samples for nutrient analysis at the most downstream transect in the sampling reach. We measured conductivity (mS/ cm), dissolved oxygen (D.O.) (mg/l), pH, and water temperature

(8C) in situ using a Hydrolab1 meter. We measured water clarity using a Secchi disk. Finally, a grab sample of approximately

250 ml of stream water was collected and preserved with concentrated sulfuric acid for nutrient analysis.

2.3. Laboratory processing

For all macroinvertebrate samples, organisms were sorted from debris using a stereo dissecting microscope (10� magnification). A grid square was randomly selected from a gridded pan and all organisms sorted from that square.

Subsequent grids were randomly selected and completely sorted until the total number of organisms was within 10% of the fixed count size. The fixed-count size was 300 organisms for most samples, but for some that were also part of a method comparison component of this study, a 500-organism sub-sample was used. All organisms were identified to the lowest practical taxonomic level, given availability of keys and condition of specimens. In addition to a fixed-count sub-sample, a large-rare sort was carried out on each sample to supplement estimates of taxonomic richness for each site with large and rare organisms that were not sorted as part of the random subsample. This two-phase processing of benthic samples has been recommended because rare taxa can be valuable for bioassessments (Cuffney et al., 1993; Courte-manch, 1996; Vinson and Hawkins, 1996).

Preserved water chemistry samples were kept on ice until transport to the laboratory facility, where they were kept refrigerated until analysis within 28 days of arrival. Subsamples for total Kjeldahl nitrogen (TKN) and total phosphorus (TP) were filtered prior to analysis. Water chemistry samples were analyzed for TP, TKN, nitrate/nitrite (NOX), and ammonia

(NH3) concentrations using segmented flow analysis.

2.4. Data analysis

2.4.1. Defining the disturbance gradient

Thefirststepindevelopingthe indexwastodefineadisturbance gradient across sites within our data set. We recognized that by the nature of the probability design we used for selecting sampling sites, itwasunlikely that wewould havesampled sites that might be considered in minimally disturbed condition

(Stoddard et al., 2006), if indeed these even exist currently for rivers in the Midwest. We therefore identified sites in least disturbed condition (LDC), as defined in Stoddard et al. (2006).

However, we were unable to delineate specific water chemistry and habitat criteria for LDC due to limited knowledge of what either historical conditions or minimally disturbed conditions were like in these rivers. Instead, we used an approach followed by others in which a disturbance gradient is characterized through ordination on abiotic variables (Zampella and Bunnell, 1998; Ferreira et al., 2005; Bressler et al., 2006), and LDC sites are identified from the ordination (Ferreira et al., 2005). We performed principal components analysis (PCA) on a set of water quality and habitat variables collected at each site that were intended to represent human disturbance. The assump-

Table 2 – Abiotic variables used in principal components analysis (PCA), along with transformations to reduce skewness, distributional characteristics (minimum, maximum, median), and Pearson correlations (N = 133) with PCA axis 1

Parameter code Description Transform Min, max value

Median value

PC 1

(l = 5.43)

DISTRV Distance to riparian vegetation (m) log10(x + 0.1) 0.00, 75.30 1.44 0.45

CANOP Canopy cover (%) 0.04, 99.17 61.90 �0.60

P_MED Cobble + coarse gravel habitat coverage in sampling zone (%)

0.0, 79.5 16.9 �0.48

P_SM Sand + fine gravel habitat cover (%) 0.0, 98.2 39.4 �0.46

P_SCM Silt + clay + muck habitat cover (%) 0.0, 90.4 26.7 0.37

P_PLANT Rootwads + undercut banks + woody debris + leaf packs + macrophytes cover (%)

0.0, 81.1 9.6 �0.84

HABSUM Sum of qualitative habitat scores

(Wilhelm et al., 2005)

18, 83 57 �0.61

COND Specific conductivity (mS/cm) 91, 1369 565 0.89

SECCHI Secchi depth (m) 0.23, 1.45 0.49 �0.79

TOTP Total phosphorus (mg/L) log10(x) 0.009, 1.017 0.093 0.78

TKN Total Kjeldahl nitrogen (mg/L) 0.262, 2.459 0.792 0.67

NH3 Ammonia (mg/L) log10(x) 0.0025, 0.5880 0.0620 0.68

NOX Nitrate + nitrite (mg/L) log10(x) 0.005, 4.408 0.580 0.51

The eigenvalue of the axis is represented by l.

e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8 317 tions of multivariate normality are relaxed in cases where the

PCA is used for exploratory purposes rather than hypothesis testing, and reducing or eliminating skewness (i.e., to a range of

�1to1) isofprimaryconcerninevaluatingvariabledistributions

(McCune and Grace, 2002, pp. 114–121). Table 2 provides a list of abiotic variables used in the PCA and thetransformation used, if any. Values flagged as below the detection limit (DL) were set at the DL before transformation or use in analysis. Sixteen variables chosen to reflect human influence were included in the PCA, including nutrient concentrations, habitat measures, a qualitative habitat index score (Wilhelm et al., 2005), and in situ water quality measurements.

The PCA was run using a correlation matrix because the variables included were on varying scales with varying units.

Only axes with an eigenvalue larger than the corresponding broken-stick eigenvalue were examined in greater detail. The broken-stick eigenvalue for the ith axis is calculated as:

lbroken-stick;i ¼ Xn i¼k li where li is the eigenvalue for the ith axis, k the axis of interest, andn is the maximumnumberofaxes (Jackson, 1993).A broken-stick eigenvalue represents the expected eigenvalue if the var-iation among components was distributed randomly and is a useful way to determine the axes worth interpreting. In exam-ining PCA axes, we evaluated the direction and strength of

Pearson correlations of component variables with axes to deter-mine the end of a given gradient reflecting more desirable conditions. We then took the 25th and 75th percentiles of the axis to designate LDC (henceforth termed reference) and stressed sites based on the site score along the PCA axis. We use the reference and stressed sites primarily to evaluate metrics with respect to responsiveness and relative reference site variability.

As a verification of the disturbance gradient, we determined land use for each site in a riparian corridor consisting of a 500 m buffer on each bank for a distance of 4 km upstream of the site.

The relationshipbetween the disturbancegradientand landuse was evaluated with Spearman rank correlations between the

PCA axis and percent forest, percent agriculture, and percent development in the riparian corridor. We compared land use and the constructed disturbance gradient to physical char-acteristics of sites (i.e., wetted width, elevation, thalweg depth, flow velocity, bankfull height, latitude) using Spearman rank correlations to evaluate the relationships between natural variation and human disturbance. Although an adjustment for longitudinal position (river mile) within a river may account for some of this natural variation within a river, this approach is most appropriately applied using only reference, or least disturbed, sites to develop the relationship. In this case, this is not possible because some rivers may not have any sites representing reference condition. Adjusting for river mile using all sites may inadvertently remove variation related to human disturbance, rendering a condition index much less useful.

However, we did examine plots of river mile with the disturbance gradient within each river to look for strong trends.

2.4.2. Metric evaluation and index development

Although all organisms were identified to the lowest practical taxonomic level, most states do not identify beyond the genus level. In an effort to make the macroinvertebrate indicator more useful to states and other potential users, we aggregated data to genus level for most taxa and family level for annelids before calculating metrics. Then we calculated a total of 97 metrics based on richness, percent of taxa, and percent of individuals within taxonomic groups, functional feeding groups (i.e., collector-gatherers, collector-filterers, scrapers, shredders, scavengers, and predators), habit preferences and behavior groups (i.e., swimmers, clingers, burrowers, spraw-lers, climbers, and skaters), and pollution tolerance (i.e., intolerant, facultative, tolerant, and the Hilsenhoff Biotic

Index (Hilsenhoff, 1987)). Designations of tolerance were based on those developed for use in Klemm et al. (2003) with the

EMAP wadeable streams database. For the purposes of calculating richness metrics, rarefaction (Hurlbert, 1971) to

300 organisms was performed on all samples to account for e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8318 differences in the subsample size for some samples. Organ-isms from the large-rare sorts were included only in calcula-tions for taxa richness metrics.

Before evaluating metrics, we randomly divided the data into calibration and validation data sets. To ensure an adequate number of reference and stressed sites in each data set, we divided the datawithineachclasssothat twothirdswereplaced into the calibration data set and one third into the validation data set. The similarity of distributions across the PCA axis was compared between calibration and validation data using a

Kolmogorov–Smirnov (K–S) test to ensure that a similar range of conditions was represented in both data sets. Only calibration data were used to evaluate metrics and select a subset for a multimetric index. We then used the validation data to verify the responsiveness and precision of the final index.

Weevaluatedmetrics with respect to several characteristics:

range, responsiveness, relative scope of impairment, and redundancy with other metrics. Wefollowed a stepwise process to reduce the number of candidate metrics remaining after each stage of evaluation. For range, we eliminated percentage metrics with a range of 10% or less, richness metrics with a range of 5 or less, and metrics with 50% or more of sites having the same value. For responsiveness, we examined box plots of reference and stressed sites and selected only those showing moderate to strong separation (i.e., neither median overlapping with the interquartile range (IQR), or 25th to 75th percentile range, of the other group, up tono overlap of IQRs at all). Relative scope of impairment (SOI) is calculated as the ratio of the range of possible impairment (values beyond the poorer quality 25th percentile) and the IQR of reference sites, a measure of variability. Relative SOI values larger than 1 are more desirable, with values less than 1 indicating too much variability among reference sites compared to the range of impairment. For simplicity, we ran Spearman rank correlations among all remaining metrics and flagged those pairs of metrics having a correlation coefficient with a magnitude of 0.75 or greater as redundant. After these steps, we grouped the remaining set of metrics by the type of metric (i.e., richness, composition, tolerance, functional feeding group, and habit). We selected metrics for the final index based largely on relative respon-siveness and variability across reference sites. In addition, particular effort was made to incorporate as many features of the assemblage as possible by including two of each type of metric while avoiding metric redundancy.

To create a multimetric index, we had to standardize or score metrics to put them all on the same scale. We encountered metrics that were related to disturbance in the positive and negative directions. We calculated continuous scores as (value � floor)/(ceiling � floor) � 100 for positive metrics, and as (ceiling � value)/(ceiling � floor) � 100 for negative metrics. The ceiling and floor were calculated as the 95th and 5th percentiles, respectively, across the entire distribution of sites. Because a probability design was used to collect the data, the entire range of conditions is assumed to be captured in our sample. In this case, use of the full distribution of sites to set thresholds for scoring makes more sense than using reference and stressed sites to set expectations for each metric because there is a relatively small number of reference sites available (Blocksom, 2003). We averaged all the metric scores to obtain the final index score.

2.4.3. Validation

To validate the index, we calculated index scores using the validation data set. Then we compared reference and stressed sites based on these scores. We also examined the correlation between the index and the PCA axis, as well as the interquartile range among reference sites. All of these tests were compared with results obtained using the calibration data set. We also compared the NMACI score within each river to longitudinal position using Spearman rank correlations and visual inspections of plots to assess the dependence of scores on river mile, based on all sites in each river.

As a secondary evaluation of the index, we performed unconstrained ordination on macroinvertebrate abundance data using non-metric multidimensional scaling (NMDS) (ter

Braak, 1995; McCune and Grace, 2002, pp. 125–142). This allowed us to summarize variability across a large number of taxa into a few dimensions and to determine how well the index and the disturbance gradient reflected this taxonomic variability. We first defined Operational Taxonomic Units (OTUs) by aggregat-ing some taxa to a higher level and dropping some that were identified above the OTU for that group, an approach intended to minimize the loss of information. For example, the vast majority of occurrences of the family Limnephilidae were identified at the family level, and only a few occurrences were identified to genus or species, so all data were aggregated to the family level. However, most occurrences of Baetidae were identified to the genus level or below, so occurrences identified to the family level were dropped and the OTU was taken to the genus level. To reduce the effects of highly abundant taxa in the data set, we transformed count data using log(x + 1) (McCune and Grace, 2002, pp. 67–79). To reduce the noise that can result from rare taxa, we dropped taxa present in fewer than 5% of samples (N = 6 or less) (McCune and Grace, 2002, pp. 75–78). We used Bray-Curtis distance in PC-ORD 4 (MjM Software, Gleneden

Beach, Oregon) and performed 40 runs using real data and 50 runs using randomized data. By identifying the number of dimensions that provided a solution with the greatest reduction in stress for real data compared to randomized data (Monte

Carlo test, McCune and Grace, 2002, pp. 131–134), we identified the appropriate number of dimensions for the ordination. We then reran the analysis using the starting configuration that resulted in the lowest final stress value for that number of dimensions.

We examined plots of the ordination axes with sites coded as reference, stressed, or other, and examined correlations of

NMDS axes with the index scores for all sites. We examined relationships of PCA axes, individual abiotic variables, metrics, and individual taxa with NMDS axes to assess correspondence among different types of data. We used two non-parametric correlation measures to assess these relation-ships: Kendall’s Tau for correlations with individual taxa counts, and Spearman rank correlations for PCA axes, individual abiotic variables, and macroinvertebrate metrics.

Spearman correlations can be interpreted similarly to Pearson correlations, but Kendall’s Tau has a slightly different interpretation. Tau is actually an estimate of the difference between the probability that two variables (e.g., NMDS axis 1 and taxon X) are in the same order (i.e., when axis scores increase, taxon X abundance increases) and the probability that they are not in the same order (i.e., one increases and the

Fig. 2 – Disturbance gradient across sites, represented by first two principal components (PC 1 and PC 2). Dotted lines represent the 25th and 75th percentile PC 1 values, delineating reference and stressed sites. Symbols indicate sites by river as follows: IL, Illinois River; MN, Minnesota

River; SC, Scioto River; ST, St. Croix River; WA, Wabash

River; WI, Wisconsin River.

e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8 319 other decreases) (Hollander and Wolfe, 1999, pp. 363–364). For example, a Tau of 0.4 indicates that the probability of the same order is 0.7 and that of different orders is 0.3, because these two probabilities must sum to 1.

2.4.4. Cumulative distribution functions

After establishing the multimetric index, we plotted the empirical cumulative distribution function (CDF) of index scores for each river separately. The CDF for a river represents an estimate of the proportion of river km with index scores at or below a given value, plotted across all possible values. As part of the sampling design, each site received a weight to indicate the number of river km represented by that site, and for this study, all sites on a given river were weighted equally.

By plotting CDFs for all rivers on a single plot, we were able to compare the relative condition across rivers qualitatively. For this comparison, confidence bands were not included in plots.

However, by calculating the proportion of river km at or below the overall median index score, along with 95% confidence intervals (CIs), we were able to make statistical comparisons among rivers. To generate data for CDF curves and CIs, we used the psurvey.analysis package (v. 1.4) for R software released by the EMAP Design Team (http://www.epa.gov/ nheerl/arm/analysispages/software.htm). Following a Bonfer-roni approach, the 95% CIs were recalculated as 99.67% CIs

(adjusted 95% CIs) to account for the 15 pairwise comparisons among rivers (adjusted a = 0.05/15 = 0.0033). Then overlap among rivers was determined using the adjusted 95% CIs.

3. Results

3.1. Defining the disturbance gradient

A human-influenced disturbance gradient was evident in the first PCA axis. Only the first axis (PC 1) was significant according to broken-stick eigenvalues (l1 = 5.60). PC 1 reflected a gradient from sites with higher canopy cover, habitat scores, Secchi depth, and percentages of plant-based substrates (see P_PLANT in Table 2 for components of this substrate) at the negative end to those with higher conductivity and nutrient concentrations at the positive end (Table 2). Thus, only PC 1 was used as a composite disturbance gradient for the purposes of identifying least disturbed and stressed sites. We identified reference sites asthosefallingbelow the 25thpercentilealongthis gradientand stressed sites as those above the 75th percentile of PC 1 (Fig. 2).

Of the 34 reference sites, 23 were from the St. Croix River and 11 from the Wisconsin River. Among stressed sites, there were 22 in the Illinois River, 8 in the Minnesota River, 1 in the Scioto

River, and 3 in the Wabash River.

Percent forest land use in the riparian corridor was most highly correlated with PC 1 (Spearman r = �0.617), followed by percent agriculture (r = 0.476). Percent development was not stronglyrelatedtoPC1(r = 0.202).Onlysiteelevationandlatitude had correlation magnitudes greater than 0.40 with either PC 1

(r = �0.64 and �0.70, respectively) or land use (percent agricul-ture: r = �0.49 and �0.70, respectively, percent forest: r = 0.40 and 0.63, respectively). Plots of river mile against PC 1 did not indicate trendsfor most rivers (Fig.3),but theuppermost siteson the Illinois River showed some trend toward less human disturbance and the Wisconsin River showed a strong relation-ship between longitudinal position and disturbance, with downstream areas having higher levels of disturbance.

3.2. Metric evaluation

Upon randomly dividing the data into calibration and validation data sets, we did not find any significant difference in the distribution of sites across the PC 1 axis (Kolmogorov–

Smirnov test: Dmax = 0.099, p = 0.94). The calibration data set contained 23 reference and 23 stressed sites, and 44 sites fitting neither category. The validation data set contained 11 each of reference and stressed sites and 21 other sites.

Only about a third of metrics passed the criteria for adequate range, responsiveness, and relative scope of impair-ment according to the criteria already described. Fifteen metrics had limited ranges of values, and of the remaining 82 metrics, 29 showed a lack of responsiveness in box plots. Of the remaining 53 metrics, 15 did not have adequate relative SOI.

All of the remaining metrics were highly correlated with at least one other metric, and five were correlated with so many other metrics that they wereexcluded from consideration for the final index (i.e., total taxa richness, insect taxa richness, percent individuals in the dominant five taxa, facultative tolerance taxa richness, and clinger taxa richness). From the remaining 38 metrics, we were able to selected 9 that represented richness, composition, functional feeding group, behavioral habit (e.g., clinger, swimmer, burrower), and tolerance but were not highly correlated with one another: Diptera taxa richness, Ephemeroptera, Plecoptera, and Trichoptera (EPT) taxa richness, % Coleoptera taxa, % oligochaete and leech taxa, % collector-filterer individuals, predator taxa richness, % http://www.epa.gov/nheerl/arm/analysispages/software.htm http://www.epa.gov/nheerl/arm/analysispages/software.htm

Fig. 3 – Relationships between river mile and the disturbance gradient (PC 1) within each river.

e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8320 burrower taxa, tolerant taxa (tolerance values greater than 6) richness, and % facultative individuals (tolerance values between 4 and 6). Each metric was selected based on responsiveness, variability in reference sites, correlations among metrics, assemblage component represented by metric, ecological meaningfulness, and ease of interpretation and calculation. All metrics selected were positively related to water quality, such that increasing values were associated with less disturbed conditions, except % oligochaete and leech taxa and % burrower taxa. Thresholds used to score metrics are provided in Table 3. The final index, the Non-wadeable Macroinvertebrate Assemblage Condition Index

(NMACI), was the average of individual metric scores. The

NMACI showed strong separation of reference and stressed sites (Fig. 4) and a strong correlation with the disturbance gradient represented by PC 1 (Pearson r = �0.82, N = 90).

3.3. Validation

The relationship between the NMACI and the disturbance gradient were similar between the validation and calibration data sets. There was again a strong separation of reference and stressed sites (Fig. 4), and the Pearson correlation with PC 1 Table 3 – Thresholds (ceiling and floor) and formulae for scori

Metric Ceiling (95th percentile)

Diptera taxa richness 24

EPT taxa richness 17

% Coleoptera taxa 13.3

% Oligochaete and leech taxa 11.5

% Collector-filterer individuals 29.6

Predator taxa richness 16

% Burrower taxa 29

Tolerant taxa richness 15

% Facultative individuals 66.2 was �0.58 (N = 43) for validation sites. Among reference sites, the mean NMACI scores were 70.9 (S.E. = 1.82, N = 23) and 63.6

(S.E. = 3.46, N = 11) for the calibration and validation data sets, respectively.

Only the Wabash River showed any strong relationship between the NMACI and river mile with a correlation of 0.62, indicating that NMACI score improved as one moved up upstream. From scatter plots, this trend was only apparent in the uppermost quarter of the river (not shown). All other rivers showed no trend and had correlation magnitudes of 0.20 or less with river mile.

A three-axis solution was selected for the NMDS, which was based on all 133 sites and 159 taxa. The final stress was relatively low for ecological data at 15.34. Axis 3 represented the largest amount of variation (38.3%), followed by axis 2

(27.5%), and axis 1 (15.4%). Ordinations of axes 2 and 3 revealed strong separation among reference, stressed, and other sites, with other sites falling between reference and stressed sites

(Fig. 5). Axis 3 mostly strongly distinguished reference and stressed sites, whereas axis 2 separated reference and stressed sites from other sites. Each of the two primary axes were strongly correlated with several taxa (Table 4), and both axes were correlated with the clam genera Pisidium and Sphaerium.

ng metrics included in the NMACI

Floor (5th percentile) Scoring formula

9 [(X � 9)/15] � 100

2 [(X � 2)/15] � 100

0 [(X)/13.3] � 100

1.8 [(11.5 � X)/9.7] � 100

0.7 [(X � 0.7)/28.9]*100

3 [(X � 3)/13] � 100

7.5 [(29 � X)/21.5] � 100

6 [(X � 6)/9] � 100

8.7 [(X � 8.7)/57.5] � 100

Fig. 4 – Comparison of distributions across site types and data sets. The box represents the interquartile range and the horizontal line in the box represents the median value. The whiskers represent the 10th and 90th percentiles of the distribution.

e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8 321 The NMACI and several component metrics were also strongly correlated with NMDS axes, particularly axis 3 (Table 5; Fig. 6).

The disturbance gradient, as represented by PC 1, was highly correlated with the third NMDS axis (Spearman rank r = �0.82) and had a strong quadratic relationship with NMDS axis 2 (Spearman r = �0.64 with (PC 1)2). Many of the same abiotic variables driving PC 1 were strongly related individu-ally to NMDS axis 3 (jSpearman rj > 0.5) (Fig. 7). Only DISTRV was linearly related to NMDS axis 2 (Spearman r = 0.530), although conductivity showed a strong quadratic relationship with NMDS axis 2 with a peak at approximately 400 mS/cm.

3.4. Cumulative distribution functions

The river-specific CDFs showed a great deal of overlap for the

Minnesota, Scioto, and Wisconsin rivers, with the St. Croix

River exhibiting the highest scores in general, and the Wabash Fig. 5 – Ordination of reference, stressed, and other sites based on abundances of macroinvertebrate taxa.

and Illinois rivers having the lowest scores (Fig. 8). However, given small sample size in the Minnesota, the level of overlap of CDFs is difficult to assess. The estimated proportion of river km falling below the overall median NMACI score of 48 varied strongly by river, although there were few significant differences among rivers (Table 6).

4. Discussion

The primary disturbance gradient observed across the rivers in this study was a contrast of high conductivity and nutrient concentrations with good water clarity and canopy cover, along with high habitat scores and plant material-based substrate. Based on this gradient, most of the St. Croix River sites were considered reference condition, and most of the

Illinois River sites were considered impaired. Much of the

Illinois River basin consists of agricultural land use, whereas agriculture only predominates in parts of the lower St. Croix

River. This contrast in land use corresponds to the location of sites along the observed disturbance gradient (PC 1 axis).

We found the St. Croix River to be the least impacted of the rivers included in this study. The St. Croix River was one of the first designated National Wild and Scenic rivers. The river is largely free flowing, and its protected status makes it the least disturbed large river in this region (Delong, 2005). Most urban development and agriculture occur only in the lower St. Croix basin where there is some concern about increasing stress associated with the Minneapolis-St. Paul metropolitan area

(Niemela et al., 2005). The St. Croix River had higher NMACI scores than any other river in this study (Fig. 8), partly attributed to higher scores for most metrics.

The Wisconsin, Minnesota, and Scioto rivers have moder-ate levels of anthropogenic disturbance and shared similar

NMACI scores (Fig. 8). The Wisconsin River is heavily dammed, particularly in the middle and upper reaches, and has more

Table 4 – Taxa with Kendall’s Tau correlations (jTauj > 0.40) between NMDS axes 2 and 3 and individual taxon abundances

Axis 2 Axis 3

Taxon Kendall’s Tau Taxon Kendall’s Tau

Brachycercus �0.46 Caecidotea 0.44

Campeloma �0.40 Cheumatopsyche 0.49

Glyptotendipes 0.59 Marstonia 0.41

Hyalella �0.48 Microtendipes 0.49

Micropsectra/Tanytarsus �0.40 Phaenopsectra/Tribelos 0.44

Pisidium �0.42 Pisidium 0.49

Sphaerium �0.48 Pseudocloeon 0.43

Stempellina �0.42 Sigara 0.47

Trichocorixa 0.52 Simuliidae 0.48

Sphaerium 0.42

Stempellinella 0.44 e c o l o g i c a l i n d i c a t o r s 9 ( 2 0 0 9 ) 3 1 3 – 3 2 8322 agriculture and other anthropogenic disturbance than the St.

Croix River. The Wisconsin, however, retains free flowing sections that are relatively undisturbed and the NMACI scores were generally higher than all rivers other than the St. Croix.

Though only eight sites were sampled, the third of the northern rivers, the Minnesota River, generally had lower

NMACI scores than either the Wisconsin or the St. Croix River, and was more similar to the Scioto River in the southern portion of the region. The basin is almost entirely used for agriculture (�95%) which has led to elevated conductivity, sedimentation, nutrient and pesticide concentrations, as well as reduced fish species diversity (Delong, 2005).

The Scioto and Wabash rivers are located in the southern portion of the region and are both tributaries to the Ohio River.

The upper Scioto River is impacted by agriculture and urbanization, while much of the lower basin remains forested.

A portion of the Big and Little Darby Creek tributaries are also designated as National Wild and Scenic rivers. The mainstem

Scioto River is generally turbid and heavily incised with higher nutrient concentrations and conductivity than other Ohio

River tributaries (White et al., 2005). The Scioto River had the highest NMACI scores of the southern rivers, but was intermediate when compared with the three northern rivers.

The Wabash River has the largest river basin of those included in this study and remains relatively free flowing along much of its length. The Wabash basin is predominantly agricultural and, like the Scioto, the channel was incised and suspended sediments, conductivity, and nutrient concentrations were relatively high. Wabash River NMACI scores were typically lower than all other rivers except the Illinois River.

Table 5 – Spearman rank correlations (N = 133) of NMACI component metrics with NMDS axes 2 and 3

Metric NMDS 2 NMDS 3

Diptera taxa richness �0.44 0.58

EPT taxa richness �0.31 0.75

% Coleoptera taxa �0.03 0.42

% Oligochaete and leech taxa 0.16 �0.61

% Collector-filterer individuals �0.30 0.72

Predator taxa richness �0.24 0.65

% Burrower taxa 0.28 �0.53

Tolerant taxa richness �0.27 0.49

% Facultative individuals �0.49 0.54 Our results suggest the Illinois River is highly impaired…