Quantifying gradients of anthropogenic stress can inform the development of sample designs, provide an important covariate in modeling relationships of response variables, identify reference and highly-disturbed sites, and provide a baseline and guidance to restoration and remediation efforts. We describe development of SumRel, a composite index of anthropogenic stress, for the U.S. and Canadian Lake Superior basin. Key elements of the project include development of high-resolution watersheds throughout the basin, summarization of the major point and non-point stressors within these watersheds, and creation of tools for scaling the watersheds and stressor summaries. SumRel was calculated at two spatial scales: for high resolution subcatchments within the Lake Superior basin (mean watershed area = 93 ha) and for coastal watersheds of Lake Superior. An assessment of subcatchments within Minnesota's St. Louis River watershed showed a correlation between the degree of disturbance, as indicated by SumRel, and impaired water quality, as evidenced by in-stream conductivity. These data and tools allow identification and visualization of reference and highly-disturbed sites at multiple spatial scales, providing decision support for individual agency and binational monitoring, assessment and restoration initiatives across the Lake Superior basin.

Introduction

Lake Superior, headwaters to the largest freshwater system in the world, faces increasing risk from human activities coupled with global climatic change. Human-induced stressors affecting Lake Superior are many, including chemical inputs from point and non-point sources, biological factors such as invasive species and the rapid spread of diseases, and physical changes ranging from shoreline alteration to effects of land use change in the watersheds. Not only do these stressors interact with one another, but they also operate under changing temperature and precipitation regimes, creating challenges to monitoring, assessment and remediation activities. In addition, the gradient of ecological condition across Lake Superior is large, ranging from the St. Louis River and Bay Area of Concern to relatively pristine waters that border sparsely populated boreal watersheds. As a result, finding the appropriate balance between protection and restoration efforts to sustain Lake Superior ecosystems presents a formidable challenge to resource managers.

One of the challenges for large lake systems is to develop monitoring and assessment programs that can (1) effectively identify time trends in habitat improvement or degradation, and (2) be coordinated across the multiple agencies responsible for habitat management (Lake Superior Binational Program 2006). A foremost issue in environmental monitoring is how to distribute resources to collect a limited number of samples that are truly representative of a target population and can detect trends in biotic and abiotic response variables. Sample designs, monitoring, and assessment can all be enhanced by understanding the gradient or range of environmental stressors impacting these units (Host et al., 2005; Danz et al., 2005). Understanding the spatial distribution and combined influence of human-related stressors is particularly challenging, in that such disturbances reflect the influence of multiple and often spatially correlated factors; these include a variety of point source discharges (e.g. U.S. Toxic Release Inventory, Canada National Pollution Release Inventory), stresses related to human populations (estimable from road and population densities), and non-point sources related to factors such as agriculture, forestry, or atmospheric deposition. This information can guide the allocation of restoration and protection efforts by (1) identifying the highly-disturbed sites in greatest need of restoration and/or mitigation, (2) characterizing reference or benchmark sites that can be used to develop restoration targets and goals (e.g. water quality standards, reference lists of plant species for a particular habitat type), (3) diagnosing the types of stress causing degradation at specific sites, which can guide decisions on remedial action necessary for protection or restoration. Also, stressors occurring upstream of a restoration project could threaten the success of that project (e.g. sediment inputs upstream of fish spawning habitat restoration work).

Multiple approaches have been used to identify environmental stressor gradients. On a global scale, (Vorosmarty et al., 2010) developed a cumulative threat index for river networks based on the distribution and intensity of 23 geospatial variables, mapped on hydrologic networks. Halpern et al. (2008) derived and applied ‘impact weights’, coefficients that describe the relative magnitude of effects, to a series of anthropogenic drivers to assess the state of global marine coastal ecosystems; they concluded that no coastal systems remain unimpacted. In the Great Lakes Environmental Indicators (GLEI) project (Niemi et al., 2007), Danz et al. (2005) assembled a suite of 207 spatial data layers classified according to seven categories of human disturbance or biophysical landscape attributes. Data were summarized using principal component analyses within and across stressor categories, followed by a cluster analysis based on the principal component scores. In a separate study to identify reference conditions (representations of minimally disturbed habitats), Host et al. (2005) derived the MaxRel metric, which consisted of the single maximum (relativized) score of any of five key stressor variables within a watershed: population density, road density, percent agriculture, percent urban, and point source density. The premise of the MaxRel approach was that, to identify the best reference conditions, all stressors should be held to their lowest levels, and that the dominant stressor was likely the one that was most limiting to the assemblages. In this paper we describe SumRel, a composite index in which component stressor values are summed to better reflect the full gradient of environmental stress. MaxRel works well for the purpose of reference area identification, but does not represent the cumulative influence of multiple stressors in more disturbed regions.

