A framework for developing ecosystem-specific nutrient criteria: Integrating biological

thresholds with predictive modeling

Patricia A. Soranno1

Department of Fisheries and Wildlife, Michigan State University, East Lansing, Michigan 48824

Kendra Spence Cheruvelil,2 R. Jan Stevenson, and Scott L. Rollins3

Department of Zoology, Michigan State University, East Lansing, Michigan 48824

Sarah W. Holden and Sylvia Heaton

Water Bureau, Michigan Department of Environmental Quality, Lansing, Michigan 48909

Eric Torng

Department of Computer Science and Engineering, Michigan State University, East Lansing, Michigan 48824

Abstract

We present a novel ecosystem-specific framework for developing nutrient criteria from biological thresholds and

predictive modeling (BTPM) and an application of this framework to lakes in Michigan, U.S. The four main

components for the BTPM framework are: (1) to predict each ecosystem’s ‘expected nutrient concentration’ in the

absence of human effects using a predictive model, (2) to identify important biological thresholds along a nutrient

gradient (i.e., biological [BIO] benchmarks), (3) to determine each ecosystem’s current nutrient concentration, and (4) to

use the above information to derive a nutrient criterion for each ecosystem using the BTPM algorithm. The BTPM

framework is extremely flexible in that it can be applied to any aquatic ecosystem type or nutrient and the four

components can be implemented in a variety of ways. Our BTPM framework has two additional features: it recognizes

that prior to human disturbance, ecosystems varied in their natural nutrient concentrations, and it incorporates risk into

the decision-making process. In the simplest scheme, a nutrient criterion is set at a BIO benchmark greater than the

expected nutrient concentration. However, to protect ecosystems more conservatively, a criterion is set at current lake

nutrient concentrations if current is less than the BIO benchmark. In our application of the BTPM framework, we

developed total phosphorus (TP) criteria for a diverse set of 374 lakes in MI. The expected lake TP concentrations in the

absence of human effects ranged from 3 mg L21 to 24 mg L21, suggesting that a single criterion approach would not be

appropriate. We found two predominant benchmarks in the biological data along the TP gradient, one for zooplankton

metrics at 8 mg L21, and one for phytoplankton metrics at 18 mg L21. We present the sequence of analyses and

decisions that could be used to apply this approach in a management context using Michigan lakes as an example.

Although the primary role that phosphorus plays in lake

eutrophication has been known for decades, setting

management goals for phosphorus has been challenging.

In the U.S., the measurement of nutrient levels in water

bodies serves as the main basis for assessing progress

towards meeting the goals of the Clean Water Act to

protect the designated uses of water bodies (USEPA 2000;

Reckhow et al. 2005). In the European Union, the Water

Framework Directive adopted in 2000 mandates that ‘good

ecological status’ be achieved in all waterbodies in

European Union countries by 2015, but provides little

guidance on methods for achieving this goal (European

Union 2000; Søndergaard et al. 2005). In the U.S., states

and tribes consistently identify excessive levels of nutrients

as a major cause of why as much as half of the surface

waters surveyed in the country do not meet water quality

standards (USEPA 2000). The U.S. Environmental Protection Agency has mandated that states either adopt the

ambient nutrient criteria developed by that agency, or

develop scientifically defensible numerical criteria for

1 Corresponding author (soranno@msu.edu).

2 Present address: Lyman Briggs College, Michigan State

University, East Lansing, Michigan 48824. 3 Present address: Environmental Studies, University of California, Santa Cruz, California 95064.

Acknowledgments

Funding for this research was provided by a Michigan

Department of Environmental Quality (MI-DEQ) grant to

R.J.S. and P.A.S., an MSU–Michigan Agricultural Experiment

Station Disciplinary grant to P.A.S., and an EPA grant to P.A.S.,

K.S.C., R.J.S., and others. We thank James Breck for providing

some of the GIS base layers, Michael Belligan, Justin Booth, and

David Lusch for GIS processing of the landscape databases, and

Katherine Webster and two anonymous reviewers for helpful

comments on earlier drafts. Jason Smith and Steve Holden

provided valuable database management help. We also thank

Carrie Scheele, Lesley Knoll, Aaron Jubar, Sherry Martin, and

Mary Bremigan for access to their databases for the biological

analyses, and Chelsea Stephens for zooplankton processing. Craig

Stow and Jim Bence provided statistical guidance, and Howard

Wandell provided help with the historical STORET water

chemistry database. Finally, we thank the MI-DEQ Nutrient

Criteria Technical committee and the MI-DEQ Nutrient Criteria

Implementation Committee for many critical and insightful

discussions during the development of this approach, and Candice

Bauer from USEPA for suggestions and critical evaluation

throughout the development of this approach.

Limnol. Oceanogr., 53(2), 2008, 773–787

E 2008, by the American Society of Limnology and Oceanography, Inc.

773
nutrients (phosphorus and nitrogen) that will protect the

designated uses of water bodies (USEPA 2000). In fact,

without nutrient criteria, a water body typically must reach

nuisance algal conditions before an action is taken, which is

often subjective at best. Despite this need, few states have

established numerical nutrient criteria on a state-wide scale,

and as of this writing, less than ten states have adopted

numeric nutrient criteria for the majority of their lakes (C.

Bauer, pers. comm.).

Several approaches for setting nutrient criteria have been

proposed, and an important assumption in each method is

that a given nutrient criterion acts as an indicator of

whether designated uses are being met (Reckhow et al.

2005). If we consider the ‘aquatic life’ designated use as

measured by biological integrity, then the different

approaches to setting nutrient criteria can be classified as

either implicit or explicit in this regard. In an implicit

approach, biological integrity is assumed to be protected at

minimal human disturbance levels defined by some human

disturbance gradient and associated nutrient value, but the

biological integrity is not measured. In an explicit

approach, changes in specific biological responses are used

as a surrogate for the designated use more explicitly by

examining how they change along a nutrient gradient either

through analytical approaches, expert judgment, or some

combination of the two (Stevenson et al. 2004; Reckhow et

al. 2005). Surprisingly, biological response gradients have

only sometimes been used for setting nutrient criteria (King

and Richardson 2003; Heiskary and Wilson 2005; Reckhow

et al. 2005). This situation may occur because many of the

relationships between lake biology and phosphorus are

linear which provide little support for establishing a nutrient criterion. However, although not always explicitly

quantified, there are many examples of nonlinear relationships between biological responses and nutrients in lakes

(Watson et al. 1997; Jeppesen et al. 2000; Downing et al.

2001; Vadeboncoeur et al. 2003) and in wetlands (Stevenson et al. 2002; King and Richardson 2003; Lougheed et al.

2007). In cases where thresholds can be clearly quantified,

a nutrient criterion can be justifiably set at or below that

threshold value (Stevenson et al. 2004).

Although there has been a wealth of research on

eutrophication, its causes, and its associated effects on

aquatic foodweb dynamics (Smith 1998), there are still few

integrated approaches available for setting nutrient criteria

in aquatic ecosystems that incorporate these insights into

a framework that can be easily applied at the state (or

national) level by agency managers and biologists. Such

a framework must obviously strike a balance between

model and analytical complexity and implementation

practicality in a management context. Management agencies are charged with managing often thousands of aquatic

ecosystems within social, political, and logistical constraints such as limited funding, limited data availability

and limited personnel for implementation. Our goals were

to develop a framework that can work within such

constraints by integrating biological thresholds and predictive modeling (BTPM) of nutrient expected conditions,

and to present a practical implementation of the framework

that can be readily adopted by management agencies. In

our approach, the biological thresholds can be viewed as an

indicator of whether designated uses are being met, which

are integrated into the overall criteria development process,

and which in different applications can take different form.

Overview of the BTPM framework—Our BTPM framework integrates nutrient modeling to predict ecosystemspecific expected condition (using any model form),

biological thresholds (i.e., nonlinear biological responses

along a nutrient gradient), and current nutrient concentrations into an algorithm to derive ecosystem-specific

nutrient criteria. Our framework has four main components (Fig. 1). The framework is flexible in that it can be

applied to any aquatic ecosystem type or nutrient, it can

incorporate any approach to predict an individual ecosystem’s expected condition or to identify thresholds, and

