# EMERGING FRONTIERS IN ECOLOGICAL STOICHIOMETRY

EDITED BY : Michelle Evans-White and James Joseph Elser PUBLISHED IN : Frontiers in Ecology and Evolution

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers.

> The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-294-7 DOI 10.3389/978-2-88963-294-7

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# EMERGING FRONTIERS IN ECOLOGICAL STOICHIOMETRY

Topic Editors:

Michelle Evans-White, University of Arkansas, United States James Joseph Elser, University of Montana, United States

Citation: Evans-White, M., Elser, J. J., eds. (2019). Emerging Frontiers in Ecological Stoichiometry. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-294-7

# Table of Contents


Kimberley D. Lemmen, Orpheus M. Butler, Thomas Koffel, Seth M. Rudman and Celia C. Symons


Ashley D. Keiser

*73 Elemental Ratios Link Environmental Change and Human Health* Rachel E. Paseka, Anika R. Bratt, Keeley L. MacNeill, Alfred Burian and Craig R. See

# Editorial: Emerging Frontiers in Ecological Stoichiometry

Michelle A. Evans-White<sup>1</sup> \*, Zoe G. Cardon<sup>2</sup> , Jennifer A. Schweitzer <sup>3</sup> , Jotaro Urabe<sup>4</sup> and James J. Elser <sup>5</sup>

*<sup>1</sup> Department of Biological Sciences, University of Arkansas, Fayetteville, AR, United States, <sup>2</sup> Ecosystems Center, Marine Biological Laboratory, Woods Hole, MA, United States, <sup>3</sup> Ecology and Evolutionary Biology, University of Tennessee, Knoxville, Knoxville, TN, United States, <sup>4</sup> Graduate School of Life Sciences, Tohoku University, Sendai, Japan, <sup>5</sup> Flathead Lake Biological Station, University of Montana, Polson, MT, United States*

Keywords: ecological stoichiometry, biological stoichiometry, global change, human health, cryosphere, eco-evo, silica, soil organic matter

**Editorial on the Research Topic**

**Emerging Frontiers in Ecological Stoichiometry**

### WHAT IS STOICHIOMETRY AND WHY IS IT RELEVANT FOR GLOBAL CHANGE?

The ecological (ES) and biological stoichiometry (BS; Elser et al., 2000) theoretical frameworks apply mass balance of energy and multiple chemical elements to fluxes at multiple scales of biological organization and integrate them into an evolutionary framework. The unprecedented rate of global change (GC) paired with future models of change that may occur unless anthropogenic impacts are mitigated (Easterling et al., 2000; Coumou and Rahmstorf, 2012; IPCC, 2018) necessitates rapid global human innovation and social change to avoid various emerging perils (Figueres et al., 2017; IPCC, 2018; UNCC, 2019). The ES/BS framework can represent a powerful tool to drive the innovation needed to address these GC perils (Urabe et al., 2010).

#### Edited and reviewed by:

*Elise Huchard, UMR5554 Institut des Sciences de l'Evolution de Montpellier (ISEM), France*

> \*Correspondence: *Michelle A. Evans-White mevanswh@gmail.com*

#### Specialty section:

*This article was submitted to Behavioral and Evolutionary Ecology, a section of the journal Frontiers in Ecology and Evolution*

> Received: *11 October 2019* Accepted: *18 November 2019* Published: *05 December 2019*

#### Citation:

*Evans-White MA, Cardon ZG, Schweitzer JA, Urabe J and Elser JJ (2019) Editorial: Emerging Frontiers in Ecological Stoichiometry. Front. Ecol. Evol. 7:463. doi: 10.3389/fevo.2019.00463*

#### WOODSTOICH 15TH AND WOODSTOCK 50TH-TWO ANNIVERSARIES OF YOUTHFUL INNOVATION

Within this perilous global context and on the 15 and 50th anniversary of the first Woodstoich Workshop (WW; August 2004) and the famous Woodstock festival, respectively, the 4th WW convened from 15 to 19 August 2019 at Flathead Lake Biological Station in Montana, USA. The WWs aim to produce a collection of publications that innovate and expand the use of ES/BS theory. The WWs take inspiration from the Woodstock festival (August 15–17 1969) where young people caught up in a time of war and turbulent politics came together in a short, energetic musical event espousing unity, peace, and love to create a lasting legacy. Three previous WW events have taken place (Hessen and Elser, 2005; Urabe et al., 2010; Sterner et al., 2015), but the sense of urgency for unity and innovative solutions to address GC issues facing the young Woodstoich participants has never been greater than in 2019 (e.g., Hagedorn et al., 2019).

# WHAT IS THE WOODSTOICH CONTRIBUTION?

Woodstoich speeds the pace of development and communication of scientific ideas; international, collaborative ties among early-career scientists (graduate students, post-docs, pre-tenured professors) with common research interests are established, and the participants' innovative ideas are transformed into a peer-reviewed publication in less than 6–7 months. Five group leaders were chosen in December 2018 based on manuscript proposals submitted with their applications. Leaders then selected team members by April 2019 based on applications that demonstrated a combination of common interest and the skill set required to complete the proposed manuscript. Groups worked on manuscripts remotely using shared, cloud-based applications for ∼4 months prior to the workshop and then, in keeping with the Woodstock concept, these innovative, young minds all came together for an intellectual jam session in the form of the 6-days workshop in August 2019. During that short time, groups networked, interacted with a few senior scientist "mentorsages" invited by group-leader consensus, and then finalized and submitted manuscripts for a pre-arranged, extremely-rapid review; a 24-h turnaround was requested from three outside reviewers recruited by Associate Editors prior to the event. Groups responded to reviews during the remainder of the workshop with the goal of obtaining provisional acceptance of manuscripts prior to departing for their home institutions.

Additionally, innovation is stimulated by bringing together young minds from diverse backgrounds whose mindsets are not locked into place by long-term tradition or professional relationships and by encouraging them to embrace the unique spirit of the workshop. Contributing to that spirit are the Three Rules of Woodstoich, about which all participants were briefed in a kickoff session. Rule 1: Obey the law. Which law? The Law of Conservation of Matter applied across multiple chemical elements. Development of stoichiometric theory has revealed significant implications of constraining ecological dynamics within such bounds, enabling insight into the feasibility envelope for eco-evolutionary systems. Rule 2: Question authority. While senior experts on ES/BS were in attendance at the event (and indeed lurking in the reviewer pool), participants were strongly encouraged to hold ES/BS paradigms up to the light and challenge them as necessary. Rule 3: Make it groovy. All were urged to follow the spirit of Woodstock, to do their best to say something different, to be bold, and to try to forge a new path of innovation on their paper's topic. Here we summarize the papers produced under this charge.

# Linking Genetic Variation and Stoichiometric Traits

Lemmen et al. examines potential importance of genetic variations in eco-evolutionary dynamics using existing data of common garden experiments for various organisms. They found that genetic variations in some of stoichiometric traits were substantial, suggesting that single species have various options in a way of mass balance for maintaining their elemental composition. These findings suggest a new direction for understanding the role of ES in eco-evolutionary dynamics.

### Ecological Stoichiometry of the Mountain Cryosphere

Ren et al. reviews nutrient and trophic stoichiometry in mountain glacial ecosystems and finds strong potential for nutrient limitation, especially by phosphorus, that is likely transmitted through the limited, but surprising, food webs of organisms living on and near glaciers. Mountain glaciers are less nutrient-limited, however, than the better-known ice sheets of Antarctica, making predictions of one landscape to another inadvisable. In a warming world with increased ice melting, nitrogen and phosphorus limitations are predicted to increase with implications for downslope organisms, including humans.

# Elemental Ratios Link Environmental Change and Human Health

Paseka et al. leaps across disciplinary boundaries to argue for the importance of a stoichiometric perspective in identifying the roots of a range of issues for human health. Building from a handful of examples in the literature where stoichiometric ratios of elements have clear effects on the health of human or animal models, Paseka et al. point to further links between ecological stoichiometry and pressing issues ranging from food quality and water security to human susceptibility to non-infectious and infectious disease.

## Stoichiometry in Soil Organic Matter Models

Buchkowski et al. presents a review and new soil organic matter (SOM) model that uses stoichiometry of plants and microorganisms to show how incorporation of mass balance information (including ratios of carbon, nitrogen, phosphorus, and sulfur) constrain SOM models to probable fluxes across landscapes, a perspective that is missing in current SOM models. A stoichiometric approach to SOM models adds boundaries to flux estimates that provides critical realism to models of carbon stability that substantially reduce uncertainty and will allow for more robust carbon accounting in a changing world.

# Silica Stoichiometry in a Large Floodplain River

In a powerful new analysis of diverse components of the Upper Mississippi River System, Carey et al. identifies dominant controls and constraints over silicon:nitrogen:phosphorus (Si:N:P) ratios that are fundamental to determining phytoplankton assemblages in fresh and marine receiving waters. Distinct drivers, from backwater residence time to north-south gradients in lithology and human land use, set the stage for shifting Si:N:P ratios, which, in turn, have the potential to result in cascading effects on food-web structure and carbon sequestration.

# CONCLUSION AND PERSPECTIVES

In 2054, Woodstoich may mark its own 50th anniversary (WW10). Likely its founders will have passed from the earthly scene by then and time will tell if the concept of Woodstoich or even of ES will persist for that long. Let's imagine that they do and that WW10 organizers do a better job than the Woodstock 50th concert organizers (who canceled the anniversary concert). What will the world of 2050 be like and how would the WW4 have to be re-written at WW10 to account for the world of 2050?

Thinking about Lemmens et al., we would ask: how much evolution will have taken place? Was it fast enough to accommodate the rapid biogeochemical and climatic changes underway? How much more will we be able to discover about that evolution given 30 more years of innovation in genomic tools? Will a re-write of Ren et al. have to account for massive global glacier loss if climate targets are not met? Or will the paper be describing the biogeochemical and ecological dynamics of reduced but stabilized glaciers in mountain regions? Of great relevance to a new version of Paseka et al. will be what emissions pathway we follow during the next 30 years. If we follow the high emissions pathway (Riahi et al., 2011), pCO<sup>2</sup> will be 550 ppm. What will this imply for the nutritional quality of crops? Will fertilization rates have been amplified in order to lower C:N and C:P and C:essential metal ratios, with implications for water quality and its associated health impacts? Soils hold tremendous potential to sequester (or release) carbon during the coming decades. Will stoichiometric approaches such as those described by Buchowski et al. have improved our ability to understand and predict those changes? Finally, under current projections and trajectories, many regions of the world in 2050 will be experiencing large shifts in precipitation regimes. Thinking of Carey et al., how will this affect weathering processes, riverfloodplain interactions, and diffuse nutrient inputs from the expanded agricultural operations of 2050? With publication of

#### REFERENCES


this Research Topic we call on scientists of 2054 to remember these efforts and keep the spirit of Woodstock and Woodstoich alive so that innovation can continue to navigate the perils of a changing planet and a changing society.

# AUTHOR CONTRIBUTIONS

ME-W and JE organized the Woodstoich IV Workshop. ME-W contributed initial editorial draft, and finalized the editorial manuscript. ZC, JS, JU, and JE revised the material and contributed working group summaries.

# FUNDING

The Woodstoich workshop was supported by the US National Science Foundation (DEB-1840408).

#### ACKNOWLEDGMENTS

We are grateful to the Woodstoich 4 participants for their efforts and good spirits during the event. We thank Joseph Vanderwall and all the staff of the Flathead Lake Biological Station for their efforts in supporting the workshop. We are especially grateful to the 15 external reviewers who returned their reviews within the requested 24-h turnaround period.

Gomis, E. Lonnoy, T. Maycock, M. Tignor, and T. Waterfield (Geneva: World Meteorological Organization), 32.


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Evans-White, Cardon, Schweitzer, Urabe and Elser. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Stoichiometric Traits Vary Widely Within Species: A Meta-Analysis of Common Garden Experiments

Kimberley D. Lemmen<sup>1</sup> \*, Orpheus M. Butler <sup>2</sup> , Thomas Koffel 3,4, Seth M. Rudman<sup>5</sup> and Celia C. Symons 6,7

<sup>1</sup> Department of Aquatic Ecology, Netherlands Institute of Ecology (NIOO-KNAW), Wageningen, Netherlands, <sup>2</sup> Smithsonian Tropical Research Institute, Ancón, Panama, <sup>3</sup> Will Keith Kellogg Biological Station, Michigan State University, Hickory Corners, MI, United States, <sup>4</sup> Program in Ecology, Evolutionary Biology and Behavior, Department of Plant Biology, Michigan State University, East Lansing, MI, United States, <sup>5</sup> Department of Biology, University of Pennsylvania, Philadelphia, PA, United States, <sup>6</sup> Department of Ecology and Evolutionary Biology, University of California, Irvine, Irvine, CA, United States, <sup>7</sup> Department of Ecology and Evolutionary Biology, University of California, Santa Cruz, Santa Cruz, CA, United States

#### Edited by:

Michelle Evans-White, University of Arkansas, United States

#### Reviewed by:

Halvor Matthew Halvorson, University of Central Arkansas, United States Puni Jeyasingh, Oklahoma State University, United States Esteban Balseiro, National University of Comahue, Argentina

\*Correspondence:

Kimberley D. Lemmen lemmen.kimberley@gmail.com

#### Specialty section:

This article was submitted to Behavioral and Evolutionary Ecology, a section of the journal Frontiers in Ecology and Evolution

> Received: 22 July 2019 Accepted: 23 August 2019 Published: 20 September 2019

#### Citation:

Lemmen KD, Butler OM, Koffel T, Rudman SM and Symons CC (2019) Stoichiometric Traits Vary Widely Within Species: A Meta-Analysis of Common Garden Experiments. Front. Ecol. Evol. 7:339. doi: 10.3389/fevo.2019.00339 "Ecological stoichiometry," a framework that focuses explicitly on the balances and flows of chemical elements within and between organisms and ecosystems, has provided crucial insights into many biological patterns and processes. Despite the proliferation of stoichiometrically-focused studies in recent decades and recognition of the potential for rapid evolution of stoichiometric traits, the prevalence of genetic variation in stoichiometric traits within species remains unclear. We compiled data from 30 published common garden studies of a broad range of taxa (including invertebrates, vertebrates, and autotrophs) to examine how genetic variation influences the acquisition, assimilation, allocation (AAA), composition, and excretion of elements. To quantify the extent of genetic variation for a given trait we calculated the absolute mean response ratio from pairwise comparisons of populations within the same common garden (820 population and 708 genotype comparisons). We observed substantial intraspecific variation of stoichiometric traits across populations and among genotypes; however, the magnitude of variation was greater in AAA traits (effect sizes of 20 and 164% for population and genotype contrasts, respectively) and excretion (effect sizes of 52 and 23%) than in content of carbon (2.1 and 3.1%) and nitrogen (4.5 and 24%). These results suggest that the content of some elements may be evolutionarily constrained relative to AAA traits that determine the processing of these elements, and that a sole focus on elemental content would underestimate the importance of intraspecific genetic variation, particularly within populations. Across many trait types the variation was greater among genotypes within a population than across populations. Finally, we compared pairs of populations from environments with different phosphorus (P) availability to pairs of populations with similar P availability. Genetic variation in the traits measured was similar regardless of the P environment from which genotypes were isolated, suggesting that differences in elemental availability across environments do not necessarily drive enhanced trait divergence. Overall, our results highlight the substantial amount of intraspecific variation in stoichiometric traits and underscore the potential importance of intraspecific variation in driving ecological and evolutionary processes.

Keywords: elemental phenotype, evolution, eco-evolutionary dynamics, phosphorus, organismal stoichiometry, nutrient excretion, intraspecific variation

# INTRODUCTION

Ecological stoichiometry (ES) is a scientific framework that views living systems as composite chemical reactions which, like all physical processes, are governed by the law of mass balance (Sterner and Elser, 2002; Sterner et al., 2015). In doing so, ES focuses explicitly on the balance and flow of elements and energy within and between organisms and ecosystems. The ES framework has applicability across a wide range of disciplines, from astrobiology to cancer research (Elser, 2003; Elser et al., 2007), and can be applied across multiple levels of biological organization and across taxa. It has also been argued that explicit consideration of an organism's "stoichiometric traits" (i.e., elemental contents and traits that influence the flux of elements between an organism and its environment) is a natural and convenient way to investigate the complex interplay between ecology and evolution (Kay et al., 2005). This is because elements can be traced through space and time, such that genetic variation in stoichiometric traits can be linked mechanistically to variation in environmental elemental availability (Matthews et al., 2011; Leal et al., 2017b).

Despite the apparent potential for stoichiometric traits to link ecology and evolution, there remain some fundamental gaps in our understanding of the genetic basis of variation in stoichiometric traits. In particular, while the composition, acquisition (A), assimilation (A), allocation (A), and excretion (E) of elements (hereafter referred to as "elemental composition" and "AAA," or "AAAE") are regarded as the defining traits of the stoichiometric phenotype (Jeyasingh et al., 2014), it remains unclear how and to what extent genetic factors actually contribute to their variation within and between natural populations of conspecific organisms. Heritable trait variation is central to evolutionary processes, and high levels of intraspecific genetic trait variation within populations can be thought of as "fuel" for adaptation (Barrett and Schluter, 2008). At the same time, intraspecific genetic variation in stoichiometric traits has been identified as an important driver of variation in local ecological processes, such as leaf litter decomposition, that can scale up to affect overall ecosystem functioning (Whitham et al., 2006; Silfver et al., 2007). Quantification of the genetic component of intraspecific variation in stoichiometric traits can, therefore, help to identify which traits have the potential for contemporary evolution, and would represent an essential step toward the development of a stoichiometrically-explicit ecoevolutionary framework.

Intraspecific genetic variation in stoichiometric traits might emerge as a direct response to nutritionally-based or other selection pressures (Leal et al., 2017b), or because the stoichiometric traits are mechanistically correlated with other physiological, life history, and morphological traits, as in the well-documented coupling between P content and individual growth rate (the Growth Rate Hypothesis or "GRH"; see Elser et al., 2000a, 2003, 2008). However, stoichiometric traits also exhibit considerable plasticity, i.e., variability in response to environmental parameters (i.e., Townsend et al., 2007; González et al., 2011). For example, plant foliar nutrient content and resorption are often strongly correlated with soil nutrient availability (Schade et al., 2003; Rejmánková, 2005), and the elemental composition of invertebrates can be highly sensitive to the elemental composition of their diets (e.g., Visanuvimol and Bertram, 2011; Zhou et al., 2018). Thus, in many prior studies that have considered intraspecific variation in stoichiometric traits, the respective genetic and environmental drivers of intraspecific variation have likely been confounded. However, as the field of ES has matured, an increasing number of studies have conducted "common garden experiments" that can be used to assess the magnitude of intraspecific genetic variation in stoichiometric traits. These experiments minimize the effects of the environmental factors by comparing genetically distinct individuals within a common "garden" environment, thus allowing for the quantification of genetic variation either within or between populations (Lynch and Walsh, 1998). Synthesis of such studies could illuminate whether genetic variation is present, the relative amount of genetic variation within and across populations, which stoichiometric traits show genetic variation, and whether elemental limitation could influence evolution in nature (Elser, 2006; Leal et al., 2017a; Rudman et al., 2019).

To determine the amount of intraspecific genetic variation in ES traits, we conducted a meta-analysis of studies that have measured ES traits on individuals reared in common gardens. These studies typically fall into two broad categories: those that compare ES traits between pairs of populations, and studies comparing ES traits of different genotypes from within a single population. From these studies we collected data on intraspecific genetic variation in ES traits, which we defined as measures of elemental content or movement, and life history data (specific growth rate, size/mass). We used these data to answer four questions about the extent of genetic variation in ES traits. (1) How much intraspecific genetic variation is present in ES traits relative to key life history traits? Here we used somatic growth rate and morphological characteristics to represent two key traits that are not explicitly stoichiometric, but have often been measured concurrently with ES traits. We predicted similar levels of genetic variation in life history traits and ES traits because all traits, at some level, have an elemental basis. (2) Do particular classes of ES traits show more genetic trait variation than others? We predicted that elemental composition would show reduced variation relative to AAAE traits, as there are well-documented constraints in elemental composition that are maintained by stoichiometric homeostasis (Sterner and Elser, 2002). (3) Is the magnitude of intraspecific variation among genotypes (i.e., within a population) similar to that observed across populations? Individuals from different populations are often reproductively isolated and likely experience a more divergent environment than genotypes from a single population (Schluter, 2000). As such, we predicted a greater magnitude of genetic trait variation between distinct populations than among genotypes from within a population. (4) How does the magnitude of intraspecific genetic variation differ when comparing populations from distinct and similar stoichiometric environments? As resource stoichiometry has been demonstrated to be an important selection pressure (Jeyasingh and Weider, 2007), we predict that populations from distinct stoichiometric environments will show greater genetic variation than populations from similar environments.

#### MATERIALS AND METHODS

#### Data Collection

We collected studies for the meta-analyses by searching the ISI Web of Science (v. 5.32) and Google Scholar databases with "title" or "abstract" search terms "('ecological stoichiometry' OR 'biological stoichiometry') AND ('intraspecific variation' OR 'evolution') AND 'common garden"'. Searches were conducted between the 15–17th of May 2019. These search terms produced 4 results in Web of Science and 197 results in Google Scholar. We screened these 201 studies for suitability based on their titles, abstracts, methods, and results. We then examined the reference lists of the seemingly-appropriate articles for additional studies that did not appear in our search results. We excluded studies that did not include data on ES traits. Thirty studies were eventually selected via this process. Importantly, this was not an exhaustive collation of the studies that might have met our criteria. This is presumably a result of our selection of search terms, which were necessarily broad in order to capture the full range of taxa, study systems, and ES traits throughout the literature without biasing our search toward particular combinations of these study foci. Thus, our collection of studies is perhaps best thought of as a large, reproducible sample of the literature with a bias toward studies conducted within the ES framework. We also note here that we explicitly excluded studies of genetic variation that was created by artificial selection (i.e., the crop and livestock literature) because the amount of genetic variation across treatments in artificial selection experiments will largely be proportional to the amount of standing variation included in the study and the strength and duration of artificial selection.

During study selection, it became clear that studies generally fell into two broad categories. Specifically, some studies reported stoichiometric trait variation within populations (i.e., variation among distinct genotypes), while others reported variation among distinct populations (i.e., variation in population trait means). We assigned each study to one of these categories (**Table 1**), and the categories were used as fixed effects variables as part of the meta-analysis.

For each study we recorded details of location, ecosystem type, focal taxa, the purported reasons for differentiation between populations and all measured stoichiometric traits. For a given trait to be included in our analyses, it was necessary that at least two studies reported common garden data about that trait. We considered traits to be explicitly "stoichiometric" if they were part of the elemental phenotype defined by Jeyasingh et al. (2014) as elemental composition, acquisition, assimilation, allocation, and excretion (CAAAE). Many organismal traits, including growth rate and body size, are often closely coupled with CAAAE traits in such ways that could justify their categorization as a "stoichiometric trait" if defined broadly. By adopting a strict definition of stoichiometric trait, we are able to qualitatively compare the intraspecific genetic variation of unambiguously stoichiometric traits with life history traits that do not explicitly involve the movement of elements.

Due to differences in how often CAAAE traits have been measured we found it necessary to subdivide or group certain traits. We considered the content of carbon (C), nitrogen (N) and phosphorus (P) and their ratios as separate traits due to the abundant reporting of elemental composition. Few studies measured acquisition, assimilation, or allocation, therefore we grouped these aspects of the elemental phenotype into a single trait category (AAA) for all elements. We also recorded somatic growth rate and measures of mass and size as "morphological traits" because they were the most commonly measured life history traits for comparison. Studies on invertebrates and vertebrates measured and reported elemental composition of whole organisms, while studies on plants used leaf tissues. In studies on plants, we included data from studies that measured either green or senesced vegetation as elemental content and data on the decomposition rate of plant material, which has modest heritability (Rodriguez-Cabal et al., 2017), was included as excretion (**Table 1**).

Means, variances, and sample sizes for traits were extracted from repositories, tables, figures, or text. Standard deviation was estimated from SE, 95% confidence intervals or quantiles using standard methods (Higgins and Green, 2008; Bland, 2015). We used DataThief III (Tummers, 2006) to extract data from figures. Where data were not provided in a format suitable for our analyses, we contacted study authors via email to request the relevant data. In cases where study designs featured more than one type of common garden (e.g., Dinh Van et al., 2013), population contrasts within different common gardens were treated as unique observations. For example, if there were two common gardens, and population A and population B were both raised in both common gardens, then there were two population contrasts included in our meta-analysis dataset, one for each common garden.

#### Effect Size Calculation

For each trait within a study, we compared the difference in genotype/population's trait values by calculating the log response ratios (LRR) for all pairwise combinations of genotypes/populations grown in the same common garden.

$$LRR = \wp\_i = \ln\left(\frac{\overline{Trait}\_1}{\overline{Trait}\_2}\right) \tag{1}$$

$$\text{variance } LRR = \left. \nu\_i = \frac{\left\| SD\_1 \right\|^2}{n\_1 \left(\overline{Train\_1}\right)^2} + \frac{\left\| SD\_2 \right\|^2}{n\_2 \left(\overline{Trait\_2}\right)^2} \tag{2}$$

Where i is the study, Trait<sup>j</sup> is the genotype/population mean, SD<sup>j</sup> the standard deviation and n<sup>j</sup> is the sample size (number of replicates), with j = 1, 2 corresponding to the two genotypes/populations being compared. We choose to use LRR because (i) it is dimensionless, allowing for comparison between studies with different methodologies and units, (ii) it is easily interpretable as it quantifies the proportionate change (Hedges et al., 1999), and (iii) the value of the LRR is not affected by non-independent samples (Noble et al., 2017). We TABLE 1 | A complete list of the studies included in the meta-analysis, with information on study type—either populations ("P") or genotypes ("G"); the number of populations/genotypes ("# Pop/Geno"), the type of ecosystem considered (Terrestrial, Freshwater or Marine); the coarse taxonomic group of the study organism (Plant, Invertebrate, or Vertebrate); the Latin name of the study organism; for populations, if they were from biogeochemically distinct environments (always "no" for genotypes); the number of common gardens in the experimental setup ("# CGs"); the list of traits that were measured in the study ("Measured Trait"), and; the trait category under which they were grouped in the meta-analysis ("Response Variable," with "SGR" standing for Somatic Growth Rate).


(Continued)

#### TABLE 1 | Continued


(Continued)


#### TABLE 1 | Continued

defined n<sup>j</sup> differently for studies comparing genotypes and those comparing populations. For genotype studies, n<sup>j</sup> was the number of individuals measured for a given genotype. For population studies, n<sup>j</sup> was the number of genetically distinct individuals measured for a given population. For example, in a sexually reproducing species n<sup>j</sup> is the number of individuals measured, but in studies of cyclical parthenogens (e.g., Daphnia) n<sup>j</sup> is the number of distinct clone lines measured. If it was not clearly indicated that genetically different individuals were included in the study, then n<sup>j</sup> was recorded as 1.

In the set of studies that compared traits of different genotypes from the same population, we were not always able to use all pairwise comparisons to calculate LRRs. This group of studies often had a larger number of genotypes to be compared than studies comparing distinct populations (maximum number of genotypes compared was 39, and populations was 9), and the use of all pairwise comparisons was too computationally demanding. Instead, if a common garden measured <5 genotypes we calculated the LRR for all pairwise combinations. However, if six or more genotypes were measured (six studies), pairwise comparisons were performed for only five genotypes (min, first quartile, median, third quartile, max). With this method we took a subset of the data to be computationally reasonable and to represent the full range of trait values. We also used a conservative approach and calculated all LRRs within a common garden using the median trait value as the denominator in Equation 1 (see **Supplementary Material**). Although the magnitude and 95% credibility interval of some mean effect estimates of these two methods differed, it did not change the interpretation of the data. In the studies that compared mean traits between populations, the maximal number of distinct populations was small enough to enable us to compute all possible pairwise comparison, even for the two studies that had more than five distinct populations.

#### Global Meta-Analysis

We were interested in the magnitude, rather than the direction of, trait differences across populations and among genotypes; therefore, we analyzed the absolute values of LRRs. We used an "analyze and transform" approach derived from Morrissey (2016) to account for the folded-normal distribution of absolute value effect sizes (Hereford et al., 2004; Kingsolver et al., 2012). First, we implemented Bayesian meta-analytic random-effects models in R (R Core Team, 2016) with the MCMCglmm package (Hadfield, 2010). We ran a separate model for each trait with study type as a fixed effect to obtain separate effect sizes for distinct population and genotype studies. Each trait may have a different level of sample error due to measurement techniques and technology. We account for measurement error in our analyses in two ways. First, we weighted the effect sizes by the inverse of the sample variance. Second, we have modeled each trait separately, allowing each trait to have a different error structure. The model structure was:

$$\hat{\wp}\_{ij} \sim \beta \,\,\,\ast \,\, X + \text{stud}\,\varphi\_i + \text{cg}\_{\vec{\jmath}|i} + m\_{\vec{\imath}\vec{\jmath}} + \text{e}\_{\vec{\imath}} \tag{3}$$

$$m\_{\vec{\imath}\vec{\jmath}} \sim N(0, M\_{\vec{\imath}\vec{\jmath}})$$

where yij is the untransformed LRR for the trait from the jth common garden (cg) which is nested in the i-th study; β is the estimate of the fixed effect; X is a dummy-coded variable representing whether the LRRs were comparing distinct populations or genotypes; study<sup>i</sup> represents the random intercept in the i-th study, cgj|<sup>i</sup> represents the random intercept for the j-th common garden in the i-th study, and eij represents residual error. Effect size estimates were weighted using mij, which accounts for the known sample variance associated with trait yij and the covariance structure of effect sizes through the use of matrices **Mij**. For each common garden j, the **Mij** matrix was constructed following Noble et al. (2017) to account for within-common-garden covariance of effect sizes due to pairwise comparisons (i.e., covariation because each population was represented in multiple contrasts). The matrix was created using the function make\_VCV\_matrix() (Noble, 2019). Every effect is represented by both a row and column in the matrix which has sample error variance on the diagonal and estimates of covariance between effect sizes on the off-diagonals. Covariance was estimated as cov k, l = r ∗ √ V<sup>k</sup> ∗ √ Vl , where k and l are the row and column ID of the within-common-garden matrix, V is the sampling variance calculated for each effect ("variance LRR" in Equation 2) and r is the correlation between effect size k and l which was assumed to be 0.5 for all offdiagonal elements.