The maps, decision tools and data from this effort will permit resource managers and decision makers across the basin to make more informed decisions on prioritizing watersheds for monitoring and restoration efforts. Specific objectives for this project were to (1) derive and map the distribution of a stressor index to quantify the amount of human disturbance in Lake Superior watersheds defined at local and regional scales; (2) identify potential reference (least subject to stress) and highly-disturbed watersheds and coastal regions within the Lake Superior basin; and (3) develop tools that allow users to visualize environmental stressors at both broad and fine spatial scales.

Methods

Watershed delineation

A high resolution delineation of watersheds for the U.S. and Canadian sides of the Lake Superior basin was derived using ArcHydro, a data model developed by ESRI (Maidment and Morehouse, 2002). Delineations were based on a 10 m digital elevation model (DEM) for the U.S. side of the Lake Superior basin, and a 20 m DEM on the Canadian side. Drainage enforcement, the process of removing depressions with no apparent outlet from the DEM, was performed using stream data from the National Hydrologic Dataset for the U.S. portion of the basin and the Water Virtual Flow Seamless Provincial Data Set for the Canadian basin. Using flow direction and flow accumulation grids derived from elevation maps, stream networks were identified based on a minimum flow accumulation threshold, i.e. the minimum drainage area needed to form a stream. Once the stream networks were delineated, flow direction was used to delineate the contributing area or sub-catchment for each stream reach between stream confluences.

The ArcHydro model maintains hydrologic continuity by assigning a unique ‘Hydro-ID’ to each subcatchment and identifying the Hydro-ID of the receiving subcatchment immediately downstream (Hollenhorst et al., 2007). These attributes are transferred to the corresponding stream reach and pour points. The Hydro-ID makes it possible to accumulate information for each stream's flow as one progresses down the drainage network, thereby allowing area-weighted means of environmental stressor values associated with each subcatchment to be accumulated down the network.

The ArcHydro procedure resulted in the identification and delineation of approximately 131,000 subcatchments in the Lake Superior basin; the average size of each subcatchment was 93 ha (230 ac). Subcatchments were combined based on their Hydro-IDs to identify coastal watersheds, i.e. watersheds emptying into Lake Superior. Approximately 7,000 Lake Superior coastal watersheds (hereafter referred to as simply ‘watersheds’) and the adjacent coastal areas that drain directly into the lake (interfluves) were identified. The GIS shapefiles for the watersheds and subcatchments can be downloaded at www.nrri.umn.edu/lsgis2 or viewed interactively at http://gisdata.nrri.umn.edu/geomoose/GLNPO.html.

SumRel: Quantifying environmental stressors

Anthropogenic stress was quantified using a suite of U.S. and Canadian spatial databases (Table 1). All data were obtained from existing, publically-available data sources with well-established and independently approved federal, state, or provincial data sources with in-house quality assurance programs. These data were selected because they (1) provide a comprehensive coverage of a broad geographic region, (2) exist at appropriate temporal and spatial scales, and (3) have recognized impacts on the structure, function, and composition of the ecological communities that comprise the basin (Danz et al., 2005; Brazner et al., 2007b; Halpern et al., 2008).

United States land cover information was derived from the National Land Cover Dataset (Vogelmann et al., 1998); Canadian land cover data obtained from the Ontario Land Cover Dataset (Spectranalysis 2004). Both land cover datasets were derived from 30 m Landsat Thematic Mapper satellite data. Hollenhorst et al. (2011) describe the procedure used to amalgamate the Canadian and U.S. classifications; the latter was developed as part of basin-wide land use change study (Hollenhorst et al., 2011). Land cover data were used to calculate percent agricultural and residential land use by area for each subcatchment. Canadian and U.S. population densities were derived from the 2001 Census of Canada and the 2000 U.S. Census, respectively. Canadian Census divisions and U.S. Census blocks were gridded to 30 m pixels and summarized by subcatchment.