these data can be integrated into an algorithm to determine

nutrient criteria in a variety of different ways. Our

framework has two additional important features. It

recognizes that prior to human disturbance, water bodies

differed in their natural nutrient concentration as a result of

their hydrogeomorphic (HGM) features, which should be

reflected in the establishment of nutrient criteria, and it

incorporates risk into the decision-making process to

choose biological thresholds. Application of this framework requires data on as many ecosystems as possible (to

obtain a wide range in HGM features). Ideally, these

ecosystems would be sampled weekly or monthly during an

important index period (such as summer) and over the

course of multiple years. However, in Michigan and many

other states, there are few lakes with monitoring data from

multiple time periods within or across summers. Our

framework can accommodate data from single time periods

or multiple time periods. Similarly, if lacking the necessary

data on relevant biological responses of interest, informaFig. 1. The four main components of an ecosystem-specific

framework for developing nutrient criteria.

774 Soranno et al.
tion from the literature can be used, preferably in

combination with sampled data. Below, we present an

application of our BTPM framework.

Application of the BTPM framework to Michigan lakes—

We applied our BTPM framework to a diverse set of 374

lakes in Michigan, in consultation with the Michigan

Department of Environmental Quality (MI-DEQ), the state

management agency responsible for setting nutrients standards for Michigan’s aquatic resources. We present results

for lake total phosphorus (TP) criteria. A key advantage of

our Michigan lake application is that these individual lake

criteria can be derived efficiently given standard and widely

available limnological monitoring data sets including variables such as: lake morphometry, lake chemistry (total

nutrients), water color, and geographic information system

(GIS) data that are becoming widely available at the state

(and national) level such as geology, ecoregion, and land use/

cover (LULC). Next we describe the main steps that we used

in this application (Fig. 2).

There were two steps to predict the lake-specific expected

TP (also called lake reference condition). First, we built

a hydrogeomorphic-land use model (HGM-LU) that

models lake TP as a function of commonly used predictors

such as HGM features and human LULC (Norvell et al.

1979; Ha ̊kanson 1996). The predictive model for this step

can take a variety of different forms; we used a multiple

regression approach of cross-sectional data on several

hundreds of lakes. Because we did not have phosphorus

loading estimates for these lakes, we used variables that are

related to phosphorus loading. In the second step, we used

the HGM-LU model to predict expected lake TP (EXP0) in

the absence of human LULC effects through backcalculation of each lake by setting the human LULC

coefficients to zero (Dodds and Oakes 2004; Baker et al.

2005). EXP0 is the expected condition for each lake defined

as the least disturbed condition given the state of today’s

landscape (Stoddard et al. 2006). We then calculated

expected condition with allowance (EXPA) by adding an

‘allowance’ to EXP0. The allowance can serve any number

of different purposes such as representing model uncertainty in the estimation of EXP0, or some low or

minimal level of allowable human disturbance to the lake

(Fig. 2). The allowance should be chosen based on the

goals of the nutrient criteria, the definition of expected

condition, ecological justification, or some estimate of

prediction uncertainty, all of which are likely to vary from

state to state. We chose to represent model uncertainty.

To identify biological thresholds, we analyzed lake

biological responses (such as zooplankton, phytoplankton

and macrophyte metrics) along a TP gradient to identify

critical thresholds where major changes in lake biology

occur (King and Richardson 2003; Lougheed et al. 2007),

which we define as BIO benchmarks (Fig. 2). These

benchmarks were used to define BIO zones, which are

defined as the TP values between two BIO benchmarks (as

well as between the lowest benchmark and the lowest TP

concentration, and the highest benchmark and the highest

TP concentration in the dataset). We determined current

TP by using existing sampling data from state databases.

The fourth component was to integrate each lake’s

expected TP, current TP, and the BIO benchmarks to

derive a nutrient criterion for each lake using the BTPM

algorithm (Fig. 2). The BTPM algorithm is made up of

four rules for setting individual lake criteria (Fig. 2): (1) If

modeled expected condition (EXPA) is greater than current

TP, then set criterion at current TP. This rule is particularly

important because it ensures that there is no further

degradation in lakes that have lower TP than is predicted

from the model. (2) If EXPA is less than current TP, but is

in the same BIO zone as current TP (but not the highest

BIO zone), then set criterion at current TP. (3) If EXPA is

less than current TP, but in the highest TP BIO zone, then

set criterion at EXPA. (4) If EXPA is in a lower TP BIO

zone than current TP, then set criterion at the higher TP

BIO benchmark for the zone EXPA is in. The integration of

the BIO zones with the EXPA recognizes that model

predictions have some uncertainty associated with them.

Our application of the BTPM framework to Michigan

lakes is based on six key assumptions that we evaluated and

developed during our criteria development process in

a work-group setting with researchers, managers, and

biologists. First, we assumed that phosphorus is the main

‘stressor’ to lakes. Second, we assumed that there is

important natural variation in TP among lakes because

of HGM features such as catchment and lake morphometry, geology, and regional scale factors such as climate,

which can be modeled. Third, we assumed that benchmarks

should be established to sustain desired levels of one or

more biological attributes of lakes which are related to

a lake’s designated use. Fourth, we assumed that important

lake biological responses should include integrative measures of lake biology from both the pelagic and littoral

zones and lake water clarity, which is related to lake

biology through phytoplankton biomass. Fifth, we assumed that human disturbance of a lake can be reasonably

approximated as the proportion of human LULC in the

lake catchment. Although LULC does not represent all of

the possible human effects on lakes, it is commonly used as

an ‘indicator’ of human effects and is the most widely

available dataset for this purpose. Finally, because lakes in

the U.S. are managed at the state level, we chose the state

as an appropriate spatial scale to build models.

Methods

Databases—Water chemistry data: For water chemistry,

we used data for 374 natural lakes ($0.20 km2) with public

access sampled by the Michigan Department of Environmental Quality (MI-DEQ) from 1975–1982 to match the

time period for our LULC data (1978–1985). We obtained

the data from the Storage and Retrieval database

(STORET). Natural lakes are those with or without a water

control structure on them, but that had an open-water

basin prior to human modification. Data were collected

from each lake on a single date during the summer

stratified season (July–September; Table 1) for TP, alkalinity, and true water color (Table 1).

HGM-LULC data: We created a GIS-based lake HGMLULC database for all 374 lakes (Table 2). Lake morEcosystem-specific nutrient criteria 775
Fig. 2. The four main components in our application of the BTPM framework for Michigan lakes. LULC is land use/cover, morph

is morphometry, and geol is geology. For the fourth component, the letters represent all of the eight possible relationships among lake

expected condition (EXPA), current condition (CUR), and the two biological response (BIO) benchmarks for the 374 Michigan lakes. The

‘% of lakes’ is the number of lakes in the Michigan lake dataset that fall within each of the eight possible relationships for the

BTPM algorithm.

776 Soranno et al.
phometry was measured as: mean depth, maximum depth,

shoreline development factor (the ratio of shoreline

perimeter divided by the circumference of a circle of the

same area), and lake basin slope [(surface area)1/2/mean

depth] (Table 2). Catchment areas were delineated as

cumulative catchments that include all upstream stream

and lake catchments. Stream length included the length of

all upstream streams directly connected to lakes. Bedrock

geology data was grouped into five categories: carbonate,

clastic, hard rock, salt, and iron (Geologic Survey Division,

MI-DEQ). Surficial geology data was grouped into five

categories: dune sand, glacial till, lacustrine, moraine, and

outwash (Michigan Natural Features Inventory and

Michigan Department of Natural Resources [MDNR]).

We obtained LULC data from the Michigan Resource

Information Service (MDNR), where urban, agriculture,

forest, upland field, and open water classes were determined from aerial photographs taken between 1978 and

1985 at a resolution of 0.025 km2.

We quantified the proportions of LULC, bedrock

geology, and surficial geology in the 500-m buffer around

each lake, which we use as an indicator of the lake’s

catchment. Although these 500-m buffers are not true

catchments, analyses on these Michigan lakes (Soranno,

unpubl. data) have shown that LULC from the 500-m

buffer is highly correlated (p , 0.001) to LULC of the

lake’s catchment (% urban, r 5 0.63,% agriculture, r 5

0.83), and these data are readily available for all lakes in

