Nitrogen loading leads to increased carbon accretion in both invaded and uninvaded coastal wetlands

Gaining a better understanding of carbon (C) dynamics across the terrestrial and aquatic landscapes has become a major research initiative in ecosystem ecology. Wetlands store a large portion of the global soil C, but are also highly dynamic ecosystems in terms of hydrology and N cycling, and are one of the most invaded habitats worldwide. The interactions between these factors are likely to determine wetland C cycling, and specifically C accretion rates. We investigated these interactions using MONDRIAN, an individualbased model simulating plant growth and competition and linking these processes to N and C cycling. We simulated the effects of different levels of (1) N loading, (2) hydroperiod, and (3) plant community (natives only vs. invasion scenarios) and their interactions on C accretion outcomes in freshwater coastal wetlands of the Great Lakes region of North America. Results showed that N loading contributed to substantial rates of C accretion by increasing NPP (net primary productivity). By mediating anaerobic conditions and slowing decomposition, hydroperiod also exerted considerable control on C accretion. Invasion success occurred with higher N loading and contributed to higher NPP, while also interacting with hydroperiod via ecosysteminternal N cycling. Invasion success by both Typha × glauca and Phragmites australis showed a strong nonlinear relationship with N loading in which an invasion threshold occurred at moderate N inputs. This threshold was in turn influenced by duration of flooding, which reduced invasion success for P. australis but not for T. × glauca. The greatest simulated C accretion rates occurred in wetlands invaded by P. australis at the highest N loading in constant anaerobic conditions. These model results suggest that while plant invasion may increase C storage in freshwater coastal wetlands, increased plant productivity (both native and invasive) due to increased N loading is the main driver of increased