A road density index (weighted km road km−2 land area) was calculated from U.S. Census TIGER line files and the Ontario Ministry of Natural Resources (OMNR) Road Segment dataset. Roads were weighted based on size, with arterials, collectors, and local roads receiving lower weights than expressways and limited access highways (typically four-lane roads; Table 2), following Forman's assumption that the impact of roads is proportional to their size and amount of traffic (Forman, 2000). The road density index was calculated as total weighted road length/subcatchment area.

United States point source data were obtained from the National Pollutant Discharge Elimination System permit system database, and includes industrial, municipal, and other facilities that discharge pollutants into U.S. waters. Canadian point source data were obtained from Environment Canada's National Pollution Release Inventory. Point sources were weighted based on the number and types of stressors potentially resulting from these sources. Stressors within the point source coverage included sewage, pathogens, PAHs, solvents, nutrients, salts, and pharmaceuticals.

We evaluated a number of normalizing transformations for each variable, including log and arcsine transformations. The use of high-resolution subcatchments resulted in a large number of zeros (i.e. absence of the stressor) for many of the variables. The best results, in terms normalizing the data, were obtained using a log10 transformation of non-zero values. All variables were log10(x+min) transformed, where min was the smallest non-zero value observed for variable x. Each transformed (x′) variable was then standardized, (x′-μx′)/σx′, where μ and σ are the true population mean and standard deviation for x′, respectively. These standardized values (x) were normalized, (x-min)/(max-min), with min and max being the minimum and maximum for all x, respectively. Finally the five standardized and normalized x values (% agriculture,% residential, population, road and point source densities) for each watershed were summed, and the basinwide set of summed values was normalized again. This procedure ultimately gave a single number–SumRel–for each watershed. Danz et al. (2005) explored numerous methods and assumptions behind combining stressor layers, including the use of weighting schemes and covariance analysis. They found that the different schemes gave results that were nearly identical to the simpler additive model. The value of SumRel ranges from 0.0–1.0, with 1.0 representing the maximum composite stress within a geographic coverage of interest. This design allows stressor scores to be calculated at any watershed scale, ranging from the high-resolution subcatchments derived from ArcHydro to coastal watersheds of Lake Superior. They can also be derived for any geographic extent, from local watersheds to an ecoregion, lake, or basin.

In a related study, the SumRel gradient was used to stratify water quality sampling sites within the St. Louis River watershed (Axler et al., 2011) Water samples were collected throughout the season in summer 2010, and analyzed for a suite of water quality variables, including temperature, pH, electrical conductivity, and dissolved oxygen. These data were used to conduct a preliminary assessment of the relationship between SumRel and a typical environmental response variable.

Results and Discussion

Mapping reference and at-risk watersheds

The SumRel scores were used to generate maps of stressor intensity and distribution for subcatchments and tributary watersheds across the Lake Superior basin. At the scale of coastal watersheds for Lake Superior, the urban areas of Duluth and Thunder Bay and south shore of Lake Superior had the highest scores, whereas sparsely-populated Canadian watersheds had low SumRel scores. Islands, including Isle Royale, the Apostle, and those in Lake Nipigon had low values for road density, population, and agriculture, and consequently, low SumRel scores (Figure 1).

Reference areas, representations of the ‘least-disturbed’ watersheds or ecosystems, are typically defined as systems within the lowest 10th to 25th percentiles of the population with respect to levels of disturbance (Davis and Simon, 1995). Figure 2 shows the distribution of the ‘tails’ of the stressor gradient–reference and highly-disturbed sites- estimated using SumRel thresholds of 10, 20, and 30 percent. The 10 percent cutoff identifies Isle Royale, the Apostle Islands, and a number of Canadian coastal interfluves as reference sites, and urban sites in Duluth and Thunder Bay as the most highly disturbed. The 20 percent cutoff included Lake Nipigon islands and coast, and much of the Canadian north shore as reference, and added several urban watersheds as most highly disturbed. The 30 percent cutoff considerably increased the land area in both reference and highly-disturbed categories, including the large St. Louis River watershed and many south shore watersheds in the most highly disturbed category. Note that the ‘most highly disturbed’ sites identified above represent the tail of the distribution for the particular geographic extent. They are not necessarily biologically degraded; such an assessment requires an analysis of the environmental stress – biological response relationships pertinent to the study area (Niemi and McDonald, 2004).