a region with a GIS coverage.

Biological data: We combined biological response

datasets from six recent studies of lakes sampled from

1998–2004 (Table 3). To maximize sample sizes, we

combined data across studies when lakes were sampled

using the same methods. Data were based on single visits to

a lake with the exception of the study by Bremigan et al.

(2005), from which we had data for six lakes and two years

of contrasting TP levels. The biological responses examined

included: phytoplankton biomass and clarity metrics,

phytoplankton community and toxin metrics, micro- and

macrozooplankton community metrics, and macrophyte

cover metrics. (For a complete list of metrics, see Tables 6

and 7 and see citations in Table 3 for further details on

methods and variables.)

Step 1: Predict lake-specific expected TP—We built the

HGM-LU model that predicts lake epilimnetic TP from

HGM variables, limnological variables that influence TP

(water color and alkalinity; Chow-Fraser and Duthie 1983;

Ha ̊kanson et al. 2005), ecoregion, and LULC (Tables 1 and

2).

To build the model, we used stepwise multiple regression

(both forward and backward) and included all predictor

variables that were significantly univariately correlated to

TP (correlation coefficient, p # 0.10). We transformed

variables as needed to meet assumptions of normality and

homogeneity of variances. We chose the best model based

on the Akaike information criterion (AIC) because it takes

into account the number of parameters in the model and

penalizes models with more parameters; models with AIC

values of 10 units lower than other models were assumed to

have the most support (Burnham and Anderson 2002). We

took several additional steps to minimize the potential

multicollinearity problems associated with using multiple

linear regression. First, we conducted both forward and

backward stepwise regression and used a conservative

significance level (p # 0.05) to retain predictor variables in

the models. Second, we quantified the level of multiTable 1. Median, minimum, and maximum values for the

lake chemistry variables. n 5 374 lakes for all variables except

alkalinity (n 5 372).

Lake chemistry Median Minimum Maximum

Total phosphorus (mg L21) 13.0 1.0 155

Alkalinity (mg L21 CaCO3) 105 1.0 225

Water color (Pt Co) 9.0 1.0 99

Table 2. Median, minimum, and maximum values for the

hydrogeomorphic (HGM) and land use/cover (LULC) features. n

5 374 for all variables. Stream length is the length of upstream

streams connected to the given lake. Precipitation is the 30-yr

annual average. Drainage area is the catchment area plus the lake

area. Shoreline development factor, all ratios, and lake basin slope

are unitless. Land use/cover percentages, geology percentages and

precipitation were calculated on the 500-m buffer surrounding

each lake.

HGM features and LULC Median Minimum Maximum

Lake and catchment morphometry

Lake area (km2) 0.61 0.20 75.8

Mean depth (m) 4.5 0.8 42.5

Maximum depth (m) 12.8 3.0 86.9

Shoreline development factor 1.8 1.0 6.3

Lake basin slope 194 46.4 1,825

Catchment area (km2) 7.1 0.21 1,948

Catchment area : lake area 10.0 0.7 2,655

Drainage area (km2) 8.4 0.44 1,949

Drainage area : lake area 11.0 1.7 2,656

Stream length (km) 1.4 0 2,805

Climate

Precipitation (mm) 835 706 1,008

Surficial geology

% dune 0 0 53

% outwash 30 0 100

% moraine 0 0 100

% exposed bedrock 0 0 100

% peat and muck 0 0 71

% lacustrine 0 0 100

% glacial till 0 0 100

Bedrock geology

% carbonate 0 0 100

% clastic 100 0 100

% hard rock 0 0 100

% salt 0 0 100

% iron 0 0 100

LULC

% agriculture 6 0 81

% urban 13 0 87

% forest 50 1 98

% upland vegetation 6 0 54

% wetland 6 0 49

% open water 1 0 17

Ecosystem-specific nutrient criteria 777
collinearity in the predictor variables in the final model by

calculating variance inflation factors and using a cutoff of

30 (Hair et al. 1995), above which predictor variables are

assumed to suffer from high multicollinearity. For the best

model, we calculated partial R2’s for each predictor

variable to assess the worth of collecting each piece of

information for practical relevance and application of the

model.

For the best HGM-LU model, we also tested whether

including an ecoregion framework explained additional

variation in lake TP. Because ecoregions are also correlated

to a variety of HGM and LULC variables (Wickham et al.

2005), we examined the residuals of the best HGM-LU

model and performed an ANOVA with ecoregion as the

predictor. We used Ecological Drainage Units (Higgins et

al. 2005) as our ecoregion because it has been shown to

partition variation in TP for Michigan lakes (Cheruvelil et

al. in press).

Using the best HGM-LU model, we then calculated

EXP0 by setting the model coefficients for agriculture and

urban LULC to zero. We calculated EXPA by adding the

allowance to EXP0, which we estimated by calculating the

upper 75th confidence interval from the HGM-LU model

to account for model uncertainty.

Step 2: Identify biological thresholds along TP gradient—

We examined BIO response variables that are commonly

sampled by agencies, are related to attributes that the

public value (e.g., algal biomass influences water clarity)

and are good indicators of ecological integrity in lakes (e.g.,

zooplankton are important component of lake pelagic

foodwebs). We quantified the relationship between each

BIO response and TP in two ways; we used linear

regression on transformed variables (where necessary) to

quantify linear relationships and regression tree analysis

(RTA) on untransformed variables to quantify nonlinear

relationships using Systat 11.0 (Systat Software, Inc.). RTA

is a recursive data partitioning algorithm that initially splits

the response variable into two subsets based on the value of

the predictor variable (TP) that maximizes the reduction in

total residual sum of squares from the parent group to the

two daughter groups (Breiman et al. 1984). Each node is

then split using the same algorithm. At each split, we

examined the ‘proportional reduction in error’ (PRE),

which is a goodness of fit statistic and similar to an r2 value.

We considered that there was scientifically supportable

evidence for a shift in the BIO response if the PRE was

$0.30 (moderate-strong) and if the PRE was greater than

the r2 from the linear model.

For these moderate-strong relationships, we quantified

uncertainty by defining the ‘TP benchmark risk’ using

a bootstrap simulation (Qian et al. 2003). We resampled the

data (with replacement) 1,000 times and calculated the

single strongest threshold for each simulated dataset. The

‘TP benchmark risk’ was defined as the cumulative

probability of a threshold based on the relative frequency

of each threshold in the distribution (King and Richardson

2003), in other words, the proportion of bootstrapped

thresholds that are equal to or lower than a given TP

threshold. The bootstrapped RTA was conducted in the R

statistical environment using the rpart and bootstrap

libraries (R Development Core Team, 2007). We then

plotted cumulative probability curves for each metric,

which described the cumulative risk of a change in the

response variable (BIO) with changing TP values (Qian et

al. 2003).

To inform our decision of choosing the BIO benchmarks, we calculated BIO response factors by first

calculating the median value of each BIO response variable

above and below a BIO threshold, and then dividing the

larger median by the smaller median. We quantified the

BIO response factor at three risk levels: 1%, 10%, and 50%.

We assumed that there was an ecologically important

change in the BIO response when the medians above and

below a threshold differed by a factor of two (although any

factor can be chosen). Because the BIO response can either

increase or decrease with increasing TP dependent on the

relationship (e.g., Secchi depth decreases and algal biomass

increases with increasing TP), we calculate the factor by

dividing the larger median by the smaller regardless of

whether the BIO response increased or decreased. We use

these BIO response factors to help identify the point at

Table 3. The biological datasets used for the biological response threshold analysis. n refers to the number of lakes in each dataset,

years refers to the year that the lakes were sampled in, and TP range is the minimum and maximum TP value in the dataset.

Biological datasets Response variables n Years TP range (mg L21)

(a) Martin and Soranno 2006 Chl a 64 2003 2.6–34.0

Secchi depth 63

(b) Cheruvelil and Soranno. In press. Macrophyte cover 53 2000–2001

2002–2003

4.4–65.6

(c) Knoll 2004 Chl a 64 4.8–64.0

Secchi depth 64

Extinction coefficient 64

Phytoplankton 56

Microzooplankton 46

(d) Scheele unpubl. data Macrozooplankton 15 2002–2003 5.3–19.3

(e) Bremigan et al. 2005 Chl a 13 1998–1999 11.8–32.3

Secchi depth 13

Macrozooplankton 13

(f) Jubar 2004 Secchi depth 15 2004 10.0–24.6

Macrozooplankton 15

778 Soranno et al.
which an important change in the BIO response has

occurred, which we define as a BIO benchmark and is one

of the inputs into the BTPM algorithm.

Step 3: Determine current TP in each lake—We used the

existing sampling data from the 374 Michigan lakes that

were sampled during one time period in the summer (see

Databases for details).

Step 4: Derive lake-specific nutrient criteria using the

BTPM algorithm—We integrated the results from step 1 to

step 3 into the BTPM algorithm to develop lake-specific

nutrient criteria (see Fig. 2 for the graphical representation

of this algorithm in step 4). The inputs to the algorithm are

the measured TP concentration for each lake (current), the

expected lake TP concentration with allowance (EXPA),

and the BIO benchmarks. We define ‘BIO zones’ as TP

values between any two adjacent BIO benchmarks, and

between a BIO benchmark and either the maximum or

minimum TP value in the dataset. The BTPM algorithm is

made up of the four rules stated earlier that defines all the

possible relationships among these variables (see Application of the BTPM framework to Michigan lakes).

Results

Step 1: Predict lake-specific expected TP—The lakes

used to build the HGM-LU model ranged in TP from

1 mg L21 to 155 mg L21, with a relatively low median of

13 mg L21 (Table 1). The lakes also ranged widely in

alkalinity and water color (Table 1). Of the many HGM

variables univariately related to lake TP (Table 4), the

strongest included common predictors of TP such as lake

depth, LULC, catchment morphometry, precipitation, and

water color.

The best HGM-LU model included variables known to

be correlated with P loading to lakes (Table 5, Fig. 3).

Although significant, outwash geology only adds a small

amount to the overall R2, therefore, if a management

agency were to use this model, they could choose whether

to retain this variable based on cost and effort of including

it. To assess the importance of ecoregions for predicting

reference TP, we performed an ANOVA with the residuals

from the best model as the response and Ecological

Drainage Units as the predictor. Because ecoregion was

not significant (p . 0.05), it was not incorporated into the

final model.

Using the best HGM-LU model, we estimated EXP0 and

EXPA for each lake. The upper TP range for current lake

condition exceeded the upper range for both EXP0 and

EXPA as expected (Fig. 4). In the absence of human

disturbance (EXP0), TP in our study lakes ranged from

3 mg L21 to 24 mg L21, demonstrating the wide range in TP

expected condition. When incorporating an allowance

(EXPA), TP expected condition ranged from 4 mg L21 to

33 mg L21.

Step 2: Identify biological thresholds along TP gradient—

There were 36 BIO response relationships with TP as the

predictor. Fourteen of them had a PRE $0.20, and seven

of those had a PRE $0.30 (moderate–strong) and a higher

fit than the r2 from linear regression models (Tables 6 and

7). We grouped these BIO responses into two ecological

classes: 1) phytoplankton and water clarity and 2)