For each model, we used non-informative priors for fixed effects and inverse-Wishart priors for random effects (Hadfield, 2010). We ran the model for 150,000 iterations, discarded the first 30,000 as burn-in and used a thinning interval of 15. We examined plots of Markov chain Monte Carlo to ensure good mixing. We generated posterior predictive distributions for each type of contrast (distinct populations or genotype comparisons) for each trait. These distributions were then transformed using a folded-normal distribution. This transformation allowed us to present the mode and 95% highest posterior density intervals (95% credible intervals) that relate to absolute effect sizes (Morrissey, 2016). When using absolute effect sizes, confidence intervals that do not cross zero are not necessarily indicative of effects that are significantly different than zero (Morrissey, 2016). Thus, we use the absolute effect size magnitude to compare each contrast of our study on the same scale.

#### Effect of Phosphorus Environment on Population Traits

We performed an analysis on a subset of data to investigate if populations from stoichiometrically distinct environments displayed a greater amount of trait variation than populations that were not identified as stoichiometrically distinct. For our stoichiometrically distinct populations we specifically looked at studies in which it was stated that resource P availability differed between the environments from which the populations originated. We focused on environmental P availability in this subanalysis because P is an essential, non-substitutable elemental resource for all life forms, and because it was the only element for which there were sufficient studies and observations to enable the analysis. Populations were considered to be from environments with similar P availability if there was no mention of nutrient resource stoichiometry in the methods or results. After filtering studies, we had limited data for this analysis, therefore we focused on the three traits which had at least 10 pairwise comparisons: P content, AAA and somatic growth rate. We used the same statistical procedure as in the global meta-analysis to determine mean effect size, now using differences in P availability (binary Yes/No) as a modifier.

#### RESULTS

From the 30 studies we found, 17 compared populations, 12 compared genotypes and one had both population and genotype comparisons (**Table 1**). The 18 studies comparing distinct populations had a total of 44 common gardens, with a range of 1–6 common gardens per study. Most of these studies were on invertebrates (n = 11) while comparatively fewer focused on vertebrates (n = 6) and plants (n = 1). The studies were overwhelmingly conducted on freshwater organisms (n = 15), with fewer terrestrial examples (n = 3). Eleven different genera were represented, with Daphnia being the most commonly studied (n = 5). Elemental composition was the ES trait measured most frequently (study n = 15) while excretion (n = 5) and AAA traits (n = 8) were less commonly measured.

The 13 studies that assessed variation among genotypes within a population had a total of 23 common gardens, with a range of 1– 4 common gardens per study. Most of the genotype studies were focused on plants (terrestrial n = 8, marine n = 1), and a few focused on freshwater invertebrates (n = 3) and vertebrates (n = 1). There were 8 genera represented, with the most common being Populus (n = 4) and Daphnia (n = 2). Similar to the population studies, the most frequently measured ES traits were elemental composition (n = 12) while excretion and AAA traits were each only measured in three studies.

Among the 30 studies included in our meta-analysis, 21 were based on populations of organisms that were obtained from ecosystems in North America. Of the nine remaining studies, eight were based on populations from north-western continental Europe, and one was based on a single species of freshwater snail from New Zealand. We found no suitable studies of tropical organisms to include in our meta-analyses.

#### How Much Intraspecific Genetic Variation Is Present in ES Traits Relative to Life History Traits?

The magnitude of genetic variation was large for both life history traits, with effect sizes of 30.9 and 84.4% for morphology and somatic growth rate across populations, and of 32.8 and 146% for the same among genotypes. Most ES traits showed a moderate or small amount of genetic variation (i.e., absolute effect size; **Figure 1**) relative to life history traits. Differences among genotypes in AAA traits was the only effect size that was

larger than life history traits (with an effect size of 164%) and the magnitude of this effect can be partially attributed to a single study that found substantial genotypic effects, particularly for P retention (Sherman et al., 2017).

# Do Particular Classes of ES Traits Show More Trait Variation Than Others?

Of the stoichiometric traits, excretion and AAA traits showed the most variation across populations and among genotypes, while elemental composition varied less (**Figure 1**). Contents of P, N and C showed decreasing levels of trait variation (**Figure 1**), with particularly low levels of variation across populations and among genotypes in both N and C. Consequently, nutrient ratios that include P show more trait variation than C:N (**Figure 1**).

# Is the Magnitude of Intraspecific Genetic Variation Similar Among Genotypes (i.e., Within a Population) to What Is Observed Across Populations?

We found a similar magnitude of trait variation among genotypes that came from the same population as we did when we compared distinct populations (**Figure 1**). For some traits, namely AAA, P content, and N content we saw more trait variation among genotypes than among populations (**Table 2**). To ensure that

of pairwise comparisons used to calculate each effect size. Somatic growth rate and morphology (measures of mass and size) effect sizes are denoted with a dashed

line to indicate that this meta-analysis is not a comprehensive review of common garden experiments that measured these traits.

TABLE 2 | Deviance information criteria (DIC) of models for each trait with and without study type (global analysis) or stoichiometric environment (sub analysis) as a categorical variable.


Positive 1DIC values indicate that the model with study type or stoichiometric environment as a fixed effect is the model that would best replicate a similar dataset. 1DIC < 2 indicates little support for the inclusion of the categorical predictor variable, 2 < 1DIC < 10 indicate some support and 1DIC > 10 indicates strong support (Burnham et al., 2002). "SGR" stands for Somatic Growth Rate. Bold values indicate at least some support for the inclusion of study type (global analysis) or stoichiometric environment (sub analysis) as a categorical variable in the model.

this result was not an artifact of statistical methodology we conducted a supplemental analysis. We calculated LRRs for genotype studies using a different, more conservative approach and still found that there was similar or greater variation in traits when comparing genotypes than when comparing populations (**Figure S1**).

## Do Populations From Different Biogeochemical Environments Show More Trait Differentiation?

Our subanalysis included a total of 17 studies; two studies that differed in resource availability, but not specifically P, were removed from this analysis. Seven studies (16 common gardens) used populations originating from environments that differed in P availability. All seven studies focused on freshwater invertebrates, and five used Daphnia spp. Populations that were from environments with different P availability did not show any more trait differentiation than those from populations without any known or measured biochemical differences (**Figure 2**). Somatic growth rate and P content showed nearly identical differentiation across populations that did and did not come from different P environments (somatic growth rate: 85 and 100% difference, P content: 9.4 and 4.8% difference), while AAA traits varied more in populations that were from similar environments (**Table 2**; **Figure 2**).

#### DISCUSSION

Our meta-analysis of common garden studies uncovered substantial intraspecific genetic variation in most ES traits. However, the magnitude of genetic variation that we observed for ES traits was generally less than what was observed for key life history traits measured in these same studies. Among ES traits, acquisition, assimilation, allocation, and excretion traits tended to vary more than elemental composition traits (i.e., the elemental content and stoichiometric ratios of organism biomass), and there tended to be markedly more variation among genotypes than across populations. Taken together, these results suggest that genetic variation in ES traits, particularly AAAE traits, within and across populations may be important for ecological and evolutionary outcomes.

## Evolution and the Magnitude of Genetic Variation in ES and Life History Traits

Our meta-analysis revealed substantial genetic intraspecific variation in ES traits across a range of taxa. These findings join a growing body of literature that has demonstrated that variation below the species level can be both substantial and biologically important (Schweitzer et al., 2004; Crutsinger et al., 2006; Whitham et al., 2006; Post et al., 2008; Violle et al., 2012; Des Roches et al., 2018). In the light of these findings it is increasingly clear that incorporating data on intraspecific variation is likely to increase understanding of ecological patterns in nature (Bolnick et al., 2011). Our comparison of the amount of genetic variation in ES traits and the genetic variation in life history traits demonstrated that nearly all types of ES traits showed lower genetic intraspecific variation than life history traits. One clear exception was genetic variation among genotypes in AAA traits, which was influenced by a large amount of variation in one study (Sherman et al., 2017). If this difference in variation between ES and life history traits is genuine, it could be inferred that less genetic variation in ES traits is maintained within populations and that ES trait means evolve slowly across populations, in comparison to life history traits. It is certainly plausible that life history traits should vary more widely than ES traits, because constraints on organism elemental composition may be strong relative to those on some life history traits, and because tradeoffs in element allocation within an organism could decouple variation in the total elemental content of an organism from variation in organism morphology and/or growth rate (Meunier et al., 2017).

The evolution of ES traits has been explicitly assessed in various ways in a small number of prior studies. Two studies have calculated the heritability of ES traits using genotypes from wild populations (H<sup>2</sup> = 0.61 [Barbour et al., 2015]; H<sup>2</sup> = 0.18–0.21 [Crutsinger et al., 2014a]) suggesting that ES traits can respond to selection when genetic variation is present (but see Hansen et al., 2011). Indeed, artificial selection experiments have shown that all livestock species appear to have genetic variation for feed efficiency and this variation often has moderate heritability (Arthur and Herd, 2005). Moreover, artificial selection experiments on ES traits tend to achieve rapid ES trait evolution in agricultural species (Neely et al., 2008; de Verdal et al., 2013; Mignon-Grasteau et al., 2017). Selection on life history traits can also lead to evolution in ES traits, for example, artificial selection on body size can produce

rapid evolution in P content (Gorokhova et al., 2002). We were able to locate two selection experiments where environmental conditions were manipulated and evolution of ES traits were tracked through time in laboratory conditions. In one of these experiments ES traits did not evolve over ∼130 generations (Declerck et al., 2015) while the second experiment, which found that N and P content evolved in carbon-limited conditions, assessed trait evolution over 50,000 generations (Turner et al., 2017). We were only able to find a single instance of tracking the evolution of ES traits of natural or naturalistic populations (Frisch et al., 2014). In this study, the comparison of genotypes hatched from Daphnia resting eggs originating from a 700-year time interval in a single lake undergoing cultural eutrophication revealed an effect of P availability on AAA traits but not elemental composition (Frisch et al., 2014). Clearly, much more work is needed to determine whether and when ES traits evolve rapidly in natural contexts. Selection experiments conducted in nature or observational studies that track the evolution of populations through time in both ES and life history traits would provide valuable empirical data on the rate of evolution of ES phenotypes with realistic population dynamics and selective landscapes. These types of approaches could improve the understanding of whether and why ES traits tend to show less genetic variation than life history traits.

# Differences in the Magnitude of Genetic Variation Among Stoichiometric Traits

One of the fundamental principles of ES theory is that the elemental composition of organisms is constrained by their basic requirements for the non-substitutable elemental resources that are allocated to fundamental biological structures and processes (Reiners, 1986; Sterner and Elser, 2002). Many of these lifeenabling structures and processes are ancient (e.g., bone and wood; Wagner and Aspenberg, 2011; Morris et al., 2018), and some are as old as life itself (e.g., RNA synthesis; Joyce, 1989). Thus, it is perhaps not surprising that we found little genetic variation in bulk elemental composition within species, because any contemporary variation is miniscule relative to that found over the entirety of evolutionary history. In other words, the magnitude of rapid evolution in elemental composition is limited relative to what is seen across species with deep evolutionary divergence.

However, there was substantially more intraspecific genetic variation present in the AAAE traits than in the elemental composition traits, and this was the case for both genotype and population comparisons (**Figure 1**). This leads to questions about constraint of genetic diversity for elemental composition relative to AAAE traits. Acquisition, assimilation, allocation, and excretion are traits which actually contribute to and maintain elemental composition within the necessary range for organismal growth, maintenance, and survival (Sterner and Elser, 2002). Thus, our results suggest that there is more genetic variability in the strength of stoichiometric homeostasis than there is in the actual bulk elemental composition of organisms. This makes sense, given that genetic variation in AAAE traits is a pre-requisite for genetic variation in elemental composition (Jeyasingh et al., 2014). Indeed, elemental composition reflects the "balance" between the acquisition, assimilation, allocation, and excretion of elements at a given point during an organism's life. An evolved change in elemental composition would, therefore, presumably entail a corresponding heritable shift in the balance among AAAE traits.

Importantly, this does not mean to say that elemental composition does not evolve. Over macroevolutionary timescales evolution has generated an enormous diversity of stoichiometric phenotypes across taxa (Sterner and Elser, 2002). Yet much of this stoichiometric diversity has emerged not necessarily because ES traits are (or were) under direct selection, but because they are mechanistically correlated with functional traits that have had direct fitness implications. For example, the high C:nutrient ratio of a woody tree is a consequence of selection for height, which imparts an obvious competitive advantage to terrestrial autotrophs. Wood provides the rigid structural material necessary for this height and has a C:N ratio 200–1000 (Levi and Cowling, 1969). Thus, as entire organisms, trees have high C:N ratios relative to other organisms, but not because C:N ratio was under direct selection per se. Another example of indirect selection on ES traits with clearer implications for rapid evolution, comes from the stoichiometric growth rate hypothesis (GRH), which asserts that the functional trait of rapid growth has a P rich signature due to the high P content of ribosomal RNA (Elser et al., 2000b; Sterner and Elser, 2002). Identifying the stoichiometric patterns associated with variation in a larger suite of functional traits will be tremendously valuable to the study of evolutionary and ecological relationships (Meunier et al., 2017).

Within the elemental composition traits, we found more variation in P content than in C and N content. This is consistent with the argument that variation in organism P content will generally be greater than variation in organism C and N content, because levels of C and N tend to be similar among major biomolecules (i.e., protein and nucleic acids) while P content varies widely among biomolecules (e.g., protein contains no P, while RNA contains around 9.2% P; Sterner and Elser, 2002). Thus, changes in proportional concentrations of major biomolecules on a cellular level can scale up to drive variation in overall organism P content. Our results provide a novel genetic perspective that supports the view that P content is more variable than C and N, at least in the case of genotype comparisons.

#### Genetic Variation in ES Traits Within Populations and Across Population Pairs

We collected data on the amount of genetic variation by comparing across population pairs and among genotypes from the same population. The extent of variation among genotypes was generally larger than variation in population trait means, even when using a conservative methodology to calculate genotypic variation (**Figure S1**). When comparing across populations, we assessed shifts in the means between pairs of populations, which does not explicitly account for variation within populations. As such, our finding that genetic variation is generally greater amongst genotypes within a population than across populations reflects that the average genetic trait difference between genotypes is larger than the average trait shift between populations. Finding lower levels of variation across populations than among genotypes indicates that there is either little differentiation across populations, or large amounts of variation among genotypes. Whether this is to be expected or not depends on where genetic variation within a species is expected to be found: within or across populations. Genomic sequencing both within and across populations has uncovered support for each, including cases where the majority of genetic variation is found within populations (Pometti et al., 2015) and cases where the majority is found among them (Bakker et al., 2006). From the latter perspective this result is surprising, given that genetic variation within a population is limited by recombination amongst genotypes, which would serve to limit divergence within a population (Jain and Bradshaw, 1966; Ehrlich and Raven, 1969). There is potential for increased variation across populations because they originated from different environments, which have the potential to act as an agent of divergent selection driving local adaptation, whereas the genotypes compared mostly originated from the same environment. Genotypes may show a high amount of variation due to fluctuations in the environment, gene flow, or any other mechanism that can enhance the maintenance of trait diversity. On the other hand, trait divergence between populations may be limited if they do not come from stoichiometrically-disparate environments (e.g., little divergent selection on ES traits) or do not show substantial local adaptation due to stoichiometric constraints.

Our data provide some ability to assess the role of environmental differences, specifically P availability, in driving divergence in ES traits across populations. Populations that come from divergent P environments and those that come from similar P environments showed no differences in the amount of genetic variation in ES traits. This suggests that natural environmental differences in elemental availability may not be a sufficiently strong agent of selection to drive rapid evolution in ES traits, which could explain the largely similar levels of variation observed across genotypes and pairs of populations. Determining whether the limited variation we observed across natural populations is due to comparatively modest environmental stoichiometric differences or other factors that could limit divergence of populations (e.g., low divergence times, connectance between populations, reduced efficacy of selection) would be valuable in addressing this disparity.

# Caveats and Limitations of the Meta-Analysis

There are several limitations in our analysis stemming from the meta-analysis framework, the nature of common gardens, and the specific designs of the studies available. First, studies that met the criteria for our meta-analysis included both plants and animals, but those that measured variation across genotypes were dominated by plants and cladocerans. Thus, inherent differences among taxonomic groups may be confounded with the genotypeand population-level comparisons, though the relative amounts of variation were similar among all taxonomic groups. Second, genotype-level variation may be overestimated relative to natural populations, because some studies included genotypes with low ecological fitness (e.g., Sherman et al., 2017) and some collected genotypes from wide geographical areas that covered a range of environments similar to the magnitude of variation we would expect across populations. Third, many of the common gardens were only conducted for a single generation, meaning that maternal and epigenetic effects could represent latent sources of trait variation (Weaver et al., 2004). Finally, a "common garden" allows for measurement of the amount of genetic variation in traits in a single environment and the chosen environmental conditions may influence the amount of trait variation displayed. For example, there was more variation in P-acquisition and assimilation under P-limited conditions relative to P-sufficient conditions (Sherman et al., 2017). However, we assumed that the common garden conditions were chosen by the original study authors to reflect ecologically-relevant environments, and that the trait variation we observed is, therefore, a reasonable representation of that found in nature.

#### Ecological Consequences of Intraspecific Variation in ES Traits and the Potential for Eco-Evolutionary Feedbacks

Over the past two decades there has been a growing recognition that intraspecific trait variation can have large effects on ecological patterns and processes (Schweitzer et al., 2004; Crutsinger et al., 2006; Whitham et al., 2006; Post et al., 2008; Harmon et al., 2009; Bassar et al., 2010; Violle et al., 2012; Des Roches et al., 2018). Intraspecific variation in ES traits has been identified as particularly likely to have strong ecological consequences because ES traits connect the organism with its environment. Our meta-analysis uncovered substantial differences in the amount of intraspecific variation across types of ES traits, with excretion and AAA traits showing considerably more intraspecific variation than most elemental ratios and measures of elemental composition. Although there have been few empirical tests of the ecological consequences of intraspecific variation in ES traits (but see Chowdhury and Jeyasingh, 2016), any evolution that shapes these traits is likely to have ecological consequences (Jeyasingh et al., 2014; El-Sabaawi et al., 2016). Future experimental work that measures the ecological consequences of genetic variation in ES traits in ecologically realistic contexts will be crucial to determining the importance of intraspecific variation in ES traits for community structure and ecosystem function.

Ideally, this future work will indicate whether rapid evolution of ES traits is likely to have ecological consequences, and whether these ecological consequences could lead to rapid evolution (i.e., an eco-evolutionary feedback). Using a strict definition of eco-evolutionary feedback (Schoener, 2011), two interacting processes would be required: ES traits would have to evolve rapidly in a population and the evolution of these ES traits would need to modify the environment in a way that significantly affects selection on that population. The potential for eco-evolutionary feedbacks mediated by the evolution of ES traits has been much discussed (Jeyasingh et al., 2014; Leal et al., 2017b; Rudman et al., 2019), and some studies have demonstrated that rapid evolution can shape elemental availability (Bassar et al., 2010; Rudman and Schluter, 2016). Yet,

#### REFERENCES


whether these differences in elemental availability feed back onto rapid evolution is wholly unknown. As such, determining the existence, magnitude, and prevalence of these eco-evolutionary feedbacks is well beyond the scope of our meta-analysis, especially as there is very little empirical evidence of ecoevolutionary feedbacks overall (Schoener, 2011; Turcotte et al., 2011; Matthews et al., 2016; Rudman et al., 2018). What little data our meta-analysis does provide to examine eco-evolutionary feedbacks based on ecological stoichiometry does not support widespread rapid evolution of ES traits, as variation across populations was similar to what was observed amongst genotypes found within populations. Moreover, we did not find evidence that divergence in P environment led to enhanced divergence in ES traits across populations, suggesting that environmental variation might need to be substantial to drive ES trait evolution. Nevertheless, our findings are mere hints at the prevalence of ecoevolutionary feedbacks mediated by ES traits and future work is certainly warranted.

### AUTHOR CONTRIBUTIONS

KL conceived of the study. KL, OB, TK, SR, and CS reviewed literature, aggregated data, and wrote and edited the manuscript. The authors were listed alphabetically after the first author.

### FUNDING

Funding was provided by the National Science Foundation (NSF) through NSF award DEB-1840408.

### ACKNOWLEDGMENTS

We would like to thank Jim Elser, Michelle Evans-White, and Joe Vanderwall for their efforts that made the Woodstoich IV workshop possible. We thank our Woodstoich IV mentors, Jotaro Urabe, and Jennifer Schweitzer for feedback and helpful discussions. We thank Daniel Noble and Vincent Fugère for valuable guidance on the statistical analysis.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00339/full#supplementary-material

genetic basis of herbivore community assembly. Funct. Ecol. 29, 995–1006. doi: 10.1111/1365-2435.12409


primrose (Oenothera biennis) from a field experiment. J. Evol. Biol. 22, 1295–1307. doi: 10.1111/j.1420-9101.2009.01747.x


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer, PJ, declared a past co-authorship with one of the authors, SR, to the handling editor.

Copyright © 2019 Lemmen, Butler, Koffel, Rudman and Symons. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Exploring Silica Stoichiometry on a Large Floodplain Riverscape

Joanna C. Carey 1,2 \* † , KathiJo Jankowski 3†, Paul Julian II 4†, Lienne R. Sethna5† , Patrick K. Thomas 6† and Jason Rohweder <sup>3</sup>

*<sup>1</sup> Division of Math and Science, Babson College, Wellesley, MA, United States, <sup>2</sup> Department of Earth and Environment, Boston University, Boston, MA, United States, <sup>3</sup> U.S. Geological Survey Upper Midwest Environmental Sciences Center, La Crosse, WI, United States, <sup>4</sup> Whitney Marine Laboratory for Marine Bioscience, University of Florida, St. Augustine, FL, United States, <sup>5</sup> O'Neill School of Public and Environmental Affairs, Indiana University, Bloomington, IN, United States, <sup>6</sup> Plankton Ecology Lab, Institute for Chemistry and Biology of the Marine Environment, University of Oldenburg, Wilhelmshaven, Germany*

#### Edited by:

*Michelle Evans-White, University of Arkansas, United States*

#### Reviewed by:

*Walter Dodds, Kansas State University, United States Michael J. Vanni, Miami University, United States Dag O. Hessen, University of Oslo, Norway*

> \*Correspondence: *Joanna C. Carey jcarey@babson.edu*

*†The preparation and creation of this manuscript was a strong collaborative effort and these authors have contributed equally to this work, with authors ordered alphabetically*

#### Specialty section:

*This article was submitted to Behavioral and Evolutionary Ecology, a section of the journal Frontiers in Ecology and Evolution*

> Received: *22 July 2019* Accepted: *29 August 2019* Published: *25 September 2019*

#### Citation:

*Carey JC, Jankowski K, Julian P II, Sethna LR, Thomas PK and Rohweder J (2019) Exploring Silica Stoichiometry on a Large Floodplain Riverscape. Front. Ecol. Evol. 7:346. doi: 10.3389/fevo.2019.00346* Freshwater ecosystems are critical zones of nutrient and carbon (C) processing along the land-sea continuum. Relative to our understanding of C, nitrogen (N), and phosphorus (P) cycling within the freshwater systems, the controls on silicon (Si) cycling and export are less understood. Understanding Si biogeochemistry and its coupled biogeochemical processing with N and P has direct implications for both freshwater and coastal ecosystems, as the amount of Si in relation to N and P exported by rivers to coastal receiving waters can determine phytoplankton species assemblages, which in turn affects C cycling and food web structure. Here we examine the relationships between dissolved Si (DSi), total nitrogen (TN), and total phosphorus (TP) concentrations, and how these relationships relate to basin land cover, lithology, and river hydrogeomorphology (i.e., aquatic areas) in the Upper Mississippi River System (UMRS) using two datasets (one from the tributaries and one from the mainstem) that span a nine-year period (2010–2018) representing >10,000 unique samples. We found significant declines in DSi concentrations, as well as Si:TP and Si:TN ratios along the north-south gradient of the mainstem UMRS across the six aquatic area types we study here. This signal was driven partially by a corresponding decline in tributary DSi inputs along this latitudinal gradient. Contrary to findings from other regions of North America, basin land cover was not an important predictor of tributary DSi concentrations, especially compared to lithology. However, Si:TN and Si:TP ratios appear to be strongly controlled by basin land cover, likely due to excess N and P loading from row-crop agriculture. DSi, and its ratio with TN and TP (i.e., Si stoichiometry), was similar across most aquatic area types, including run-of-river impoundments and the main channel, suggesting similar processes affecting DSi, TN, and TP concentrations in these reaches. However, backwater lakes had lower DSi and TN concentrations compared to the other aquatic area types, highlighting the importance of water residence time and nutrient uptake in controlling Si stoichiometry in inland waters. Together, our results show that rivers are not simple pipes for Si, but rather the complexity in watershed characteristics, hydrology, and biological uptake results in dynamic Si stoichiometry along the river continuum.

Keywords: silica, river, residence time, nutrients, land cover, Mississippi river, ecological stoichiometry, biogeochemistry

"Letting the days go by, let the water hold me. . . Letting the days go by, water flowing. . . Water dissolving and water removing Into the blue again into silent water Once in a lifetime. . . " –The Talking Heads

#### INTRODUCTION

The concentrations, ratios, and fluxes of nutrients in freshwater systems fundamentally connect terrestrial and marine systems and are tightly linked with biodiversity and ecosystem functioning (Smith et al., 1999; Rabalais et al., 2002; Elser et al., 2007; Maranger et al., 2018). Increasing silicon (Si) supply from continental mineral weathering gave rise to the diatoms in the Cenozoic era (Cermeño et al., 2015), which now produce nearly half the oxygen we inhale, substantially altering the global carbon (C) cycle (Berner et al., 1983; Brady and Carroll, 1994). Diatoms are highly productive autotrophs in the oceans, accounting for ∼50% of oceanic net primary productivity (Rousseaux and Gregg, 2014). Unlike other phytoplankton, marine diatoms generally require as much Si as nitrogen (N) on a molar basis, meaning the ratio of Si:N in receiving waters has direct implications for phytoplankton species assemblages and C sequestration rates (Officer and Ryther, 1980; Beusen et al., 2009; Frings et al., 2016). Compared to marine species, freshwater diatoms have more variable stoichiometric nutrient demands, often requiring greater Si relative to N and P (Conley et al., 1988). Thus, the functioning of aquatic ecosystems and the stability of our climate are linked to the processing and export of Si, and its ratio with critical nutrients.

A central tenet of stoichiometric theory (Redfield, 1934, 1958) is based on the consistency of organismal elemental composition, which led to the understanding of the conditions limiting phytoplankton growth relative to C, N, and phosphorus (P). While Redfield initially focused on C:N:P, a "Redfield-like" ratio was developed to include Si in order to better understand the conditions limiting diatom growth (Brzezinski, 1985). Despite the importance of Si and its ratio with other nutrients (i.e., Si stoichiometry) in controlling aquatic productivity, we know less about the complex processes controlling Si exports along the land-sea continuum relative to C, N, and P. For example, Maranger et al. (2018) recently examined the global controls on riverine C:N:P processes along the river continuum, expanding upon our appreciation for rivers as "processors" rather than "pipes" (Cole et al., 2007), as C is processed during its transit to the oceans (Tranvik et al., 2009; Hotchkiss et al., 2015). However, we still lack a holistic understanding regarding the physical and biological factors that shape Si:N:P stoichiometry along the river continuum.

Chemical weathering is the primary source of Si to the biosphere, as Si is the second most abundant element in the lithosphere (Exley, 1998). In turn, variation in geochemical factors, such as lithology and climate, largely control baseline availability of Si in a given system (Bluth and Kump, 1994). Rivers transport DSi, the bioavailable form of Si, from terrestrial to aquatic systems in large quantities, exporting >80% of annual DSi inputs to the global oceans (Tréguer and De La Rocha, 2013). Relative to Si, human activities (e.g., fertilizer runoff, urban wastewater treatment, fossil fuel combustion) have increased the export of N and P to downstream receiving waters. In the presence of excess N and P, diatoms can bloom until available Si is depleted. Such Si-limited waters can then favor the growth of non-siliceous, potentially toxic algae (Officer and Ryther, 1980; Smayda, 1990; Conley et al., 1993; Justic et al., 1995; Anderson ´ et al., 2002; Parsons and Dortch, 2002; Garnier et al., 2010) and can also affect the toxicity of diatom communities (Pan et al., 1996; Bargu et al., 2016). In turn, altered ratios of nutrients (N, P, Si) relative to the Redfield ratio (16N:1P:16Si) have raised serious concerns about anthropogenic effects on diatom-based food chains in freshwater and marine ecosystems (Officer and Ryther, 1980; Conley et al., 1993; Turner et al., 1998; Anderson et al., 2002). Globally, this type of stoichiometric imbalance in N, P, and Si has been implicated in influencing toxic and harmful algal blooms in both fresh and marine systems, including the Baltic Sea, North Sea, Great Lakes, and the California Current Ecosystem (Schelske and Stoermer, 1971; Anderson et al., 2009; Garnier et al., 2010; Davidson et al., 2012; Bargu et al., 2016). For example, agricultural activities in the Mississippi River basin, which accounts for ∼40% of the land area in the continental United States, have altered nutrient ratios of the Mississippi River and its exports to the Gulf of Mexico over the last several decades (Turner et al., 1998; Rabalais et al., 2002; Royer, 2019). These changes in stoichiometry are likely to compound already severe environmental stressors, such as massive hypoxic dead zones in the Gulf of Mexico (NOAA; Rabalais et al., 2002) caused by increases in N and P loads.

In addition to the impact of excess N and P on the relative availability of Si in aquatic systems, human activities can directly alter global Si cycling through several means. First, land use change can directly alter DSi export rates from terrestrial to aquatic systems (Conley et al., 2008; Struyf and Conley, 2012; Carey and Fulweiler, 2013). Terrestrial plants have the capacity to release and retain Si on terrestrial landscapes. Plant growth can stimulate Si mineral dissolution, releasing DSi into soil solution (Drever, 1994). Terrestrial plants also sequester large quantities of Si within their tissue (Conley, 2002; Hodson et al., 2005) (84 ± 29 Tmol yr−<sup>1</sup> ; Carey and Fulweiler, 2012a). Therefore, shifts in watershed land use can alter Si retention and export from terrestrial landscapes (Cooke and Leishman, 2011; Cornelis et al., 2011; Carey and Fulweiler, 2016). In temperate and tropical landscapes, urbanization, deforestation, and agricultural land uses have been shown to alter fluvial Si exports (Conley et al., 2008; Struyf et al., 2010; Carey and Fulweiler, 2012b, 2013; Vandevenne et al., 2012; Chen et al., 2014), groundwater DSi concentrations (Maguire and Fulweiler, 2016, 2019) and soil Si cycling (Clymans et al., 2011; Vander Linden and Delvaux, 2019). In addition, river damming alters Si cycling in aquatic systems by increasing water residence time and stimulating diatom blooms and Si burial behind impoundments (Humborg et al., 1997, 2000, 2006). Because Si is recycled slower than N and P, waters directly downstream of dams are often Si depleted (van Bennekom and Salomons, 1980). More subtle hydrogeomorphic variation in dynamic river-floodplain systems also has the capacity to alter Si processing by creating heterogeneous patterns in sedimentation and uptake (Junk et al., 1989; Houser and Richardson, 2010).