At finer spatial scales, the hydrologic continuity within the ArcHydro subcatchments provides the ability to ‘accumulate’ stressors’ influence along the length of a stream network. This revealed cases in which nonmonotonic longitudinal patterns of stress occur. For example, in the large St. Louis River watershed, SumRel values were high in the upper reaches of the watershed, which include the extensively developed minelands of Minnesota's Iron Range (Figure 3). Lands downstream of the Iron Range comprise extensive forest and wetlands, with little agriculture or population development. Consequently, the cumulative stressor values downstream became ‘diluted’ as the contributing land area expands to include these more rural subcatchments. Accumulated stressor scores then increased with the inclusion of areas subject to increasing urbanization in the lower reaches of the St. Louis River.

The delineation of these fine-scale and non-linear patterns of environmental stress along the stream network can provide guidance to local watershed managers for targeting remediation efforts or identifying candidate areas for protection or restoration. For example, the numerous highly-disturbed watersheds in the northern portion of the St. Louis River watershed (Figure 3) are likely important to local water quality, and candidates for restoration activities. Similarly, the heavily-forested subcatchments along the St. Louis River likely contribute to improved water quality and are important areas for protection and use of Best Management Practices. The urbanized subcatchments adjacent to the St. Louis River estuary have direct impacts on Lake Superior, making them high-priority areas for remediation activities. Simulations of the expected effects of remedial activities in priority areas could be used to determine the relative benefit (i.e. reduction in Sumrel value) of investing in different types of restoration projects.

SumRel also provides a means of stratifying samples to capture the gradient of conditions across a watershed. To evaluate the degree to which the watershed attributes summarized by SumRel compare with field-collected data, we regressed SumRel scores against in-stream electrical conductivity (EC25, a surrogate for salt concentrations) measured in summer 2010 at 32 locations within the St. Louis River watershed (Figure 3). Sample sites were chosen to span the gradient of previously-determined SumRel scores (Axler et al., 2011). There was a significant positive relationship (R2 = 0.65, p < 0.01) between SumRel and EC25 (Figure 4). While there is a wide range of biological and water quality environmental indicators, each varying in their response to environmental stressors, this preliminary analysis demonstrates that the stressors summarized by SumRel do relate to measures of environmental quality.

Identifying the appropriate spatial extent of particular classes of stress is one of the key issues in using a stressor index to address watershed management issues. Across the Great Lakes basin, watersheds of lake Erie and Ontario basins have much higher stressor scores than those of the Lake Superior basin, confounding a basinwide interpretation of reference condition (Danz et al., 2005). Using the GLEI composite stressor index (Danz et al., 2005), (Brazner et al., 2007a) found that more variance in the scores of fish, bird and vegetation indicators could be accounted for by considering the Great Lake in which a sample was collected (Superior, Erie, etc.), rather than by accounting for the ecoregion or wetland type. The Lake Superior basin itself has such a broad range of watershed conditions that the appropriate balance of protection, restoration and remediation efforts will vary between its urban and boreal watersheds in a much different manner than if Sumrel scores are calculated for the entire Great Lakes basin. The ability to scale both the grain (minimum watershed size) and geographic extent of the stressor gradient provides managers with flexibility needed to hone in on local, regional or basinwide management issues.

The SumRel and the GLEI principal component approaches represent two extremes in complexity. The GLEI gradient was a ‘one-off’ analysis originally designed to guide site selection as part of a basin-wide indicator development study (Danz et al., 2005); the large number of input data layers (207) and effort involved in reconstructing them would make it impractical to repeat this approach at a future time. In contrast, SumRel is based on a relatively few data layers that are institutionally maintained and publically available. Consequently, the component layers can be re-assessed at time intervals associated with landscape monitoring, typically 5–10 years (Bourgeau-Chaves et al., 2008). Although site-specific SumRel scores from different time periods are not directly comparable (because of scaling), the ranking of watersheds relative to one another within time periods gradient are comparable. Understanding temporal changes in the intensity of individual stressors at a site or the relative ranking of watersheds would allow managers to identify areas at risk of degradation and prioritize restoration activities accordingly.

Visualizing the stressor gradient