zooplankton. Figure 5 shows the raw data for each

moderate-strong BIO response relationship and the location of the RTA threshold(s). This figure also depicts the

cumulative frequency plot of the 1,000 bootstrapped

thresholds. The interquartile ranges of the bootstrapped

thresholds ranged from 5 mg L21 to 37 mg L21 TP for

phytoplankton and clarity metrics, and from 6 mg L21 to

13 mg L21 TP for zooplankton metrics (Table 7). Because

the thresholds calculated using RTA were high relative to

the range of bootstrapped thresholds (Fig. 5 and Table 7),

we suggest caution in using RTA values for benchmark

Table 4. Univariate correlations between hydrogeomorphic

(HGM), land use/cover (LULC), and water chemistry predictor

variables and TP. Units and descriptions are the same as for

Table 2. Correlations were performed on transformed variables

when necessary. Correlation coefficients that are significant at

alpha #0.10 are in bold and were included in the predictive

modeling step. n 5 374 for all variables, except for alkalinity (n

5 372).

Predictor variables r p

Lake and catchment morphometry

Lake area (km2) 20.050 0.336

Mean depth (m) 20.292 ,0.001

Maximum depth (m) 20.301 ,0.001

Shoreline development factor 0.044 0.394

Lake basin slope 0.199 ,0.001

Catchment area (km2) 0.107 0.038

Catchment area : lake area 0.160 0.002

Drainage area (km2) 0.096 0.063

Drainage area : lake area 0.093 0.073

Stream length (m) 0.101 0.051

Climate

Precipitation (mm) 0.218 ,0.001

Surficial geology

% dune 0.095 0.066

% outwash 20.097 0.062

% moraine 0.107 0.039

% exposed bedrock 20.032 0.531

% peat and muck 0.053 0.310

% lacustrine 20.062 0.233

% glacial till 0.052 0.313

Bedrock geology

% carbonate 20.104 0.044

% clastic 0.171 0.001

% hard rock 20.077 0.137

% salt 0.013 0.799

% iron 20.083 0.109

LULC

% agriculture 0.287 ,0.001

% urban 0.148 0.004

% forest 20.318 ,0.001

% upland vegetation 0.093 0.072

% wetland 0.007 0.894

% open water 20.058 0.265

Water chemistry

Alkalinity 0.084 0.131

Water color 0.410 ,0.001

Ecosystem-specific nutrient criteria 779
identification without consideration of the uncertainty or

risk associated with those values. We defined risk as the

proportion of bootstrapped thresholds that were equal to

or lower than a given TP threshold, which we incorporated

into the next step that identifies the BIO benchmarks.

To choose the BIO benchmarks for the BTPM algorithm,

we used our estimates of risk in the seven moderately-strong

BIO threshold relationships at three levels (1%, 10%, and

50%) and the BIO response factors (Table 8). For phytoplankton and water clarity metrics, medians for both

chlorophyll a (Chl a) and phytoplankton biomass differed

by more than a factor of two (i.e., biologically important) at

the 1% risk level, which was exacerbated at higher risk

values. For extinction coefficient, the medians differed by

a factor of two at the 10% risk level. In contrast, medians for

Secchi depth never differed by a factor of two. This suggests

that although the thresholds derived for Secchi depth were

statistically significant, they were not biologically meaningful. In addition, the PRE for Secchi depth was the lowest of

all of the moderate-strong relationships (0.32, Table 7), and

in fact, was very similar to the r2 for the linear regression

Table 5. Multiple regression results for the best hydrogeomorphic-land use model (HGM-LU) for predicting lake TP (R2 5

20.39, p , 0.001, F 5, 368 5 46.6). SE refers to the parameter

estimate standard errors and partial R2’s provide the contribution

of each variable to the overall model.

Predictor variables

Parameter

estimate SE Partial R2

Constant 1.80 0.138

Mean depth 20.24 0.050 0.05

Surficial geology

(% outwash) 20.18 0.065 0.02

Water color 0.36 0.037 0.18

% agriculture 1.10 0.146 0.13

% urban 1.24 0.193 0.10

Fig. 3. Predicted versus observed (current) TP for the best

hydrogeomorphic-land use model. See Table 5 for model description.

Fig. 4. Frequency histograms of lake TP for (A) current TP,

(B) the expected TP with allowance (EXPA), and (C) the expected TP without allowance (EXP0). The x-axis ends at

75 mg L21, thus two lakes are not shown at 118 mg L21 and

155 mg L21 for (A).

780 Soranno et al.
(0.29). Because the nonlinear relationship was not well

supported by the data, thresholds for Secchi depth did not

contribute towards setting criteria in Michigan lakes. For the

zooplankton metrics, the difference in the medians above

and below the benchmarks was biologically important for

proportional biomass of Daphnia and zooplankton biomass

at the 1% risk level. Based on these analyses, we identified

two BIO benchmarks, one associated with phytoplankton

metrics, and the other associated with zooplankton metrics.

BIO benchmarks must be chosen with care to minimize

the risk of losing biological integrity in an ecosystem by

setting a benchmark either too high or too low. Although

choosing the lowest possible BIO benchmark may seem

protective and conservative, this may not actually be the