The Upper Mississippi River System (UMRS) provides a useful template to evaluate the effects of land use and river landscape complexity on large river Si processing and stoichiometry. The northern reaches of the UMRS contain a diversity of side channels and backwater lakes that vary in their degree of connectivity to main channel flow, whereas the river becomes increasingly disconnected from its floodplain and faster-flowing in its lower reaches (Wilcox, 1993; De Jager et al., 2018). These shifts along the longitudinal continuum result in shifts in processing rates (Richardson et al., 2004) and phytoplankton communities (Manier, 2014) that likely shape the way Si moves through the system.

Here, we evaluated the extent and the origins of variation and coupling in large river DSi, total N (TN) and total P (TP) stoichiometry using nearly a decade (2010–2018) of empirical data from a spatially extensive dataset that spans six reaches of the mainstem UMRS and the mouths of 23 tributaries. We investigated the role of land use, lithology, and in-situ ecological controls in altering external inputs of Si:TN:TP from the 23 tributaries. We then examined how mainstem river hydrogeomorphology altered internal patterns and coupling of Si:TN:TP along the ∼1,900 km of the mainstem UMRS. We hypothesized that:

H1: Basin land use and lithology are important controls on stoichiometric inputs of Si:TN:TP from the tributaries to the mainstem;

H2: Si:TN:TP varies longitudinally along the mainstem UMRS due to variable nutrient inputs and internal processing, and;

H3: Si:TN:TP varies across aquatic area types due to differences in water residence time and other physical characteristics.

#### MATERIALS AND METHODS

#### Study System

The UMRS is defined as the commercially navigable reaches of the Upper Mississippi and Illinois Rivers (**Figure 1**). Flowing ∼1,900 km, the river and its tributaries encompass a substantial gradient in land use, from forested headwaters in the north to increasingly intensive agriculture as the river winds its way south. We used limnology data collected by personnel of the Long-term Resource Monitoring element (LTRM) of the US Army Corps of Engineers' Upper Mississippi River Restoration Program (UMRR). These data were collected from 2010 to 2018 in six reaches of the Upper Mississippi River System (UMRS), representing 490,000 km<sup>2</sup> (**Table 1** and **Figure 1**). Five of the six reaches are "navigation pools," defined as the reach of river between two consecutive locks and dams (**Supplemental Information**).

The UMRS and its floodplain contain a diversity of aquatic areas that vary in depth, aquatic vegetation, water residence time, and hydraulic connectivity to the main river flow (De Jager et al., 2018). We examine six aquatic areas types within the mainstem of the UMRS, including the main channel, side channels, impounded areas, contiguous backwater lakes, riverine lakes (e.g., Lake Pepin), and isolated backwater lakes. Isolated backwater lakes are connected to the main channel only at extremely high flows. These aquatic areas are primarily differentiated based on their depth and velocity characteristics (**Table 2**), which we used as an indicator of relative water residence time. The navigation impoundments in the system are created by low-head, run-of river dams built to maintain sufficient depth in the river for navigation during low flow and were designed to have little impact on discharge or water level during high flow conditions (Sparks et al., 1998). Despite the large-scale shifts in land use and hydrologic character of the system, parts of the UMRS remain a complex floodplain riverscape that contains a diversity of aquatic areas that vary widely in water residence time, biogeochemical processing rates, and freshwater communities (Wilcox, 1993; Richardson et al., 2004; Manier, 2014; De Jager et al., 2018).

#### Data Collection