The scale issues and sheer volumes of data available for systems as large, complex, and multijurisdictional as the Lake Superior basin present a challenge to making intelligent resource management decisions. To aid in interpreting the SumRel index and the distribution of stress across the Lake Superior basin, we developed an online interactive stressor viewer (ISV) application to allow users to explore and summarize stressor data The application uses GeoMOOSE, an open-source JavaScript framework, to allow a user to view basemaps (low and high resolution watersheds, streams, shorelines), along with individual and composite stressors (Figure 5). The user can select various layers to display on the map including streams, tributary or fine scale watersheds, as well as the SumRel index values for roads, development, population and other attributes. The ISV also provides a unique analysis tool – selecting an individual catchment retrieves information on the magnitude of stressors associated with that particular catchment. The tool identifies associated upstream and downstream catchments, summarizes stressors and generates a simple visualization comparing the magnitude of each stressor relative to other catchments within the Lake Superior basin (Figure 5).

Conclusions

The SumRel stressor metric identified in this study achieves several important objectives relevant to sample design, interpretations of ecological indicator data, reference area identification and long-term monitoring and assessment. The delineation of highly-resolved subcatchments and visualization tools enables managers and decision-makers to identify: (1) specific tributaries that may account for disturbances in the coastal and nearshore zone of the lakes, (2) specific locations within the tributary, and (3) specific stressor types that may potentially result in impairments to that part of the river system. Identifying the location and magnitude of point and nonpoint source stressors permits identification of reference or highly-disturbed conditions, which in turn can inform the process of prioritizing restoration and protection efforts. Furthermore, identification of ‘least impacted’ areas within a watershed can serve as a benchmark for restoration efforts. Lastly, the development of tools to identify the stress gradient over a user-specified region can inform the design of future monitoring, assessment and restoration programs.

Acknowledgements

Ms. Jane Reed created the website design for ‘Explore Lake Superior,’ and Mr. Gerald Sjerven served as website administrator. This work was funded in part by a grant from the U.S. EPA Great Lakes National Program Office, Grant number: GL00E28801-0. Although this research has been funded wholly or in part by the U.S. EPA, it has not been subjected to the agency's required peer and policy review and therefore does not necessarily reflect the views of the agency and no official endorsement should be inferred. Dr. Richard Axler generously shared water quality data collected in the St. Louis River watershed; funding for this analysis was provided by the Minnesota Pollution Control Agency. The integration of the U.S. and Canadian land use classifications was funded by grants from the National Fish and Wildlife Foundation, and Environment Canada through the Lake Erie Lakewide Area Management Plan to the University of Windsor and the Lake Erie Millennium Network. This is contribution number 517 from the Center for Water and the Environment, Natural Resources Research Institute, University of Minnesota-Duluth.

The text of this article is only available as a PDF.

References