case once the benchmarks are incorporated into the BTPM

framework. For example, consider what happens in a twobenchmark algorithm with 20 mg L21 TP as the upper BIO

benchmark and either 10 mg L21 or 12 mg L21 TP as the

lower benchmark. If we use 12 mg L21 as the lower BIO

benchmark, lakes with EXPA of 11 mg L21 will be managed

to maintain a benchmark of 12 mg L21. Although choosing

10 mg L21 as the lower BIO benchmark seems as if it would

be more protective, doing so means that lakes with EXPA

of 11 mg L21 will be managed for the upper BIO

benchmark of 20 mg L21. As a result, choosing excessively

low or high benchmarks may not be prudent. Therefore, we

chose intermediate benchmark values, using risk levels

from either 1% or 10% (depending on the BIO response

factor) for thresholds that were associated with an

important ecological change. For phytoplankton metrics,

we chose 18 mg L21 (range from 16.0 mg L21 to

20.9 mg L21); for zooplankton metrics, we chose 8 mg L21

(range from 7.5 to mg L21 9.1 mg L21).

Step 3: Determine current TP in each lake—We show the

histogram of the current TP data in Fig. 4.

Table 6. Results of the regression tree analysis (RTA) of the relationship between TP and the untransformed biological responses

(BIO) where the fit statistic, the proportional reduction in variance (RTA-PRE), is ,0.30; and r2 and p values for the linear regression

between TP and the BIO responses (transformed where necessary). For the dataset description, letter designation, and citation, see

Table 3. n is the number of lakes in each analysis. For relationships that had more than one threshold, the values in bold are the

thresholds that occur first, with the greatest reduction in deviance. Toxin is microcystin toxin.

Dataset Biological responses

TP range

(mg L21) n

TP threshold(s)

(mg L21) RTA-PRE

Linear

regression

r2 p

c Phytoplankton community and toxins

Microcystis biomass (mg L21) 4.8–64.0 56 34.6 0.11 0.06 0.07

Toxin concentration per Chl a (mg mg21) 4.8–64.0 64 8.1, 9.8 0.21 0.00 0.91

Environmental concentration of toxin (mg L21) 4.8–64.0 64 34.6 0.25 0.17 0.001

Toxin per microcystis biomass (mg mg21) 4.8–64.0 56 34.6 0.17 0.06 0.08

Cyanobacteria biomass (mg L21) 4.8–21.0 51 — 0 0.03 0.26

Microzooplankton

Ciliate biomass (mg L21) 4.8–64.0 54 — 0 0.04 0.15

Rotifer biomass (mg L21) 4.8–44.0 61 20.7 0.16 0.02 0.23