We used two datasets to conduct our analysis (retrieved from https://umesc.usgs.gov/ltrm-home.html). The first dataset (hereafter referred to as "mainstem" data) consists of 216 unique sampling events representing >9,800 individual samples from seasonal stratified random sampling of six aquatic area types along the mainstem of the UMRS (Soballe and Fischer, 2004) (see **Supplemental Information**). The second dataset consists of >2,500 distinct sampling events collected ∼monthly (14 times per year) from fixed sampling locations at the mouth of 23 tributaries, just upstream of their confluence with the mainstem UMRS.

All measurements and water samples were taken 0.2 m below the surface of the water, or 0.2 m below the submerged base of the ice when present (i.e., during winter sampling). Both water speed and direction were measured at the time of sampling using a Hach model FH950.0 Portable Velocity Meter (Hach, Inc. Loveland, CO) and a standard magnetic compass. At sites where flow velocity was below detection of the meter (<2 cm s−<sup>1</sup> ), velocity was reported as half the detection limit of the instrument, or 1 cm s−<sup>1</sup> . Water samples were returned to the lab immediately after sampling for the determination of turbidity, total nutrients (e.g., TN, TP), dissolved silica (DSi), and chlorophyll a (CHL) concentrations. DSi samples were filtered in the field (0.45µm membrane) and frozen until analysis. TN and TP were preserved in the field with concentrated H2SO4, transported on ice, and refrigerated until analysis. DSi, TN, and TP were determined colorimetrically following standard methods (American Public Health Association, 2005). CHL was determined fluorometrically (CHLF) for all sites, and spectrophotometrically (CHLS) for a randomly chosen subset of sites (random selection of 10% of the sites sampled) (see **Supplemental Information**). Water column CHL values represent the majority of algal biomass in the system; phytoplankton are abundant in many areas of the river due to the open nature of the channels that provide ample light availability in the upper water column (Gardner et al., 2019). Light typically does not reach the benthos due to

high turbidity and depth, thereby limiting benthic production in channels. Thus, while we recognize benthic algal production could be important in off channel areas with higher water clarity (i.e., backwaters), this dataset did not include benthic algal biomass data. We used Streamstats (https://water.usgs. gov/osw/streamstats/ss\_documentation.html) and Geographic Information Systems (GIS) to delineate watersheds and calculate basin characteristics for the tributaries, respectively.

Land use data for each tributary basin were generated from the 2011 National Land Cover Dataset (Homer et al., 2015). Primary bedrock lithology coverage data for each tributary basin were generated from the combined geologic map data for the conterminous US (Schweitzer, 2011).

#### Data Analysis

We used TN and TP as a proxy for all available N and P in the system. TN and TP include the dissolved inorganic portion and the organic portion of the elements. Because the organic portion of N and P can rapidly remineralize, the use of TN and TP captures the dynamic N and P pool available in the system (Downing, 1997; Guildford and Hecky, 2000; Wetzel, 2001; Dodds, 2003; Bergström, 2010). The stoichiometric



*Standard error of mean given in parentheses for nutrient values and ratios. Catchment area, study reach area, mean depth, median discharge, and median retention time from Houser et al. (2010). NA* = *data not available.*

TABLE 2 | Average characteristics of aquatic area types across all reaches.


*TURB, turbidity; VSS, volatile suspended solids; CHL, chlorophyll a concentration; Temp, temperature. Standard error given in parentheses.*

ratios we present here reflect the particulate and dissolved fractions of N and P, but only the dissolved fraction of Si. Because we did not measure amorphous Si in our samples, our ratios are likely underestimates of Si availability relative to N and P in the UMRS, as amorphous Si, which represents 16% of total river non-mineral Si fluxes globally (Conley, 1997), is potentially a biologically important component of river Si availability (Cornelis et al., 2011).

We analyzed TP, TN, and DSi concentrations and molar ratios of Si:TP, Si:TN, and TN:TP for each site within tributaries and mainstem sites from 2010 to 2018. We used multiple linear regression to determine the role of land use and primary lithology in explaining DSi concentrations and stoichiometric ratios in tributaries entering the mainstem UMRS. To create appropriate regression models and meet the assumptions of normality, we first log transformed the response variable (DSi, Si:TN, Si:TP) and used stepwise variance inflation factors (VIFs) to remove co-varying predictors (e.g., to address correlations between land use and lithology). In stepwise VIF ("car" package; Fox et al., 2012), we set a threshold of VIF (in this case 10), and sequentially removed predictors with the highest VIFs over the threshold until all model predictors had VIFs below the threshold. Cultivated crops, wetlands, and dolostone had elevated VIFs (>10) and were removed as predictors from the regression. Upon removing covarying predictors, we then used Stepwise Akaike Information Criteria (AIC; Burnham and Anderson, 2002) ("MASS" package; Ripley et al., 2013) to determine the most appropriate predictors for the regression models of stoichiometric ratios and nutrient concentrations. We developed these models for dependent variables of interest (i.e., TN, TP, DSi, Si:TN, Si:TP) with land use and lithology predictors together, and then again separately in order to determine the relative roles of both land use and lithology in explaining the nutrient stoichiometry in the tributaries. We also used Spearman's rank sum correlation between latitude and tributary basin characteristics (i.e., percent land cover or primary lithology) to determine if trends in basin coverage existed across the study area.

We evaluated differences in DSi concentrations, Si:TP and Si:TN across river reaches and among aquatic areas along the mainstem UMRS in three ways. First, to evaluate changes in stoichiometry along the north-south gradient of the river, we used a linear model that regressed average nutrient concentrations and ratios with river mile in both the mainstem and tributaries (tributary river mile is where the tributaries enter the mainstem). Second, we compared DSi and its ratio with TN and TP across reaches and aquatic areas using mixed effects models with effects of river reach, aquatic area, and their interaction to a model with only an intercept. These models included a random effect of year to account for potential yearspecific variability. We compared these models using AICc and report marginal (R2m, fixed effects) and conditional (R<sup>2</sup> c, full model) R 2 values for each model (Barton and Barton, 2013; Nakagawa and Schielzeth, 2013; Nakagawa et al., 2017; MuMIn package). Subsequent pairwise comparisons among aquatic areas and river reaches were done for the "best" model (i.e., lowest AICc) for each dependent variable (DSi, Si:TP, and Si:TN) using the "emmeans" package (Lenth, 2018). Third, we examined stoichiometric relationships by comparing DSi to TP (TP x Si) and DSi to TN (TN x Si) within each aquatic area by evaluating power law slopes and intercepts using standardized major axis (SMA) regression ("smart" package; Warton et al., 2006). Unlike standard regression techniques, SMA regression allows us to evaluate how DSi relates to TP or TN by assessing best fit line between two variables, rather than using one variable to predict another (Cleveland and Liptzin, 2007; Sterner et al., 2008; Julian et al., 2019). In the SMA framework, a slope different from the absolute value of one indicates the variables are independent and do not change proportionally with each other. If the slope is not statistically different from the absolute value of one, then the variables exhibit proportional changes and have more constrained stoichiometry between nutrients. In addition, the sign of the slope provides information on how nutrients change relative to each other; positive slopes (not significantly different from one) indicate close positive coupling and near constant ratios across samples, whereas negative slopes indicate negative coupling and highly variable ratios among samples as one nutrient increases and the other decreases. In addition to the slope test, we did pairwise comparisons of the slopes within aquatic areas of each reach (Warton and Weber, 2002; Warton et al., 2006). The coefficient of determination (SMA R<sup>2</sup> ) from each model was also used to assess the degree to which the variables were coupled, where high R 2 values indicate higher relative coupling.

We then evaluated how multiple limnological variables [i.e., depth, velocity, temperature, CHL, TN, TP, turbidity, and volatile suspended solids (VSS)] related to DSi concentrations, Si:TP, and Si:TN in two ways. First, we used a principal components analysis (PCA) to look at covariation of these drivers across reaches and aquatic areas using the "vegan" package (Oksanen et al., 2018). We used the mean decadal value for each reachaquatic area combination (n = 23). We considered variables to be significant in the PCA if they had a correlation coefficient >0.7 with either axis. Second, in order to assess how much variation was explained by these variables in different aquatic area types, we used all data in a mixed model framework with turbidity, VSS, temperature and CHL as fixed effects and random effects of "year" and "reach" on the intercept to account for variation resulting from year-to-year changes or longitudinal differences

among reaches not captured by the fixed effects. Sample size therefore varied by aquatic area type (**Table S5**). We did not include TN, TP in these analyses because of their inherent correlation with each ratio. Since we ran this analysis specifically by aquatic area, we did not include velocity or depth as predictors in the model. We used the lme4 package for the mixed model analysis (Bates et al., 2007). We report marginal and conditional R 2 values for these models as above and compare models using AICc (Burnham and Anderson, 2002). All statistical operations were performed using base R unless otherwise noted (Ver 3.5.2, R Foundation for Statistical Computing, Vienna Austria) and the critical level of significance was set at α = 0.05.

#### RESULTS

#### Spatial and Seasonal Variation in Si:TN:TP Stoichiometry

During the 2010–2018 period of record within the mainstem UMRS, DSi concentrations, Si:TN, and Si:TP ranged widely (**Table 1**). The percentage of dissolved inorganic nutrients varied between 32 and 77% of the N pool (dissolved inorganic N to TN), and 16–31% of the P pool (soluble Reactive P to TP) (**Table S1**). In general, we found that variation across space in Si:TN and Si:TP ratios was more closely linked to TN and TP concentrations than to variation in DSi concentrations (**Figure 2**). These patterns are likely due to differences in the magnitude of variation in each nutrient's concentration; we observed higher variation [as the coefficient of variation (CV)] across tributaries in TN (59.3) and TP (52.9) compared to DSi (20.5) (**Table S2**). Moreover, TN and TP were also more variable than DSi in the mainstem, but the CV for TN was dramatically lower in the mainstem (24.2) than the tributaries (59.3) (**Table S2**). Consequently, there was greater variation in the Si:TN ratio across tributaries (CV = 61.9) than across mainstem sites (CV = 32.7), highlighting the relatively high heterogeneity of the UMRS tributary riverscape compared to then main channel, which integrates these diverse sources.

We did not observe any trends over time in DSi, TN, TP concentrations or their ratios across our study period (2010– 2018). However, we observed similar seasonal patterns in DSi and Si:TN between the mainstem and tributaries, with the highest concentrations and ratios in winter and late summer and declines in the spring and fall. Si:TP followed a distinct seasonal pattern from DSi and Si:TN; it was highest in the tributaries in the late fall through early winter and lowest in the spring and summer. Mainstem Si:TP ratios follow a similar pattern to the tributaries, with ratios also declining in the late fall (**Figure S1**). Across the 23 UMRS tributaries, basin land cover and lithology were highly variable (**Table S3**). Developed land cover ranged from 4.1 to 47.1%, forest cover ranged from 1.3 to 57.2%, wetlands ranged from 0.2 to 24.4, pasture/hay ranged from 2.9 to 27.8%, and cultivated cropland ranged from 10.4 to 81.6%. Open water accounted for <5% of basin coverage across all tributaries. Primary lithology showed even larger ranges, with basins having 0 to >75% coverage of each individual primary lithology (dolostone, limestone, shale, or sandstone), and other lithology (including schist, gneiss, granite, mafic, basalt, and clay) making up 0–47% of basins (**Table S3**).

#### H1: Land Use and Lithological Controls on Si Stoichiometry

The most appropriate regression models varied among various nutrients and ratios, and the role of land use and lithology differed in controlling tributary stoichiometry. Regression models that included both land use and lithology predictors together indicated tributary Si:TP ratios were explained by developed land cover and shale lithology (adj. R <sup>2</sup> = 0.42, df = 20, p < 0.01) (**Table 3**). Tributary Si:TN ratios were better explained by basin characteristics, specifically forest, developed, shale, and sandstone coverage (df = 18, p < 0.01), which explained 74% of Si:TN observations (adj. R <sup>2</sup> = 0.74). With regards to concentrations (rather than ratios), stepwise AIC identified only lithological, rather than land use, parameters as important model parameters for DSi concentrations (limestone and shale; df = 20, p < 0.05), whereas land use parameters were important predictors of TP and TN. Specifically, TP concentrations were best explained by developed, forest, and shale coverage (adj. R <sup>2</sup> = 0.38, df = 19, p < 0.01), and TN concentrations were best explained by developed, forest, shale, and sandstone basin coverage (adj. R <sup>2</sup> = 0.78, df = 17, p < 0.01).

In order to determine the relative roles of watershed land use vs. lithology on tributary Si:TN:TP stoichiometry, we also ran regression models with land use and lithology predictors separately. Together, we found Si:TP ratios were better explained by lithology alone than land cover (adj R <sup>2</sup> = 0.18 vs. 0.30 for land use vs. lithology, respectively), whereas Si:TN ratios were better explained by land cover (adj R <sup>2</sup> = 0.72 vs. 0.35 for land use vs. lithology, respectively). Regressions using both land use and lithology together indicate developed and shale coverage is important for both Si:TP and Si:TN, and forest and sandstone are additional predictors for Si:TN.

The four primary lithology cover classes represented in tributary basins were all significantly correlated with latitude. Limestone and shale follow a north-to-south gradient with an inverse correlation between percent cover and latitude (R 2 = 0.19, p < 0.05, n = 23 for both lithology types). Dolostone and sandstone follow the north-to-south latitudinal gradient but are positively correlated with latitude [R <sup>2</sup> = 0.17, p = 0.05 and R <sup>2</sup> = 0.21 p < 0.05 (n = 23), respectively]. Land use was not significantly correlated with latitude along the UMRS.

We observed varying amounts of stoichiometric coupling within the tributaries (**Figure 3**). The majority (52%) of tributaries showed relatively uncoupled Si:TP ratios, with the slopes of TP x DSi concentrations significantly differing from one. DSi to TP relationships vary from weakly to strongly coupled with R 2 values ranging from <0.001 to 0.82 (**Table 4**). Conversely to Si:TP, the majority of Si:TN ratios (61%) displayed coupled relationships, as indicated by slopes of TN vs. DSi not different than one. DSi to TN relationships vary from weakly to strongly coupled relationships with R 2 values ranging from 0.008


TABLE 3 | Multiple linear regression model results describing ability of land use and lithology to describe DSi, TN, and TP concentrations, as well as Si:TN and Si:TP ratios.

*Model types were chosen using Stepwise AIC. Models include land use and lithology predictors combined and individually. Bold indicates significant (p*<*0.05) predictors.*

to 0.528 (**Table 4**). In instances where slopes were significantly different from one, both TN and TP changed at faster rates relative to DSi.

#### H2: Si:TN:TP Varies Longitudinally Down the River

We observed significant changes in Si stoichiometry down the longitudinal continuum of the river, across both tributaries and the mainstem. With regards to the tributaries, average DSi concentrations and Si:TP ratios significantly (p < 0.01) decreased moving south across the watershed from northern Wisconsin and Minnesota to southern Illinois, whereas Si:TN did not change with latitude (**Figure 4**). The declines in tributary Si:TP ratios were a result of significant declines in DSi concentrations (R 2 = 0.45, p < 0.01) and increases in TP concentrations (R <sup>2</sup> = 0.38, p < 0.01) as we move south across the study region. The longitudinal TABLE 4 | Standardized major axis regressions results for the comparison of silicon (Si) to total phosphorus (TP) and Si to total nitrogen (TN) for each major tributary entering the Upper Mississippi River System.


patterns observed in the tributaries were likely driven by the significant shifts in lithology across the study basin, especially given lithology was such a strong predictor of DSi concentrations and stoichiometry. Differences in river size were not a likely driver of these longitudinal shifts, as we found no significant relationships between tributary stoichiometry and basin size (R 2 < 0.05 for all constituents).

In the mainstem, we also observed significant longitudinal shifts in Si stoichiometry; DSi concentrations, Si:TN, and Si:TP ratios all declined significantly (p < 0.01) from north to south (**Table 1** and **Figures 4**, **5**). The decline was steepest for DSi and Si:TP, although it was also statistically significant for Si:TN (**Figure 4**). Downstream declines in mainstem Si:TN and Si:TP were functions of decreasing DSi and increasing TN and TP concentrations (p < 0.01 for both TN and TP). The consistent decline in mainstem Si:TN resulted in ratios below Redfield downstream of Reach 3, whereas Si:TP shifted below Redfield prominently only in Reach 6 (**Figure 5**). DSi and TN were generally more closely coupled than DSi and TP (higher R 2 values; **Tables 5** and **6**), and in both cases the coupling was stronger in the more northern reaches.

Longitudinal differences in mainstem river conditions were captured along first PC axis (PC1), and explained the majority of variation in DSi concentrations and ratios (58%; **Figure 6**). These longitudinal differences in DSi concentrations, Si:TN, and Si:TP were inversely correlated with turbidity, temperature, and VSS (**Figure 6** and **Table S4**). Regression models showed CHL, VSS, and temperature were the most important predictors of main channel DSi concentration, Si:TP and Si:TN patterns (**Tables S5A,B**). However, these factors did not explain a great deal of variation in any case (Si—R2m = 0.20; Si:TN—R2m = 0.13; Si:TP—R2m = 0.17). The R<sup>2</sup> c of the full models, which reflect the random effect of the reach, explained more variation than these limnological factors alone. This suggests these limnological changes along the river continuum are not the dominant influence on longitudinal changes we observed in the mainstem (R<sup>2</sup> c: DSi = 0.53; Si:TN = 0.56, and Si:TP = 0.63).

#### H3: Hydrogeomorphology Explains Variation in Si:TN:TP Across Aquatic Areas

Aquatic areas were differentiated most strongly by depth and velocity (**Table 2** and **Figure 6**). Average depth was highest in the main channel (6.5 ± 0.04 m) and shallowest in the isolated backwater (0.91 ± 0.04 m). Average velocity was highest in the main and side channels (0.48 and 0.52 ± 0.01 m s−<sup>1</sup> , respectively) and lowest in the riverine lake and isolated backwater areas (<0.01 m s−<sup>1</sup> ). Turbidity, VSS, and temperature were generally highest in the main and side channels, whereas CHL was highest in backwaters and riverine lake areas (31.1 ± 0.39 and 27.1 ± 0.66 µg L−<sup>1</sup> , respectively; **Figure 6**). TN was significantly higher in main and side channels than contiguous and isolated backwaters and TP was highest in channels and lowest in the impounded area (**Table 2** and **Figure 6**).

Although differences among reaches were greater than among aquatic areas, DSi, Si:TN, and Si:TP varied among aquatic areas within the mainstem UMRS, indicating the importance of river

TABLE 5 | Standardized major axis regressions results for the comparison of dissolved silicon (DSi) to total phosphorus (TP) for each aquatic area and reach along the mainstem Upper Mississippi River System.


*Slope and intercept values were compared between river reaches within each aquatic area (where sufficient data is available).*

hydrogeomorphology in influencing stoichiometric patterns (**Figures 5**, **6**, **Table 2**, and **Table S4**). Pairwise comparisons showed DSi concentrations were similar among aquatic areas, except isolated backwaters, which had lower DSi than all other aquatic area types (**Table 2**, **Figure 7**, and **Table S6**). Si:TN was higher in the main channel, side channels and impounded areas than in contiguous and isolated backwaters (**Table 2** and **Figure 7**). Si:TP was highest in the riverine lake, followed by the impounded areas, main and side channels, and isolated backwaters (**Table 2** and **Figure 7**). Overall, DSi was generally more similar among aquatic areas than Si:TN or Si:TP.

Although most of the variation in DSi, Si:TN, and Si:TP was explained by longitudinal differences (PC1; **Figure 6**), DSi and Si:TN also differed across the aquatic area gradient (PC2 in **Figure 6** and **Table S4**). The "lateral" or residence time gradient of the river was captured by the second PC axis (PC2), which explained 25.3% of the variation in river conditions. PC2 shows the lateral gradient (i.e., aquatic areas) of the river was best described by differences in depth, velocity, and CHL concentration (**Figure 6** and **Table S4**). DSi was negatively correlated with CHL, and Si:TN was positively correlated with depth and velocity along PC2 (**Table S4**). Further, when we evaluated controls on DSi, Si:TN, and Si:TP in specific aquatic area types using regression models, CHL, VSS, turbidity consistently explained the most variation. However, these factors were able to explain more variation in DSi concentrations in the mainstem (R2m = 0.20) and impounded areas for Si:TN and Si:TP (R2m = 0.34 and 0.39, respectively) compared to contiguous backwaters, where little explanatory power of these factors was observed (R2m = 0.04, 0.08, and 0.15, for DSi concentrations, Si:TN, and Si:TP, respectively).

Across aquatic areas, we observed a relatively high degree of coupling in Si stoichiometry. That is, changes in one nutrient (i.e., DSi) corresponded with proportional changes in the other (i.e., TN or TP). Specifically, 83 and 96% of aquatic areas displayed relatively strong coupling of Si to TN and Si to TP, respectively, defined as having slopes not different than one (**Tables 5**, **6** and **Figure 5**). In general, Si was more coupled with TN, compared to TP, as indicated by tighter coefficients of determination (i.e., R 2 ) between TN × Si compared to TP × Si (**Tables 5**, **6**). However, Si stoichiometry showed weaker coupling in several aquatic areas, especially those with lower residence time. For example, connected backwaters showed weaker coupling of Si × TP (**Figure 5**), with R <sup>2</sup> < 0.1 in all cases (n = 5) and absolute slope values significantly different than one. In these cases, we observe an abundance of DSi relative to TP (**Table 5**). We observed variability across aquatic areas in Si stoichiometry in relation to Redfield ratios, with some areas showing a stronger degree of Si limitation. For example, main and side channels were more Si-limited with regard to N compared to other aquatic area types (**Figure 5**).

TABLE 6 | Standardized major axis regressions results for the comparison of silicon (Si) to total nitrogen (TN) for each aquatic area and reach along the mainstem Upper Mississippi River System.


*Slope and intercept values were compared between river reaches within each aquatic area (where sufficient data is available).*

# DISCUSSION

The stoichiometric relationships we observed among DSi, TP, and TN in the UMRS integrate both external inputs and internal processing of these nutrients (Maranger et al., 2018). We hypothesized basin land use and lithology would alter tributary inputs into the mainstem UMRS (H1), Si:TN:TP would vary longitudinally across the basin (H2), and variation in river hydrogeomorphology would result in heterogeneous Si stoichiometry among across aquatic area types (H3). We found significant downstream declines in average DSi concentrations, and Si:TN and Si:TP in the UMRS, and altered stoichiometric ratios associated with lower residence time. We explored potential causes for this observed variability, highlighting shifts in tributary inputs and the limnological complexity of the river (e.g., altered residence time, suspended sediment, and phytoplankton dynamics).

### Effects of Basin Land Use and Lithology on Tributary Inputs

We hypothesized basin land use and lithology would be strong predictors of river Si stoichiometry. Contrary to work in other temperate regions, our results indicate basin land cover was not a critical driver of tributary DSi concentrations. In fact, we found lithology, rather than land use, was a much stronger predictor of DSi concentrations in UMRS tributaries (**Table 4**). We suspect the strong signal of lithology here is due largely to the wide range of lithology relative to the differences in land use. For example, the New England (USA) watersheds, where land use was shown to be a strong control of river DSi exports, were one to three orders of magnitude smaller, and exhibited less variety in lithology and larger shifts in land use, compared to the tributary watersheds we study here (Carey and Fulweiler, 2012b).

However, basin land cover appeared to be an important driver of tributary ratios (Si:TN and Si:TP) (**Table 3**). Specifically, we observed a significant decline in Si:TN with increasing watershed agricultural land use. These patterns were driven by the significantly positive relationship between cultivated crop lands and TN concentrations, rather than altered DSi concentrations with cultivated land area, which is consistent with recent findings from the same region (Downing et al., 2016). Despite the significant relationship between basin cultivated crops coverage and Si:TN ratios in the tributaries, our regressions indicated forest and developed coverage were more important drivers of tributary nutrient stoichiometry compared to cultivated crop land, likely due to co-variation of cultivated crops with other land uses (see **Supplemental Information**). Thus, shifts in primary lithology across the longitudinal gradient of the system appear to be an important factor driving declining DSi concentrations at the mouths of tributaries southward across the region.

The fact that N is more responsive to basin land use, specifically cultivated crop area, is not surprising; the direct application of fertilizers on cultivated crops are known to impact N and, to a lesser extent, P loads, to the UMRS (Raymond et al.,

FIGURE 6 | Principal component analysis (PCA) showing variables that differentiate reaches (PC1) and aquatic areas (PC2). Length and direction of arrows indicates how strongly each variable loads along each axis. Variables shown have correlation coefficients of >0.60 with at least one axis. Points are colored by their reach and text label indicates aquatic area. Main, main channel; Side, side channel; Backwater, contiguous backwater; RL, riverine lake; Impound, impounded; and Iso-BW, isolated backwater lake. Variables: TN, total N; TP, total P; CHL, chlorophyll *a* concentration; TEMP, water temperature; VSS, volatile suspended solids; TURB, turbidity (NTU); Depth, depth; VEL, current speed.

2012; Robertson et al., 2014; Van Meter et al., 2018). Interestingly, recent work on tributaries in this basin showed increasing trends in N loading, but declining trends in P loading from agricultural watersheds (Kreiling and Houser, 2016). It is possible that the same processes that reduced P movement off the landscape (e.g., reduced erosion) may have also influenced the transport of Si, resulting in a non-significant relationship of DSi with land use in recent years.

#### Variation in Si:TN:TP Along the Length of the UMRS

We hypothesized we would observe varying Si stoichiometry along the longitudinal length of the river due to variable nutrient loading and processing along the UMRS. We observed significant declines in DSi concentrations, and stoichiometric ratios (Si:TN and Si:TP) in the mainstem across all aquatic area types from north to south across the study region. These patterns appear to reflect both changes in inputs and in-river processing. We saw a similar decline in external Si inputs from the tributaries, which similarly changed from north to south. In addition, increased runoff from human activities supplement concentrations of TN and TP along the mainstem (Kreiling and Houser, 2016). Other sources of N, such as from N deposition, also have influenced these patterns, given its similar longitudinal distribution (NADP, 2016). We can also attribute some of the longitudinal changes in DSi to downstream increases in CHL, temperature, and organic suspended material (**Table S5A** and **Figure 6**), suggesting instream uptake of Si by aquatic plants and diatoms play a role in shaping Si stoichiometry (Johnson and Hagerty, 2008; Manier, 2014; Crawford et al., 2016). Our results confirm previous work that has shown higher loads of TN relative to DSi and TP to the Gulf of Mexico, creating a more Si and P-limited system for marine phytoplankton (Houser et al., 2010; Royer, 2019).

### Effects of River Hydrogeomorphology on Variation in Si:TN:TP

We hypothesized Si stoichiometry would vary across aquatic area types due to differences in water residence time and other physical characteristics. Water residence time is a key driver of Si cycling in aquatic systems. Higher residence time can facilitate diatom blooms and subsequent Si burial and retention (van Bennekom and Salomons, 1980; Humborg et al., 1997, 2000, 2006; Downing et al., 2016). We found decreasing DSi concentrations with increasing residence time among aquatic areas. For example, contiguous and isolated backwaters, which had the highest relative residence times, displayed the lowest DSi concentrations compared to other aquatic areas (152 ± 1.4 and 126 ± 7.1µM, respectively vs. 158–195µM; **Table 2**). These areas also have high abundances of macrophytes, epiphytic algae, and higher water clarity that could facilitate benthic algal growth (Johnson and Hagerty, 2008; Kreiling and Houser, 2016). The uptake of DSi by primary producers other than phytoplankton in these environments likely reflects the combined role of multiple autotrophs, but their relative roles in Si cycling in rivers are not yet well-understood (Schoelynck et al., 2010, 2012).

The relationships of Si:TN and Si:TP varied more across aquatic area types compared to DSi concentrations alone (**Figure 7**), suggesting more dynamic processing of N and P compared to DSi. Both isolated backwaters and backwater lakes had statistically higher Si:TN than other aquatic area types. These differences are due to lower TN concentrations in isolated backwaters and backwater lakes, potentially due to increased denitrification rates in these areas and their lower connectivity to the main channel. Denitrification rates are likely higher in these systems as a result of higher sediment organic matter and lower sediment oxygen availability (Richardson et al., 2004; Downing et al., 2016). In contrast, Si:TP ratios were more variable along the residence time gradient (**Table 2**), although they were lowest in isolated backwater lakes, which was driven by the dramatically lower DSi concentrations in these aquatic areas. While anoxic conditions associated with the higher water residence times can foster P release from sediments (Wetzel, 2001; Houser and Richardson, 2010), we do not see that signal in our data. Rather, we observe the highest TP concentrations in main and side channel aquatic areas, which were the most turbid and have the shortest residence times. We suggest that the elevated TP concentrations may be due to P sorption to suspended particulates (Froelich, 1988). These same factors are likely responsible for the differences in coupling we observed across aquatic areas, where connected backwaters again behaved differently than other aquatic area types.

Consistent with prior studies from the region (Downing et al., 2016), we found no significant difference in DSi concentrations or Si:TN and Si:TP ratios between low head impoundments and main channel reaches of the UMRS (**Table 2** and **Figure 7**).

The lack of observed DSi depletion within these run-of-river impoundments is likely due to the relatively unaltered water residence time compared to most other aquatic areas (**Table 2**). An additional reason for the lack of DSi depletion behind the impoundments we study here may be due to the potential dominance of cyanobacteria, rather than diatoms, in the system (Manier, 2014). In the impoundments studied by Downing et al. (2016) in the same region, diatoms made up <10% of the phytoplankton community, with cyanobacteria dominating. In cases where non-siliceous algae dominate over diatoms, DSi uptake and burial is limited, resulting in a lack of change in Si stoichiometry in these systems.

#### Implications for Coastal Receiving Waters

The implications of nutrient concentrations and their ratios supplied by freshwater rivers to coastal receiving waters can range from minor shifts in the community composition of diatoms to large blooms of toxic algae and extensive hypoxia, all of which have direct implications for marine ecology and C cycling. Early work by Smayda (1990) suggested decreases in the Si:N ratio below the canonical Redfield 1:1 ratio (Redfield, 1963; Brzezinski, 1985) would shift phytoplankton communities entirely away from diatoms and would promote nuisance blooms dominated by haptophytes (e.g., Phaeocystis) and dinoflagellates (e.g., Alexandrium). As such, comparing the UMRS stoichiometry to the Redfield ratio can indicate the potential for nutrient limitation of algal communities in receiving waters. Stoichiometric ratios in the UMRS basin varied above and below the Redfield ratio (1:1) for Si:TN, dependent on season. During the winter and spring (January through June), average Si:TN ratios are consistently below Redfield, indicating possible Si-limitation relative to N in the system. Conversely, Si:TP ratios are consistently above Redfield throughout the year. Declines in DSi concentrations in the spring and fall could be attributed to increased discharge and dilution, as well as biological uptake and assimilation of DSi in terrestrial plants and diatoms (Humborg et al., 2006; Carey and Fulweiler, 2013). We also observed variation in ratios spatially in the UMRS, with increasing Si limitation relative to TN and TP as we move south across the study region. Downstream of Reach 3, Si becomes limiting relative to TN, and Si approaches limitation relative to TP at the southern end of our study region (Reach 5) (**Figure 4**). On average the main channel is DSi limited across all reaches relative to N, likely due to high N inputs.

Although the Redfield ratio is useful for comparison, it was developed for marine algae, and freshwater diatoms can require up to an order of magnitude greater Si than marine species (Conley et al., 1988). Moreover, recent work suggests extreme flexibility in the Si:N ratio of diatoms; some genera, such as Pseudo-nitzchia, can vary by 15-fold in their Si quota and thus, form blooms under highly divergent Si:N ratios (Pan et al., 1996; Davidson et al., 2012). While Si:N ratios may not strictly determine the proportion of diatoms in a community to the extent previously believed, low Si:N ratios can favor weakly silicified diatoms, which can have significant ecosystem-wide effects. For example, weaker diatom silicification can increase grazing efficiency (Liu et al., 2016; Zhang et al., 2017), potentially affecting trophic transfer efficiency, food web dynamics, and ecosystem services. In addition, lower Si content can result in slower sinking rates, with direct implications for the rate of vertical C export and bottom-water hypoxic zone formation (Turner et al., 1998). Finally, factors affecting diatom growth rates (e.g., light or temperature) also likely affect cellular Si concentration since the rate of Si uptake is dependent upon the rate of cellular growth (Taylor, 1985). In turn, any driver of global change that alters diatom growth rates may alter diatom Si content, independent from the stoichiometric ratios of supplied nutrients. Therefore, even small deviations from balanced stoichiometry as a result of riverine nutrient processing, in combination with other facets of global change, could exacerbate shifts toward potentially harmful algal communities and alter ecological functioning.

In order to evaluate whether the ratios in the UMRS are representative of nutrient exports to coastal receiving waters from the entire basin, we compared our ratios from our southernmost study site (Reach 5) to those at the mouth of the Mississippi River as it enters the Gulf of Mexico (Belle Chasse station: USGS 07374525) for the 2010–2018 period. We found very similar Si stoichiometry between the two stations, with both Si:TN and Si:TP showing only slight changes at the mouth (0.83 and 15.3, respectively) compared to Reach 5 (0.79 and 19.6, respectively). This indicates the majority of external drivers and internal modulators of nutrients relevant for coastal waters may occur in the UMRS portion of the basin. Although the complex hydrology of the UMRS may differentially influence the loads of DSi, TN, and TP (Smits et al., 2019), we only considered concentrations in our analyses because of the difficulty in calculating loads across diverse hydrogeomorphic areas. Our results agree with recent work by Royer (2019), which shows that although the UMRS represents 15% of the total Mississippi-Atchafalaya River Basin, it explains the majority (57%) of the variation in N loading at the Gulf of Mexico. This highlights the importance of studying Si:TN:TP in the UMRS, as this stoichiometric signal is reflected in the Si stoichiometry entering the Gulf of Mexico.

In general, the values we observed show the UMRS has lower relative Si:TN, but higher relative Si:TP compared to rivers worldwide. Compared to values from the GEMS-GLORI database, which represents data from 126 rivers on six continents, UMRS Si:TN ratios fall within the 18–32nd percentile, and Si:TP values fall within the 50–97th percentile of the values in the database (Garnier et al., 2010) (see **Supplemental Information**). This indicates that the UMRS has relatively high TN concentrations compared to other systems, and the combined processes that shape Si:TN stoichiometry act in concert to deliver a relatively Si- and P-limited nutrient supply to downstream systems.

Together, our results show rivers are not simple pipes for Si, but rather the complexity in watershed characteristics, hydrology, and biological uptake result in dynamic Si cycling along the river continuum. Shifting Si:TN:TP ratios are driven by basin land cover and lithology and local limnological conditions, such as water residence time, turbidity, and inorganic nutrient concentrations. Although there is a great deal of work on Si stoichiometry on rivers globally, our manuscript is unique in evaluating Si variability as a function of river hydrogeomorphology and across the river continuum. The conceptual framework of this paper translates well to understanding the drivers of variable Si stoichiometry in large rivers globally. This enhanced understanding of the causes of variation in Si:TN:TP stoichiometry can inform management decisions by shedding light on the relative extent to which these ratios are driven by natural vs. anthropogenic factors.

# DATA AVAILABILITY

All data are publicly available for download at the following https://umesc.usgs.gov/ltrm-home.html.

# AUTHOR CONTRIBUTIONS

The preparation and creation of this manuscript was a strong collaborative effort and JC, KJ, PJ, LS and PT contributed equally to the effort, and are ordered alphabetically. Specifically, JC, KJ, PJ, LS and PT conceived the idea, performed the statistical analyses and wrote the paper. JR performed spatial analyses and created all maps.

# FUNDING

This paper was a product of the workshop Woodstoich 2019. Financial support was provided by National Science Foundation (DEB-1840408). KJ was funded by the US Army Corps of Engineers' Upper Mississippi River Restoration Program as were many of the datasets used in our analyses.

# ACKNOWLEDGMENTS

We would like to thank Kate Henderson, Dr. Todd Royer, Dr. Zoe Cardon, and the peer reviewers for constructive comments on earlier versions of this manuscript. We thank Drs. Jim Elser and Michelle A. Evans-White for organizing Woodstoich4 (keeping it groovy). We acknowledge University of Montana's Flathead Lake Biological Station, where we spent an intensive week finalizing this manuscript. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00346/full#supplementary-material

#### REFERENCES


Barton, K., and Barton, M. K. (2013). Package 'MuMIn.' Version 1, 18.


mixed-effects models revisited and expanded. J. R. Soc. Interface 14:20170213. doi: 10.1098/rsif.2017.0213


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Carey, Jankowski, Julian, Sethna, Thomas and Rohweder. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Edited by:

Michelle Evans-White, University of Arkansas, United States

#### Reviewed by:

Rana El-Sabaawi, University of Victoria, Canada Edward Hall, Colorado State University, United States James Cotner, University of Minnesota Twin Cities, United States

\*Correspondence:

Scott Hotaling scott.hotaling@wsu.edu

†These authors have contributed equally to this work

‡Present address: Isabella A. Oleksy, Cary Institute of Ecosystem Studies, Millbrook, NY, United States

#### §ORCID:

Ze Ren orcid.org/0000-0002-1473-2697 Nicolas Martyniuk orcid.org/0000-0002-2423-7040 Isabella A. Oleksy orcid.org/0000-0003-2572-5457 Anshuman Swain orcid.org/0000-0002-9180-2222 Scott Hotaling orcid.org/0000-0002-5965-0986

#### Specialty section:

This article was submitted to Behavioral and Evolutionary Ecology, a section of the journal Frontiers in Ecology and Evolution

Received: 22 July 2019 Accepted: 11 September 2019 Published: 26 September 2019

#### Citation:

Ren Z, Martyniuk N, Oleksy IA, Swain A and Hotaling S (2019) Ecological Stoichiometry of the Mountain Cryosphere. Front. Ecol. Evol. 7:360. doi: 10.3389/fevo.2019.00360

# Ecological Stoichiometry of the Mountain Cryosphere

Ze Ren1†§, Nicolas Martyniuk 2†§, Isabella A. Oleksy 3†‡§, Anshuman Swain4†§ and Scott Hotaling<sup>5</sup> \* §

<sup>1</sup> Flathead Lake Biological Station, University of Montana, Polson, MT, United States, <sup>2</sup> Laboratorio de Limnología, INIBIOMA, CONICET—Universidad Nacional del Comahue, Bariloche, Argentina, <sup>3</sup> Natural Resource Ecology Laboratory, Colorado State University, Fort Collins, CO, United States, <sup>4</sup> Department of Biology, University of Maryland, College Park, MD, United States, <sup>5</sup> School of Biological Sciences, Washington State University, Pullman, WA, United States

Roughly 10% of the Earth's surface is permanently covered by glaciers and ice sheets and in mountain ecosystems, this proportion of ice cover is often even higher. From an ecological perspective, ice-dominated ecosystems place harsh controls on life including cold temperature, limited nutrient availability, and often prolonged darkness due to snow cover for much of the year. Despite these limitations, glaciers, and perennial snowfields support diverse, primarily microbial communities, though macroinvertebrates and vertebrates are also present. The availability and mass balance of key elements [(carbon (C), nitrogen (N), phosphorous (P)] are known to influence the population dynamics of organisms, and ultimately shape the structure and function of ecosystems worldwide. While considerable attention has been devoted to patterns of biodiversity in mountain cryosphere-influenced ecosystems, the ecological stoichiometry of these habitats has received much less attention. Understanding this emerging research arena is particularly pressing in light of the rapid recession of glaciers and perennial snowfields worldwide. In this review, we synthesize existing knowledge of ecological stoichiometry, nutrient availability, and food webs in the mountain cryosphere (specifically glaciers and perennial snowfields). We use this synthesis to develop more general understanding of nutrient origins, distributions, and trophic interactions in these imperiled ecosystems. We focus our efforts on three major habitats: glacier surfaces (supraglacial), the area beneath glaciers (subglacial), and adjacent downstream habitats (i.e., glacierfed streams and lakes). We compare nutrient availability in these habitats to comparable habitats on continental ice sheets (e.g., Greenland and Antarctica) and show that, in general, nutrient levels are substantially different between the two. We also discuss how ongoing climate warming will alter nutrient and trophic dynamics in mountain glacier-influenced ecosystems. We conclude by highlighting the pressing need for studies to understand spatial and temporal stoichiometric variation in the mountain cryosphere, ideally with direct comparisons to continental ice sheets, before these imperiled habitats vanish completely.

Keywords: alpine, glacier biology, nutrient dynamics, C:N, supraglacial, subglacial, glacier-fed lakes, glacier-fed streams

**40**

# INTRODUCTION

Approximately 75% of the freshwater on Earth is frozen at any given time (Jain, 2014; Talalay et al., 2014). This global "cryosphere" includes continental ice sheets, mountain glaciers, snowfields, permafrost, and sea ice. Beyond extreme cold and frequent subzero temperatures, frozen environments place additional harsh controls on life including a lack of available water, high ultraviolet radiation, and highly dynamic, but generally limited, nutrient supplies (Anesio and Laybourn-Parry, 2012). Nevertheless, the cryosphere supports diverse biological communities on the surface of glaciers and continental ice sheets (Anesio et al., 2009; Hotaling et al., 2017a), beneath them (Hamilton et al., 2013), in their meltwater (Hotaling et al., 2019a), in permafrost (Jansson and Tas, 2014), and in sea ice (Boetius et al., 2015). Even though these biotic communities are dominated by unicellular microbial life (e.g., bacteria and algae: Boetius et al., 2015; Anesio et al., 2017; Hotaling et al., 2017a), macroinvertebrates (e.g., ice worms and rotifers; Shain et al., 2016; Hotaling et al., 2019b) and vertebrates (e.g., gray-crowned rosy finches; Rosvold, 2016) are also represented. However, rising global temperatures are driving widespread recession of the cryosphere (Lyon et al., 2009; Notz and Stroeve, 2016; Roe et al., 2017) with profound implications for biodiversity (Hotaling et al., 2017b), human populations (Pritchard, 2017), and ecological feedbacks (e.g., trophic interactions and biogeochemical cycles, Hood et al., 2009; Sommaruga, 2015; Milner et al., 2017). Therefore, understanding the interplay between resource availability and biological communities across cryosphere-associated habitats, as well as how climate change may alter them, is a pressing research need.

The field of ecological stoichiometry, which seeks to understand the balance of energy and chemical elements in ecological interactions and processes (Sterner and Elser, 2002), provides a valuable framework for understanding how ecosystems function. The availability and mass balance of key elements [carbon (C), nitrogen (N), and phosphorous (P)] place crucial controls on the population dynamics of organisms (Cross et al., 2005), ultimately regulating the structure, function, and processes of ecosystems (Sterner and Elser, 2002; Elser et al., 2007). On mountain glaciers and snowfields, solar radiation, atmospheric CO2, allochthonous (e.g., arthropod fallout, Edwards, 1987; black carbon, Skiles et al., 2018), and autochthonous (e.g., cryoconite hole metabolism, Anesio et al., 2009; snow algae primary productivity, Hamilton and Havig, 2018) inputs all interact to shape elemental balances and the structure of food webs. Furthermore, glacier environments support higher trophic levels (Tynen, 1970; Kohshima, 1984; De Smet and Van Rompu, 1994; Kikuchi, 1994), highlighting the key role that stoichiometric dynamics likely play in shaping ice-associated food webs. Spatial links also exist among cryospheric habitats; for instance, primary production and atmospheric deposition occuring on the surface of glaciers (supraglacial zone) influences stoichiometric processes beneath them in the subglacial zone. Glacier-associated biological and biogeochemical processes ultimately affect downstream lakes and streams as far as the world's oceans (Hotaling et al., 2017a). Thus, ecological stoichiometry can provide mechanistic insight into connections among different environments, interactions between trophic levels, and can act as a means to link the influences of glacier melt dynamics on ecosystem structure and function in cryosphere-associated environments and beyond.

In this review, we provide the first synthesis of ecological stoichiometry in the mountain cryosphere. We identify mountains (and by proxy, mountain ecosystems) in the same way as a related review (Moser et al., 2019), which used the definition provided by Körner et al. (2017): mountainous areas exhibit more than 200 m of elevation difference within a 2.5′ grid cell, irrespective of elevation. Our motivation for focusing on the mountain cryosphere was 3-fold. First, much is known of the general stoichiometry of continental ice sheets (i.e., Greenland and Antarctica), however there has been considerably less focus on similar mountain ecosystems and thus, it is unclear how comparable the two environments are. Second, mountain glaciers and snowfields harbor extensive biological diversity, including multiple trophic levels from microbes to birds (e.g., Anesio and Laybourn-Parry, 2012; Hotaling et al., 2019a). Third, many human populations rely on mountain glacier meltwater for agriculture, hydropower, and drinking water (Milner et al., 2017), highlighting the potential relevance of ecological processes occurring in headwaters beyond adjacent biotic communities.

Within the mountain cryosphere, we focus on three distinct, but interconnected, habitats: the supraglacial zone, subglacial zone, and downstream glacier-fed streams and lakes (**Figure 1**). The supraglacial zone encompasses the interface between surface ice/snow and atmospheric conditions. The subglacial zone comprises the area where glacier ice interacts directly with bedrock and meltwater. And, glacier-fed streams and lakes include the downstream habitats which are directly adjacent to glaciers (**Figure 1**). By considering the mountain cryosphere in a stoichiometric framework, we clarify key environment-organism interactions, including the effects of nutrients on psychrophiles (organisms capable of living in extremely low temperatures) and their reciprocal effects on biogeochemical processes (e.g., respiration, nitrification) and habitat characteristics (e.g., snow albedo, light attenuation). For each habitat, we address two stoichiometric questions: (1) What is the origin, availability, and variation (spatial and temporal) of C, N, and P? What dynamics of stoichiometric ratios exist? (2) To what extent are complex food webs (i.e., multiple trophic levels) present? What does stoichiometry mean for their trophic interactions and biogeochemical cycling? As part of this, we explicitly test the degree to which nutrient concentrations in the mountain cryosphere are comparable to those measured for continental ice sheets to better understand how insight gained from ice sheets can be translated to mountain systems. We conclude by discussing how contemporary climate change, and particularly rapid glacier decline, will impact these ecosystems. As the cryosphere continues to change, a stoichiometric understanding will allow for refined ecological predictions of nutrient dynamics and consequences for trophic interactions.

NH<sup>+</sup> 4 , ammonium; CH4, methane; NO<sup>−</sup> 3 , nitrate. Figure is not drawn to scale.

# LIVING ON ICE: THE SUPRAGLACIAL ZONE

During melt, considerable primary production occurs in the supraglacial zone, primarily due to snow algal blooms (e.g., Chlamydomonas nivalis and cyanobacteria, Anesio et al., 2017), in spite of high levels of ultraviolet radiation (Morgan-Kiss et al., 2006; Yallop et al., 2012) and often limited nutrient availability (Hawkings et al., 2016; Wadham et al., 2016). This biomass production has key implications for local biota (e.g., heterotrophs), hydrologically connected habitats downstream, and even glacier albedo (Ganey et al., 2017; **Figure 1**). While surface ice and snow represent the bulk of available habitat in the supraglacial zone, and thus are where most of the biomass on mountain glaciers resides, the supraglacial zone is highly heterogeneous (e.g., with respect to debris cover) and habitat types within it (e.g., cryoconite holes) must also be considered (**Figure 1**).

#### Nutrients

Supraglacial environments are typically nutrient poor (**Figure 2**; **Table S1**) and available nutrients are often present in organic forms (e.g., residues of microbial and/or plant cells, decayed organic matter, microbial exudates; Antony et al., 2017). For reference eutrophic lakes typically have total N concentrations of 650–1,200 µg L−<sup>1</sup> and total P concentrations of 30– 100 µg L−<sup>1</sup> (Dodds and Whiles, 2010). Mountain glaciers typically only contain a fraction of that level of nutrient concentrations (**Figure 2**). Surface ice and snow in particular exhibit low nutrient concentrations when compared to other microhabitats (e.g., cryoconite holes; **Figure 2**). Through C fixation, supraglacial primary production promotes the accumulation of autochthonous organic C (Musilova et al., 2017; Williamson et al., 2019), thereby driving C and N cycling (Stibal et al., 2012; Havig and Hamilton, 2019). Due to near constant melt during periods of warm temperature and/or high solar radiation, any C that is fixed on the supraglacial surface and not consumed by local heterotrophs (e.g., bacteria, fungi, ice worms) is exported to subglacial and downstream habitats. Thus, metabolic activity of supraglacial organisms may alter the environmental availability of nutrients and shape trophic relationships in hydrologically connected habitats.

The presence and growth of supraglacial microbes, and in particular snow algae (**Figure 1**), can reduce local albedo and induce additional melting (Ganey et al., 2017). Since meltwater itself is often a limited resource in the supraglacial zone, the release of additional water can locally stimulate further microbial growth, which drives more melt in a bio-albedo feedback loop

(Anesio et al., 2017). Melting can also promote the formation of cryoconite holes in the supraglacial zone. Cryoconite holes form when particles with a lower albedo than the surrounding area induce locally elevated melt rates which in turn form small melt depressions (< 10 cm deep and < 1 m wide)—"cryoconite holes"—on the glacier surface (Fountain et al., 2004; **Figure 1**). These unique microhabitats are hotspots of microbial activity (Anesio and Laybourn-Parry, 2012) but also typically harbor larger organisms (e.g., De Smet and Van Rompu, 1994).

Overall, average dissolved organic carbon (DOC) concentrations are significantly higher in mountain glacier ice and snow than the same habitats on continental ice sheets (**Figure 2**; **Table S2**). While the mechanism(s) underlying this pattern remain unclear, it likely exists because mountain glaciers are generally closer to terrestrial sources of C vs. continental ice sheets. Related factors may also include higher ambient temperatures (and thus more meltwater) as well as more primary production on mountain glaciers. In general, supraglacial ecosystems accumulate organic matter through in situ primary production as well as deposition from terrestrial and anthropogenic sources (Hood et al., 2009; Singer et al., 2012; Stibal et al., 2012; **Figure 2**). Both autochthonous and allochthonous dissolved organic matter (DOM) are actively transformed by microbial communities through degradation and synthesis (Antony et al., 2017). Through this microbial processing, low quality DOM (i.e., high C:N, C:P) can be processed and "upgraded" (converted to lower C:N and/or C:P), thereby becoming more bioavailable to heterotrophs both within and downstream of the supraglacial zone (Musilova et al., 2017). Due to the tight recycling of DOM in the supraglacial zone, DOM in downstream systems can still be dominated by allochthonous C sources derived from aerial and terrestrial deposition (Stubbins et al., 2012). However, recent work from the outflow of surface and rock glaciers in the western United States demonstrated that exported DOM is still readily consumed by bacteria (Fegel et al., 2019). This suggests that even among mountain glacier ecosystems, there is substantial variation in exported DOC quality and quantity.

The diversity and activity of supraglacial microorganisms also exert a major influence on the cycling of N and P (Tranter et al., 2004; Hodson et al., 2005). Supraglacial meltwaters are low in nutrients (Säwström et al., 2007a; Grzesiak et al., 2015; **Figure 2**) and bacterial activity on glaciers is largely P-limited (Mindl et al., 2007; Säwström et al., 2007b; Stibal et al., 2008). This is reasonable given that P concentrations in snow, ice, and supraglacial streams are typically very low with an average of 0.005, 0.043, and 0.002 mg L−<sup>1</sup> , respectively (**Figure 2**). Cryoconite holes are an exception and exhibit higher P concentrations (0.31 mg L−<sup>1</sup> ), likely due to leakage of nutrients from algal cells, cell lysis, and general breakdown of organic matter (**Figure 2**). Thus, P demand likely exceeds supply in the supraglacial zone (except in cryoconite holes) because atmospheric deposition and microbial rock weathering alone cannot support existing needs (Mueller et al., 2001; Stibal et al., 2009). However, in addition to terrestrial inputs, debris cover on glaciers, and weathering from atmospheric deposition or rock fall can also supply P (Modenutti et al., 2018). These additional P sources can greatly alter nutrient dynamics in glaciers, and therefore productivity. However, since concentrations of P are much less frequently reported, spatial variation in P and the degree to which additional P inputs alter supraglacial productivity on local, regional, and global scales remains unclear.

Supraglacial microbial communities cycle N through the oxidation of ammonium to nitrate, chemical decomposition of organic matter (Hodson et al., 2005; Wynn et al., 2007), and nitrogen fixation (Telling et al., 2011). Terrestrial debris and anthropogenic N deposition are significant sources of reactive N to glacier surfaces worldwide (Anderson et al., 2017; Havig and Hamilton, 2019); this is particularly true in mountainous areas located near centers of industrial and/or agricultural activity and may partially explain the large differences in nitrate (NO<sup>−</sup> 3 ) and ammonium (NH<sup>+</sup> 4 ) concentrations between mountain glaciers and continental ice sheets (**Figure 2B**). For example, rates of total N deposition in the Rocky Mountains (∼3–5 kg N ha−<sup>1</sup> yr−<sup>1</sup> ; NRSP-3, 2019) and the southern Andes (∼8.2 kg N ha−<sup>1</sup> yr−<sup>1</sup> ; Godoy et al., 2003) are nearly an order of magnitude higher than Greenland (∼0.5–1 kg N ha−<sup>1</sup> yr−<sup>1</sup> ; AMAP, 2006). Since deposition is the most important factor influencing supraglacial N (Tranter et al., 1993; Hodson et al., 2005), the highest N concentrations are typically observed in snow and ice (**Figure 2**). Most N fixation occurs in the ablation zone of glaciers, where more ice mass is lost due to melting, evaporation, and related processes than accumulates via new snowfall. The zone of N fixation progressively migrates upslope during melt seasons as temperatures rise, solar radiation increases, and lower elevation N reserves are taken up (Telling et al., 2011). Due to high microbial metabolism, cryoconite holes act as sinks for dissolved inorganic N, and typically exhibit the lowest levels of dissolved inorganic N (i.e., ammonium and nitrate) in the supraglacial zone (**Figure 2**). Breakdown of organic matter, either from allochthonous deposition or historical autochthonous production, within cryoconite holes is likely an important additional source of recycled N to microbial communities (Stibal et al., 2012).

#### Trophic Interactions

To a greater extent than continental ice sheets, mountain glaciers support multi-trophic food webs which often include microbial diversity, macroinvertebrates, and, in some instances, vertebrates. The structure of supraglacial food webs varies geographically and the reason for this remains largely unknown (Zawierucha et al., 2015; Hotaling et al., 2019a). Recently, the gut microbiome of the Patagonian dragon (Andiperla willinki), a ∼2.5 cm stonefly which lives on the surface of Patagonian glaciers, revealed glacier-specific community modifications (Murakami et al., 2018). Indeed, bacterial taxa associated with glacier ice (e.g., Polaromonas and Rhodoferax) were present and are speculated to contribute to both host nutrition as well as the recycling and breakdown of organic matter where A. willinki is present. The gut bacteria-host relationship between A. willinki and glacier microbiota highlights the potential for direct trophic relationships between glacier-endemic bacteria and eukaryotes (Murakami et al., 2018).

In North America, millions of glacier ice worms (Mesenchytraeus solifugus) inhabit coastal glaciers from Oregon to Alaska, USA (Hotaling et al., 2019a). Ice worms feed on glacier bacteria and snow algae (Murakami et al., 2015) and are, in turn, heavily predated upon by a variety of highelevation nesting birds, including Gray-crowned Rosy finches (Leucosticte tephrocotis, Hotaling et al., 2019a). Thus, nutrient limitations occurring at lower trophic levels have the potential to be transferred from microbial communities to vertebrates. This could be particularly true on mountain glaciers because stoichiometric flexibility may be reduced at low temperatures (Godwin and Cotner, 2015). The potential for nutrient variation at the base of the food web to alter higher level trophic ecology, however, is not specific to glaciers in North or South America. In addition to A. willinki and M. solifugus, at least three other species of macroinvertebrates live on glaciers, including another ice worm in Tibet (Sinenchytraeus glacialis, Liang, 1979) and two chironomid midges, one in the Himalayas (Diamesa kohshima, Kohshima, 1984) and another in New Zealand (Zealandochlus latipalpis, Boothroyd and Cranston, 1999). Vertebrates feeding on glacier-endemic invertebrates may provide nutrient subsidies to supraglacial zones, potentially driving ecological stoichiometry and food web structure in a manner similar to the role of bird guano on remote islands (e.g., Anderson and Polis, 1999; Vizzini et al., 2016). Collectively, the presence of these unique, understudied organisms highlights the need for global, multitrophic studies of ecological stoichiometry in mountain glacier ecosystems to better understand how nutrient inputs (and limitations) flow through these vanishing habitats.

## LIVING BELOW ICE: THE SUBGLACIAL ZONE

Though perpetually dark and long assumed to contain little to no life, an extensive microbial community has been documented below mountain glaciers and continental ice sheets around the world (Hamilton et al., 2013; Anesio et al., 2017; Hotaling et al., 2017a; **Figure 1**). The subglacial zone exists at the ice-bedrock interface and includes any habitat that is perpetually covered by glacier ice. However, current understanding of biodiversity, nutrient availability, and associated stoichiometric composition of life in these habitats remains limited compared to supraglacial and downstream environments. This lack of understanding is likely due in large part to the extreme challenges associated with sampling habitats below tens to hundreds of meters of ice in rugged terrain (Anesio et al., 2017). To our knowledge, all studies of microbial ecology and ecological stoichiometry in mountain subglacial environments have only included samples collected near glacier snouts where subglacial runoff and sediments (e.g., within ice caves) are accessible rather than coring through ice to the bedrock itself.

#### Nutrients

During the melt season, water flows over and within the surface of glaciers, inducing a dynamic drainage network that moves water from the glacier surface to the base where much of it flows into the subglacial zone, and eventually into headwater streams and lakes (Fountain et al., 2005). The subglacial zone is rich in biogeochemical activity (Hamilton et al., 2013) that mobilizes nutrients and organic matter of both supraglacial and subglacial origin. Through this processing, microbially produced DOC is exported downstream, ultimately contributing to cycling of C in downstream and adjacent ecosystems (e.g., glacier forefields and headwater streams; Anesio et al., 2009, 2010; Edwards et al., 2014). Unlike continental ice sheets, the steep gradients of mountain systems generally preclude the development of subglacial lakes (though they can occur, e.g., Capps et al., 2010). However, substantial water is still present at the bedrock-ice interface which provides key habitat and resources for an active, diverse subglacial microbiome (Hamilton et al., 2013). Indeed, flowing water is the primary driver of nutrient movement both within glaciers (e.g., from the surface to subglacial habitats, Hotaling et al., 2017a) and into glacier outflows (e.g., delivery of labile C to glacier-influenced streams, Hood et al., 2009). Subglacial microbial communities may have even acted as refugia for microbial biodiversity during glacial periods due to the relative stability of subglacial habitats compared to their surface counterparts paired with a continual exposure to fresh mineral surfaces due to the grinding of bedrock by glaciers (Hodson et al., 2008).

An array of energy sources sustain life beneath glaciers (Hotaling et al., 2017a) including the aforementioned bedrock grinding and supraglacial inputs (e.g., primary productivity) flowing to the subsurface, as well as in some instances, geothermal energy (Hodson et al., 2008; Boyd et al., 2014; Telling et al., 2015). There is also considerable evidence for chemolithoautotrophic primary productivity (e.g., via methane production pathways; Boyd et al., 2010; Hamilton et al., 2013). In polar and subpolar regions, chemolithoautotrophs can fix several micrograms of C per square meter per day (Christner et al., 2014). Though C fixation by chemolithoautotrophs has not been estimated for any mountain glacier ecosystem, their abundance below mountain glaciers (e.g., Hamilton et al., 2013) suggests this activity likely rises to non-negligible levels with direct implications for food webs both near and far. In the Canadian Rockies, subglacial primary production is likely sustained in large part by the oxidation of pyrite and nitrification (Boyd et al., 2011, 2014). However, the full scope of bedrock lithologies in mountain ecosystems remains unsampled, making it difficult to generalize beyond specific study regions and local geology (Hotaling et al., 2017a). For N, microbially mediated nitrification and denitrification occur in subglacial environments (e.g., Boyd et al., 2010) and similar to C, its export during times of glacial melt affects downstream communities (see below; Hodson et al., 2010).

#### Trophic Interactions

Subglacial food webs appear to generally be sustained by microbially-mediated chemical weathering of bedrock through sulfide oxidation (e.g., Wadham et al., 2010; Boyd et al., 2014) or carbonic acid weathering (e.g., Havig and Hamilton, 2019). This chemical weathering releases organic C (Wadham et al., 2004) as well as N and P (Hodson et al., 2005) from the bedrock, thereby providing crucial nutrients. However, the degree to which nutrient limitations affect subglacial food webs, whether chronically or seasonally, remains largely unknown. For instance, bacterial and archaeal sediment communities beneath a midlatitude, temperate Canadian glacier, with a ∼210:1 C:N of dissolved organic matter, appear N limited (Boyd et al., 2011). A particulate organic C:N of ∼137:1 beneath the same Canadian glacier is elevated relative to comparable geologies (Ingall et al., 1993), whether or not this elevated C:N ratio is the product of contemporary or historical processes is not known (Boyd et al., 2011). Furthermore, the degree to which nutrients influence higher trophic levels (e.g., fungi, the likely largest organisms in subglacial ecosystems; Hotaling et al., 2017a) below glaciers has not been investigated.

## LIVING DOWNSTREAM OF ICE: GLACIER-FED STREAMS AND LAKES

Glacier-fed streams and lakes reflect the integration of both upstream and in situ processes (**Figure 1**; Brittain and Milner, 2001; Mindl et al., 2007; Robinson et al., 2016; Hotaling et al., 2017b; Ren and Gao, 2019). Melting glaciers supply key water, sediment, and nutrients to aquatic ecosystems, determining their optical properties, resource availability, and the structure and composition of biological communities (Laspoumaderes et al., 2013; Martyniuk et al., 2014; Rose et al., 2014; Hotaling et al., 2017b; Ren et al., 2017a; **Figures 4**, **5**). In general, glacier-fed streams and lakes are characterized by considerable sediment input from upstream grinding of bedrock, cold temperatures (e.g., < 10◦ even in summer), and dynamic water levels. Climate change is accelerating glacier retreat (Masiokas et al., 2008; Zemp et al., 2015) and will alter downstream ecosystems from biogeochemical cycling to elemental flow through food webs (Baron et al., 2009; Slemmons and Saros, 2012; Fell et al., 2017; **Figures 4**, **5**).

#### Glacier-Fed Streams Nutrients

Glaciers export organic C, inorganic and organic N, and soluble P to adjacent stream environments (Fegel et al., 2016; Milner et al., 2017; Colombo et al., 2019). In headwaters, glaciers are a major source of DOC (Hood et al., 2009; Milner et al., 2017; Hemingway et al., 2019), which is assimilated into multiple levels of stream food webs (Fellman et al., 2015). Concentrations of DOC in glacier-fed streams vary widely from near zero to more than 3 mg L−<sup>1</sup> (Zah et al., 2001; Hood et al., 2009; Wilhelm et al., 2013; Martyniuk et al., 2014; Robinson et al., 2016; Hemingway et al., 2019) and are generally comparable to values observed for snow and ice (**Figures 2**, **3**). Previously ice-locked DOC is particularly valuable to glacier-fed stream food webs as it is highly bioavailable for heterotrophic consumption (i.e., composed of > 50% bioavailable DOM; Hood et al., 2009, 2015; Singer et al., 2012; Fegel et al., 2019).

Three substantial pools of N exist near glacier-fed streams and lakes: glacier ice, subglacial sediment, and nearby soils (Malard et al., 2000; Hood and Scott, 2008; Robinson et al., 2016). In glacier-fed streams, nitrate concentrations range from 0.08 to 0.22 mg L−<sup>1</sup> , which is lower than supraglacial concentrations, potentially reflecting uptake upon export (Rinke et al., 2001; Robinson et al., 2016; **Figure 3**). While soluble reactive P is often below the detection limit in glacier-fed streams (i.e., < 2 µg L −1 ; Rinke et al., 2001; Robinson et al., 2016), depending on the underlying geology, dissolved and particulate P can be quite high. For example, Patagonian glacier-fed streams, with basaltic, granitic, pyritisic, and silicic metamorphic bedrock, have up to ∼40 µg of dissolved P L−<sup>1</sup> (Martyniuk et al., 2014; Miserendino et al., 2018). It is not clear, however, how bioavailable particulate P released from bedrock weathering is though very small grain sizes of P can be readily consumed by microbial life (Smith et al., 1977; Brahney et al., 2015a). Generally, high concentrations of DOC and N but limited P drive correspondingly high ratios of C:P and N:P in glacier-fed streams (Rinke et al., 2001; Robinson et al., 2016; **Figure 4**).

#### Trophic Interactions

Glacier-fed stream biofilms are generally limited by both N and P (**Figure 4**), but they exhibit stoichiometric flexibility and are capable of producing biomass with high C:P and N:P, which can in turn drive P limitation of consumers (e.g., macroinvertebrates, fish; **Figure 4**). Experimental addition of P increases microbial richness in glacier-fed streams (Kohler et al., 2016). Combined N and P additions stimulate biofilm primary production and shift microbial community composition (Robinson et al., 2003; Kohler et al., 2016; **Figure 4**). However, hydrological characteristics (e.g., flow and proportion of glacier influence) may override the influence of nutrient additions in some areas, such as the Swiss Alps (Rinke et al., 2001) and the Tian Shan Mountains of central Asia (Ren et al., 2017b). Other physical factors (e.g., suspended solids) can enhance C fixation by algae, thereby increasing C:P and lowering food quality (Martyniuk et al., 2014). In lightlimited streams with high suspended solids, this reduction in food quality (e.g., higher biofilm C:nutrient) may consequently lead to elemental imbalances between primary producers and invertebrate consumers (Martyniuk et al., 2019).

#### Glacier-Fed Lakes Nutrients

In glacierized regions, mountain lakes are typically fed either directly by glacier meltwater or indirectly by glacier-fed streams. Either way, glaciers drive lake ecosystem processes through the

blue), as well as water from mountain glacier-fed streams and lakes. No data were available for phosphate in subglacial systems. Numbers in red along the x-axis indicate sample sizes for each group. The data for this figure, including details of the studies represented, are provided in Table S1.

regulation of nutrients, temperature, and contaminant inputs (Saros et al., 2010; Slemmons and Saros, 2012; Rose et al., 2014). However, unlike nearby streams, mountain lakes are more likely to be fed by multiple hydrological sources than streams. To this end, an important distinction must be made between meltwater from glaciers and perennial snowfields. Glacierderived meltwater stems from a permanent, moving body of ice (i.e., the glacier) and is often rich with suspended sediments and may reflect longer term accumulation of nutrients (e.g., N from atmospheric deposition). Perennial snowfields encompass any other permanent bodies of snow. In general, glacier-derived meltwater exhibits higher nitrate concentrations than snowmelt (Hodson, 2006; Wynn et al., 2007; Saros et al., 2010), which in turn leads to more nitrate in primarily glacier-fed vs. snow-fed lakes (Saros et al., 2010; Slemmons and Saros, 2012; Williams et al., 2016; Warner et al., 2017). For example, Williams et al.

(2016) reported average nitrate concentrations of 13 µg L−<sup>1</sup> in glacier-fed lakes vs. 5 µg L−<sup>1</sup> in snow-fed lakes of the northern Rocky Mountains, USA. In the southern Rocky Mountains, nitrate concentrations are remarkably higher in glacier-fed lakes (> 50 µg L−<sup>1</sup> ) compared to nearby snow-fed lakes (< 15 µg L−<sup>1</sup> ; Saros et al., 2010; Slemmons and Saros, 2012). In western North America, rock and debris-covered glaciers are an additional, generally overlooked, resource subsidy to headwater aquatic ecosystems. In addition to providing another source of meltwater, rock glacier outflows are substantially higher in total N and other solutes vs. surface glaciers (Fegel et al., 2016).

In the absence of significant anthropogenic N inputs (pollution and atmospheric deposition) and meltwater input from glaciers, mountain lakes are typically low in N (Bergstrom and Jansson, 2006; Elser et al., 2009). Phytoplankton in mountain lakes are thus generally N-limited and increased N input can stimulate their growth and easily shift them to P-limitation (Elser et al., 2009, 2010). Meltwater from glaciers and snowfields provide additional N to mountain lakes, alleviating N-limitation and escalating P-limitation of phytoplankton. In addition to providing a source of N (Saros et al., 2010; Slemmons and Saros, 2012; Williams et al., 2016), glacier meltwater carries large amounts of suspended P-rich silt and clay, which is mainly generated by comminution (i.e., grinding) of the underlying bedrock, and supplies important, albeit less bioavailable, P to glacier-fed lakes (Hodson et al., 2004, 2005). In practice, nutrient limitation experiments have demonstrated that phytoplankton communities are limited by P in glacier-fed lakes while colimitation of P and N occurs in snow-fed lakes (Slemmons and Saros, 2012). While total phytoplankton biomass is not consistently higher in glacier- vs. snow-fed lakes (Saros et al., 2010; Slemmons and Saros, 2012), phytoplankton communities in glacier-fed lakes tend to be more similar to one another than their snow-fed counterparts (Warner et al., 2017). At highertrophic levels, heterotrophic production appears primarily limited by P in glacier-fed lakes (Mindl et al., 2007).

Suspended sediments in glacial meltwater affect the optical properties of the lake water column by intensifying light attenuation via absorption and reflectance (Gallegos et al., 2008; Laspoumaderes et al., 2013; Sommaruga, 2015), inducing a shallower photosynthetic zone, reduced photosynthetically active radiation (PAR), and ultimately reduced photosynthetic rates (Modenutti et al., 2000; Rose et al., 2014; **Figure 5**). According to the light:nutrient hypothesis, the variation of PAR intensity and phosphorus availability (light:P) determine phytoplankton stoichiometry (Sterner et al., 1997; Huisman et al., 2004). Under high light:P, phytoplankton are severely P-limited and exhibit high biomass C:P, while low light:P results in low biomass C:P. The light:nutrient hypothesis is supported in glacier-fed lakes. Phytoplankton tend to have low C:P in turbid water (higher contribution of glacier meltwater) while high C:P ratios are associated with clearer water (i.e., less contribution of glacier meltwater; Laspoumaderes et al., 2013, 2017; **Figure 5**).

#### Trophic Interactions

Both nutrient availability and physical characteristics of glacier-fed lakes drive trophic interactions (**Figure 5**). Nutrient availability ultimately governs phytoplankton C:P (an indicator of food quality for zooplankton) and zooplankton taxonomic groups differ significantly in their energetic P-demand which results in decoupled community responses to altered food quality (Elser et al., 1996; Elser and Urabe, 1999; Sterner and Elser, 2002; Urabe et al., 2002). For example, Daphnia require relatively high P food sources and exhibit consistent, low C:P body content. In contrast, some copepods exhibit more variable, but generally high C:P ratios. Thus, with declining food quality (high C:P), P-rich species are less competitive than P-poor species (Sterner and Elser, 2002). Consequently, as glaciers recede, increasing phytoplankton C:P could alter grazer communities by increasing the abundance of P-poor copepods (e.g., Boeckella gracillipes, Sommaruga, 2015; Laspoumaderes et al., 2017) but decreasing P-rich Daphnia (e.g., Daphia commutate; Laspoumaderes et al., 2013, 2017). Temperature also affects food quality and may intensify P limitation in P-rich consumers of glacier-fed lakes (Laspoumaderes et al., 2013). According to the stoichiometric growth rate hypothesis, high temperature accelerates consumer growth rates, leading to elevated demand for P in ribosomal RNA to support rapid growth (Main et al., 1997; Acharya et al., 2004). However, the growth rate hypothesis has not been tested in glacier-fed lakes. Although an important component of lake ecosystems (Vadeboncoeur et al., 2002), sediment stoichiometry and benthic communities are also largely overlooked (but see Lepori and Robin, 2014; Oleksy, 2019). The framework of ecological stoichiometry provides a means for testing these hypotheses to develop new understanding of how nutrient availability, elemental ratios, and light availability alter ecosystem structure and function in glacier-fed lakes (e.g., the predictions shown in **Figure 5**).

# CLIMATE CHANGE IMPLICATIONS

Climate change is dramatically altering mountain landscapes. At high elevations, atmospheric temperatures are rising up to three times more quickly than the global average (Nogués-Bravo et al., 2007). Rapid warming is driving the most obvious physical change in mountain ecosystems worldwide: the ongoing recession of glaciers and perennial snowfields (e.g., Zemp et al., 2015). Considerable attention has been devoted to understanding the effects of cryosphere decline on the ecology of mountain ecosystems, and this work has been summarized in recent syntheses and reviews, focused on streams (Hotaling et al., 2017b), lakes (Moser et al., 2019), and terrestrial plants (Vitasse et al., 2018). Generally, the story is clear; climate-induced recession of the mountain cryosphere will alter hydrological regimes, nutrient fluxes, and aquatic biogeochemistry (Huss et al., 2017) and induce global biodiversity loss across taxonomic scales from genetic diversity (e.g., Finn et al., 2013; Jordan et al., 2016) to species (e.g., Giersch et al., 2017) and communities (e.g., Fell et al., 2018; Hotaling et al., 2019a). Efforts to predict biodiversity loss in response to a fading cryosphere, while valuable, represent an end-member focus (e.g., the loss of a species) which overlooks an intermediary discussion of how factors that underlie biotic responses (e.g., nutrient limitations) may be affected.

As we seek to predict how the ecological stoichiometry of cryosphere-influenced headwater lakes and streams will change, additional factors must be considered. First, even though atmospheric temperatures are rising more quickly at high elevations than almost anywhere on Earth (Nogués-Bravo et al., 2007), mountain lakes may actually be warming more slowly than those in low elevations (Christianson et al., 2019). The same is likely also true for headwater streams. This lag in aquatic temperature change as they relate to atmospheric temperatures is likely due to buffering from dwindling glaciers and perennial snowfields (Zhang et al., 2014; Zemp et al., 2015). However, this buffering effect may soon be lost, and in combination with earlier melting of snowpack and longer ice-free stretches in headwaters (Preston et al., 2016), increased productivity in glacier-influenced lakes and streams is almost sure to follow. Second, nutrient deposition legacies and spatial variation in deposition rates matter. Globally, deposition rates of N, P, and dust are geographically variable (Mahowald et al., 2008) and are known to alter aquatic stoichiometry even in remote locations (Mladenov et al., 2011; Brahney et al., 2015b). For instance, historical N deposition leads to P limitation in lakes (Elser et al., 2009), thus N deposition legacies must be considered when trying to anticipate how N:P stoichiometry might change.

#### Nutrients

Warming temperatures, and altered precipitation regimes (e.g., more winter precipitation falling as rain rather than snow, Knowles et al., 2006), will directly influence the supraglacial zone and lead to declines in glacier volume, shorter periods of seasonal snow cover, and increased outflows in the near-term that will dwindle as ice sources fade (Huss and Hock, 2018). A longer "growing season" for microbial communities in the supraglacial zone paired with rising atmospheric CO2, which will stimulate snow algae growth (Hamilton and Havig, 2018; **Figure 1**), may yield a concurrent increase in autochthonous organic C production by primary producers (e.g., algae and cyanobacteria). Such a rise in primary production is likely to escalate existing glacial-melt feedback loop by decreasing surface ice albedo (e.g., Ganey et al., 2017). While N deposition rates are decreasing or stabilizing in some regions (Engardt et al., 2017; Yu et al., 2019), N emissions are rising globally and resulting deposition will increase N subsidies in remote mountainous areas with reactive N (Holtgrieve et al., 2011; Battye et al., 2017; Milner et al., 2017). In general, DOC and N concentrations will increase in downstream environments as glaciers recede due to ice-locked pools of both nutrients being released while P concentrations will decline as the loss of ice mass reduces the erosive power of the glacier (Hood and Scott, 2008; Hood and Berner, 2009; Hood et al., 2015; **Figure 4**).

In the short term, increased N and DOC inputs from glacial meltwaters may alleviate N limitation of primary producers and heterotrophic microorganisms in biofilms, increasing gross primary production and ecosystem respiration in glacier-fed streams (Uehlinger et al., 2010; Cauvy-Fraunié et al., 2016; Kohler et al., 2016; **Figure 4**). However, while DOC in glacier-fed streams may initially increase due to accelerated glacier melting, the lability (i.e., potential rate of turnover) will substantially decrease as the proportion of ice-locked DOC declines relative to less labile terrestrially derived DOC (Singer et al., 2012; Milner et al., 2017; Hemingway et al., 2019). Glacier-fed lakes and streams should continue to be P-limited in the near term (**Figures 4**, **5**); however, as glaciers completely recede, downstream environments may become co-limited by both N and P (Saros et al., 2010; Slemmons and Saros, 2012; **Figures 4**, **5**). Ultimately, it is clear that future climate conditions will favor increased in situ C fixation in the mountain cryosphere. Assuming no corresponding increases in N or P, via atmospheric deposition or otherwise, supraglacial communities will become increasingly N- and P-limited as climate change proceeds. However, heterogeneity in N and P deposition rates make general predictions difficult.

Changes in downstream ecosystems induced by atmospheric warming will not be solely dependent on upstream deglaciation. Indeed, warmer temperatures will increase chemical weathering of bedrock minerals (e.g., calcite and apatite) which will subsequently alter water chemistry, particularly via increased N and P (Heath and Baron, 2014; Price et al., 2017). Rock glaciers, which are likely to be less affected by atmospheric warming and will therefore persist on the landscape longer than their surface counterparts (Knight et al., 2019), release solute-rich outflows into headwaters and labile C that can fuel heterotrophic production (Fegel et al., 2016, 2019). Thus, the ratio of surface:rock glacier meltwater in headwaters is likely to decline with solute concentrations in available water shifting toward levels typical of rock glaciers occurring simultaneously.

As climate change proceeds, treeline dynamics will also alter ecological stoichiometry in mountain headwaters (Martyniuk et al., 2016). Generally speaking, lakes and streams surrounded by terrestrial vegetation have higher water column C:P than those above treeline (Stenzel et al., 2017). In addition to changing quality and quantity of C melting from glaciers, lakes above treeline are typically net autotrophic, however, as treelines advance, they will become more influenced by terrestrial allochthonous inputs (Rose et al., 2015). This will fundamentally alter how C is cycled, since presently DOM is mainly of autochthonous origin, which drives a greater proportion of highly labile, easily respired organic C (Kortelainen et al., 2013). Shifts in DOC quality coupled with climate-driven warming may therefore alter algal-bacterial interactions, which will ultimately cascade through the trophic food web and shift community compositions (González-Olalla et al., 2018). Thus, as treeline encroachment occurs, C and P inputs will increase for lakes and streams (Kopàcek et al., 2011 ˇ ), however, a corresponding increase in N uptake potential will occur in the surrounding landscape which will influence the amount of N export.

#### Trophic Interactions

Changing nutrient dynamics in the supraglacial zone under climate warming have the potential to alter linked habitats through changes in the amount and quality of nutrients being exported. Through the biology-albedo feedback loop described earlier, an increase in microbial growth may induce both elevated C fixation on the glacier surface and greater export of nutrients to the subglacial zone as well as glacier-fed streams and lakes. However, temporal dynamics of microbial communities on glaciers and snowfields—and particularly the interplay between primary producers (e.g., cyanobacteria) and consumers (e.g., fungi, invertebrates)—are largely unknown. Evidence from cryoconite hole bacterial communities of the Italian Alps suggests that early season melting gives rise to a wave of primary productivity (e.g., cyanobacteria) which seeds a later season heterotrophic community (Pittino et al., 2018). The ratio of producers to consumers throughout the melt season, and specifically that of resource production and consumption, paired with potential albedo variation among groups (e.g., if snow algae reduce albedo to a greater degree than fungi), are important avenues for future study. This is particularly true as we seek to move from observational studies of supraglacial biological activity to a more predictive framework that also includes the ecological stoichiometry of adjacent downstream habitats.

With limited knowledge of biological and chemical diversity of subglacial habitats on global scales, including the degree to which they are nutrient-limited in space and time, it is difficult to make stoichiometric predictions of their future. As discussed above, climate change will drive an increase in meltwater flowing from the surface to the subglacial zone with a concomitant increase in nutrients (particularly organic C) in tow. However, given that subglacial communities appear N-limited (Boyd et al., 2011), an increase in available C may not dramatically alter existing stoichiometric ratios in these habitats. Furthermore, in the shortterm, glacial sediment loads will increase, leading to increases in algal biomass and trends toward increased C:P ratios and thus lower food quality for higher trophic levels (Martyniuk et al., 2014). From a hydrological perspective, extended melt seasons at the glacier surface may translate to more stable, and open, flow paths within the glacier, perhaps including a reduction in anoxic conditions at the glacier-bedrock interface and elevated potential for direct nutrient transport. Because so little is known of trophic ecology at higher levels (e.g., ice worms, birds) on mountain glaciers, it is also difficult to predict how their interactions will be altered in the future. However, for high-elevation birds (e.g., Gray-crowned Rosy finches) which appear to derive substantial food resources from ice worms (S.H., pers. obs.), loss of glacier habitat may induce a major diet shift across portions of their range, fundamentally altering their present-day stoichiometry.

# CONCLUSIONS AND FUTURE DIRECTIONS

Although frozen environments place severe constraints on life, the mountain cryosphere harbors complex, dynamic biological communities which often experience some degree of nutrient limitation. In this review, we synthesized the ecological stoichiometry of three major components of the mountain cryosphere—supraglacial, subglacial, and adjacent downstream habitats—which vary in their nutrient inputs and availability. These habitats have received varying levels of research attention to date.

#### Habitats

Supraglacial environments are generally nutrient-limited however the vast majority of our understanding of these habitats stems from continental ice sheets (e.g., Greenland, Antarctica), which differ from mountain glaciers in many ways, including geomorphology, nutrient availability (**Figure 2**), and biological diversity (e.g., presence of higher trophic levels). Chemical and biological differences between mountain glaciers and continental ice sheets may be due in large part to differing proximity to terrestrial and anthropogenic nutrient sources (e.g., N deposition). Moreover, chemical contrasts between mountain glaciers and continental ice sheets suggest that the former cannot simply serve as an analog for the latter, and that mountain glaciers must be treated as unique components of the global cryosphere. However, given the spatial variation of samples included in our analysis, we consider our conclusions to be intriguing but largely preliminary as unforeseen biases (e.g., in sampling locations) may exist in the data being compared.

Below mountain glaciers in the subglacial zone, much less is known of ecological stoichiometry beyond a handful of studies (e.g., Boyd et al., 2011) which point to high water C:N driving N limitation in sediment communities. However, due to sampling challenges (i.e., extremely difficult access), the nature of biogeochemical cycling, the general rules governing subglacial life, and the nutrient limitations on subglacial food webs are largely unknown. Similarly, the englacial environment—i.e., the portion of the glacier below the supraglacial but above the subglacial—remains a black box with presumably little to no biodiversity present (Hotaling et al., 2017a). However, future efforts to generate ice cores in mountain glaciers from the supraglacial zone to bedrock would provide an excellent opportunity to explore patterns of diversity, nutrient availability, and stoichiometry along a spatial continuum in these enigmatic habitats. This type of effort may also provide new insights on a challenge that has long puzzled glacier biogeochemists: why is N export from glaciers high even in regions with relatively low rates of N deposition (Slemmons et al., 2013)?

Like supraglacial environments, the ecological stoichiometry of glacier-fed lakes and streams has received relatively intense study. As the most influential factor mediating biological processes in these habitats, it is perhaps unsurprising that glacier coverage largely drives nutrient availability, elemental ratios, and light availability. However, the biological and biogeochemical connections between upstream (supraglacial and subglacial) and downstream environments remain unclear, with an overwhelming stoichiometric focus on the effects of N loads and turbidity. To this end, future efforts should focus on the relative ratios of N:P, overall ecosystem function, and temperature-nutrient interactions, as these avenues are likely to be particularly important to understanding algal community structure and productivity (Oleksy, 2019).

# Looking Ahead

Glaciers and perennial snowfields will continue to disappear in the decades to come (Zemp et al., 2015), particularly at low latitudes (Hall and Fagre, 2003). As the mountain cryosphere fades, the entire biome they support including most, if not all, of the existing genetic, species, and stoichiometric diversity will also be lost. However, while some habitats will be lost entirely (e.g., supraglacial and subglacial), others will persist in fundamentally altered form (e.g., high-elevation lakes). In the short-term, glacier recession will generally lead to elevated levels of C and N, and reduced P, in headwaters. Eventually, however, downstream environments will likely transition to co-limitation by N and P as the large store of N currently present in mountain glaciers is depleted.

In this review, we identified gaps in stoichiometric knowledge of the mountain cryosphere, as well as what appear to be key, systemic differences in nutrient concentrations between mountain glaciers and continental ice sheets (**Figure 2**). These differences are important because they highlight the need for targeted studies of mountain glacier ecosystems since for the most part, they indicate that patterns of ecological stoichiometry in mountains cannot be inferred from existing ice sheet-focused research efforts. Furthermore, given the potential for substantial variation in nutrient concentrations in mountain cryosphere-influenced habitats (**Figure 3**), more comparative studies that explicitly integrate spatial and temporal sampling will shed light on this variation and refine understanding of the mechanism(s) underlying it. Multi-trophic stoichiometric perspectives are also needed. For instance, in North America, it is currently unknown how nutrient dynamics that shape elemental ratios at the microbial scale translate to heterotrophic fungi, macroinvertebrate consumers (e.g., ice worms), and ultimately birds. Finally, refined understanding of ecological stoichiometry in the mountain cryosphere can also improve general understanding of stoichiometric principles. Indeed, with nutrient limitations that vary in space and time paired with dramatic differences in turbidity and light availability across aquatic habitats, the mountain cryosphere and the habitats it directly influences provide a natural laboratory for testing fundamental stoichiometric hypotheses at environmental extremes.

#### REFERENCES


# AUTHOR CONTRIBUTIONS

ZR conceived of the review. ZR, NM, IO, AS, and SH wrote the manuscript. All authors read and approved the final version.

#### ACKNOWLEDGMENTS

This manuscript stemmed from the fourth Woodstoich meeting at Flathead Lake Biological Station (FLBS) and was supported by the National Science Foundation (#DEB-1840408). We thank FLBS staff for their logistical support. Trinity Hamilton, Jen Schweitzer, and three reviewers provided feedback that improved the manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00360/full#supplementary-material


dissolved organic matter in clear alpine lakes. Nat. Commun. 2:405. doi: 10.1038/ncomms1411


physicochemical conditions in glacier-fed streams. Freshw. Sci. 36, 479–490. doi: 10.1086/693133


affects ecosystem structure and process. Am. Nat. 150, 663–684. doi: 10.1086/ 286088


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer, EH, declared a shared affiliation, though no collaboration, with one of the authors, IO, to the handling Editor.

Copyright © 2019 Ren, Martyniuk, Oleksy, Swain and Hotaling. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Constraining Carbon and Nutrient Flows in Soil With Ecological Stoichiometry

Robert W. Buchkowski <sup>1</sup> \*, Alanna N. Shaw<sup>2</sup> , Debjani Sihi <sup>3</sup> , Gabriel R. Smith<sup>4</sup> and Ashley D. Keiser 5,6

*<sup>1</sup> Yale School of Forestry & Environmental Studies, New Haven, CT, United States, <sup>2</sup> Department of Ecosystem and Conservation Sciences, W.A. Franke College of Forestry & Conservation, University of Montana, Missoula, MT, United States, <sup>3</sup> Environmental Sciences Division and Climate Change Science Institute, Oak Ridge National Laboratory, Oak Ridge, TN, United States, <sup>4</sup> Department of Biology, School of Humanities and Sciences, Stanford University, Stanford, CA, United States, <sup>5</sup> Stockbridge School of Agriculture, University of Massachusetts Amherst, Amherst, MA, United States, <sup>6</sup> Department of Ecology, Evolution and Behavior, University of Minnesota, St. Paul, MN, United States*

#### Edited by:

*James Joseph Elser, University of Montana, United States*

#### Reviewed by:

*Joshua Schimel, University of California, Santa Barbara, United States Chengrong Chen, Griffith University, Australia Michael Kaspari, University of Oklahoma, United States*

> \*Correspondence: *Robert W. Buchkowski robbuchko@gmail.com*

#### Specialty section:

*This article was submitted to Behavioral and Evolutionary Ecology, a section of the journal Frontiers in Ecology and Evolution*

> Received: *22 July 2019* Accepted: *24 September 2019* Published: *10 October 2019*

#### Citation:

*Buchkowski RW, Shaw AN, Sihi D, Smith GR and Keiser AD (2019) Constraining Carbon and Nutrient Flows in Soil With Ecological Stoichiometry. Front. Ecol. Evol. 7:382. doi: 10.3389/fevo.2019.00382* Soil organic matter (SOM) is central to soil carbon (C) storage and terrestrial nutrient cycling. New data have upended the traditional model of stabilization, which held that stable SOM was mostly made of undecomposed plant molecules. We now know that microbial by-products and dead cells comprise unexpectedly large amounts of stable SOM because they can become attached to mineral surfaces or physically protected within soil aggregates. SOM models have been built to incorporate the microbial to mineral stabilization of organic matter, but now face a new challenge of accurately capturing microbial productivity and metabolism. Explicitly representing stoichiometry, the relative nutrient requirements for growth and maintenance of organisms, could provide a way forward. Stoichiometry limits SOM formation and turnover in nature, but important nutrients like nitrogen (N), phosphorus (P), and sulfur (S) are often missing from the new generation of SOM models. In this synthesis, we seek to facilitate the addition of these nutrients to SOM models by (1) reviewing the stoichiometric bias the tendency to favor one element over another—of four key processes in the new framework of SOM cycling and (2) applying this knowledge to build a stoichiometrically explicit budget of C, N, P, and S flow through the major SOM pools. By quantifying the role of stoichiometry in SOM cycling, we discover that constraining the C:N:P:S ratio of microorganisms and SOM to specific values reduces uncertainty in C and nutrient flow as effectively as using microbial C use efficiency (CUE) parameters. We find that the value of additional constraints on stoichiometry vs. CUE varies across ecosystems, depending on how precise the available data are for that ecosystem and which biogeochemical pathways are present. Moreover, because CUE summarizes many different processes, stoichiometric measurements of key soil pools are likely to be more robust when extrapolated from soil incubations to plot or biome scale estimates. Our results suggest that measuring SOM stoichiometry should be a priority for future empirical work and that the inclusion of new nutrients in SOM models may be an effective way to improve precision.

Keywords: soil organic matter, ecological stoichiometry, carbon, nutrients, litter decomposition, soil organisms, stabilization, mineral associated organic matter

# INTRODUCTION

Why does some soil organic matter (SOM) persist for hundreds of years while other SOM turns over quickly, rapidly losing its carbon (C), and nutrients? Throughout the twentieth century, most scientists believed that stable SOM was composed primarily of plant compounds that persisted in soil because their complex chemical structures resisted microbial degradation (Brady and Weil, 2007). Now, new evidence indicates that more than 50% of stable SOM may instead be made of chemically simple compounds incorporated into microbial residue and dead microbial cells (i.e., necromass; Simpson et al., 2007; Ludwig et al., 2015). Microbial residues and necromass can become attached to mineral surfaces or trapped within soil aggregates, rendering them inaccessible to decomposition (Prescott, 2010; Bradford et al., 2013; Cotrufo et al., 2013; Kallenbach et al., 2015; Lehmann and Kleber, 2015; Jilling et al., 2018). These new findings represent a major shift in our understanding of SOM formation, and have spurred the development of a new generation of SOM models. These models outperform those without microbial and mineral pathways (Wieder et al., 2013; Robertson et al., 2019).

However, to represent the influence of microbial productivity and mineral sorption accurately, models need to include not only a new set of pathways, but also the factors that control those pathways. One of the most important determinants of microbial productivity is the degree to which the relative elemental composition, or stoichiometry, of microbial biomass matches the substrate on which it grows (Sterner and Elser, 2002; Schimel and Weintraub, 2003; Cleveland and Liptzin, 2007; Buchkowski et al., 2015). A closer match enables greater thermodynamic efficiency and reduces the potential for nutrient limitation (Sterner and Elser, 2002; Schimel and Weintraub, 2003). The stoichiometries of plant inputs, SOM, and microbial biomass are known to affect C and nutrient flows in soil (Griffiths et al., 2012; Tipping et al., 2016) but remain absent from many SOM models (Allison et al., 2010; Wieder et al., 2014).

SOM models developed in the past decade often contain no more than two elements (Moorhead et al., 2012; Abramoff et al., 2017; Sulman et al., 2019) because analytical tractability decreases sharply with the inclusion of each additional elemental cycle. Instead, new models designed to assess the microbial components of SOM often rely on substrate use efficiency parameters, which succinctly represent the proportion of nutrient uptake that can be converted into organism biomass (Allison et al., 2010; Manzoni et al., 2012). In fact, substrate use efficiency is a summary of a number of factors such as resource chemical quality, stoichiometry, and climatic conditions. Modeling work has repeatedly demonstrated that adding the underlying mechanisms behind substrate use efficiency can improve our predictions about SOM cycling (Schimel and Weintraub, 2003; Wieder et al., 2013; Abramoff et al., 2017; Sulman et al., 2019).

To facilitate the development of stoichiometrically comprehensive SOM models capable of better prediction and accuracy, we (1) qualitatively review any changes in stoichiometry that occur during the flow of organic matter into soil, and (2) use this information to build a stoichiometrically explicit SOM budget, which we refer to as the Linear SOM Description (LSD). The LSD is a static representation of the SOM cycle, including microbial and mineral pathways as well as flows of C, nitrogen (N), phosphorus (P), and sulfur (S). Rather than assessing changes in flow rates over time, as in dynamic SOM models like CENTURY (Parton et al., 1987), LSD allows us to calculate the possible range of possible flow rates in the system (van Oevelen et al., 2010). Just as a set of pool measurements can be used to build a whole-ecosystem C and nutrient budget, LSD shows us how data on substrate use efficiency and the stoichiometry of different pools changes the range of flows that mass balance can allow. Using this approach, we can quantify how effectively new information reduces uncertainty (i.e., the range between minimum and maximum possible flows) in the system.

The LSD can be used to efficiently analyze many model structures and parameter sets. When selecting model structures, we considered recent conceptual advances about SOM in the context of stoichiometry, so we could determine when stoichiometry can be an effective tool for incorporating field data into numerical models (sensu Blankinship et al., 2018). We hypothesized that the specific stoichiometric constraints on microbes and mineral-associated organic matter (MAOM) would reduce model uncertainty as effectively as carbon use efficiency (CUE), especially when we included the stoichiometry of many elements (i.e., C, N, P, and S).

We also parameterized versions of the LSD for 27 terrestrial ecosystem × biome combinations with open-source data (Bond-Lamberty and Thomson, 2014; Iversen et al., 2017; Qiao et al., 2019). Our goal was to ground the analysis of model structures and parameters by using case studies from real ecosystems. We hypothesized that the difference in the information available for these systems combined with the difference in structure across them would reverse whether stoichiometric or substrate use efficiency measurements were better at constraining flows. By addressing these hypotheses, we (1) demonstrate the value of stoichiometric data for improving our understanding of SOM cycling, (2) provide upper and lower bounds for feasible element flows across different ecosystem types, which can be used to evaluate the performance of predictive SOM models, and (3) identify high-value targets for future theoretical and empirical research efforts.

#### STOICHIOMETRY AND SOM FORMATION

The structure of the LSD reflects our current understanding of SOM cycling. Following the new generation of SOM models, we focused on microbial productivity and mineral sorption as primary pathways of SOM stabilization. We also included other processes, like physical protection, that are important contributors to SOM formation, transformation, and turnover (Six et al., 2002).

In this section, we provide a qualitative overview of stoichiometry in four important SOM pathways represented in the LSD: litter losses through abiotic processes and soil animals, belowground inputs of organic matter, microbial productivity and CUE, and the stabilization of MAOM. Other processes important for SOM formation, transformation, and loss have been reviewed elsewhere (e.g., stoichiometry and leaf litter decomposition by soil microorganisms; Manzoni et al., 2010). **Figure 1** provides a simplified visualization of the pools and flows reviewed here (see **Figure S1** for a diagram of the entire system in the LSD). In each subsection, we identify both how a particular process is dependent on stoichiometry and how that process itself can alter stoichiometry.

## Aboveground Litter Fragmentation and Decomposition

Plants contribute large pulses of C, N, P, and S into soil ecosystems through litterfall (Schlesinger and Bernhardt, 2013). Leaf litter can be incorporated into the soil through abiotic and biotic pathways, each of which has the potential to alter the stoichiometry of SOM inputs (**Figure 1**). Here, we focus on abiotic pathways of litter losses and processing of litter by soil animals.

#### Abiotic Litter Loss

Litter losses through abiotic processes can change stoichiometry or be stoichiometrically neutral. For example, physical fragmentation caused by freeze-thaw or animal movement isn't associated with stoichiometric changes (Fahey et al., 2013). Because fragmented litter is more easily incorporated into the particulate organic matter (POM) pool, this process helps to transfer higher C:N:P:S compounds from litter into SOM. In contrast, photodegradation preferentially degrades lignin (Austin and Ballare, 2010; Brandt et al., 2010; Austin et al., 2016). Limited data suggest that changes in litter stoichiometry across the first year of decomposition are similar between photodegradation and degradation by microbial exoenzymes (Frouz et al., 2011; Wang et al., 2015). Photodegradation is not, however, independent of microbial decomposition; photodegraded litter is more accessible to microbial decomposers, which facilitates subsequent transformation and incorporation into microbial biomass or leaching into dissolved organic matter once it is broken into small enough molecules (DOM; Kaiser and Kalbitz, 2012; Cotrufo et al., 2015; Austin et al., 2016).

FIGURE 1 | A framework for the formation and loss of soil organic matter. Dashed orange arrows represent pathways with a stoichiometric bias, while solid blue arrows represent pathways with stoichiometry matching the source pools. This is not a model diagram and for readability does not show that most pools have variable chemical quality, multiple species, and multiple elements. The diagram also does not show inputs (root exudates, carbon to symbionts, nitrogen fixation, leaf, wood, and root litter) or outputs (leaching and respiration). See the model diagram in the supplemental (Figure S1; Table S1).

DOM lost from leaf litter can be incorporated into the MAOM pool as it moves into the soil profile. Once adsorbed, DOM can desorb, partially as a function of redox-sensitive elements, and undergo processing by microbes (Fahey et al., 2011; Kaiser and Kalbitz, 2012). As a result, DOM stoichiometry in upper soil profiles more closely matches the original litter, while DOM in lower soil profiles has a lower C:N because of microbial processing (Cotrufo et al., 2015). In some ecosystems, the DOM pathway can account for losses of up to a third of litter C (Cotrufo et al., 2015). Carbon and nutrients in litter that are not photodegraded or leached out are instead decomposed by soil animals and microorganisms.

#### Litter Consumption by Soil Animals

Through consumption and digestion, soil animals alter the stoichiometry of leaf litter entering the soil (Kautz et al., 2002; David, 2014; Frouz et al., 2015). However, their effects are extremely variable because soil animals have a wide range of stoichiometric requirements (Teuben and Verhoef, 1992; Martinson et al., 2008). For example, soil macrofauna, like isopods and millipedes, usually lower the C:N ratios of their food when converting them to feces (Kautz et al., 2002; Bastow, 2011; Frouz et al., 2015). However, soil mesofauna, like collembola, can raise the C:N ratio of leaf litter passing through their guts (Teuben and Verhoef, 1992). Data on the assimilation of P and S by soil animals are sparse, but studies on isopods and collembola indicate that gut passage tends to increase bulk C:P and C:S ratios, as well as the amount of P and S in dissolved forms (Morgan and Mitchell, 1987; Teuben and Verhoef, 1992). Some animals may reduce C: nutrient ratios because easily degradable C-rich compounds in the leaf litter are decomposed and assimilated in their guts (David, 2014). If this is true, understanding the impact of soil animals on SOM stoichiometry will require data on chemical quality as well as stoichiometry (David, 2014). These generalizations are based on limited data sets because the nutrient ratios of soil faunal feces, especially those of smaller-bodied organisms, are rarely measured (Osler and Sommerkorn, 2007). While the role of soil animal gut processing is partially dependent on the animal's nutrient requirements, it is also driven by the activity of litter-degrading microorganisms.

#### Belowground Inputs of Organic Matter

Our understanding of SOM inputs that originate belowground is poorly constrained compared to that of aboveground litter transformations (Sokol et al., 2019). Since a sizable portion of terrestrial ecosystem productivity is allocated belowground (Gill and Finzi, 2016), this represents a major knowledge gap hindering our ability to model SOM formation. Here, we discuss two major pathways by which belowground organic matter inputs can be incorporated into SOM: roots and mycorrhizal fungi.

#### Roots

A majority of the decomposition literature is devoted to leaf litter (Sokol and Bradford, 2019), despite the disproportionate importance of fine root turnover and the marked difference between foliar and root chemistry. Indeed, root-derived C may be more important to SOM formation than aboveground C inputs by a factor of 2.4–2.7 (Crow et al., 2009; Kong and Six, 2010; Clemmensen et al., 2013; Sokol et al., 2019). Roots also have a longer mean residence time (Rasse et al., 2005) and are more likely to be incorporated into stable SOM. Like annual inputs of leaf litter, fine root litter chemical quality strongly influences initial (<1 year) litter decomposition dynamics (Adair et al., 2008; Cotrufo et al., 2013; See et al., 2019).

In the first stage of decomposition, fine roots decompose slowly (See et al., 2019). This is a function of chemical recalcitrance, which summarizes not only stoichiometry but also the composition of chemical constituents in organic matter. Therefore, the impact of chemical quality and stoichiometry in roots matches trends for leaf litter (Aerts, 1997; Hobbie, 2005; Adair et al., 2008; Harmon et al., 2009; Keiser and Bradford, 2017). However, while foliar litter inputs tend to be nutrient poor, with a global average C:N:P of 3007:45:1 (McGroddy et al., 2004), fine roots (<2 mm) are relatively nutrient rich, with a global average C:N:P of 522:12:1 (Gordon and Jackson, 2000). Fine root litter maintains the C:N of fine roots, but has higher C:P because, unlike leaf litter, only modest P resorption has been identified in fine roots (Gordon and Jackson, 2000; Yuan et al., 2011).

The finest roots have extremely slow decomposition rates despite their low C:N ratios (∼21) and high turnover rates (Sun et al., 2018). The chemical structure of roots is important for determining decomposition rates, because microbes need to depolymerize these compounds to metabolize them. Organic matter with less favorable stoichiometry (e.g., high C:N ratio) and high glucose content, for example, is more easily degraded than organic matter with favorable stoichiometry (e.g., lower C:N ratio) and a high portion of lignin and suberin. Therefore, stoichiometry is one indicator of the chemical quality and decomposition of fine roots, but cannot be totally separated from differences in other dimensions of chemical quality (Sugai and Schimel, 1993). As in leaf litter, where higher proportions of structural (i.e., lignin) and inhibitory compounds decrease microbial decomposition and the assimilation of C and nutrients (Manzoni et al., 2010; Sinsabaugh et al., 2013), many first order roots contain low energy C compounds that limit their decomposition rates, and may contribute to C limitation in soil microbes (Mooshammer et al., 2014; Sun et al., 2018).

Root litter is an important belowground input, but several other classes of root-derived C also make large contributions to SOM (Crow et al., 2009; Kong and Six, 2010; Clemmensen et al., 2013; Sokol et al., 2019). These inputs include several classes of root exudates, including those that contribute to priming (Kuzyakov et al., 2000; Kuzyakov, 2010; Di Lonardo et al., 2017), desorption (Keiluweit et al., 2015), defense (Jung et al., 2012), and symbiont signaling/support (Badri and Vivanco, 2009; Jung et al., 2012).

Root exudation is an important component of belowground C inputs because it can promote the decomposition of SOM compounds that require more energy to degrade or don't have a stoichiometric ratio that matches microbial demand (Brzostek et al., 2013; Drake et al., 2013; Di Lonardo et al., 2017). Such root exudates may promote microbial N mining in some circumstances by providing the energy required to decompose N rich SOM (Dijkstra et al., 2013). By providing lower C:N inputs to mine SOM, root exudates are an indirect pathway that allows soil microorganisms to overcome C limitation (Drake et al., 2013). Exudation into rhizosphere soils occurs alongside another, more direct plant-to-microbe C exchange involving symbiotic microorganisms that support plant nutrient acquisition.

#### Symbionts

Mycorrhizal fungi and symbiotic N-fixing bacteria are soil microorganisms that provide nutrients to plants in exchange for photosynthetic C (Smith and Read, 2008; Menge et al., 2017). The activities of these specialized symbioses can alter soil stoichiometry. For example, by facilitating plant nutrient uptake from soil, mycorrhizal fungi could increase SOM C:nutrient ratios (Smith and Read, 2008). In contrast, by incorporating new N into the plant-soil system, N-fixers could decrease SOM C:N ratios.

Almost all plants rely on root-associated mycorrhizal fungi to obtain the soil nutrients needed for growth (Smith and Read, 2008). Most of these plants associate with one, or sometimes both, of the two most common types of mycorrhizal fungi: ectomycorrhizal and arbuscular mycorrhizal. Due to their geographic ubiquity and importance to plant nutrition, mycorrhizal fungi constitute a major pathway by which plant C enters the soil (Leake et al., 2004) and by which soil N, P, and S exit the soil, mediating changes in SOM stoichiometry (Allen and Shachar-Hill, 2009; Orwin et al., 2011; Rosling et al., 2016). Crucially, these two mycorrhizal types are functionally different, resulting in systematic variation in nutrient cycling between ecosystems where most plants associate with one or the other (Phillips et al., 2013).

In general, ectomycorrhizal fungi promote plant uptake of organic nutrients from soil more than arbuscular mycorrhizal fungi (Smith and Smith, 2011; Phillips et al., 2013; Hodge and Storer, 2014; Shah et al., 2016), which can allow ectomycorrhizal plants to increase SOM C:nutrient ratios (Orwin et al., 2011). Leaf litter from ectomycorrhizal plants also decomposes more slowly than leaf litter from arbuscular mycorrhizal plants (Cornelissen et al., 2001), reducing litter nutrient accessibility. These aboveand belowground traits interact, enabling ectomycorrhizal plants to thrive in environments where slow mineralization limits nutrient availability (Steidinger et al., 2019), and, potentially, to also reinforce these conditions (Lin et al., 2017; Fernandez et al., 2019; Smith and Wan, 2019). Together, these dynamics contribute to correlations between global ectomycorrhizal plant abundance and higher SOM C:N ratios (Averill et al., 2014).

Mycorrhizal fungi also form extensive hyphal networks in soil that can directly contribute to SOM (Leake et al., 2004). In a Mediterranean poplar plantation, C derived from the turnover of mycorrhizal hyphae exceeded that of fine roots and leaf litter, forming the majority of total SOM inputs (Godbold et al., 2006). In boreal forests, a large proportion of primary productivity is allocated belowground (Gill and Finzi, 2016), and much of the SOM in these environments is derived from roots and mycorrhizal fungal biomass (Clemmensen et al., 2013). Fungal biomass has a higher nutrient content than plant roots, so inputs from mycorrhizal mycelia are likely to lower SOM C:nutrient ratios (Zhang and Elser, 2017).

In contrast to mycorrhizal fungi, which help plants acquire nutrients already present in soil, symbiotic N fixers introduce new N into the system (Menge et al., 2017). The majority of N fixation is performed by plant-associated symbionts, so these symbionts likely have the largest impact on SOM stoichiometry (Cleveland et al., 1999). Most N fixed by plant symbionts enters SOM after first traveling through the host plant; returning to the soil as root exudates and litter deposited above- or belowground (Sulman et al., 2019). Some of the fixed N, however, can enter the soil system directly through death of the symbionts themselves. This flow of nutrients may impact the stoichiometry of inputs into SOM (Brophy and Heichel, 1989).

### Microbial Productivity and Carbon Use Efficiency

Microbial metabolism controls the speed at which organic matter is decomposed. Microbial CUE measures the fate of metabolites as the proportion of total C uptake that is allocated to biomass (Manzoni et al., 2012). Microbial CUE is also highly variable (Qiao et al., 2019), affected by the chemical composition of substrate (Sugai and Schimel, 1993; Frey et al., 2013), as well as by the match between resource stoichiometry and microbial stoichiometry (Sinsabaugh et al., 2013, 2016). Carbon use efficiency determines how efficiently substrate inputs are transformed into stable, microbially-derived SOM. As such, even small changes in CUE can have important downstream impacts on SOM formation (Six et al., 2006). Because CUE succinctly captures the relationship between organic matter inputs and the SOM precursors produced by microorganisms, it appears as an important parameter in the new generation of microbe-focused SOM models (Allison et al., 2010; Wieder et al., 2014; Sihi et al., 2016).

### Stabilization of Mineral-Associated Organic Matter

Soil organic matter can be stabilized by association with mineral surfaces. The strength and nature of organo-mineral interactions varies with edaphic conditions like soil mineralogy and texture (Dungait et al., 2012; Lehmann and Kleber, 2015; Kallenbach et al., 2016; Jilling et al., 2018). Soil pH also plays an important role in SOM stabilization. Organo-metal complexes are associated with aluminum- and iron-dominated clay minerals in acidic systems and with calcium in alkaline systems (Rasmussen et al., 2018). The MAOM pool can persist over long time scales when physically protected within soil aggregates (Tipping et al., 2016), even if the organic compounds that are attached to mineral surfaces would otherwise be prone to microbial degradation. Predicting the stoichiometry of MAOM thus requires an understanding of which compounds are adsorbed onto mineral surfaces under what conditions.

The quantity, stability, chemistry, and distribution of MAOM in an ecosystem is determined by the relative dominance of biotic vs. abiotic pathways to sorption, whether DOM is adsorbed directly to mineral surfaces or after it is processed by soil microbes (Mikutta et al., 2019; Sokol et al., 2019). Because microbes tend to outcompete minerals for organic matter (Fischer et al., 2010), the biotic pathway increases in prevalence with microbial biomass such that sorption of microbiallyprocessed compounds dominates in surface soil horizons and rhizosphere soils (Sokol et al., 2019). In the rhizosphere, sugars and amino acids exuded by roots can release SOM from mineral associations, allowing for microbial processing and turnover of previously adsorbed compounds (Keiluweit et al., 2015). However, the high-quality substrates exuded into the rhizosphere also promote the production of large quantities of microbial necromass, which can be readily converted into stable MAOM and aggregates (Knicker, 2011; Cotrufo et al., 2013; Schrumpf et al., 2013; Craig et al., 2018). The balance of decomposition and formation of MAOM in the rhizosphere appears to produce a net increase in stable MAOM, likely because root exudates provide plenty of C and nutrients for sorption (Sokol et al., 2019).

The two classes of biomolecules most susceptible to sorption on mineral surfaces are proteins (peptide compounds) and inositol phosphates (Sollins et al., 2006; Celi and Barberis, 2007; Kleber et al., 2007; Tipping et al., 2016; Newcomb et al., 2017), both of which can introduce strong stoichiometric biases. The accumulation of N-rich proteins lowers MAOM C:N while the accumulation of inositol phosphates strongly (and variably) lowers MAOM C:P (Celi and Barberis, 2007; Tipping et al., 2016; Jilling et al., 2018). Inositol phosphates are the most common phosphate monoesters, making them the most abundant, highest-affinity organic P compounds in soils (Turner et al., 2002). Though they interact with exchange sites via their phosphate groups much like inorganic phosphates (Goldberg and Sposito, 1985), each inositol phosphate can have one to six phosphate groups. Many minerals, especially iron oxides, have a much higher affinity for these organic P compounds than their inorganic counterparts (Celi et al., 1999; Celi and Barberis, 2007). However, the extent of MAOM P enrichment resulting from the sorption of inositol phosphates (Adams et al., 2018) depends on the minerals in a given soil (Celi and Barberis, 2007). The flow of organic matter into the MAOM pool is thus subject to a variable stoichiometric filter, which generally concentrates nutrients relative to C and specifically concentrates N and/or P according to soil mineralogy.

#### THE LINEAR SOM DESCRIPTION

The LSD is a static description of the soil system (van Oevelen et al., 2010; Liang et al., 2011) with pools and processes included in recent dynamic SOM models as well as several new additions identified in our review (Wieder et al., 2014; Robertson et al., 2019; Sulman et al., 2019). It contains leaf/wood litter, root litter, soil animals, POM, DOM, saprotrophic microorganisms, symbiotic microorganisms, MAOM, and inorganic nutrients, as well as pools of POM and MAOM that are physically protected (**Figure 1**; **Figure S1**). We track the possible flows of four elements (C, N, P, and S) through these pools. We chose to include N, P, and S in addition to C because they are often in high demand relative to their supply in soil (Sterner and Elser, 2002; Brady and Weil, 2007). However, rather than tracking the mass of C, N, P, and S in these pools, we instead use pool stoichiometry to constrain the ranges of possible C and nutrient flows between them.

Most SOM models use a dynamic approach (e.g., Jenkinson et al., 1987; Parton et al., 1987; Wang et al., 2013; Sulman et al., 2014; Wieder et al., 2014; Abramoff et al., 2018), where pool sizes are explicitly modeled so that flows or fluxes between them can be quantified (Equation 1). We will use the term flow to describe the rate of carbon and nutrient movement between pools (van Oevelen et al., 2010), but note that the term flux is used interchangeably in the ecosystem literature (Chapin et al., 2011). In a linear dynamic model, flows (y; mass • time−<sup>1</sup> ) are calculated by multiplying a rate constant (k; time−<sup>1</sup> ) by the pool size (X; mass):

$$y^\* = kX \tag{1}$$

In contrast to dynamic linear models like CENTURY (Parton et al., 1987) or Roth-C (Jenkinson et al., 1987), the LSD is always at equilibrium, so neither rate constants (k) nor pool sizes (X) are explicitly modeled over time. The LSD does not include equations like Equation 1. Instead, equations setting mass balance, stoichiometry, and substrate use efficiency relationships are used to constrain flows (see Model Example; Equation 2). By constraining possible flows, the LSD indirectly determines the reasonable combinations of rate constants and equilibrium pool sizes (Equation 1). The flows in the LSD are comparable to those that would be reached at equilibrium in a linear model such as CENTURY (Parton et al., 1987) or a non-linear SOM model such as MIMICS or the Millennial model (Wieder et al., 2014; Abramoff et al., 2018).

Our approach has three major benefits. First, a linear system can effectively account for the interaction between flows that have stoichiometric constraints placed upon them (van Oevelen et al., 2010; Yang et al., 2017). Second, linear systems at equilibrium can be solved in high dimensions, allowing us to evaluate many structural and parameter changes (van Oevelen et al., 2010). Third, a linear system built with minimal assumptions yields the set of equilibria obtainable from structurally analogous nonlinear models (Stevens, 2009). An equilibrium approach does not simulate or predict how systems evolve over time, which means that the LSD does not replicate the efforts of new nonlinear models that include many of the same processes (Wieder et al., 2014; Abramoff et al., 2018; Sulman et al., 2019). Instead, it is useful for comparing equilibrium scenarios and evaluating what new data would provide the most effective constraints on equilibrium soil C and nutrient cycling.

#### Model Constraints

We imposed a few initial constraints to capture the biologically feasible range of state space in the absence of stoichiometric information. The baseline model conditions were:

1. System inputs, defined as leaf/wood litter, root litter, plant root exudates to DOM, and plant C to microbial symbionts, are constrained between 180 and 2,500 g C m−<sup>2</sup> year−<sup>1</sup> (Chapin et al., 2011).


Constraint 1 places the model output within a realistic range of elemental flows for ecosystems. Constraints 2 and 3 establish a reasonable model structure for soil elemental cycling including the important components discussed in our review. Finally, constraints 4–7 prevent infinite cycling between pools that exchange nutrients (e.g., microorganisms and inorganic N). These constraints are not necessary in models like CENTURY that mass balance pools, but are necessary in systems like the LSD, which mass balances flows (van Oevelen et al., 2010).

The pools with different chemical qualities defined in constraint 2 capture important mechanisms linking organisms to soil elemental cycling. They also help confirm that the importance of stoichiometry is not nullified by adding a categorical difference in chemical quality to our accounting. Many SOM models include some measure of quality, but two general strategies predominate: ontogeny models track the quality change of each unit of organic matter as it is processed, effectively treating quality as a continuous variable (Bosatta and Agren, 1991; Moore et al., 2004). Other models divide organic matter into discrete quality pools (Parton et al., 1998; Adair et al., 2008; Wieder et al., 2014). Because tractability is one of our primary goals, we chose to follow the latter approach. We include two pools with different levels of chemical recalcitrance, which are distinguished by their stoichiometry and the CUE with which they can be used (Moore et al., 2004; David, 2014). We recognize there are limitations in defining chemical recalcitrance purely on stoichiometry and differences in CUE, as opposed to a metric which also includes constituent compounds like lignin:N. Future expansions of the LSD could include many more biochemical categories, allowing it to track molecular differences in CUE that extend beyond stoichiometry (Yang et al., 2017).

Analytical models of SOM decomposition use different strategies to constrain the cycling of nutrients through microbial biomass. Many models combine the stoichiometric imbalance between organic matter, microorganisms, and their exoenzymes with flexible microbial CUE (Manzoni and Porporato, 2009; Moorhead et al., 2012, 2013; Manzoni, 2017) and overflow respiration (Schimel and Weintraub, 2003). Others represent stoichiometric biases in SOM decomposition by shifting microbial community composition or microbial community homeostasis (Sinsabaugh and Shah, 2012; Sinsabaugh et al., 2013, 2015; Waring et al., 2013; Warton et al., 2015; Hartman et al., 2017). Both of these modeling approaches demonstrate that stoichiometric constraints, soil community composition, and substrate use efficiencies are intertwined. In the LSD, we impose stoichiometry and substrate use efficiency as constraints that set a range of possible values (e.g., minimum to maximum of CUE or C:N ratio for bacteria). Since these constraints act in similar ways, we can identify the relative ability of either stoichiometry or substrate use efficiency to constrain the range of possible flows.

#### Model Example

In **Figure 2**, we illustrate how the LSD works by following C and N partitioning through ectomycorrhizal fungi. This example demonstrates how the flow changes as new constraints (i.e., substrate use efficiency or stoichiometric data) are added. It also highlights that the LSD does not currently incorporate empirical data on the flows between specific SOM pools or track the dynamics of SOM pool sizes, but instead uses C and nutrient input rates plus stoichiometry and substrate use efficiency to derive the range of possible flows.

Total C taken up by the ectomycorrhizal fungi (EYc) is equal to losses, since the model is mass balanced and static. Carbon lost from ectomycorrhizal fungi can exit the system through respiration (YEc), or can be partitioned into the MAOM pool (YMc + YJc), the POM pool (YPc + YQc), or the DOM pool (YDc + YUc), each with two levels of chemical quality:

$$EY\_C = YE\_C + YM\_C + YI\_C + YP\_C + YQ\_C + YD\_C + YU\_C \quad \text{(2)}$$

When applicable, CUE constrains the proportion of C uptake (EYc) that ectomycorrhizal fungi respire (YEc). Respiration is equal to 1-CUE multiplied by the total quantity of C taken up (EYc):

$$YE\_C = (1 - CUE) \times EY\_C \tag{3}$$

In the example given by **Figure 2A**, the quantity of C taken up by ectomycorrhizal fungi (EYc) is between 0 and 10 g C m−<sup>2</sup> y −1 . If neither CUE nor stoichiometric constraints are imposed, the C partitioned to respiration (YEc), DOM (YDc+YUc), MAOM (YMc +YJc), or POM (YPc + YQc) all have the same potential size range, which is also 0–10 g C m−<sup>2</sup> y −1 . However, all pathways must sum to inputs (Equation 2), which means that as one flow grows, all others must shrink.

If we define CUE as 0.5 (**Figure 2B**), the size of the respiration flow, YEc, becomes 0.5 <sup>∗</sup> EYc, which constrains it between 0 and 5 g C m−<sup>2</sup> y −1 . Because the flows are mass balanced (Equation 3), the other three flows must now be between 0 and 5 g C m−<sup>2</sup> y −1 . In this way, including CUE reduces the uncertainty, i.e., the range between minimum and maximum flows, in the system.

stoichiometric limitations, since the ECM pool must have ten units of C for every unit of N. The Lost N flow (i.e., transfer to plants) now has a maximum flow size of 7, because the system must remain mass balanced.

The LSD can also be constrained by imposing stoichiometric ratios, such as a fixed C:N ratio, as follows (van Oevelen et al., 2010):

$$\begin{aligned} EY\_C - YE\_C &= Y\_{CN} \left[ LY\_N + RY\_N + VY\_N + WY\_N + DY\_N \right.\\ &+ UY\_N + PY\_N + QY\_N + IY\_N - YI\_N - YE\_N \right] \end{aligned} \tag{4}$$

Here, net C gain by the ectomycorrhizal fungi (EYc–YEc) is equal to the C:N ratio of the ectomycorrhizal pool (YCN) multiplied by its net N gain. Net N gain is calculated by adding up N acquired from leaf/wood litter (LY<sup>N</sup> + RYN), root litter (VY<sup>N</sup> + WYN), DOM (DY<sup>N</sup> + UYN), POM (PY<sup>N</sup> + QYN), and inorganic N (IYN), and subtracting inorganic N losses through mineralization (YIN), and N provided to plants (YEN; Equation 4). In the example given in **Figure 2C**, no substrate use efficiency or stoichiometric constraints are imposed and the net N gain is constrained between 0 and 8 g N m−<sup>2</sup> y −1 . This means that without a fixed C:N ratio, the range of flows to POM, MAOM, and DOM as well as the quantity of N lost can be between 0 and 8 g N m−<sup>2</sup> y −1 . In contrast, if we define the C:N ratio of ectomycorrhizal fungi as 10 (**Figure 2D**), the potential flows for N lost to POM, MAOM, and DOM have a maximum of only 1 g N m−<sup>2</sup> y −1 . This change occurs because ectomycorrhizal fungi mineralize or send to plants up to 7 g N m−<sup>2</sup> y −1 to maintain their C:N ratio.

With each additional stoichiometric relationship defined in the model (i.e., going from C:N to C:N:P), a further equation of the same form as Equations 3 or 4 is added to represent the partitioning of each element. Each additional equation is likely to shrink the difference between the minimum and maximum flows, thus reducing uncertainty. However, the impact of each constraint on the entire LSD network is not obvious a priori because not all constraints will be binding. For example, the C:N constraint we described does not reduce uncertainty in C flow, because there is enough N to span the entire range of available C (i.e., N is limiting; **Figure 2**). Consequently, we can use the LSD as a tool to determine if and when stoichiometry and substrate use efficiency matter in the network of ∼400 flows that we consider (**Figure S1**).

#### Application of the LSD

We used the LSD to consider three fundamental stoichiometric biases: soil organism homeostasis (stoichiometric demands), mineral sorption biases, and abiotic litter loss biases (i.e., leaching through DOM). We then evaluated how three classes of empirical information constrained the state-space in which our linear system could exist. First, we considered the effects of broad stoichiometric ranges (e.g., animal C:N ratio is between 4 and 25; Sterner and Elser, 2002). Second, we considered specific stoichiometric constraints, where the exact C:N:P:S ratios of each process are known (e.g., animal C:N ratio is 6; Equation 4). Finally, we considered the case where substrate use efficiency is a set as at exact value (Equation 3). For each, we assessed the effect of these data on model uncertainty, defined as the range between the minimum and maximum flows possible under the prescribed constraints.

We also quantified the value of stoichiometry and substrate use efficiency information across 27 different biome × ecosystem type combinations. These systems have different animal and microbial communities, as well as different mineral properties, albeit with varying levels of data to constrain them (Chapin et al., 2011; Bond-Lamberty and Thomson, 2014; Adams et al., 2018; Steidinger et al., 2019). We first confirmed that the LSD was producing reasonable flows by comparing its litter inputs and heterotrophic respiration estimates to values from the SRDB (version 3.0) dataset (Bond-Lamberty and Thomson, 2014; **Figure S6**). Then, we used the LSD to examine how stoichiometric vs. substrate use efficiency information constrained the range of possible C and nutrient flows into stabilized pools across ecosystem types.

# Solving the LSD and Quantifying Uncertainty

We solved the LSD using inverse linear programming in R (van Oevelen et al., 2010; R Core Team, 2014; see section Model Example for a case study on one pool). We explored the importance of stoichiometry and substrate use efficiency in the LSD across differences in structure, parameterization, and type of stoichiometric information available, as detailed above. Structures varied the number of microbial pools or the presence of different stabilization pathways. Modifying parameters allowed us to confirm that our results were general across possible parameter ranges, and modifying the pools that receive constraints allowed us to assess the impact of stoichiometric or substrate use efficiency information for different pools.

We present these data at the level of overall uncertainty to test whether changes in structural conditions or parameters influence the value of stoichiometric information. Here, overall uncertainty means the range of possible flows across the entire linear system. To calculate a value of overall uncertainty (Equation 6), we solved for the individual flows that maximize (ymax) and minimize (ymin) the total flow throughout the linear system. Uncertainty is the overall difference between all i flows in the maximum and minimum cases for each element.

$$\text{Uncertainty} = \sum\_{i}^{N} |\nu\_{\text{max},i} - \nu\_{\text{min},i}| \tag{5}$$

Given that the LSD is static and does not explicitly account for pool sizes and rate constants, our uncertainty estimates are not directly comparable with those derived from dynamic SOM (and carbon cycle) models, which do not always operate at equilibrium (Sihi et al., 2018; Sulman et al., 2018). Instead, the magnitude of uncertainty in the LSD indicates the precision of each flow at equilibrium. Although this metric may differ from the general measure of uncertainty used within the soil (and terrestrial) ecology community, inverse linear programming has been extensively used to solve network flow problems and their associated uncertainties in community ecology, metabolic networks, engineering, nutrition, and transport planning (Dantzig, 1963; van Oevelen et al., 2010; Yang et al., 2017; Goldford et al., 2018).

# RESULTS

The general model scenario constrained soil processes within ranges set by the described structural conditions (see section Model Constraints) and by mass balance. The uncertainty calculated for the base model was at least two orders of magnitude higher than the system inputs (i.e., 180–2,500 g C m−<sup>2</sup> y −1 ; assumption 1). Uncertainty was higher than the system inputs because many flows recycled C and nutrients among pools in the soil system (This does not occur in the one pool example provided; **Figure 2**).

Specific stoichiometric constraints, such as an exact C:P ratio for bacteria, had a much more important impact on uncertainty than broad stoichiometric ranges (**Figure 3A**). Adding broad constraints on the relationship between C, N, P, and S encompassing all potential minimum to maximum ratios had a negligible impact on uncertainty in C (**Figure 3A**), but reduced uncertainty in N, P, and S flows whenever that nutrient was included (**Figure S2**). We tested the system using specific stoichiometric constraints (fixed C:N:P:S ratios across all pools) and found that they reduced the uncertainty in C flow by 20% when we included only C:N and by 80% when we included at least C:N:P (**Figure 3A**).

We then tried including substrate use efficiency for C, N, P, and S rather than imposing stoichiometry (**Figure 3B**). Information on CUE for microorganisms and soil animals was effective at reducing the uncertainty in C flows (**Figure 3B**) but the substrate use efficiency of other nutrients (i.e., N, P, and S) did not provide additional benefits (**Figure 3B**; **Figure S2**). Unlike specific stoichiometric information, which became more useful with additional elements, most of the value of substrate use efficiency data came from CUE (**Figure 3B**; **Figure S2**). The general trends in uncertainty across stoichiometry and substrate use efficiency were robust to changes in the parameter values (**Figure 3**: error bars). These results show that stoichiometry, especially precise C:N:P ratios, is on par with substrate use efficiency at constraining C flows (**Figure 3**).

Modifying the model structure (see section Model Example) altered the overall uncertainty and, in certain cases, altered the value of stoichiometric and substrate use efficiency information. Removing differences in chemical quality and microbial pools dramatically reduced the overall uncertainty, while removing MAOM pools and physical protection had a more modest effect

FIGURE 3 | A summary of C cycle uncertainty in the LSD using different amounts of stoichiometry and substrate use efficiency information. (A) The scenarios represent an increasing number of links between element cycles, ranging from entirely independent to linkage between "C:N," "C:N:P," and "C:N:P:S." The general scenario has the widest possible bounds, while still avoiding infinite recycling between pools. Broad scenarios are constrained by a wide range of potential ratios between C and the indicated nutrients, while specific scenarios are constrained to a point estimate of that ratio. (B) Scenarios with CUE, NUE, PUE, and SUE are constrained to a point estimate of C, N, P, and S use efficiency, respectively. Uncertainty is shown for the general model containing all pools and information along with error bars that span the 25th−75th quantiles of 100 random parameters draws within ± 20% of the base value. Uncertainty is calculated as the sum of the absolute difference in flows between the minimum and maximum flow scenarios. See Figure S2 for N, P, and S.

(**Figure 4A**: inset). The fact that differences in chemical quality, microbial biomass, and MAOM increased uncertainty in the model is consistent with our new understanding of them as central to SOM formation and loss (Lehmann and Kleber, 2015). When microbial pools were taken out of the budget, substrate use

efficiency information had no effect on uncertainty even though animals were still included (**Figure 4B**). Removing chemical quality from the budget also made substrate use efficiency information less useful. Similarly, stoichiometric ratios were much less important without differences in chemical quality,

MAOM pools, or saprotrophic microbial pools (**Figure 4A**). Since recent empirical findings indicate that such interactions drive the formation of large pools of stable SOM (Lehmann and Kleber, 2015), our results establish stoichiometry as a necessary component of SOM models that focus on microbial and mineral SOM formation. These results highlight an important interaction between stoichiometry and the structural assumptions used to build predictive models of SOM formation and loss.

Changing the organisms or pools that received stoichiometric or substrate use efficiency constraints altered their impact on uncertainty. For reducing C flow uncertainty, stoichiometric and substrate use efficiency constraints on saprotrophic microorganisms (i.e., bacteria and fungi) were the most important (**Figure 5**). Stoichiometric constraints on MAOM also reduced uncertainty, but were much less impactful (**Figure 5A**). For N flow, constraints on symbionts were also important. For P and S, constraints on MAOM stoichiometry reduced uncertainty (**Figure S4**). Overall, specific stoichiometric data for microbial C:N:P ratio was the most useful for reducing model uncertainty (by > 50%).

Nutrient flows responded differently to a reduction in the possible range of pool stoichiometry and organism CUE across ecosystem types. The differences reflect both the pathways present in a given ecosystem (analogous to structural changes) and starting uncertainty in the potential range of CUE or C:N:P:S for each ecosystem type (**Table S2**). For example, we found that stoichiometric information is much more effective at constraining flows in tropical systems, where uncertainty in stoichiometric ratios is high and all three plant symbionts can occur (**Figure 6**; Steidinger et al., 2019). In contrast, stoichiometry and CUE could both reduce uncertainty in temperate systems, with CUE being particularly effective given the wide range of CUE estimates available for temperate forests (**Figure 6**; Qiao et al., 2019). In general, specific stoichiometric information was more effective than substrate use efficiency at constraining the flows of all elemental cycles (**Figure 6**). Microbial and animal CUE was most effective at constraining C (**Figure 6**). Even though temperate forests have more data with which to estimate CUE, the trend in **Figure 6** held when we made the range of initial CUE values the same for temperate and tropical systems.

In general, our analysis indicated that acquiring more precise stoichiometric data will reduce uncertainty in soil nutrient cycling, confirming our hypothesis and the conclusions of our qualitative review (**Figure 6**; **Figure S5**). We also found that a budgeting approach, like the LSD, can quantify the usefulness of new stoichiometry and CUE data for a particular system.

### CONCLUSIONS AND FUTURE DIRECTIONS

Our review traced several of the many stoichiometric changes that can occur along the path separating soil inputs from SOM. With the LSD, we quantified whether these stoichiometric relationships could constrain uncertainty in the size of SOM flows at equilibrium. We discovered that adding stoichiometric information narrowed the range of possible flows, and that the inclusion of elements other than C and N substantially improved the precision of overall flow patterns. In fact, stoichiometric data allowed us to constrain all elemental cycles, improving our representation not only of nutrient flows but also of C flows.

The new generation of SOM models generally use CUE parameters, either fixed or as a function of abiotic conditions,

to constrain the flow of C inputs into microbially-derived SOM pools (Allison et al., 2010; Wieder et al., 2014). Though appealing for its simplicity, this approach cannot capture the dependency of CUE on resource availability and often relies on empirical point measurements, which can obscure the variability inherent in CUE (Manzoni et al., 2012). Fixed parameters also do not capture the changes in CUE at different levels of mineral nutrient availability, meaning that models relying on prescribed CUE parameters instead make implicit stoichiometric assumptions. By contrast, stoichiometrically explicit models allow this component of CUE to emerge dynamically as a function of nutrient availability. Our findings demonstrate that a stoichiometric approach can effectively reduce uncertainty in C flow to SOM, justifying the addition of multiple elements into a dynamic accounting of SOM (Schimel and Weintraub, 2003; Buchkowski et al., 2015, 2017; Sulman et al., 2019). Parameterizing such a model would require stoichiometric measurement of important pools, but this is likely to be both more cost-effective and more consistent than measuring substrate use efficiencies at a large scale. Our results also show that such SOM models could be effectively constrained with open-source C and nutrient data on soil and plant inputs available from environmental research and observation networks (Weintraub et al., 2019).

Measurement of SOM pool stoichiometry is a clear opportunity for high-value empirical work. For instance, our review and analysis strongly support several recent theoretical advances (Lehmann and Kleber, 2015; Jilling et al., 2018) highlighting the importance of the microbial to MAOM transition as a hub of C and nutrient exchange in the SOM system. Unfortunately, data on the nutrient content of the MAOM pool remain sparse, especially outside of agricultural and temperate forest ecosystems and for ratios other than C:N. While flows into this pool are predicted to be nutrient-enriched, the quality and extent of this enrichment is highly dependent upon soil characteristics and processes, such as sorption, that are non-linear and require data across substrate × mineral combinations to parameterize SOM models.

Despite the fact that stoichiometry is more consistent than CUE (Cleveland and Liptzin, 2007), the stoichiometry of MAOM and other large soil pools can be sensitive to disturbance, and may therefore be vulnerable to global change. For example, fire (e.g., Butler et al., 2019), warming (e.g., Sihi et al., 2019), and nutrient addition (e.g., Crowther et al., 2019) have been known to influence the stoichiometry of soil, invertebrate, saprotrophic microorganisms, mycorrhizal fungi, and enzyme kinetic activities, which in turn can affect SOM pools and flows. Likewise, MAOM stoichiometry is likely to change with soil pH and associated organo-metal complexation with sesquioxides (Fe and Al oxides) and exchangeable Ca (Rasmussen et al., 2018). Adding these stoichiometric mechanisms into dynamic SOM models may help them replicate larger-scale changes in soil stoichiometry and resolve global patterns.

By synthesizing recent advances in SOM cycling and developing a simple, stoichiometrically explicit budget, we were able to quantify the size of organic matter flows in soil and identify which SOM flows were well-constrained by new data across 27 different ecosystem × biome combinations. Information not only on soil organic C and N but also on P and S was central to achieving this goal, which indicates that collecting data on the distribution of these soil nutrients is an effective way to strengthen our capacity to predict SOM stocks and flows. With our work, we have provided a foundation for incorporating these new stoichiometric data into dynamic SOM models, where they have the potential to significantly enhance accuracy. Together, empirical and theoretical advances like these can allow us to move away from equilibrium and toward a better understanding of the non-linear dynamics that characterize our changing world.

#### DATA AVAILABILITY STATEMENT

The model code and data necessary to recreate our analyses is available on Zenodo at https://doi.org/10.5281/zenodo.3464585.

#### REFERENCES


#### AUTHOR CONTRIBUTIONS

RB initiated the project. RB and DS ran the linear model analyses. All authors defined the scope of the project, synthesized the available data, and wrote and revised the manuscript.

#### FUNDING

The work of RB was supported by the Yale School of Forestry & Environmental Studies. The work of AS was supported by NSF grant DEB-1556643. The work of DS was funded by the U.S. Department of Energy Office of Biological and Environmental Research through the ORNL TES SFA program and an Early Career Award. ORNL is managed by the University of Tennessee-Battelle, LLC, under contract DE-AC05–00OR22725 with the US DOE. The work of GS was funded by a US NSF Graduate Research Fellowship. The work of AK was funded by a NatureNet Science Fellowship.

#### ACKNOWLEDGMENTS

The ideas in this manuscript were synthesized at the Woodstoich 4 workshop, funded by the National Science Foundation under the award DEB-1840408. The authors are grateful to Jim Elser, Michelle A. Evans-White, Jen Schweitzer, Eli Fenichel, Craig See, Paul Julian, Joe Vanderwall, and all the participants in the Woodstoich 4 workshop for their valuable advice and feedback. The final paper was substantially improved by invaluable comments from Josh Schimel and two reviewers.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00382/full#supplementary-material


plant inputs form stable soil organic matter? Glob. Chang. Biol. 19, 988–995. doi: 10.1111/gcb.12113


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a shared affiliation, though no collaboration, with one of the author AS.

Copyright © 2019 Buchkowski, Shaw, Sihi, Smith and Keiser. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Elemental Ratios Link Environmental Change and Human Health

Rachel E. Paseka<sup>1</sup> \*, Anika R. Bratt <sup>2</sup> , Keeley L. MacNeill <sup>3</sup> , Alfred Burian<sup>4</sup> and Craig R. See<sup>1</sup>

<sup>1</sup> Department of Ecology, Evolution, and Behavior, University of Minnesota, St. Paul, MN, United States, <sup>2</sup> Environmental Studies Department, Davidson College, Davidson, NC, United States, <sup>3</sup> Department of Forestry, Ecosystems, and Society, Oregon State University, Corvallis, OR, United States, <sup>4</sup> Environmental Sustainability Research Centre, University of Derby, Derby, United Kingdom

Humans have fundamentally altered the cycling of multiple elements on a global scale. These changes impact the structure and function of terrestrial and aquatic ecosystems, with many implications for human health. Most prior studies linking biogeochemical changes to human health have evaluated the effects of single elements in isolation. However, the relative availability of multiple elements often determines the biological impact of shifts in the concentration of a single element. The balance of multiple elements is the focus of ecological stoichiometry, which highlights the importance of elemental ratios in biological function across all systems and scales of organization. Consequently, ecological stoichiometry is a promising framework to inform research on the links between global changes to elemental cycles and human health. We synthesize evidence that elemental ratios link global change with human health through biological processes occurring at two scales: in the environment (natural ecosystems and food systems) and within the human body. Elemental ratios in the environment impact the key ecosystem processes of productivity and biodiversity, both of which contribute to the production of food, toxins, allergens, and parasites. Elemental ratios in diet impact processes within the human body, including the function and interactions of the immune system, parasites, and the non-pathogenic microbiome. Collectively, these stoichiometric effects contribute to a wide range of non-infectious and infectious diseases. By describing stoichiometric mechanisms linking global change, ecological processes, and human health, we hope to inspire future empirical and theoretical research on this theme.

Keywords: ecological stoichiometry, biological stoichiometry, public health, global change, food security, water quality, nutrition, disease

# INTRODUCTION

Understanding the effects of global change on human health is a major challenge for the twenty-first century (Leaf, 1989; Whitmee et al., 2015). Anthropogenic activities drive substantial shifts in biogeochemical cycles in ecological systems throughout the world (Peñuelas et al., 2013), many of which can impact human health. For example, application of nitrogen (N) fertilizer directly improves crop yields and reduces malnutrition but indirectly contributes to the prevalence and severity of multiple infectious and non-infectious diseases (Townsend et al., 2003). Similarly, increases in atmospheric carbon dioxide (CO2) concentrations driven by fossil fuel combustion reduce the nutritional quality of food crops and other plants, with myriad human health outcomes (Loladze, 2002).

#### Edited by:

James Joseph Elser, University of Montana, United States

#### Reviewed by:

Robert Ptacnik, Wasser Cluster Lunz, Austria Irakli Loladze, Bryan College of Health Sciences, United States Krista Capps, University of Georgia, United States

> \*Correspondence: Rachel E. Paseka rachelpaseka@gmail.com

#### Specialty section:

This article was submitted to Behavioral and Evolutionary Ecology, a section of the journal Frontiers in Ecology and Evolution

> Received: 21 July 2019 Accepted: 20 September 2019 Published: 10 October 2019

#### Citation:

Paseka RE, Bratt AR, MacNeill KL, Burian A and See CR (2019) Elemental Ratios Link Environmental Change and Human Health. Front. Ecol. Evol. 7:378. doi: 10.3389/fevo.2019.00378

Biogeochemical effects on human health are typically studied from a single-element perspective. However, the cycling of carbon (C), N, and other elements relevant to health are often linked through biogeochemical mechanisms and shared anthropogenic drivers (Peñuelas et al., 2012), and the ratios of these elements in nature are fundamental to biological function (Sterner and Elser, 2002). Therefore, we need a multielemental approach to better understand the mechanisms linking anthropogenic effects on biogeochemical cycles with human health. The balance of multiple elements in nature is the focus of ecological stoichiometry (ES), a research framework that highlights the importance of elemental ratios in biological function across all systems and scales of organization (Sterner and Elser, 2002).

Stoichiometric theory has been applied to many fundamental research topics in ecology and evolution (Hessen et al., 2013; Van de Waal et al., 2018). However, ES has been used relatively infrequently in socio-ecological research (but see Ptacnik et al., 2005; Cease et al., 2015) and has been applied to human health in just a few studies (**Table 1**). Given the central importance of elemental ratios to biological function, we suggest that ES is a powerful tool to study links between global change and multiple dimensions of human health.

Here we outline how stoichiometric effects on ecological processes (hereafter, "stoichiometric effects") link anthropogenic activities, ecological function, and human health. We outline stoichiometric effects at two major scales: in the environment (natural ecosystems and food systems) and within the human body. In the environment, anthropogenic activities including fossil fuel combustion, fertilizer application, land use change, and industrial production drive changes to elemental ratios in air, water, and soil. The alteration of these ratios impacts ecosystem productivity and biodiversity, both of which contribute to the production of food, toxins, allergens, and parasites. Stoichiometric constraints on the quantity and quality of food production have implications for food security, and elemental ratios in human diet impact immune function, metabolism, parasite function, and the microbiome. Collectively, these stoichiometric effects contribute to a wide range of noninfectious and infectious diseases (**Figure 1**).

By synthesizing evidence that elemental ratios link global change and human health, we hope to inspire future stoichiometric research on this theme. To ensure that this synthesis is accessible to a broad audience, including readers who lack a background in ES, we summarize fundamental aspects of stoichiometric theory (**Figure 2**) and begin by describing links between human actions, elemental ratios in the environment, and ecological processes that impact human health.

# A BRIEF PRIMER ON GLOBAL CHANGE AND ENVIRONMENTAL STOICHIOMETRY

All life on earth is composed of ∼25 elements (Kaspari and Powers, 2016), each of which has a unique biogeochemical cycle (Schlesinger and Bernhardt, 2013). Human activity has dramatically changed the absolute and relative availabilities of these elements around the world, which can limit the growth and function of organisms in both terrestrial and aquatic environments.

Fossil fuel combustion has increased the concentration of CO<sup>2</sup> in the environment globally, as well as other elements found in fossil fuels such as N and sulfur (S) (Smith et al., 2001). Fertilizer application has increased the availability of N and phosphorus (P) in croplands and adjacent landscapes, including aquatic ecosystems. Commercial fertilizer often includes additional elements, such as potassium (K), calcium (Ca), magnesium (Mg), and S that can also be transported via dust, erosion, or water (Kaspari and Powers, 2016). Biomass burning (to clear land) releases C from plant material back to the atmosphere, but it can also enrich local environments with trace elements, especially K (Sardans and Peñuelas, 2015) and zinc (Zn, Echalar et al., 1995). Industrial production has released trace elements into the environment, such as mercury (Hg) and arsenic (As, Adriano, 1986). These changes in the availability of elements in the environment occur at different spatial scales. For example, the elevation of atmospheric C occurs globally, while N enrichment is a more localized or regional phenomenon (Galloway et al., 2008).

Changes to biogeochemical cycles in aquatic ecosystems are complicated by interactions with soils, plants, and microbes along hydrologic pathways. Fossil fuel combustion has led to increased concentrations of dissolved CO<sup>2</sup> in aquatic ecosystems, especially oceans, thereby increasing the availability of C relative to other elements (Caldeira and Wickett, 2003). Terrestrial inputs increase the availability of N and P in freshwaters, but P is preferentially retained by soils relative to N. Nitrogen deposition can increase the P-limitation of aquatic ecosystems due to increasing N:P (Elser et al., 2009) but analysis of rivers across the United States shows that the change in N:P is not unidirectional and is influenced by abiotic factors (Dodds and Smith, 2016). Silica (Si) is also increasing in streams and rivers due to increased erosion, but it can be retained by dams (Humborg et al., 1997, 2006; Carey and Fulweiler, 2012).

These shifts in the availability of multiple elements are not occurring in synchrony across space and time, leading to changes in elemental ratios in the environment. We present a non-exhaustive list of human-driven changes to environmental stoichiometry (**Table 2**) and briefly discuss their effects on ecosystem productivity and patterns in biodiversity, both of which are linked to multiple dimensions of human health.

#### Ecosystem Productivity Primary Productivity

The absolute and relative concentrations of elements in the environment commonly limit the primary productivity of both aquatic and terrestrial ecosystems (Elser et al., 2007b). While traditional perspectives on this topic often focus on a single limiting element (sensu Liebig, 1855), recent studies highlight the importance of co-limitation, the simultaneous limitation of primary productivity by multiple elements, across aquatic and terrestrial ecosystems (Elser et al., 2007b; Harpole et al., 2011; Fay

#### TABLE 1 | Examples of prior research on ecological stoichiometry in a human health context.


FIGURE 1 | A stoichiometric framework linking global change with human health. Stoichiometric effects (ecological processes mediated by elemental ratios, dashed lines) occur at two scales: in the external environment (to the left of the human figure) and within the human body (gray circle to the right of the human represents internal processes, especially diet-mediated processes within the gut). Global environmental changes (icons represent fossil fuel combustion, land use change, fertilizer application, and industrial pollution) alter the absolute (concentration) and relative (ratios) availability of elements in the environment (1). Environmental stoichiometry impacts the key ecosystem processes of biodiversity and productivity (2). Both of these processes contribute to toxin, allergen, and parasite production (3), all of which can impact human health directly by causing non-infectious and infectious disease (4). Stoichiometric effects on biodiversity and productivity also determine food production in both managed (e.g., agricultural) and unmanaged (e.g., marine and freshwater) systems (3). The quantity and quality of food production can impact human health directly (4), but these impacts may be mediated by stoichiometric effects occurring within the human body (5). Elemental ratios in human diet impact the function of parasites, the non-pathogenic microbiome, and the immune system (6), which collectively influence both non-infectious and infectious disease (7).

et al., 2015). Primary productivity is often synergistically limited by N and P, such that additions of N and P together increase primary productivity beyond the additive effects expected based on individual inputs of these elements (Elser et al., 2007b; Harpole et al., 2011). Given that life is constructed from over 25 elements, it is also crucial that we expand our view of elemental limitation and co-limitation to include elements beyond N and P (Kaspari and Powers, 2016). The relative availability of elements is also important for productivity in managed systems, such as agriculture, where fertilizer application aims to provide the correct balance of N, P, and other elements to maximize crop yield (Van der Velde et al., 2014).

FIGURE 2 | A primer of ecological stoichiometry and an overview of its application to human health.

TABLE 2 | Human drivers of altered elemental ratios in the environment and ecological consequences relevant to human health.


Global change symbols indicate drivers (fossil fuel combustion, fertilizer application, land and water alteration, industrial production). Numbers refer to publications listed at the bottom of the table. These examples highlight how changes to the cycling of a single element impact elemental ratios in the environment.

Loladze (2002); Sardans and Peñuelas (2012); Bezemer and Jones (1998); Rusterholz and Erhardt (1998); Lake and Hughes (1999); Thompson et al. (1993); Coakley et al. (1999); Selin (2009); Ptacnik et al. (2005); Van de Waal et al. (2009); Van de Waal et al. (2014); Bobbink et al. (1998); Cisneros and Godfrey (2001); Peñuelas et al. (2013); Tilman (1996); Tilman (1997); Bobbink et al. (2010); Stevens et al. (2010); Payne et al. (2013); Maavara et al. (2015); Flynn (2002); Granéli and Johansson (2003); <sup>23</sup>Pan et al. (1996); Peñuelas et al. (2012); Gilmour et al. (1992); Humborg et al. (2006); Humborg et al. (1997); Carey and Fulweiler (2012); Parsons and Dortch (2002); <sup>30</sup>Zhang et al. (2017); Bjerregaard et al. (2011).

#### Consumer Productivity

Elemental ratios in primary producers influence how energy and elements move through food webs to fuel consumer productivity at higher trophic levels (Sterner and Elser, 2002). Mismatches between the ratios of elements required for consumer growth and the elemental ratios of their resources lead to elemental imbalance, which can constrain productivity at any trophic level (Boersma et al., 2008). The effects of elemental imbalance on trophic transfer may be exacerbated when interacting organisms differ in their degree of stoichiometric homeostasis (maintaining consistent organismal stoichiometry despite changes to resource stoichiometry) vs. plasticity (shifting organismal stoichiometry to reflect changes in resource stoichiometry) (Meunier et al., 2014). Algae and terrestrial plants show much higher stoichiometric plasticity than herbivores and other consumers, which can result in substantial elemental imbalance and reduce consumer growth efficiency (Boersma, 2000, Demott et al., 2010). Therefore, changes to elemental ratios in the environment can propagate through food webs to impact productivity at all trophic levels.

#### Biodiversity

The relative availability of elements in the environment can influence the local composition and richness of ecological communities. Primary producers vary among species in growth response to environmental C:N:P, influencing community composition across ecosystems (Agren et al., 2012). For example, both legumes and N-fixing cyanobacteria have a competitive advantage in low N:P environments (Schindler, 1977; Vitousek and Walker, 1989). Stoichiometric effects on the competitive dominance of primary producers may also impact community composition at higher trophic levels (Singer and Battin, 2007), although stoichiometry is often overlooked as a potential driver of biodiversity patterns.

Stoichiometric effects may also contribute to global declines in species richness. For example, N deposition reduces species diversity in grasslands on a continental scale (Stevens et al., 2010). This is likely driven by both non-stoichiometric (acidification, Stevens et al., 2010) and stoichiometric effects (shifting N:P ratios altering competition; Olde Venterink et al., 2003), which are difficult to separate. The global increase in N availability (primarily driven by N-deposition) without proportionate increases in other elements may contribute to large-scale declines in species richness (sensu Tilman et al., 1982).

# FOOD SECURITY AND WATER QUALITY

Two key aims of the United Nations Sustainable Development Goals are to achieve food security and to ensure access to clean water for all people by 2030 (UN General Assembly, 2015). Food security incorporates both food quantity and food quality. Currently, 820 million people suffer from acute food shortage, and ∼2 billion experience chronic deficiency of essential nutrients (FAO, 2019). Additionally, 785 million people lack access to potable water, including 144 million who rely on untreated surface water (WHO and UNICEF, 2019). Elemental ratios in air, water, and soil impact both food security and water quality (e.g., Loladze, 2002, Mueller et al., 2012). We summarize stoichiometric effects on crop production, rangeland and pasture-based production of meat and dairy, water quality, and fisheries. Collectively, these stoichiometric effects on food security and water quality propagate to impact both noninfectious and infectious disease.

# Food Crop Production

Between 1.2–1.5 billion hectares are used for crop production globally (O'Mara, 2012), and yields have increased consistently from 877 million metric tons in 1961 to 2,351 million tons in 2007 (Tester and Langridge, 2010; Foley et al., 2011). However, ongoing human population growth and changes in consumption patterns necessitate continued increases in crop production, with predicted requirements of more than 4,000 million metric tons in 2050 (Tester and Langridge, 2010; Foley et al., 2011). Given the potential for primary productivity to be limited by the availability of multiple elements, understanding the stoichiometry of food crop production under global change is essential to achieving these production goals.

For example, P-limitation of crops is a global concern that may be exacerbated by stoichiometric effects. In Sub-Saharan Africa and Eastern Europe, soil P availability is a major factor limiting productivity (e.g., 75% of soils in Sub-Saharan Africa are P-deficient; Cordell et al., 2009; Mueller et al., 2012). However, the global availability of mineable P is limited, which will increase prices of P fertilizers (Van Vuuren et al., 2010), exacerbate agricultural P-limitation, and threaten global food production (Cordell et al., 2009). Human-induced increases in environmental C and N availability outpace increases in P availability, intensifying C:N:P imbalances globally (Peñuelas et al., 2013; Carnicer et al., 2015).

Global shifts in the cycling of C and N impact food security through additional stoichiometric effects on food quantity and quality. Increases in atmospheric CO<sup>2</sup> concentrations tend to increase crop yields (Ainsworth and Long, 2005; Kimball, 2016), with positive impacts on food quantity. However, crops produced under elevated CO<sup>2</sup> exhibit increased C:element ratios (e.g., C:N, P, Se, Fe, or Zn), lowering nutritional quality (Loladze, 2002). Nitrogen deposition or over-fertilization can leach K from soils and increase the likelihood of plant K-limitation (Sardans and Peñuelas, 2015), which can alter plant biomass allocation and increase drought vulnerability (Cakmak et al., 1994). Additionally, the availability of many elements in soil (including N, P, K, Zn, and others) can impact the prevalence and severity of infectious disease in crops through several mechanisms, with substantial consequences for yields (Dordas, 2008; Veresoglou et al., 2013). These effects have not previously been studied from a stoichiometric perspective, although elemental ratios mediate many infectious disease processes (Sanders and Taylor, 2018). Overall, human-induced shifts in the ratios of multiple elements in the environment lead to changes in both the quantity and quality of food crop production.

## Rangeland and Pasture-Based Production of Meat and Dairy

Grasslands account for ∼60% of agricultural land world-wide and support global meat and dairy production (O'Mara, 2012). Like crop production, the productivity of grasslands depends on both the relative and absolute availability of elements in the environment. Elevated atmospheric CO<sup>2</sup> and N deposition generally increase grassland productivity (Thornton et al., 2009; Hatfield and Prueger, 2011; Chapman et al., 2012), but the effects of CO<sup>2</sup> also contribute to increases in plant C:N that represent decreased forage quality for livestock (Craine et al., 2017; Augustine et al., 2018).

# Water Quality

The relative balance of N, P, and other elements in aquatic ecosystems are major drivers of algal bloom dynamics, and anthropogenic inputs of these elements threaten the quality of surface water as it pertains to human health. N:P ratios have long been viewed as drivers of algal blooms (Schindler, 1977), and more recent work has shown that the relative availability of dissolved elements can influence water quality in a variety of ways. For example, low N:P ratios provide a competitive advantage for N-fixing cyanobacteria (Mazur-Marzec et al., 2006), which can form harmful algal blooms (HABs). Elemental ratios in water impact both the abundance of harmful algae or flagellates and the amount of toxins they produce (increasing toxin production in low C:N, low C:P, and high N:P waters; Townsend et al., 2003; Ptacnik et al., 2005, Van de Waal et al., 2009; Anderson et al., 2012). Toxins from HABs have threatened drinking water supplies of entire cities (e.g., Toledo, USA and Wuxi, China; Liu et al., 2011; Steffen et al., 2017), and the specific effects of HABs on human health are discussed in the non-infectious disease section below.

#### Fisheries

Fish and shellfish (hereafter, "fish") represent ∼5% of global protein consumption (FAO, 2000). However, fish provide a major source of essential fatty acids, vitamins, and other nutrients, thereby playing an important role in global food quality (Golden et al., 2016). The threshold of sustainable harvests in natural systems is strongly linked to primary productivity and trophic transfer efficiencies through food webs (Brander, 2007; Boyce et al., 2010), both of which are subject to stoichiometric constraints (Winder et al., 2009). Eutrophication resulting from anthropogenic changes to biogeochemical cycles can alter fish production in both marine and freshwater systems, and HABs can lead to fish kills (Camargo and Alonso, 2006; Anderson et al., 2012). Eutrophication also impacts the trophic transfer of fatty acids through food webs, reducing the nutritional value of fish for human consumption (Taipale et al., 2016).

# NON-INFECTIOUS DISEASE

Elemental ratios link global change with non-infectious disease (encompassing a broad suite of physical and mental health outcomes) at two scales: through the effects of nutrition on processes occurring within the human body and through environmental processes occurring outside the body (**Figure 3**).

## Nutritional Stoichiometry and Non-infectious Disease

Traditional views of human nutrition focus on how imbalances between dietary macromolecules and nutritional requirements impact health, including a broad range of non-infectious diseases. Ecological stoichiometry provides a parallel perspective on this imbalance, with the focus on elemental ratios rather than macromolecules. While ES is a reductionist approach that does not capture the full complexity of nutritional biochemistry, elemental ratios provide a common currency to link food production in the environment with diet quality and its consequences for health.

#### Dietary Stoichiometry and Physical Health

Perhaps the clearest stoichiometric link between global change and human dietary health is evidence that elevated atmospheric CO<sup>2</sup> leads to higher C:element ratios in plants, including staple crops (Loladze, 2002, 2014). This reduction in crop quality means that consumption of staple crops in the future will correspond to higher caloric intake relative to nutritional value (as protein, micronutrient, and vitamin content). This change may contribute to caloric overconsumption, a problem that already drives many major public health crises, including global diabetes and obesity epidemics (James, 2008; Zimmet et al., 2010). Increasing C:element ratios in crops may also contribute to a broad range of health effects caused by mineral and vitamin deficiencies. Along with decreases in Ca, K, Mg, Fe, Zn, and Cu concentrations (Loladze, 2014), elevated CO<sup>2</sup> leads to declines in many essential vitamins in rice (Zhu et al., 2018). Among the B vitamins which decline under elevated CO2, the more N-rich forms (such as folate and thiamine) show the greatest reductions in concentration (Zhu et al., 2018).

Changes to elemental ratios in diet can also impact human health through overconsumption of micronutrients and interactions among elements that impact their absorption within the gut. For example, most people in the United States consume far more than the recommended daily allowance of P, and overconsumption of P is associated with increased mortality rates (Chang et al., 2014). However, P content is not required on the labels of packaged foods, making consumption difficult to monitor (Calvo et al., 2014). The overconsumption of P may be particularly important given the high concentrations of readily absorbed inorganic P present in many processed foods (Ritz et al., 2012), which are absorbed much more efficiently (80–100%) in the gut than organically bound P forms (40– 60%; Calvo et al., 2014). Stoichiometric effects may underlie some of the observed links between P overconsumption and negative health outcomes. For example, the ratio of dietary P:Ca is important to maintaining normal physiological function and bone health, independently of the consumption of P and Ca individually (Brot et al., 1999; Ito et al., 2011; Lee et al., 2014). Similar stoichiometric effects on other dietary micronutrients are also possible. For example, interactions among Fe, Zn, and Cu impact their absorption and bioavailability within the gut (Sandström, 2001). Overall, these examples suggest that a stoichiometric perspective can help to link global changes in food production, human diet quality, and the effects of nutrition on non-infectious disease.

Diet stoichiometry may also have consequences for cancer due to the elemental requirements of tumor growth. One particularly interesting prior application of ES to human health tested the growth rate hypothesis in tumor dynamics. The growth rate hypothesis (GRH) states that fast growing cells have lower N:P requirements than slower growing cells due to the high P content of ribosomal RNA (Elser et al., 1996, 2003). Rapidly growing lung and colon tumors had 2-fold higher P content and lower N content than normal tissue (although this was not the case for kidney or liver tumors), suggesting some degree of P limitation in cancer cells (Elser et al., 2006, 2007a). While these studies did not specifically test the effects of dietary P content on tumor dynamics, the ultimate source of P used by cancer cells is human diet. High concentrations of inorganic P in food lead to increased tumor proliferation and growth in mice (Jin et al., 2009), which may be due to the P demand of rapidly-dividing cancer cells. Additional research on the GRH in the context of cancer should test the effects of elemental ratios in diet on tumor dynamics, which is especially relevant given the many effects of global change on food quality.

Stoichiometric effects on nutrition and related health outcomes may be most severe among low-income populations (within or among nations) because people in more affluent populations have a greater ability to make dietary choices. For populations obtaining the majority of their nutrients from

staple crops, increasing C:element ratios are more likely to lead to diseases stemming from malnutrition. Similarly, as lowincome populations continue to increase their consumption of processed foods (Monteiro et al., 2013), the stoichiometric mismatch between calories and micronutrients is likely to increase (Chopra et al., 2002).

#### Dietary Stoichiometry and Mental Health

The inefficiency of current medical treatments for common mental illnesses such as major depression, along with emerging evidence from longitudinal datasets, has led to a renewed interest in the field of nutritional psychiatry (Sarris et al., 2015). Recent meta-analyses have shown clear links between poor diet and cognitive impairment, dementia, Alzheimer's disease, and especially depression (Psaltopoulou et al., 2013; Lai et al., 2014; Jacka et al., 2017). Adults with mood disorders have been shown to increase functioning when consuming diets high in micronutrients (Davison and Kaplan, 2012; Jacka et al., 2017), especially Mg (Black et al., 2015; Tarleton and Littenberg, 2015). It is fair to assume that decreased nutritional quality, such as increased C:Mg in staple crops, is likely to have pronounced effects on mental health. Therefore, any stoichiometric shifts in diet quality that are relevant to physical health outcomes should also be evaluated in the context of mental health.

#### Microbiome Mediation of Diet-Health Links

The functional composition of the human gut microbiome is inextricably linked with the physiological processes governing many aspects of health. Microbial community composition within the gut has been linked to a broad array of non-infectious diseases, including cancer (Schwabe and Jobin, 2013), obesity (Mathur and Barlow, 2015), diabetes (Barlow et al., 2015), allergies (McKenzie et al., 2017), depression and anxiety disorders (Foster and McVey Neufeld, 2013), autism (Mulle et al., 2013), and a suite of neurological disorders (Tremlett et al., 2017).

Resource stoichiometry impacts the composition and function of microbial communities across a broad suite of environments (Cherif and Loreau, 2007; Hibbing et al., 2010; Larsen et al., 2019), but this has not yet been explored in the human gut. Both dietary fiber intake (Hamaker and Tuncil, 2014) and dietary N (Holmes et al., 2017) have clear effects on gut microbial community composition and function, suggesting that dietary C:element ratios play a role in shaping the microbiome and its effects on health. Additionally, the GRH predicts that fastergrowing microbes should dominate in environments with greater relative P availability (Elser et al., 2003), but this has not been tested in the context of human dietary P intake and the microbiome. Future studies are necessary to determine the extent to which elemental ratios in human diet impact the functional composition of the microbiome and how these effects propagate to human health.

#### Environmental Stoichiometry and Non-infectious Disease

Elemental ratios also link non-infectious disease with global change through ecological processes that occur in terrestrial and aquatic ecosystems (**Figure 3**). Stoichiometric effects on the productivity and biodiversity of these systems impact human health through HABs, pollen production, and toxin exposure.

#### Harmful Algal Blooms

As previously discussed, elemental ratios in water (especially N:P) can drive both the abundance of phytoplankton involved in HABs and the amount of toxins they produce (Van de Waal et al., 2009; Anderson et al., 2012). Toxins from HABs can have both direct and indirect impacts on human health, including endocrine, neurological, and digestive effects (Paerl and Otten, 2013). Depending on the causative organism, exposure to HABs toxins can cause nausea, abdominal pain, diarrhea, respiratory distress, chills, fever, memory loss, seizures, paralysis, and coma (Backer, 1995).

#### Pollen Production

Elemental ratios in the environment are important drivers of plant productivity (as discussed previously), and stoichiometric effects on pollen production may impact human health through both allergies and effects on pollinator populations that are essential to the production of some food crops. Soil N and P and atmospheric CO<sup>2</sup> contribute to pollen production rate and pollen grain size (Lau and Stephenson, 1993, 1994; Lau et al., 1995; Perez-Moreno and Read, 2001), which may exacerbate pollen allergies. Although pollen production has not been specifically linked to elemental ratios in the environment, the frequency of co-limitation in terrestrial ecosystems suggests that stoichiometric effects on plant productivity may underlie pollen production. Stoichiometric effects on biodiversity loss may also contribute to pollen allergies; biodiversity loss has been linked to increased prevalence of allergies and inflammatory diseases in urban settings (Hanski et al., 2012).

Biogeochemical shifts can also impact pollen quality, with implications for pollinator health and food crop production. For example, elevated atmospheric CO<sup>2</sup> led to increased pollen C:N in both historical and experimental analyses, which corresponds to a decline in pollen protein content (Ziska et al., 2016). Pollen protein is essential to bee diet, so this stoichiometric shift in pollen quality may have negative effects on bee health and the pollination of food crops that are integral to human diet.

#### Toxin Exposure

While a large portion of ES literature focuses on ratios of C, N, and P, the ratios of many other elements may impact human health through exposure to toxins. For example, selenium (Se) concentrations in water impact the amount of Hg transferred through food webs (e.g., Sørmo et al., 2011; Walters et al., 2015), thus mediating human exposure to Hg via fish consumption. Following Hg consumption, the amount of Se and Hg circulating through and released from the body also depends on their relative concentrations (Drasch et al., 1996), suggesting that Hg:Se ratios mediate the toxic effects of Hg.

The health effects of arsenic (As), another toxic element that is globally prevalent and linked to a wide range of noninfectious diseases (Amini et al., 2008; Naujokas et al., 2013), are also mediated by its relative balance with other elements. For example, As (in the form of arsenate) has the same shape as phosphate, a molecule essential for all life (Finnegan and Chen, 2012). Consequently, As can be taken into cells in place of P and decouple processes that require high amounts of P, like the formation of adenosine triphosphate (ATP, the main currency of energy in humans). The toxic effects of As are greater when P is low since increased demand for P leads to greater As acquisition (Rodriguez Castro et al., 2015). Arsenic uptake by organisms (and potential exposure to humans) is further complicated by the stoichiometry of additional elements. For example, As uptake by freshwater microbial communities is controlled by N:P ratios (rather than P concentration alone, MacNeill, 2019), and Si:As ratios in soil determine the degree of As accumulation in rice (Zhang et al., 2017). Overall, these examples demonstrate that complex stoichiometric effects can mediate the human health outcomes associated with exposure to toxic elements.

### INFECTIOUS DISEASE

Anthropogenic changes to N and P cycles are associated with increased prevalence of many infectious diseases in human hosts (McKenzie and Townsend, 2007; Johnson et al., 2010). Although previous research has used stoichiometric theory to link biogeochemical cycles with parasite and pathogen infection in a range of non-human hosts (Aalto et al., 2015; Sanders and Taylor, 2018), this topic has not been explored in the context of human health. We outline stoichiometric effects linking infectious disease in humans with global change at two scales: through the effects of nutrition on parasite-host interactions within the human host and through environmental processes outside the human host (**Figure 4**). Hereafter, we use "parasite" to refer to a broad suite of pathogenic organisms (e.g., viruses, protozoa, bacteria, and helminths) that infect human and nonhuman hosts.

### Nutritional Stoichiometry Shapes Parasite-Host Interactions and Disease Development

The effects of human nutrition on infectious disease have received extensive consideration (Field et al., 2002; Cunningham-Rundles et al., 2005), though not yet from a stoichiometric perspective. Within a human host, we can view parasites, the nonpathogenic microbiome, and human immune cells as consumers that require specific resource quality to function optimally (Smith and Holt, 1996). The collective function of the consumers in this scenario determines the symptomatic response observed as infectious disease. Ecological stoichiometry provides a common currency to link elemental ratios in human diet with nutritional impacts on infectious disease development through several mechanisms (**Figure 4**).

#### Effects of Host Nutrition on Parasite Function

Parasites are ecological consumers whose survival, growth, and reproduction are limited by the quality of resources provided by their hosts (Smith, 2007). From a stoichiometric perspective, the rates of these physiological processes may depend on elemental ratios originating from host diet. Several core concepts from stoichiometric theory may contribute to our understanding of links between parasite function and human diet: the GRH, threshold elemental ratios, and stoichiometricallymediated interactions among species.

The GRH provides a promising avenue to link human diet quality with parasite function. The rapid growth rates of parasites relative to their hosts suggests that parasitic taxa should be limited by C:P or N:P within hosts. Experimental evidence from other disease systems supports this prediction: viral production increased in algal hosts grown at low C:P (Clasen and Elser, 2007), and bacterial infection rates increased when zooplankton hosts were fed a low C:P diet (Frost et al., 2008). While the GRH has not been tested for parasites infecting human hosts, understanding links between dietary C:P or N:P and parasite growth are especially relevant given the overconsumption of P in many human diets (Razzaque, 2011).

However, parasites vary widely in elemental ratios across species (Paseka and Grunberg, 2019), so the limitation of parasite growth rate by P is unlikely to be universal. For example, fungal parasites infecting cyanobacteria had the highest replication rate when hosts were grown at high N:P, suggesting that N is more limiting to these parasites than P (Frenken et al., 2017). The concept of threshold elemental ratios (TERs) provides a framework for predicting how nutritionally-limited consumers respond to variation in resource quality (Boersma and Elser, 2006). In contrast to the GRH, TERs are not specific to P but could instead be applied to any elemental ratio. A consumer's TER represents an "ideal" resource ratio where consumer growth is equally limited by both elements under consideration. The cost of processing excess elements is predicted to create a humpshaped response of growth rate to diet quality around the TER (Boersma and Elser, 2006). TERs have not been explored as a means to study stoichiometric limitation of parasite growth, but this concept could aid in making specific predictions about how parasites will respond to shifts in human diet.

Finally, ES may also provide a framework for predicting how human diet influences the outcome of interspecific interactions among parasites and the non-pathogenic microbiome. These organisms interact through complex mechanisms, and the composition and function of the microbiome can limit or promote infectious disease (Baümler and Sperandio, 2016). As described earlier, we lack critical data on potential effects of elemental ratios in host diet on the structure and function of the human microbiome and its effects on disease.

#### Nutritional Controls on the Immune System

There is ample evidence that human diet regulates infectious disease through immune function (Field et al., 2002; Cunningham-Rundles et al., 2005), though this has not been studied in the context of ES. From a stoichiometric perspective, human immune cells can be viewed as consumers that function optimally at specific elemental ratios. For this reason, applying the concept of TERs to immune cells (and comparing host immune TERs to parasite TERs) could help to draw links between human diet quality, optimal immune function, and defense against infectious disease.

### Environmental Stoichiometry Mediates Human Infection Risk

Elemental ratios also link infectious disease with global change through processes occurring in terrestrial and aquatic ecosystems that mediate human exposure to environmentally transmitted parasites (**Figure 4**). Stoichiometric effects on ecosystem productivity shape the population dynamics of consumers,

development of infectious disease (3).

stoichiometry may impact the function of parasites, the microbiome, the immune system, and their interactions. Collectively, these functions determine the

including parasites that require development outside the human host. We use "environmental parasites" to refer to species that replicate during free-living stages in the environment and those that replicate within vector or intermediate host species.

#### Population Dynamics of Parasites, Intermediate Hosts, and Vectors Impact Human Risk of Exposure

Stoichiometric effects on ecosystem productivity will impact parasites that have free-living stages in the environment (i.e., those that live independently of a host or vector). For example, Vibrio cholerae, the causative agent of cholera, lives as a heterotrophic bacterium in aquatic environments when outside the human host and is therefore impacted by environmental resource availability (Cottingham et al., 2007). Associations between algal blooms and cholera outbreaks have been observed for decades (Epstein, 1993; Colwell, 1996), yet it remains unclear what mechanisms link elemental availability, primary production, and V. cholerae population dynamics. Elemental ratios in aquatic ecosystems may regulate V. cholerae directly by providing resources required for growth or indirectly by promoting the growth of attachment surfaces (phytoplankton, zooplankton, and macrophytes) on which the bacteria aggregate during blooms (Cottingham et al., 2007). Elemental ratios constrain the growth of nonpathogenic bacteria in aquatic ecosystems and mediate algal bloom dynamics (Sterner and Elser, 2002), suggesting that stoichiometric effects may impact the population dynamics of V. cholerae and other human parasites with free-living environmental stages.

The effects of elemental ratios on ecosystem productivity can also constrain the distribution and population dynamics of vectors or intermediate hosts, thereby limiting the populations of their associated parasite species. For example, elemental ratios in aquatic ecosystems impact the distribution and population dynamics of larval mosquitoes (Murrell et al., 2011; Yee et al., 2015), which in turn determine the populations of adult mosquitoes that transmit a wide array of viruses and protozoan parasites.

In addition to the number of vectors or hosts in the environment, the quality of these organisms as resources for parasites may also be mediated by stoichiometric effects. Like parasites within human hosts, environmental parasites may be nutrient limited within intermediate hosts or vectors. Evidence from several disease systems indicates that shifts in the stoichiometry of the environment can cascade to alter the elemental composition or other aspects of host and vector physiology, thus altering the resources available to parasites and mediating transmission risk (Sanders and Taylor, 2018). For example, mosquitoes that consumed high quality diets as larvae matured into adults that were higher in %N, which corresponded to lower rates of Zika virus infection and transmission potential (Paige et al., 2019). While the mechanism behind this shift is unknown, it may reflect elemental requirements of mosquito immune function. These studies on mosquito vectors are a rare demonstration of the importance of elemental ratios for vector population dynamics, traits, and transmission of a human virus, and this topic represents an important direction for future research in additional disease systems. Like within-host disease dynamics, we suggest that TERs and the GRH are core stoichiometric concepts that should be used to study the effects of elemental ratios on environmental parasites, intermediate hosts, and vectors.

#### Ecological Communities Mediate Parasite Transmission to Human Hosts

In addition to the effects of environmental stoichiometry on the direct interactions between environmental parasites and their intermediate hosts and vectors, broader properties of ecological communities also have important implications for infectious disease. Stoichiometric effects on species richness and composition may mediate human infection risk through the dilution effect, a well-supported, inverse relationship between infectious disease and biodiversity (Keesing et al., 2006; Ostfeld and Keesing, 2012). There is robust support for the dilution effect across a broad range of infectious disease systems, including human diseases (Civitello et al., 2015). The effects of environmental stoichiometry on human infectious disease as mediated by community richness and composition have not been studied, though results of a grassland experiment suggest that Nfertilization can reduce the strength of the dilution effect and lead to greater levels of fungal disease in plant communities (Liu et al., 2017).

Overall, stoichiometric constraints on ecosystem productivity and biodiversity may mediate the rate at which humans encounter infective parasite stages produced in the environment. Future research on how elemental ratios cascade through food webs to impact parasite production and transmission will provide new opportunities to link global change with infectious disease risk.

# DISCUSSION

# Limitations of a Stoichiometric Approach to Human Health

Ecological stoichiometry is a powerful framework for studying the effects of global change on human health because it enhances our understanding of how elements and energy are transferred across scales of biological organization. However, one limitation to the application of ES is that organisms (including humans) cannot assimilate all molecular forms of elemental nutrients. For example, humans require N in complex molecules because we are only able to synthesize half of the amino acids required for physiological function. Such nuance is lost when considering only the C:N ratios (or N content alone) of human diet. Nevertheless, a stoichiometric approach can often complement more traditional considerations of dietary macromolecules. For example, the elevation of C:N content in rice corresponds to decreased concentrations of N-rich B vitamins (Zhu et al., 2018). In this case, stoichiometric and macromolecular perspectives provide complementary insight into links between global change and human nutrition.

However, some biogeochemical effects of global change on human health can be effectively characterized using a single-element approach. For example, nitrate leaching due to over-fertilization has numerous human health implications (Townsend et al., 2003). The human health effects of nitrate in groundwater are largely a function of N alone, not N:P or another stoichiometric interaction. Similarly, the health effects of toxic elements such as Hg are better expressed as Hg concentration in the body than as an C:Hg ratio. However, the ratios of Hg with other elements (such as Se:Hg) influences human health risk through stoichiometric processes both in the environment and within the human body. While a reductionist, stoichiometric approach has limitations, the interdisciplinary, cross-scale power of ES also holds great potential to inform human health research.

#### Stoichiometric Insight for Human Health Under Global Change

We have outlined many potential mechanisms through which elemental ratios link changing biogeochemical cycles with food security and water quality, as well as a wide array of infectious and non-infectious diseases in humans. However, we note that stoichiometrically-explicit data on this topic are generally scarce. We advocate for future research to harness the power of ES by testing some of its core tenets in the context of human health. This review is not an exhaustive exploration of these ideas, the current literature, or potential future directions. Instead, we hope that this work will inspire future research on the application of stoichiometric theory to human health under global change.

#### REFERENCES


#### AUTHOR CONTRIBUTIONS

RP conceived the original idea for this review and coordinated manuscript completion. All authors (RP, AB, KM, AB, and CS) contributed to development and writing. Author order is random after first.

#### FUNDING

Woodstoich 4 was supported by the National Science Foundation DEB-1840408. RP was was supported by the National Socio-Environmental Synthesis Center (SESYNC) under funding received from the National Science Foundation DBI-1639145.

#### ACKNOWLEDGMENTS

This manuscript is a product of Woodstoich 4, and we thank Jim Elser and Michelle Evans-White for organizing a groovy and productive event for early career researchers. Thank you to our three reviewers, whose comments and suggestions led to substantial improvements on this manuscript. We thank Zoe Cardon for her sage advice and enthusiasm throughout Woodstoich. Thank you to artists at the Noun Project (https:// thenounproject.com; **Supplementary Information 1**) for icons used in tables and figures and to Steven Ross Davidson for help assembling figures.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2019.00378/full#supplementary-material


of a toxic cyanobacterium and its fungal parasite. Front. Microbiol. 8, 1–11. doi: 10.3389/fmicb.2017.01015


food web, Grand Canyon, USA. Environ. Toxicol. Chem. 34, 2385–2394. doi: 10.1002/etc.3077


**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Paseka, Bratt, MacNeill, Burian and See. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.