IntroductIon
Gaining a better understanding of carbon (C) dynamics across the landscape has become one of the major research initiatives in ecosystem ecology due to the critical role of C in global climate change (Shaver et al. 2000, Davidson andJanssens 2006). Wetlands are key habitats that regulate C and nutrient flows through the landscape and loss to the atmosphere because of their position at the interface between terrestrial and aquatic zones (McClain et al. 2003). As a result, when wetlands are flooded for extended periods, anaerobic conditions can reduce decomposition September 2016 v Volume 7(9) v Article e01459 2 v www.esajournals.org MARTINA ET AL. rates causing C fixed in high-productivity wetlands to have a long residence time as inundated litter and soil (Holden 2005, Reddy andDelaune 2008). This combination of high productivity with low decomposition rates has made inland and coastal wetlands significant reservoirs of C, with freshwater wetlands storing 20-25% of the world's soil carbon while occupying only 4-6% of the global land surface (Mitra et al. 2005, Hopkinson et al. 2012.
While research on these processes has mostly focused on high latitudes (Roulet 2000), temperate wetlands are also important sinks of C (Euliss et al. 2006) and furthermore are often more susceptible to anthropogenic influences. Wetland C dynamics as part of inland C budgets have also been identified as a key point of uncertainty (Regnier et al. 2013). It is therefore important to understand how the main drivers of C accretion interact in these more anthropogenic-influenced landscapes, such as those of the Laurentian Great Lakes region of North America. We use term "C accretion" or "C accretion rate" as the accumulation (on a yearly basis) of the sum of the major pools of organic C, including living biomass, litter, muck (a highly organic, sapric soil surface layer), and mineral soil organic matter (MSOM; organic matter that occurs within mineraldominated soil layers). We find accretion to be a more useful term than C sequestration, which usually refers to the long-term storage of C in resistant soil pools (Lal 2004) because large pools of litter and muck layers can accumulate and may not be recalcitrant, but simply inundated, thus severely slowing decomposition (González-Alcaraz et al. 2012. Therefore, understanding how these four different C pools are affected by biotic and abiotic factors gives us a more mechanistic understanding of wetland C dynamics.
Elevated nitrogen (N) inputs into wetlands in the Great Lakes region in recent decades have likely resulted in significant increases in community NPP (net primary productivity). Before strong anthropogenic influence, lakeshores in this region mainly comprised low nutrient systems where native vegetation was adapted to oligotrophic conditions. Due to the widespread use of agricultural fertilizer, combustion of fossil fuels, and cultivation of N-fixing legumes (Vitousek et al. 1997, Holland et al. 2005, Han et al. 2009), N influx to wetlands through atmospheric deposition, groundwater flow, and surface water runoff has significantly increased relative to preindustrial conditions (Mitsch 1992, Morrice et al. 2004, Galloway et al. 2008. Along with causing an increase in community NPP, the increase in N loading has likely altered community composition through plant invasions into nutrient-rich systems (Farrer and Goldberg 2009, Tuchman et al. 2009, Currie et al. 2014. Wetlands in general are often highly invaded ecosystems because their placement on the landscape makes them sinks of water runoff, nutrients, and plant propagules; they are also prone to disturbance (including flooding) that facilitates invasive plant success (Davis et al. 2000, Zedler and Kercher 2004, Eschtruth and Battles 2009. In our study region, invasions by aggressive plant species, such as Phragmites australis (Cav.) Steud. and Typha × glauca Godr. (hereafter Phragmites and Typha, respectively), have drastically changed the plant community composition in many inland and coastal wetlands (Zedler and Kercher 2004). Phragmites and Typha are both large-stature clonal graminoids that positively respond to N enrichment (Woo andZedler 2002, Rickey andAnderson 2004).
As is well known, hydroperiod (the degree and duration of flooding) strongly controls wetland C accretion by mediating aerobic or anaerobic conditions. Hydroperiod has also been strongly influenced by humans over the past century. Wetlands were drained in the Midwest starting in the nineteenth century to accommodate farming in wetlands across the region (Mitsch and Gosselink 2000). Humans continue to alter wetland hydroperiod directly by diking and draining and indirectly through upstream hydrologic manipulation and climate change (Mitsch andGosselink 2000, Angel andKunkel 2010). Climate change is predicted to further alter wetland hydroperiod in our study region (Hartmann 1990). Therefore, it is critical to understand how hydroperiod interacts with other drivers of wetland C accretion across this region.
Wetland plant invasions can result in many negative consequences to local biodiversity and habitat quality (Spyreas et al. 2010. Less well studied in wetlands, invasion may alter the manner in which hydroperiod affects rates of C accretion. Introduced wetland v www.esajournals.org species can drastically differ from native species in a number of key plant traits, such as maximum height, tissue chemistry, and growth rate (Chapin and Eviner 2003, Bourgeau-Chavez et al. 2012, Currie et al. 2014). These differences in plant traits can feedback to enhance C influx into wetlands if productivity of the invasive species is greater than that of natives. Invasion can further influence C accretion rates through altered decomposition if invasive species differ from natives in litter chemistry (Chapin and Eviner 2003, Ehrenfeld 2003, Eviner 2004. Empirical results on the consequences of plant invasion on C accretion have been mixed and attributed to differences in plant traits. Cheng et al. (2006) showed that when Spartina alterniflora invaded native sedge tidal wetlands in China, the greater rooting depth of the invasive greatly increased organic C in the top 60 cm of soil. Conversely, even though Agropyron cristatum invasion into native grasslands doubled belowground productivity, there was no increase in soil C content because the invader's belowground biomass was more labile than natives (Macdougall and Wilson 2011). These examples illustrate the importance of not only knowing how traits differ among natives and invaders, but also which ecosystem C pools (e.g., aboveground litter or root litter and other detrital pools) are affected by invasion over time. It is also key to understand how dynamics in ecosystem C pools are likely to interact with wetland hydroperiod and N loading.
Here, we examine the manner in which the well-known relationship between flooding and wetland C accretion is affected by N loading and large-plant invasions in Great Lakes coastal wetlands. We expect strong interactions between hydroperiod and N loading because the lowering of decomposition rates associated with anaerobic conditions can slow the cycling of N held in undecomposed organic matter (Scholz 2011). We suggest this slowing of N cycling could decrease plant productivity and subsequent C accretion under oligotrophic conditions, but may have minimal effects under high N loading where ample influx of N is available for plant growth. The C fixed by the macrophyte plant community provides a major influx of autochthonous C in most wetlands (Wetzel 2006). Understanding how plant community dynamics, including the interaction with N loading and hydroperiod, affect C accretion will enable us to better understand the mechanisms behind wetland C sequestration, as well as where and when it is likely to occur. Elevated N loading increases NPP (LeBauer and Treseder 2008) and thus likely has a direct and large impact on wetland C accretion. Indirectly, N loading is known to influence invasion rates in wetlands Kercher 2004, Currie et al. 2014) and this interaction may further alter the outcome of increased N loading on wetland C dynamics. Besides their interactions, we are also interested in exploring the relative magnitudes of effect of N loading, hydroperiod, and plant invasion on C accretion. As described above, empirical evidence exists for the influence of each driver and some of the underlying mechanisms, but their relative importance is less well known.
The interactive effects and feedbacks among N loading, plant invasions, and hydroperiod are complex, making them difficult to disentangle through empirical work alone. Detailed, multifactor empirical data would be needed over multiple points in time (Fukami 2010). Mechanistic models provide an alternative tool to understand the complex and likely nonlinear relationships among these drivers. Models allow us to manipulate hydrology, N loading, and invasive plant traits in a precise and controlled manner not possible in the field. In this study, we present and apply an enhanced version of a mechanistic community-ecosystem model, MONDRIAN (Modes of Nonlinear Dynamics, Resource Interactions, And Nutrient cycling; Currie et al. 2014), that incorporates dynamic water levels and the anaerobic slowing of decomposition in submerged litter, muck, and MSOM. The effects of hydroperiod on decomposition and C accretion are fully integrated with N cycling and species competition in the model, allowing us to examine how N loading, plant invasions, and hydroperiod interact with control rates of C accretion in Great Lakes coastal wetlands. We expect our general results to be applicable to a variety of temperate wetlands.
We asked the following specific questions: 1. How does variation in N loading affect wetland community NPP and rates of C accretion, as mediated by aerobic and anaerobic conditions that affect both C mineralization and N cycling feedbacks?

MaterIals and Methods
MONDRIAN is an individual-based model that spans several major levels of ecological organization, from individual plant physiology to ecosystem function, and is formulated through a set of algorithms in an object-oriented programming language (Visual Basic.Net). MONDRIAN was fully described by Currie et al. (2014) so we start with only a brief description of the original model, followed by more detail on additions for the research described here. We used MONDRIAN to model a 52.5 × 52.5 cm area consisting of 49 grid cells, each 7.5 × 7.5 cm in area. Plant competition takes place in these grid cells, along with most C and N cycling, but plants can grow clonally across grid cells. At the individual level, MONDRIAN simulates up to thousands of ramets per square meter, modeling both internal source-sink C and N translocation within each plant and explicit spatial size-symmetric competition for available N, which leads to heterogeneous N availability. At the population level, plants can produce new ramets from rhizomes if they have enough C and N to create a daughter ramet; this C and N demand connects resource competition among individuals to population dynamics in a heterogeneous environment. Mortality can lead to the loss (and conversion to litter) of individual ramets or whole genets. Emergent community dynamics include species coexistence and competitive exclusion, biodiversity changes over time, and both successful and unsuccessful plant invasions. These dynamics arise from competition among neighboring plant individuals for resources and population-level expansion and mortality among up to four species simulated together. Ecosystem processes are a function of individual, population, and community-level dynamics and arise from plant growth, population fluctuations, and community composition shifts, along with externally driven N inputs. For a further description of C and N cycling in MONDRIAN, including controls on decomposition, decomposition feedbacks on N mineralization, plant growth and uptake of N, and the emergent feedback linking increased invader success to ramped up ecosystem-internal N cycling mediated by litter see Currie et al. (2014).
Although the basic features of MONDRIAN have been previously described (Currie et al. 2014), the model is undergoing further development. Processes in MONDRIAN were augmented in two important ways to conduct the research described here. First, daily fluctuation in water level and the effects of anaerobic conditions on decomposition rates were added. Second, light competition among individual stems was added to the existing N competition.
User-controlled daily fluctuation in water level was added to the model and integrated with model processes to affect ecosystem C and N cycling. A new "muck" pool was added to the model ( Fig. 1) to represent the accumulation of highly organic, sapric soil that often develops in productive wetlands where litter decomposition is slowed by inundation. In the model, muck sits physically below the litter layer and atop the surface of wetland sediments; organic matter within the sediments is part of the mineral soil organic matter (MSOM) pool. As muck mass accumulates, its upper surface rises vertically in MONDRIAN, at a rate that depends on its bulk density. On the timescale of years to decades, plant bases rise with the surface of the muck so that aboveground stems remain above the muck surface and rhizomes can grow within the muck. This vertical accumulation of muck can result in "terrestrialization" if the muck layer rises above the water surface. In anaerobic conditions at high N loading (see below for treatment description), the depth of the muck reached 10-15 cm, which was below the water surface; thus, "terrestrialization" did not occur during these simulation runs. In previous research on wetland plant invasions and N cycling with MONDRIAN (Currie et al. 2014), the model did not include a muck pool; water level and the development of anaerobic conditions under inundation were not explicitly addressed. These features have been added to the v www.esajournals.org MARTINA ET AL. model for the present analysis. To include a delay in the onset of anaerobic conditions following inundation (Reddy and Delaune 2008), a 5-d trailing average in water level is calculated. All detrital pools (or proportions thereof), including aboveand belowground litter, muck, and MSOM pools lying below the level of the 5-d trailing average in water level are considered anaerobic.
At the ecosystem level, C and N flow starts in living tissue where C is fixed through photosynthesis and N uptake occurs in the roots. This C and N enter the litter pools (above-or belowground) after tissue senescence (Fig. 1). These fluxes (living tissue C and N flux to litter C and N pools) are not directly affected by anaerobic conditions, although anaerobic conditions limit the availability of inorganic N (due to decreased decomposition), which limits plant growth; hydroperiod thus has a realistic effect on plant growth via N mineralization. During decomposition, a portion of C and N in the aboveground litter pool is transferred to the muck pool. Decomposing belowground litter likewise is transferred to muck, MSOM, or a combination of the two depending on its depth relative to the muck-MSOM interface. A small portion of the muck pool also transfers to the MSOM pool each day ( Fig. 1), representing bioturbation and particle eluviation. For each day that any detrital pool (or portion thereof) is anaerobic, decomposition in that pool is slowed by a multiplicative modifier (0.2, Reddy and Delaune 2008). Thus, floods enhance C and N accretion in detritus, while slowing the release of both C and N from detrital pools via mineralization (Fig. 1). Taken together, MONDRIAN incorporates hydroperiod feedbacks on soil moisture and productivity; anaerobic conditions allow muck to accumulate, which raises plant level, which can decrease anaerobic conditions if terrestrialization occurs (effectively decreasing soil moisture) resulting in increased productivity by increased N mineralization.  Currie et al. 2014, with muck C and N pools added). Plant pools of C and N are specific to individuals; grid cell pools are specific to each spatially explicit cell (7.5 × 7.5 cm) within the modeled area, each containing numerous individual plants and allowing heterogeneous nutrient depletion and light availability; the regional nutrient pool is a single pool across the entire modeled area, akin to a pool of standing water. Asterisks indicate C and N fluxes that are influenced by hydroperiod. C flows are internally simulated in g·C·m −2 ·d −1 , N flows in g·N·m −2 ·d −1 .

MARTINA ET AL.
The enhanced MONDRIAN model also now includes light competition by calculating shading from neighboring plants and its effect on the growth rate of each individual. Because light is especially limiting in highly productive eutrophic wetlands (Güsewell and Edwards 1999), this enhancement allowed us to confidently simulate community interactions under higher levels of N input than those used by Currie et al. (2014). Testing of MONDRIAN confirmed that shading effects on growth rates became particularly important at higher rates of N input and NPP. Light availability is calculated in 10-cm vertical segments separately in each spatially explicit grid cell (7.5 × 7.5 cm). The shading calculation uses species-specific light extinction curves applied to the plant biomass (stem + foliar) present on a daily basis, by species, in each vertical segment of each grid cell. Plant height is determined using species-specific biomass-height allometric equations obtained from our own field data ( Table 1). The effect of shading on each individual stem is simplified by the light environment at a fixed proportion of its height, a species-specific parameter that represents the typical vertical distribution of photosynthetic tissue for the species (J. Knops and H. Hager, unpublished data; Table 1). Growth rate is then scaled back using a Michaelis-Menten equation of relative growth rate as a function of light availability based on species-specific data on photosynthetic rate as a function of irradiance (Knops and Hager, unpublished data).
We parameterized MONDRIAN using realistic species parameters, rather than hypothetical species traits as used by Currie et al. (2014). Three native species (Eleocharis palustris (L.), Juncus balticus (Willd.), and Schoenoplectus acutus (Bigelow) A. Love & D. Love) and two invasive species (Phragmites and Typha) were parameterized and used in the in silico experiments we report here. These native and invasive species commonly occur in Great Lakes coastal wetlands. All plant species were parameterized using multiple values found in both the literature and our own unpublished data collected from the Great Lakes region (Table 2). If multiple values were found for a species within the Great Lakes region, an average was used for that trait.
We conducted sets of contrasting simulation runs, each lasting 45 yr, with a fully factorial design of N loading, hydroperiod, and plant community scenarios. The three plant community scenarios were natives only (three-species community), an established native community invaded by Typha, and an established native community invaded by Phragmites. In all community scenarios, natives were randomly distributed into the modeling area in four cohorts of 65 genets in years 1, 3, 5, and 7 and had stabilized in terms of NPP and density by year 15. In the invasion scenarios, two cohorts of 15 invader genets each were introduced at random locations in years 15 and 20. After initial introduction of a species, one individual genet was randomly added to the modeling area per year to represent natural colonization. The seven levels of N input ranged from 0.86 to 30 g·N·m −2 ·yr −1 and were constant throughout each model run. The lowest N input represents present-day rain-fed N deposition in northern Michigan (wet + dry inorganic N deposition plus atmospheric organic N deposition) (Neff et al. 2002, NADP 2009), and the highest N Notes: Canopy architecture is modeled using a polynomial-shading curve (% full light = Ax 2 + Bx + C) based on biomass. Native species and Typha used the same equation for light extinction (with different height ranges) because of their similar leaf-stem architecture. Phragmites used a distinct light extinction equation because of its dramatically different leaf-stem architecture compared with the other parameterized species. The weighted photosynthetic tissue height parameter was used to model individual responses to shading and is expressed as a proportion of the height of each individual.
† Biomass-height allometric equation, where height is in meters and biomass is in g dry mass: height = A × biomass B .
The three hydrologic regimes were as follows: (1) always aerobic (water level constant at 15 cm below the MSOM surface, i.e., −15 cm); (2) always anaerobic (constant water level 30 cm above the MSOM surface, i.e., +30 cm); and (3) sinusoidal fluctuation in the water level of ±15 cm about the MSOM surface with an annual period. In the fluctuating scenario, the wet period occurred in spring and early summer, while the dry period occurred in late summer and fall similar to fluctuations seen in Great Lakes coastal wetlands. Simulations of smaller water-level fluctuations (±5 cm) showed comparable results to the ±15 cm scenario and are not presented here for simplicity. We selected these three hydroperiod scenarios to represent possible water levels found in coastal wetlands in Michigan. While wetlands closer to the coast likely fluctuate similar to our ±15 cm scenario, wetlands slightly further from the coast can have a less fluctuating hydroperiod over a year and/or can go through periods of flooding or drying, comparable to our anaerobic and aerobic hydroperiod treatment endpoints (Wilcox et al. 2002). It should be also noted that in MONDRIAN water level has no direct effect on plant survival, although this should not affect the realism of our results because water levels above 30 cm are usually needed to negatively affect growth of established vegetation (Waters and Shay 1990, 1992, van der Valk 2000. The factorial design of three plant community scenarios, seven levels of N loading, and three hydrologic regimes produced 63 combinations of model settings. Each combination was replicated three times (with stochastic differences both in initial plant distributions and spatial movements during clonal reproduction) for a total of 189 model runs (model run = one 45-yr simulation). In all model runs, our key dependent variables stabilized by 30-40 yr and so for all statistical tests and figures, the average of the last 5 yr (years 41-45) of each model run was used. Total NPP, invader proportion of community NPP, and C accretion were analyzed as dependent variables using a three-factor ANOVA with community scenario, N loading, and hydroperiod as main factors and all two-way and threeway interactions. Magnitude of effect differences among the three main drivers of C accretion were determined by comparing difference in means and percentage change in the most extreme treatment levels and by calculating η 2 for each main driver. η 2 is the proportion of total variance attributed to an effect and was calculated as the sum of squares of an effect divided by the total sum of squares (similar to a partial R 2 ).

results
The range of NPP (aboveground and total), litter mass, and C accretion rates produced in our in silico experiments were comparable to values found in the literature, as well as our field data from temperate wetlands in Michigan (Table 3).

Effects of N loading on NPP and plant invasions
Total community NPP was highly sensitive to the amount of N loading as expected (Fig. 2,  Table 4). NPP showed a saturating response to Table 2. Species-specific model parameters for the three native species and two invasive species used in simulating Great Lakes coastal wetlands. K-constant of litter refers to the first-order decomposition constant (K) for litter (Currie et al. 2014  increasing N inputs beginning at ~15 g·N·m −2 ·yr −1 , resulting from light competition and shading, both within and between plant species, in the model. While total community NPP changed smoothly along the N gradient regardless of invasion scenarios, NPP of both invasive species exhibited a steep, nonlinear threshold in invasion success (Fig. 3). At low N loading (<5 g·N·m −2 ·yr −1 ), neither Phragmites nor Typha dominated over established native communities (had greater than 50% percentage of NPP), although both could persist in a native community. At high N loading (≥15 g·N·m −2 ·yr −1 ), each invasive species was almost or completely dominant (90-100% of NPP).
Generally, successful Phragmites invasion resulted in an increase in total community NPP compared with uninvaded native communities, with a few exceptions. For example, under aerobic conditions at an N input of 9 g·N·m −2 ·yr −1 , Phragmites was moderately invasive but had no effect on the total community NPP relative to the uninvaded native community at the same N level (compare Figs. 2 and 3). However, at an N loading of 15 g·N·m −2 ·yr −1 , Phragmites increased dominance to 100% and caused a substantial increase in total community NPP relative to both the uninvaded native community and the Typhainvaded community at the same N level ( Fig. 2; species × N loading interaction; Table 4). Overall, successful Typha invasion did not significantly affect total community NPP compared with the uninvaded native community. This difference in influence between Typha and Phragmites on total community NPP was likely due in part to differences in maximum size and canopy architecture as represented in MONDRIAN.

C accretion rates in simulated wetlands
Ecosystem C accretion rates depended significantly on all three factors tested (plant community scenario, N input, and hydroperiod), as well as on all two-way and three-way interactions among these factors (Table 4). We focus first on the effects of N loading and plant invasions before considering the direct and indirect effects of hydroperiod.
A key finding in our results was that C accretion rates increased with N loading regardless of other treatments (Fig. 4) and that N loading provided the strongest overall control on C accretion. While invasion success was also driven by N loading, its effects on C accretion were both smaller and more complex (nonlinear) than those of N loading alone. Invasion of both Phragmites and Typha increased C accretion rates over that of the native community at low N loading, but not at medium N loading (Fig. 3). At high N loading, invasion of both Phragmites and Typha increased C accretion rates over that of natives in aerobic conditions but not when conditions were seasonally or continually anaerobic.
Differing community compositions, including the identity of the invasive species, drove ecosystem C accretion through changes in different organic C pools (live tissue, litter, muck, and MSOM; Table 4). For example, the Typha-invaded community caused the greatest change in the summed above-and belowground litter C pools, which was up to 3× that of the native-only community and Phragmites invasion scenarios (Fig. 5). Phragmites, on the other hand, increased C accretion at the high end of the N gradient through a combination of higher live biomass and muck C accretion (especially under anaerobic conditions, Fig. 5). In contrast to this high muck C accretion under Phragmites, Typha invasions had the lowest muck C accretion rates across gradients of N loading and hydroperiod. MSOM C accretion, which increased with N loading, was more similar among communities, although was generally higher in Phragmites invasion scenarios, especially  Vaccaro et al. (2009), and E. Farrer, D. Goldberg, and K. Elgersma (unpublished data) C accretion rates (g·C·m −2 ·yr −1 ) −7.90-573 20.0-500 Rabenhorst (1995), Reddy and Delaune (2008), and Bernal and Mitsch (2012) Notes: NPP, net primary productivity. References cited in the table refer to observed values given for comparison.
September 2016 v Volume 7(9) v Article e01459 9 v www.esajournals.org at high N loading (Fig. 5). Compared with the native-only and Typha-invaded communities at high N loading, the greater rates of C accretion produced by Phragmites-invaded communities in living biomass, muck, and MSOM were consistent with the higher Phragmites NPP (Fig. 2).

Interactive effects of hydroperiod
Hydroperiod had little effect on community total NPP (Fig. 2), but had both direct and indirect effects on ecosystem C accretion through the slowing of decomposition and N mineralization under anaerobic conditions (Fig. 4). Interestingly, the N loading value at which the threshold of invasion success occurred depended on species, hydroperiod, and their interaction (all Ps < 0.001; Table 4, Fig. 3). This indicates that ecosystem-level N cycling feedbacks in the model (Currie et al. 2014) were strongly mediated by the hydroperiod and its effects on detrital accretion and decomposition. Under aerobic conditions, Phragmites was able to invade at a lower N loading threshold than Typha (due to greater availability of N mineralized from its litter), while under anaerobic conditions both species had a higher and similar N threshold for invasion (Fig. 3).
Under high N loading, when conditions were seasonally or continually anaerobic, only Phragmites invasion, not Typha invasion, increased C accretion rates compared with native-only community scenarios. In anaerobic conditions, the percentage increase in C accretion rates for Phragmites invasion compared with natives only was greater at low N loading (69%) than at the highest N loading (12%), although the absolute increase was more at highest N loading (highest N loading: 59.3 g·C·m −2 ·yr −1 , lowest N loading: 30.5 g·C·m −2 ·yr −1 ) (Fig. 4). As expected, the simulated rate of C accretion was always greater in constant anaerobic conditions compared with all other hydrologic scenarios (significant main effect of hydroperiod). This is a direct effect of slowed C mineralization under anaerobic conditions.
The effect of hydrologic regime on C accretion rates differed among C pools. Muck C accretion was the most sensitive to anaerobic conditions and was up to four times higher in anaerobic conditions than any other hydrologic regime. Variability in hydrologic regime (variable ± 15 scenario) seemed to lower muck accretion relative to anaerobic conditions (Fig. 5). Conversely, MSOM C accretion rates under variable conditions were more similar to anaerobic conditions, likely due to a significant proportion of MSOM C being inundated in anaerobic and variable hydrologic regimes. N loading had less effect on MSOM C accretion rates in aerobic conditions (Fig. 5). Living biomass and litter C accretion rates were higher in aerobic and anaerobic conditions than the variable hydrology scenario across the N loading gradient (Fig. 5).

Magnitude of effect for the drivers of C accretion
We compared the main drivers of ecosystem C accretion rates using values of the drivers occurring at extremes, but still likely under field conditions (invaded vs. uninvaded communities, N inputs ranging from oligotrophic (atmospheric deposition only) to highly eutrophic (representing agricultural and urban runoff), aerobic vs. anaerobic). The validity of these comparisons is supported by the similarity of the ranges of simulated C accretion rates to those found in the field (Table 3). Increasing N loading from 0.86 to 30 g·N·m −2 ·yr −1 increased C accretion by 958% (mean difference = 290.4 g·C·m −2 ·yr −1 ), and shifting conditions from constant aerobic to anaerobic increased C accretion by 220% (mean difference = 190.2 g·C·m −2 ·yr −1 ). It should be noted that the hydroperiod treatments in this study refer to effects on near-surface soil (top 15 cm), while all sediments below 15 cm were always considered anaerobic in this wetland model. In systems where aerobic conditions can penetrate below this 15 cm depth or where the water level can rise >30 cm above the sediment surface, it is likely that hydroperiod would have a much larger effect size on N cycling and plant survival. Hydrologic regime and N loading also explained a large proportion of variation in C accretion (η 2 = 0.654 and 0.262 and for N loading and hydroperiod, respectively). Invasion by Typha or Phragmites, however, only increases C accretion by 7% (mean difference = 11.0 g·C·m −2 ·yr −1 ) and 15% (mean difference = 23.5 g·C·m −2 ·yr −1 ), respectively, and explained a small proportion of the total variance in C accretion (η 2 = 0.005). Thus, the majority of effect on C accretion rates in our simulations came from N loading, and secondarily from hydrologic regime. At the same time, the effects of invasion and the species of invader on rates of C accretion were detectable and statistically significant.

C accretion in simulated coastal wetlands
Hydroperiod, plant community, N loading, and their interactions all influenced C accretion, illustrating the importance of examining these drivers simultaneously to understand C accretion rates in Great Lakes coastal wetlands. The nonlinear nature of invasion success, community NPP, and C accretion rates across the N loading gradient and their dependence on hydrology shows the importance of using models that allow nonlinear relationships to arise, as many ecological phenomena are expected to be nonlinear (May 1986, Turner 2005. In our simulations, we found that N loading had the greatest effect on C accretion, which consisted of a 958% increase in C accretion from a N inflow of 0.86-30 g·N·m −2 ·yr −1 . It should be noted that a potential reason a 34fold increase in N loading only resulted in a 9.58fold increase in C accretion is N retention, which is low in these systems (~17% for invasion scenarios, Currie et al. 2014), so most of the N flowing into the wetland is lost. Fig. 3. The proportion of total community NPP (above-and belowground dry mass) attributed to the invader species (either Typha or Phragmites) across the N loading gradient (mean ± SE). Results here are averaged over the last 5 yr and across replicates as in Fig. 2. Panel labeled Typha inv. represents Typha invasion scenarios and panel labeled Phragmites inv. represents Phragmites invasion scenarios. Different lines in each panel represent different hydroperiods. Native wetlands under the aerobic hydroperiod at the lowest N loading had the lowest C accretion rates (slightly negative), while the largest C accretion rates were seen in anaerobic wetlands dominated by Phragmites at the highest N loading 30 yr after initial invasion. This large magnitude of change should not be surprising given that native Great Lakes plant communities without anthropogenic N have very low NPP and not much of a litter layer (Angeloni et al. 2006, Tuchman et al. 2009), while the introduction of highly productive invasives at high N greatly increases NPP (Ehrenfeld 2003 and thus ecosystem C accretion. As explained in Currie et al. (2014), we used an ecosystem N inflow of 1.5 g·N·m −2 ·yr −1 as a baseline condition (ecosystem C and N pools are in equilibrium). In this set of model runs, we expanded the range of N inflows a little lower (to 0.86 g·N·m −2 ·yr −1 ). The negative C accretion is small and only occurs in aerobic conditions at the lowest N inflow, likely caused by a small loss of sediment C under aerobic conditions compared with this baseline condition. Additionally, while the low NPP and C accretion rates were lower than normally observed in the literature (Table 3), those studies did not include aerobic sites with extremely low N loading, where low NPP and C accretion rates might occur. There is evidence that sandy coastal wetlands receiving low N inputs in the Great Lakes region accumulate very little soil organic matter (~0.4% soil C) (Elgersma et al. 2015).
Although we expected invasion by Typha and Phragmites would increase C accretion rates at high N loading levels, we did not expect the observed increase at the low end of the N gradient. This result is especially surprising because Phragmites and Typha accounted for less than 20% of the community NPP and did not increase total community NPP compared with the native community at low N loading. The increase in C accretion rates at low N levels is likely due to differences in plant traits between natives and invaders (van Kleunen et al. 2010), including maximum size of individual plants by species (Table 2). Although NPP did not differ between natives with or without invaders at low N loading, living biomass C pools were larger for invasives (Fig. 5) owing to larger rhizomes (belowground C storage; Table 2). Additionally, Typha litter decomposes much more slowly than that of the native species (Table 2), thereby allowing litter to build up over multiple years. Overall, these differences in plant traits increased C accretion under invasion scenarios at low N loading and illustrate the importance of considering how plant traits, even from a subdominant species, can mediate ecosystem function (Chapin and Eviner 2003, Eviner 2004, Elgersma and Ehrenfeld 2010. Although plant invasion significantly increased C accretion at the low and high ends of the N loading gradient, plant community had the smallest effect size (η 2 = 0.005) of the three drivers and can only be considered a small influence on C accretion rates. This result is likely influenced by the ability of the native community, especially S. acutus, in our simulations to increase NPP under high rates of N input, much like the larger-sized invaders Phragmites and Typha. This is an interesting result because most eutrophic wetlands in the Great Lakes region are invaded or in the process of invasion; uninvaded native communities virtually do not exist at high N loading in this region Kercher 2004, Lishawa et al. 2010). Eutrophication is therefore an important driver of C accretion because it increases C accretion in all wetland communities (including natives), and it is a strong driver for invasion (Currie et al. 2014), which then further increases C accretion. Because these two factors (invasion and N loading) virtually always co-occur in natural systems, it is difficult to tease apart their separate effects in field studies and virtually impossible to compare native and invaded communities at high N loading in the field.
The important role hydrology plays in controlling organic matter accumulation in wetlands has been emphasized previously using mathematical models of northern latitude peatlands (Hilbert et al. 2000). Our study focused on midlatitude temperate wetlands and found hydroperiod exerted a strong effect on C accretion, although not as strong as variability in N loading. Higher-latitude wetlands are generally oligotrophic, while in present-day mid-latitude temperate wetlands greater variation in N loading occurs (Bedford et al. 1999). The potential importance of hydrology on C accretion was strongly inferred in a study of a semiarid grassland that investigated the effects of nitrogen addition on ecosystem carbon pools. Zeng et al. (2010) showed that after 5 yr of N addition (20 g·N·m −2 ·yr −1 ) aboveground C pools (living biomass and litter) significantly increased, but soil C pools did not. This lack of response of the soil C pools likely has to do with the aerobic nature (and thus relatively high decomposition rates) of semiarid grassland soils. In addition, plants responded to the N addition by increasing aboveground biomass while reducing belowground biomass, resulting in almost no change to total biomass (Zeng et al. 2010). Conversely, it has been shown that N addition increases both aboveground and belowground biomass of the most problematic invasives in the Great Lakes area (Woo and Zedler 2002, Rickey and Anderson 2004, Martina and von Ende 2012. Similar to the findings of Zeng et al. (2010), González-Alcaraz et al. (2012) found that Phragmites-invaded wetlands increased aboveground, belowground, and litter C pools, but failed to increase soil C pools at dry sites.
Flooding conditions also can shift C losses in wetlands from CO 2 to CH 4 , a shift that we did not explicitly model but is relevant to the broader context of climate change due to the greater warming potential of CH 4 vs. CO 2 . Therefore, while flooding might be viewed as beneficial because it leads to greater C accretion, it is important to understand how flooding favors methanogenesis over other anaerobic processes in any particular wetland ecosystem (Reddy and Delaune 2008).

Ecosystem pools responsible for C accretion
In flooded, anaerobic conditions, the development of the muck layer greatly increased ecosystem C accretion rates in our simulations, especially at high N loading where NPP (and thus litter production) was high. Under these conditions, C accretion was three to four times that of the other hydrologic regimes in our results. The important role of muck in C accretion is consistent with field data showing high muck accumulation in highly productive (usually invaded) wetlands that are flooded for the majority of the year (Angeloni et al. 2006, Fickbohm and Zhu 2006, Bernal and Mitsch 2012. In anaerobic conditions, our simulated rates of muck C accretion followed a nonlinear trend along the N loading gradient that mirrored the living biomass C pool, suggesting that the muck C pool was responding to changes in plant community NPP (Fig. 5). Conversely, v www.esajournals.org MARTINA ET AL.
MSOM C followed a more linear relationship along the N loading gradient, suggesting that MSOM C was more buffered from changes in plant community NPP.
While both Phragmites and Typha were successful invaders over some part of the N gradient, their trait differences resulted in distinct consequences to total community NPP and C accretion rates. The two most notable differences between Phragmites and Typha in our simulation runs were Phragmites' size and stature (greater maximum size [above-and belowground] and leaf and stem architecture [ Table 1]) and Typha's lower litter quality ( Table 2). The greater maximum size (and different canopy architecture) of Phragmites likely allowed it to increase community NPP at the high end of the N loading gradient to a greater degree than either Typha or natives. The greater NPP subsequently increased the muck C pool and overall ecosystem C accretion. The lower decomposition rate constant of Typha's litter resulted in a much greater increase in the litter C pool along the N loading gradient compared with Phragmites or the native community alone. The large quantity of litter buildup in Typha-invaded wetlands has been observed before and shown to increase its invasion success Goldberg 2009, Larkin et al. 2012). In our simulations, litter buildup allowed Typha to increase C accretion rates at low and high N loading in aerobic conditions. The greater litter buildup in Typha-invaded wetlands reduced N cycling by locking N up in the litter layer, reducing the positive feedback of invasions on N cycling seen previously (Currie et al. 2014). Litter buildup also slowed C cycling because C locked up in the litter layer reduced C entering the muck or MSOM C pool. These differences illustrate the importance of studying how plant traits of invasive species, such as litter decay rates, affect ecosystem-level processes (Godoy et al. 2010, van Kleunen et al. 2010, Steers et al. 2011) even in two species that share many plant traits (e.g., clonal, highly productive) and might otherwise be considered ecologically interchangeable.

Comparison to previous wetland ecosystem models and limitations
Previous models have explored how ecosystem dynamics influence plant community composition and wetland C storage. Tian et al. (2010) used the Wetland Ecosystem Model (WEM) to determine the effects of prescribed fires on the dominance of the invasive Typha domingensis in a P-enriched wetland in the Everglades (Florida, USA). WEM was used to better understand the complex P dynamics that occur postburn. Through model simulations, air temperature and hydrologic conditions were determined to be driving factors influencing postfire vegetation change. The biogeochemical model Wetland-DNDC (Li et al. 2004) was used to predict how changing environmental conditions affects carbon and hydrologic dynamics in forested wetlands in Florida and Minnesota (USA). In addition, a modified version of Wetland-DNDC was used to determine the effects of management practices on C sequestration and trace gas emissions in forested wetlands (Cui et al. 2005).
Similar to MONDRIAN, C dynamics in Wetland-DNDC are controlled by physiological plant factors, plant C pools, turnover rates, and environmental factors. MONDRIAN differs from Wetland-DNDC and WEM because it is an individual-based model that allows for resource competition among thousands of individuals from different species (Currie et al. 2014). The process-driven individual-based nature of MONDRIAN allowed us to detect patterns of invasion in different hydrologic settings across a N loading gradient. Without the ability to model community dynamics within an individualbased framework, we would not have been able to link complex population/community dynamics to their effects on ecosystem C cycling. The role of vegetation and plant functional groups on C dynamics and alternative stable equilibria have also been demonstrated using mathematical models of northern peatlands (Frolking et al. 2001, Pastor et al. 2002. While these mathematical models are useful, they too lack the individual and process-based nature of MONDRIAN. MONDRIAN does not explicitly model denitrification or methanogenesis, and thus, those losses of N and C, respectively, from the modeled wetland are not completely accounted for. Instead, because of the well-established complexity of modeling and measuring denitrification (Boyer et al. 2006, Groffman et al. 2006, we modeled simplified daily loss of N from the combined effects of hydrologic flushing and denitrification that we fixed as a constant proportion of v www.esajournals.org MARTINA ET AL. soil available N on a daily basis. These two mechanisms of losses are thus implicitly included but not distinguished or modeled with mechanistic detail. In reality, denitrification would be expected to be more dominant in anaerobic conditions and could account for a significant loss of N from the ecosystem (Reddy and Delaune 2008). Additionally, it has been shown that both Phragmites and Typha can increase denitrification rates (Findlay et al. 2003, Lishawa et al. 2014 that would further affect the dynamic loss of N by denitrification. The lack of these effects means our simulations likely underestimate the negative impact of anaerobic conditions on the invasion success of Phragmites because they would further remove available N from the wetland. The lack of the inclusion on methanogenesis in MONDRIAN is likely less of a concern on the accuracy of our C accretion rates because CH 4 production is usually much less significant in the C cycle relative to CO 2 losses (Whiting and Chanton 2001). Future research could incorporate more mechanistic processes controlling both denitrification and methanogenesis in future versions of MONDRIAN.

conclusIons
Understanding the relative importance of the main drivers of C accretion in wetlands and how they interact is critical to our ability to predict how current and future land use and climate change are likely to influence C sequestration in temperate wetlands. We found that N loading played a significant and perhaps dominant role in controlling rates of C accretion in simulated coastal wetlands and that elevated N inputs from human activities are likely to result in significantly increased rates of wetland C storage. At the same time, N inputs showed complex interactions with hydroperiod and plant invasions (including species identity), which also had significant effects on rates of C accretion. While plant invasions did increase C accretion rates at both low and high N loading, in our set of model scenarios the statistical effect size of invasions on C accretion after a 30-yr period was far less than the effects of N loading or hydroperiod. We interpret these results as evidence that while invasions are likely increasing C accretion on the landscape (Liao et al. 2008), N loading is the greater driver of C accretion. Uninvaded native communities could sequester close to the same amount of C as their invaded counterparts even at high N loading, if they were theoretically able to remain uninvaded. Using a community-ecosystem model that could simulate uninvaded eutrophic communities in direct comparison with eutrophic invaded communities, this study allowed us to tease apart these drivers in a unique way. This shows a strength of using ecosystem models to study complex interactions, similar to previous modeling studies of forested wetlands (Li et al. 2004, Cui et al. 2005. The strong correlation between N loading and invasion in our simulations could potentially lead observational studies to attribute an increase in C storage in invaded wetlands more to the invasive plants than the primary driver, N loading. At the same time, modeling studies are only one facet of increasing our understanding of complex community-ecosystem interactions and should be tested and validated with large-scale observational and experimental studies. acknowledgMents This research was funded by the NASA ROSES program (grant NNX11AC72G), the University of Michigan, and the University of Michigan Water Center (grant N017148). We thank Xin Xu for early contributions to the light competition submodel and figure creation. We thank Emily Farrer, Heather Hager, and Johannes Knops for providing data used for species-specific traits and parameterizing the light competition submodel. JM, WC, DG, and KE conceived and designed the study; JM and WC wrote the model code, tested the model, and conducted simulations; DG, KE, and JM conducted field research; JM, WC, DG, and KE analyzed results and wrote the paper. lIterature cIted