Axler, R. P., Brady, V. J., Breneman, D. H. and Johnson, L. B.
2011
. “
Assessment of stream habitat, biological condition, and water quality of the St Louis River watershed
”.
Duluth, MN
:
Natural Resources Research Institute
.
NRRI/TR-2011/22
Bourgeau-Chaves, L., Lopez, R. D., Trebitz, A., Hollenhorst, T. P., Host, G. E., Huberty, B., Gauthier, R. L. and Hummer, J.
2008
. “
Landscape-based indicators
”. In
Great Lakes Wetland Coastal Monitoring Plan
, Edited by: Burton, T. A., Brazner, J. C., Ciborowski, J. J. H., Grabas, G. P., Hummer, J., Schneider, J. and Uzarski, D. G.
143
171
.
Ann Arbor
:
Great Lakes Commission
.
Brazner, J. C., Danz, N. P., Niemi, G. J., Regal, R. R., Trebitz, A. S., Howe, R. W., Hanowski, J. M., Johnson, L. B., Ciborowski, J. J.H., Johnston, C. A., Reavie, E. D. and Brady, V. J.
2007a
.
Evaluation of geographic, geomorphic, and human influences on Great Lakes wetland indicators. A multi-assemblage approach
.
Ecological Indicators
,
7
:
610
635
.
Brazner, J. C., Danz, N. P., Trebitz, A. S., Niemi, G. J., Regal, R. R., Hollenhorst, T. P., Host, G. E., Reavie, E. D., Brown, T. N., Hanowski, J. M., Johnston, C. A., Johnson, L. B., Howe, R. W. and Ciborowski, J. J.H.
2007b
.
Responsiveness of Great Lakes Wetland Indicators to Human Disturbances at Multiple Spatial Scales: A Multi-Assemblage Assessment
.
Journal of Great Lakes Research
,
33
(
3
):
42
66
.
Danz, N. P., Niemi, G. J., Regal, R. R., Hollenhorst, T. P., Johnson, L. B., Hanowski, J. M., Axler, R. P., Ciborowski, J. J.H., Hrabik, T., Brady, V. J., Kelly, J. R., Morrice, J. A., Brazner, J. C., Howe, R. W., Johnston, C. A. and Host, G. E.
2005
.
Integrated gradients of anthropogenic stress in the U.S. Great Lakes basin
.
Environmental Management
,
39
:
631
647
.
Davis, W. S. and Simon, T. P.
1995
. “
Biological Assessment and Criteria: Tools for Water Resource Planning and Decision Making
”.
Boca Raton
:
Lewis Publishers
.
Forman, R. T.T.
2000
.
Estimate of the area affected ecologically by the road system in the United States
.
Conservation Biology
,
14
(
1
):
31
35
.
Halpern, B. S., Walbridge, S., Selkoe, K. A., Kappel, C. V., Micheli, F., D’Agrosa, C., Bruno, J. F., Casey, K. S., Ebert, C., Fox, H. E., Fujita, R., Heinemann, D., Lenihan, H. S., Madin, E. M.P., Perry, M. T., Selig, E. R., Spalding, M., Steneck, R. and Watson, R.
2008
.
A Global Map of Human Impact on Marine Ecosystems
.
Science
,
319
(
5865
):
948
952
.
Hollenhorst, T. P., Brown, T. N., Johnson, L. B., Ciborowski, J. J.H. and Host, G. E.
2007
.
Methods for generating multi-scale watershed delineations for indicator development in Great Lakes coastal ecosystems
.
Journal of Great Lakes Research
,
33
(
3
):
13
26
.
Hollenhorst, T. P., Johnson, L. B. and Ciborowski, J. J.H.
2011
.
Monitoring land cover change in the Lake Superior basin
.
Aquatic and Ecosystem Health and Management Society
,
14
(
4
):
433
442
.
Host, G. E., Schuldt, J., Johnson, L. B., Hollenhorst, T. and Richards, C.
2005
.
Use of GIS and remotely-sensed data for a priori identification of reference areas for Great Lakes coastal ecosystems
.
International Journal of Remote Sensing
,
26
:
5325
5342
.
Lake Superior Binational Program.
2006
. “
Chapter 3: Ecosystem Goals, Objectives, Indicators, Monitoring, and Beneficial Use Impairments
”. Vol.
2006
,
1
25
.
Lake Superior Lakewide Management Plan (LaMP)
.
Maidment, D. R. and Morehouse, S.
2002
.
ArcHydro: GIS for Water Resources
,
Redlands, CA
:
ESRI Press
.
Niemi, G. J. and McDonald, M. E.
2004
.
Applications of ecological indicators
.
Annual Review of Ecology and Systematics
,
35
:
89
111
.
Niemi, G. J., Kelly, J. R. and Danz, N. P.
2007
.
Environmental Indicators for the Coastal Region of the North American Great Lakes: Introduction and Prospectus
.
Journal of Great Lakes Research
,
33
(
Supplement 3
):
1
12
.
Spectranalysis, I.
2004
. “
Introduction to the Ontario Land Cover Database, 2nd edition: Outline of production methodology and description of 27 land cover classes
”.
Ontario, CA
Report to the Ontario Ministry of Natural Resources edition
Vogelmann, J. E., Sohl, T. L., Campbell, P. V. and Shaw, D. M.
1998
.
Regional Land Cover Characterization Using Landsat Thematic Mapper Data and Ancillary Data Sources
.
Environmental Monitoring and Assessment
,
51
(
1–2
):
415
428
.
Vorosmarty, C. J., McIntyre, P. B., Gessner, M. O., Dudgeon, D., Prusevich, A., Green, P., Glidden, S., Bunn, S. E., Sullivan, C. A., Liermann, C. R. and Davies, P. M.
2010
.
Global threats to human water security and river biodiversity
.
Nature
,
467
(
7315
):
555
561
.