Rotifer species number (#) 4.8–21.0 57 — 0 0.02 0.26

d, e, f Macrozooplankton community

Bosmina spp. mean length (mm) 5.3–25.6 31 19.3 0.07 0.01 0.68

Ceriodaphnia spp. mean length (mm) 5.3–18.2 20 12.3 0.23 0.18 0.06

Daphnia spp. mean length (mm) 5.3–32.3 42 11.2, 12.7 0.25 0.02 0.44

d, e Copepod biomass (mg L21) 5.3–32.3 28 11, 13.4 0.27 0.27 0.004

Daphnia pulicaria total biomass (mg L21) 5.3–25.6 17 11.2 0.16 0.02 0.63

Bosmina spp. total biomass (mg L21) 5.3–25.6 18 14.0 0.25 0.11 0.18

b Macrophyte cover (%, unless otherwise indicated)

Macrophyte edge (perimeter : area, m21) 4.4–65.6 53 8.3, 22.7, 28.4 0.27 0.00 0.90

Lake cover 4.4–65.6 53 8.3, 22.7, 28.4 0.18 0.00 0.69

Lake dense cover 4.4–65.6 53 8.4, 18.7 0.14 0.06 0.74

Lake emergent cover 4.4–65.6 53 19.9 0.08 0.03 0.23

Lake floating cover 4.4–65.6 53 — 0 0.01 0.50

Lake submersed cover 4.4–65.6 53 8.3, 22.7 0.14 0.01 0.44

Lake EWM cover 4.4–65.6 53 12.1 0.10 0.07 0.05

Lake dense EWM1 cover 4.4–65.6 53 28.4 0.14 0.07 0.05

Littoral cover 4.4–65.6 53 19.9 0.10 0.05 0.10

Littoral dense cover 4.4–65.6 53 8.4, 10.8 0.13 0.00 0.97

Littoral emergent cover 4.4–65.6 53 10.8 0.06 0.03 0.19

Littoral floating cover 4.4–65.6 53 — 0 0.00 0.71

Littoral submersed cover 4.4–65.6 53 19.9 0.11 0.07 0.06

Littoral EWM cover 4.4–65.6 53 12.2 0.08 0.05 0.11

Littoral dense EWM cover 4.4–65.6 53 — 0 0.05 0.11

1 EWM is the non-native Eurasian watermilfoil (Myriophyllum spicatum L.)

Ecosystem-specific nutrient criteria 781
Step 4: Derive lake-specific nutrient criteria using the

BTPM algorithm—We entered the following into the BTPM

algorithm: two BIO benchmarks of 8 mg L21 and

18 mg L21, expected lake TP with allowance for each lake

(EXPA), and current lake TP (CUR) for each lake. Figure 2

shows the graphical representation of this algorithm

showing all possible relationships (A–H) among the inputs.

For the 374 Michigan lakes, we calculated how many lakes

fell within each category (Fig. 2). For the 23% of lakes for

which there was a BIO benchmark between the CUR and

EXPA, criteria were set at those benchmarks (B and D). For

the 26% of lakes that did not exceed a BIO benchmark,

criteria were set at CUR (A and C). The only situation in

which criteria were set at EXPA was when both CUR and

EXPA were greater than the largest BIO benchmark and

CUR was greater than EXPA (E), which only occurred in 3%

of lakes in our database. For 48% of the lakes, CUR TP was

less than EXPA TP (F–H). However, recall that we added an

‘allowance’ to each lake’s EXP0, including those lakes with

relatively low human disturbance, and consequently low TP.

Therefore, EXPA is likely to be higher than CUR for many

minimally disturbed lakes that we have in our dataset. To

accommodate this situation and to prevent degradation, one

of the rules in the algorithm states that under these

conditions, criteria are set at CUR rather than EXPA (F–H).

Discussion

We have presented a framework to develop ecosystemspecific nutrient criteria in a scientifically-defensible way

that is flexible enough to incorporate many different

context-specific situations that might occur in different

political jurisdictions or geographic regions. We have also

presented an application of our framework to set TP

criteria in a diverse set of Michigan lakes. We will first

discuss implications of the BTBM framework and general

topics that must be addressed in any application of this

framework using our application to Michigan lakes as

a concrete example, and then we discuss some features and

issues that arose from our specific application of the BTPM

framework to Michigan lakes.

Implications and general features of the BTPM Framework—One critical advantage of our BTPM framework,

compared to other approaches for setting nutrient criteria,

is that our framework explicitly addresses several interdependent questions that should be answered when

setting nutrient criteria: What is the expected condition?

What is the effect of nutrients on biological responses?

What level of nutrients protects biological integrity? And,

what is a reasonable level of protection? Our BTPM

framework addresses these questions by: (1) integrating

multiple pieces of information for setting nutrient criteria

(expected conditions, BIO response thresholds, and current

nutrient concentrations) into a single algorithm, (2)

allowing nutrient criteria to be developed on an ecosystem-specific basis, instead of a single criterion for all

ecosystems within a given geographic region, (3) incorporating uncertainty into the benchmark-identification process which allows policy makers to make risk-level

decisions based on scientific evidence, and (4) allowing

for updating as additional data become available through

refinements of the predictive model or the BIO benchmark

identification procedure.

Although we have presented a fairly detailed description

of a framework to set nutrient criteria using Michigan lakes

as an example, each application of this framework will

likely result in a different set of decisions and justifications

based on available data. In addition, there are still many

decisions that need to be made at the implementation stage.

For example, in our case, lake managers must still decide

the lake sampling period for assessment, the management

actions for lakes that exceed its nutrient criterion, how

‘close’ a lake has to be to its criterion to be considered

attaining or nonattaining (i.e., how much measurement

error is acceptable), and how to regulate surface water

discharges of phosphorus to ensure that lakes meet the

nutrient standards. These decisions will vary region to

Table 7. Results of the regression tree analysis (RTA) of the relationship between TP and the untransformed biological responses

(BIO) where the fit statistic, the proportional reduction in variance (RTA-PRE), is $0.30; and r2 values for the linear regression between

TP and the BIO responses (transformed where necessary); all linear regressions are significant at p # 0.05, unless noted otherwise. For the

dataset description, letter designation, and citation, see Table 3. n is the number of lakes in each analysis. For relationships that had more

than one threshold, the values in bold are the thresholds that occur first, with the greatest reduction in deviance. The IQ-range is the

interquartile range of the bootstrapped thresholds.

Dataset Biological responses

TP range

(mg L21) n

TP threshold(s)

(mg L21)

RTA -

PRE

Linear

regression r2

IQ-range in

bootstrapped

thresholds

Phytoplankton and clarity

a, c, e Chl a (mg L21) 2.6–64.0 140 34.1 0.59 0.42 16–37

c Phytoplankton biomass (mg L21) 4.8–64.0 56 26.7 0.78 0.44 16–31

c Extinction coefficient (m21) 4.8–64.0 64 26.7 0.71 0.47 12–37

a, c, e, f Secchi depth (m) 2.6–64.0 153 5.5, 11.4, 22.0 0.32 0.29 5–21

Zooplankton

d, e, f Cladoceran mean length (mm) 5.3–32.3 43 11.2, 15.3 0.43 ns 6–11

d, e Proportion Daphnia biomass 5.3–32.3 28 11.7, 16.5 0.64 0.53 9–12

Zooplankton biomass (mg L21) 5.3–32.3 28 13.0, 15.6 0.30 0.21 8–13

782 Soranno et al.
region and are beyond the scope of our study here, but are

necessary to address when applying this or any other

framework.

Another challenging issue when developing nutrient

criteria is determining an acceptable allowance (i.e.,

deriving EXPA) for the expected nutrient concentrations.

In our application, TP criteria were set higher than

‘reference condition’ (EXP0) through the incorporation of

an allowance term. However, lakes should be managed to

maintain resilience in response to potentially high interannual variability in nonpoint P loading (Carpenter and

Cottingham 1997) and resulting high variability in lake TP

Fig. 5. (A–D) Biological (BIO) responses for phytoplankton and clarity metrics versus TP for the variables with moderate-strong

regression tree analysis (RTA) fits ($0.30). Each data point represents a single lake. The plotted line is the cumulative frequency of the

1,000 bootstrapped thresholds (right-hand y-axis, proportion). The arrows refer to the thresholds identified from the individual RTA

analysis. Where there is more than one threshold identified, the solid arrow is the first threshold and dotted arrows are subsequent

thresholds. n is the number of lakes in each dataset. (E–G) Same as above except for zooplankton metrics.

Ecosystem-specific nutrient criteria 783
and Chl a (Knowlton and Jones 2006). There is evidence

from modeling studies of single ecosystems that suggests

that we should take a precautionary approach in light of

the large uncertainties in limnological relationships, stochastic interannual P inputs, and lags in implementing

policy (Carpenter et al. 1999). An important limitation to

our application of the BTPM framework to Michigan lakes

is the lack of quantification of interannual P variability

associated with CUR and EXPA, which can influence

resulting management recommendations (Stow et al. 1997).

Unfortunately, quantification of prediction uncertainty

associated with interannual variability requires a large

amount of data typically not available for substantial

numbers of lakes. Future applications of the BTPM

framework should incorporate these important factors

where data are available (Knowlton and Jones 2006) or

develop novel approaches, such as examining lake variance

as an indicator of thresholds (and regime shifts) in

ecosystems (Brock and Carpenter 2006) or taking a Bayesian approach (Borsuk et al. 2004; Lamon and Stow 2004).

Finally, an important component in establishing any

standards for management is the quantification and

communication of uncertainty (Harris et al. 2003). In the

BTPM framework, uncertainty can be quantified for each

of the four major steps. However, there are tradeoffs

among complete modeling of uncertainty, implementation

ease, communication between managers and the public,

and data availability. In our application of the BTPM

framework for Michigan lakes, we incorporated uncertainty in the first two steps. In our predictive HGM-LU model,

we incorporated model error by using the 75th confidence

interval for EXPA. For the BIO benchmarks, we quantified

uncertainty (defined as risk) using a bootstrap simulation

(King and Richardson 2003). Although we did not formally

propagate errors through to the final step of assigning each

nutrient criterion, we found that the BTPM algorithm is

somewhat robust for uncertainties from the previous two

steps, particularly the prediction of EXPA. For example, as

long as the predicted value of EXPA is in the same BIO

zone as its true value, then the nutrient criterion (and

consequently, the management action) remains the same.

We chose the simplest case in our Michigan lakes example

where we treat both the BIO benchmarks and the expected

conditions as single point values for demonstration

purposes and due to lack of data. To handle more complex

methods of quantifying uncertainty including Bayesian

analyses, one could treat both the BIO benchmarks and

expected conditions as ranges of values with probabilities

associated with those ranges. In sum, the BTPM framework can be adapted to incorporate uncertainty at any level

of complexity.

Application of the BTPM framework to Michigan lakes—

In our application of the BTPM framework to Michigan

lakes, we used an HGM-LU model to predict TP in

Michigan lakes, and we generated BIO benchmarks by

searching for nonlinear relationships of commonly measured biological responses with TP. Next, we discuss issues

that arose during our application of these two specific steps

in more detail, pointing out limitations of our application,

Table 8. The median of the biological (BIO) responses above and below a given TP threshold and the BIO response factor (the larger median divided by the smaller

median). n is the number of lakes below the threshold (top) and the number of lakes above the threshold (bottom). Units of BIO response metrics are the same as Table 7.

Resp. is response, and daph. biom. is Daphnia biomass. Values in bold have factors

.2.0, which we define as biologically important. The medians and factors are calculated at

three risk levels (1

%, 10

%, and 50

%).

BIOresponsemetric

1% risk: 10

% risk: 50

% risk:

TP

threshold Median below Median above

BIO

resp.

factor n

TP

threshold Median below Median above

BIO

resp.

factor n

TP

threshold Median below Median above BIO resp. factor n

Phytoplankton and clarity 110 132 137

Chl a 16.0 2.82 6.68 2.4 30 25.8 3.00 20.20 6.7 8 34.3 3.00 51.30 17.1 3

Phytoplank- 45 50 52

ton biomass 16.4 37.85 121.4 3.2 11 20.4 38.10 500.6 13.1 6 30.4 38.90 658.8 16.9 4

Extinction 42 58 60

coefficient 11.9

20.24

20.33 1.4 22 20.9

20.24

20.78 3.3 6 33.1

20.24

21.12 4.7 4

Secchi 10 18 75

depth 5.1 4.55 3.40 1.3 143 5.8 4.40 3.38 1.3 135 11.3 3.90 3.00 1.3 78

Zooplankton

Cladoceran 2 4 7

mean length 6.1 1.31 0.78 1.7 41 7.6 1.27 0.78 1.6 39 11.1 1.21 0.76 1.6 36

Proportion 5 6 9

daph. biom 9.1 0.48 0.01 95.6 23 11.1 0.49 0.01 97.4 22 11.6 0.48 0.00 119.5 19

Zooplankton 4 5 13

biomass 7.5 39.50 132.4 3.4 24 8.4 39.60 138.6 3.5 23 12.8 85.80 142.8 1.7 15

784 Soranno et al.
topics that need further study, and implications for other

applications of this framework.

For modeling lake TP, we took a hydrogeomorphic

perspective that recognizes that lake nutrients are a function

of both local and regional hydrogeomorphic features; and

we used data readily available to state agencies. Predictive

ability would be improved if we had sufficient data to use

well-established process-based models that predict lake TP

concentrations from loading rates and lake hydrology, but

this is rarely the case for large numbers of lakes that need to

be managed. For example, MI has 6,595 lakes .0.04 km2,

most of which are public and must be managed by the state,

often without any lake-specific data available. Thus,

building models from GIS databases that are available

statewide has obvious appeal and warrants further research

to improve predictive ability.

We found that lake HGM predictors captured variation

in TP better than the regional predictor ecoregion. This

result agrees with previous studies that have shown that

ecoregions alone do not capture TP variation well for

streams (Dodds and Oakes 2004, Wickham et al. 2005) or

lakes (Cheruvelil et al. in press), most likely because

ecoregion ignores important local, site-level features such

as lake depth, catchment area or finer-scaled patterns in

surficial geology. Furthermore, the ability of ecoregion to

explain variation in TP is even lower once the effect of

LULC is removed (Wickham et al. 2005; Cheruvelil et al. in

press). Nevertheless, ecoregions should still be considered

when developing a HGM-LU model because ecoregions are

useful for a variety of management applications (Stoddard

2004), continue to be used by state and federal agencies,

and may become more relevant for larger spatial extents

than studied here (Cheruvelil et al. in press). In addition,

ecoregion is likely to become important for study areas

with more heterogeneous landscapes than Michigan

(Heiskary and Wilson 2005) or for other variables besides

TP (such as total nitrogen) (Cheruvelil et al. in press).

One way to test our HGM-LU model is to compare

EXP0 of our lakes with expected TP from presettlement

times obtained from paleolimnological reconstruction of

TP using diatom communities. In a study of northeastern

U.S. lakes, Dixit et al. (1999) found that 47% of their lakes

had TP ,10 mg L21, and 49% were between 10 mg L21 and

30 mg L21. For our study lakes, EXP0 was ,10 mg L21 TP

for 66% of our lakes, and the remaining lakes had TP

,30 mg L21 . Although lakes in northeastern U.S. are not

identical to our study lakes, they share many similar HGM

features, including climate. Paleolimnological reconstruction of historical TP levels for a limited number of lakes

within the region under study is a good supplement to our

approach and can be integrated into our framework.

An important challenge in setting nutrient criteria is the

selection of BIO responses to use for identifying benchmarks. The phytoplankton metrics that we used for the

Michigan lakes analysis have obvious relevance to lake

management because these metrics measure some component of algal biomass which subsequently influences water

clarity. Cladoceran mean length has also been shown to be

strongly related to Secchi transparency because cladocerans

favor algal particles that have high light attenuation

(Stemberger and Miller 2003). Similarly, the proportion

of daphnid biomass is relevant given the dominance of

daphnids in filtering particles, their importance in the

zooplankton community, and their critical role in pelagic

lake foodwebs as a source of food for planktivorous and

juvenile fishes. Depending on data availability, there may

be additional, highly relevant biological responses that

have threshold responses to changes in TP for other

management units.

Our results suggest that nonlinear relationships of BIO

responses with TP may be more common than previously

thought, especially when examined in the context of the risk

of nuisance conditions, such as cyanobacteria blooms

(Downing et al. 2001). Visual inspection of published

relationships suggests the existence of nonlinear relationships, even though they have not been quantified (Jeppesen

et al. 2000). For example, several taxonomic groups of

phytoplankton exhibit nonlinear relationships to increasing

TP, where the thresholds appear to occur between

10 mg L21 and 30 mg L21 TP (Watson et al. 1997), or in

the case of % cyanobacteria biomass at ,30 ug L21 and

,100 ug L21 (Downing et al. 2001). In another study, the

proportion of picophytoplankton declined dramatically

above a TP value of ,8 mg L21 (Watson et al. 1992).

Finally, Vadeboncoeur et al. (2003) examined the proportion of total lake primary production due to periphyton

versus pelagic algae, and visual inspection of their data

suggests a breakpoint at ,10 mg L21 TP, above which the

proportion of total primary production due to benthic

production dramatically declines. Having relatively high

benthic production could be viewed as a desirable state

because it signifies that there is enough light for production

to occur on sediment surfaces in the littoral zone,

facilitating both periphyton and aquatic plant growth and

consequently littoral zone foodwebs. Evidence such as the

above suggests that it will be fruitful to reexamine and

quantify many existing datasets for nonlinear responses of

lake biological and ecosystem dynamics along a TP

gradient.

These nonlinear responses are important to consider in

a management context as they are likely points where

there is an abrupt change in an ecosystem, or where small

changes in the driver produces a large response, which is

conceptually linked to the idea of resilience (Groffman et

al. 2006). An assumption that is often implicit in nutrient

criteria development is that during years of unusually

high phosphorus loads due to climate and anthropogenic

variability, lake biology and phosphorus dynamics may

respond, but not irreversibly so, as might occur with

a regime shift (Carpenter 2003). Thus, during these times,

the lake would not shift to a ‘new state’ beyond which it

would likely not return, as occurs in regime shifts

between clear and turbid states in lakes (Genkai-Kato

and Carpenter 2005). If threshold responses are in fact

tied to regime shifts, we should manage lakes in such

a way to avoid them because reversal is either difficult or

impossible over relevant time scales (Carpenter 2003).

Unfortunately, phosphorus is just one component that

pushes a lake to another regime; factors such as food web

structure, lake physical modifications, climate change,

Ecosystem-specific nutrient criteria 785
and non-native species introductions also determine

where a lake is in relation to regime shifts. Therefore,

all of these factors need to be managed concurrently with

phosphorus and provide further support for a precautionary approach to setting nutrient criteria to ensure that the

criteria correctly represent the designated uses they are

intended to protect.

References

BAKER, E. A., K. E. WEHRLY, P. W. SEELBACH, L. WANG, M. J.

WILEY, AND T. SIMON. 2005. A multimetric assessment of

stream condition in the Northern Lakes and Forests

Ecoregion using spatially explicit statistical modeling

and regional normalization. Trans. Am. Fish. Soc. 134:

697–710.

BORSUK, M. E., C. A. STOW, AND K. H. RECKHOW. 2004. A

Bayesian network of eutrophication models for synthesis,

prediction, and uncertainty analysis. Ecol. Model. 173:

219–239.

BREIMAN, L., J. H. FRIEDMAN, R. OLSHEN, AND C. J. STONE. 1984.

Classification and regression trees. Wadsworth International

Group.

BREMIGAN, M. T., S. M. HANSON, P. A. SORANNO, K. S.

CHERUVELIL, AND R. D. VALLEY. 2005. Aquatic vegetation,

largemouth bass and water quality responses to low-dose

fluridone two years post-treatment. J. Aquat. Plant Manage.

43: 65–75.

BROCK, W. A., AND S. R. CARPENTER. 2006. Variance as a leading

indicator of regime shift in ecosystem services. Ecol. Soc. 11: 9

[accessed on October 12, 2006]. Available from http://www.

ecologyandsociety.org/vol11/iss2/art9/

BURNHAM, K. P., AND D. R. ANDERSON. 2002. Model selection and

multimodel inference, 2nd ed. Springer-Verlag.

CARPENTER, S. R. 2003. Regime shifts in lake ecosystems: Pattern

and variation. V. 15. Excellence in Ecology Series. Ecology

Institute.

———, AND K. L. COTTINGHAM. 1997. Resilience and restoration

of lakes. Conservation Ecology 1: 2 [accessed on October 12,

2006]. Available from http://www.consecol.org/vol1/iss1/art2.

———, D. LUDWIG, AND W. A. BROCK. 1999. Management of

eutrophication for lakes subject to potentially irreversible

change. Ecol. Appl. 9: 751–771.

CHERUVELIL, K. S., AND P. A. SORANNO. In press. Relationships

between lake macrophyte cover and landscape features.

Aquatic Botany doi:10.1016/j.aquabot.2007.10.005.

———, ———, M. T. BREMIGAN, T. WAGNER, AND S. L. MARTIN.

In press. Grouping lakes for water quality assessment and

monitoring: The roles of regionalization and spatial scale.

Environ. Manage.

CHOW-FRASER, P., AND D. H. DUTHIE. 1983. An interpretation of

phosphorus loadings in dystrophic lakes. Arch. Hydrobiol.

97: 109–121.

DIXIT, S. S., J. P. SMOL, D. F. CHARLES, R. M. HUGHES, S. G.

PAULSEN, AND G. B. COLLINS. 1999. Assessing water quality

changes in the lakes of the northeastern United States using

sediment. Can. J. Fish. Aquat. Sci. 56: 131–152.

DODDS, W. K., AND R. M. OAKES. 2004. A technique for

establishing reference nutrient concentrations across watersheds affected by humans. Limnol. Oceanogr. Methods 2:

333–341.

DOWNING, J. A., S. B. WATSON, AND E. MCCAULEY. 2001.

Predicting Cyanobacteria dominance in lakes. Can. J. Fish.

Aquat. Sci. 58: 1905–1908.

EUROPEAN UNION. 2000. Directive 2000/60/EC of the European

Parliament and of the council establishing a framework for

the community action in the field of water policy. European

Commission, Off. J. Eur. Commun. L327 (2000) 1.

GENKAI-KATO, M., AND S. R. CARPENTER. 2005. Eutrophication

due to phosphorus recycling in relation to lake morphometry,

temperature, and macrophytes. Ecology 86: 210–219.

GROFFMAN, P. M., AND oTHERS. 2006. Ecological thresholds: The

key to successful environmental management or an important

concept with no practical application. Ecosystems 9: 1–13.

HAIR, J. F., R. E. ANDERSON, R. L. TATHAM, AND W. C. BLACK.

1995. Multivariate data analysis with readings. Prentice Hall.

HA ̊ KANSON, L. 1996. Predicting important lake habitat variables

from maps using modern modeling tools. Can. J. Fish. Aquat.

Sci. 53(Suppl. 1): 364–382.

———, A. C. BLENCKNER, B. HELLSTROM, AND S. S. HELLSTROM.

2005. The influence of calcium on the chlorophyll–phosphorus

relationship and lake Secchi depths. Hydrobiologia 537:

111–123.

HARRIS, G. P., AND oTHERS. 2003. The role of models in ecosystem

management, p. 300–306. In C. D. Canham, J. J. Cole and W.

K. Lauenroth, [eds.], Models in ecosystem science. Princeton

University Press.

HEISKARY, S. A., AND C. B. WILSON. 2005. Minnesota lake water

quality assessment report: Developing nutrient criteria. 3rd

ed. Water Assessment and Environmental Information

Section Environmental Analysis & Outcomes Division

Watershed Section Regional Division. Minnesota Pollution

Control Agency. September.

HIGGINS, J. V., M. T. BRYER, M. L. KHOURY, AND T. W.

FITZHUGH. 2005. A freshwater classification approach for

biodiversity conservation planning. Conserv. Biol. 19:

432–445.

JEPPESEN, E., J. P. JENSEN, M. SØNDERGAARD, T. LAURIDSEN, AND

F. LANDKILDEHUS. 2000. Trophic structure, species richness

and biodiversity in Danish lakes: Changes along a phosphorus

gradient. Freshw. Biol. 45: 201–218.

JUBAR, A. K. 2004. Quantifying effects of residential lakeshore

development on littoral fishes and habitat: Toward a framework for lake ecosystem conservation. M.S. thesis. Michigan

State University.

KING, R. S., AND C. J. RICHARDSON. 2003. Integrating bioassessment and ecological risk assessment: An approach to

developing numerical water-quality criteria. Environ. Manage. 31: 795–809.

KNOLL, L. B. 2004. Understanding zebra mussel (Dreissena

polymorpha) impacts on autotrophic and heterotrophic

plankton of inland lake ecosystems. M.S. Thesis. Michigan

State University.

KNOWLTON, M. F., AND J. R. JONES. 2006. Natural variability in

lakes and reservoirs should be recognized in setting nutrient

criteria. Lake and Reserv. Manage. 22: 161–166.

LAMON, E. C., AND C. A. STOW. 2004. Bayesian methods for

regional-scale eutrophication models. Water Res. 38: 2764–

2774.

LOUGHEED, V. L., C. A. PARKER, AND R. J. STEVENSON. 2007. Using

nonlinear responses of multiple taxonomic groups to establish

criteria indicative of wetland biological condition. Wetlands

27: 96–109.

MARTIN, S. L., AND P. A. SORANNO. 2006. Defining lake landscape

position: relationships to hydrologic connectivity and landscape features. Limnol. Oceanogr. 51: 801–814.

NORVELL, W. A., C. R. FRINK, AND D. E. HILL. 1979. Phosphorus

in Connecticut lakes predicted by land use. Proc. Natl. Acad.

Sci. 76: 5426–5429.

786 Soranno et al.
QIAN, S. S., R. S. KING, AND C. J. RICHARDSON. 2003. Two

statistical methods for the detection of environmental thresholds. Ecol. Model. 166: 87–97.

R DEVELOPMENT CORE TEAM. 2007. A language and environment

for statistical computing. R Foundation for Statistical

Computing.

RECKHOW, K. H., AND oTHERS. 2005. A predictive approach to

nutrient criteria. Environ. Sci. Technol. 39: 2913–2919.

SMITH, V. H. 1998. Cultural eutrophication of inland, estuarine,

and coastal waters, pp. 7–49. In M. L. Pace, and P. M.

Groffman [eds.], Successes, limitations, and frontiers in

ecosystem science. Springer-Verlag.

SØNDERGAARD, M., E. JEPPESEN, J. P. JENSEN, AND S. LILDAL

AMSINCK. 2005. Water framework directive: Ecological

classification of Danish lakes. J. Appl. Ecol. 42: 616–629.

STEMBERGER, R. S., AND E. K. MILLER. 2003. Cladoceran body

length and Secchi disk transparency in northeastern U.S.

lakes. Can. J. Fish. Aquat. Sci. 60: 1477–1486.

STEVENSON, R. J., AND oTHERS. 2004. Interpreting results

of ecological assessments, p. 85–111. In M. T. Barbour, S.

B. Norton, H. R. Preston and K. W. Thornton [eds.],

Ecological assessment of aquatic resources: Linking science to

decision-making. Society of Environmental Toxicology and

Chemistry.

———, Y. PAN, AND P. VAITHIYANATHAN. 2002. Ecological

assessment and indicator development. in wetlands: the case

of algae in the Everglades, USA. Verh. Int. Ver. Theor.

Angew. Limnol. 28: 1248–1252.

STODDARD, J. L. 2004. Use of ecological regions in aquatic

assessments of ecological condition. Environ. Manage. 34:

S61–S70.

———, D. P. LARSEN, C. P. HAWKINS, R. K. JOHNSON, AND R. H.

NORRIS. 2006. Setting expectations for the ecological condition of streams: the concept of reference condition. Ecol.

Appl. 16: 1267–1276.

STOW, C. A., S. R. CARPENTER, AND R. C. LATHROP. 1997. A

Bayesian observation error model to predict cyanobacterial

biovolume from spring total phosphorus in Lake Mendota,

Wisconsin. Can. J. Fish. Aquat. Sci. 54: 464–473.

U.S. ENVIRONMENTAL PROTECTION AGENCY [USEPA]. 2000.

Ambient water quality criteria recommendations information

supporting the development of state and tribal nutrient

criteria: Lakes and reservoirs in nutrient ecoregion VI. Office

of Water. Washington, D.C.

VADEBONCOEUR, Y., E. JEPPESEN, M. J. VANDER ZANDEN, H.-H.

SCHIERUP, K. CHRISTOFFERSEN, AND D. M. LODGE. 2003. From

Greenland to green lakes: Cultural eutrophication and the loss

of benthic pathways in lakes. Limnol. Oceanogr. 48: 1408–1418.

WATSON, S. B., E. MCCAULEY, AND J. A. DOWNING. 1997. Patterns in

phytoplankton taxonomic composition across temperate lakes

of differing nutrient status. Limnol. Oceanogr. 42: 487–495.

———, ———, AND ———. 1992. Sigmoid relationships between

phosphorus, algal biomass, and algal community structure.

Can. J. Fish. Aquat. Sci. 49: 2605–2610.

WICKHAM, J. D., K. H. RIITTERS, T. G. WADE, AND K. B. JONES.

2005. Evaluating the relative roles of ecological regions and

land-cover composition for guiding establishment of nutrient

criteria. Landscape Ecol. 20: 791–798.

Received 23 October 2006

Accepted 23 September 2007

Amended 29 October 2007

Ecosystem-specific nutrient criteria 787