Stoichiometric multitrophic networks reveal significance of land-sea interaction to ecosystem function in a subtropical nutrient-poor bight, South Africa.

Nearshore marine ecosystems can benefit from their interaction with adjacent ecosystems, especially if they alleviate nutrient limitations in nutrient poor areas. This was the case in our oligo- to mesotrophic study area, the KwaZulu-Natal Bight on the South African subtropical east coast, which is bordered by the Agulhas current. We built stoichiometric, multitrophic ecosystem networks depicting biomass and material flows of carbon, nitrogen and phosphorus in three subsystems of the bight. The networks were analysed to investigate whether the southern, middle and northern bight function similarly in terms of their productivity, transfer efficiency between trophic levels, material cycling, and nutrient limitations. The middle region of the bight was clearly influenced by nutrient additions from the Thukela River, as it had the highest ecosystem productivity, lower transfer efficiencies and degree of cycling. Most nodes in the networks were limited by phosphorus, followed by nitrogen. The middle region adjacent to the Thukela River showed a lower proportion of P limitation especially in summer. Interestingly, there was a clear distinction in sensitivities to nutrient limitations between lower and higher trophic level organisms. This is a reflection of their discrepant nutrient turnover times that are either higher, or lower, than that of the systems, and which might provide a balance to the system through this antagonistic influence. Furthermore, by tracking the stoichiometry through entire food webs it appeared how important the role of lower trophic level organisms was to regulate stoichiometry to more suitable ratios for higher trophic level requirements. Although we gained good insight into the behaviour of the three subsystems in the KZN Bight and the role of terrestrial influence on their functioning, a merged approach of incorporating data on metabolic constraints derived from experiments could further improve the representativeness of multitrophic stoichiometric ecosystem networks.

Introduction

Coastal areas, including estuaries, wetlands and the nearshore are among the most productive habitats on Earth. These habitats do not function in isolation, but are interconnected through chemical and biological interchange between them [1,2]. This means that realms may depend on each other for their sustained productivity [3], and an understanding of the interconnectivity's influence on ecosystem functioning is therefore important to their management and conservation [4,5]. How this interconnectivity of habitats is reflected in energy and nutrient pathways and affects multitrophic ecosystem function has hardly been explored [4,6].

Riverine input to the nearshore is one form of habitat connectivity that has impacts in coastal oceans. River plumes, for instance, are a distinct feature of many coastlines. Kang et al. [7] estimated that the world's largest 19 rivers have a mean annual plume area that is > 14% of the world's continental shelf area. Although not all coasts feature large rivers, many coastal areas have experienced negative impacts from riverine runoff due to nutrient outflow, consequent blooms followed by low class="Chemical">oxygen conditions that may lead to the formation of dead zones with wider reaching influences on food webs [8,9]. Concentrations of nutrients, class="Chemical">pan class="Chemical">chlorophyll-a and anthropogenically derived chemicals (e.g. pesticides) are often found in higher concentrations within plumes when present [10]. However, on oligo- to mesotrophic coasts, the terrigenous runoff may actually play a highly significant role in stimulating nearshore production [11], noticeable especially in areas featuring strong seasonal signals for river flow.

The functioning of an ecosystem, from a nutrient point of view, depends not only on its sources, but how the system's organisms consume, transform, recycle and excrete nutrients, and use them for growth and reproduction [12,13]. Within an ecosystem, class="Chemical">carbon [C], class="Chemical">pan class="Chemical">nitrogen [N] and phosphorus [P] are important macronutrients in processes across various levels of organisation from the cellular to the ecosystem level. For instance, the growth of individual organisms can be limited by low nutrient concentrations which can affect population dynamics [14] and interspecific interactions [15,16] within communities. Nutrient limitations can also affect key ecosystem processes such as nutrient cycling [17-20] and ecosystem structure [21]. Nutrient limitation at the levels of organisms, populations and ecosystems could thus extend to fisheries yields, especially that of top predators affecting fisheries and food security in nutrient-poor coastal regions in the tropics and subtropics. In areas where riverine outflow is a major nutrient source, it has been shown to increase production of lower trophic levels and that of target and forage fish species [22-24]. Productivity in the Gulf of Aden decreased after reduced river flow following the construction of the Aswan Dam on the Nile River [25], and increased in the region following higher nutrient delivery [26]. Such autochthonous sources to nutrient-poor systems can be important, but may also be responsible for altering elemental composition and stoichiometry, of autotrophs and heterotrophs [27].

Freshpan class="Chemical">water systems [28,29] have received more attention in comclass="Chemical">parison to marine ecosystems in terms of the stoichiometry of class="Chemical">poclass="Chemical">pulations and communities. In the latter, most studies are restricted to class="Chemical">phytoclass="Chemical">plankton or organic matter stoichiometry [30-32], and other investigations are often confined to single grouclass="Chemical">pings [e.g. benthos [33], Crustacea [34] or Teleostei [35]. Few ecosystem-level studies have incorclass="Chemical">porated stoichiometric data [20,36-39], and to our knowledge only one of these ensured adherence to node and flow stoichiometry during multitroclass="Chemical">phic network construction and mass balancing [20]. A more comclass="Chemical">prehensive understanding of nutrient content and multitroclass="Chemical">phic stoichiometry in an ecosystem can be class="Chemical">particularly useful for interclass="Chemical">preting the structure and functioning of food webs. Our study area, the KwaZulu-Natal Bight (KZN Bight) on the subtroclass="Chemical">pical east coast of South Africa, relies on seasonal or intermittent oceanic and riverine nutrient sources. Previous studies have investigated various asclass="Chemical">pects of the class="Chemical">physical and biological oceanograclass="Chemical">phy of the region, as well as effects of riverine inclass="Chemical">put from the Thukela river on selected class="Chemical">pan class="Species">biota in the KZN Bight [40,41].

Considering issues such as increased river pan class="Chemical">water abstraction for agriculture and industry, and climate change affecting rainfall and food security, it is imclass="Chemical">portant to understand the role of riverine nutrient inclass="Chemical">put to the class="Chemical">productivity and function of coastal marine ecosystems. In this study, a class="Chemical">particular emclass="Chemical">phasis was class="Chemical">placed on the sclass="Chemical">patial variability of various abiotic comclass="Chemical">ponents and biotic communities in three areas of the KZN Bight, which reclass="Chemical">present its northern, middle and southern regions. We aimed to investigate nearshore multi-nutrient food webs with class="Chemical">particular reference to riverine influence on the heterogeneity of the bight ecosystem in terms of its class="Chemical">productivity, nutrient transfer efficiency, stoichiometry and nutrient limitations.

Methods

Study areas

The KZN Bight extends from Duclass="Chemical">rban in the south (29.8587° S, 31.0218° E) to Caclass="Chemical">pe St. Lucia in the north (28.5151° S, 32.3995° E) and to the edge of the Agulhas current at the 200m isobath (Fig 1) [42]. Three subsystems were chosen for analyses based on samclass="Chemical">pling sites during the African class="Chemical">pan class="Species">Coelacanth Ecosystem Programme (ACEP II) research cruises in summer (February) and winter (August) 2010 [43]. These were based on hypothesised locations of the three major nutrient sources to the bight [southern eddy, central river outflow, northern upwelling]. The “Durban Eddy” subsystem (DE) is located at the 200m isobath in the southern area of the bight (Fig 1). The “Thukela River Mouth” subsystem (TM) is located between the 30m and 40m isobaths in the central area of the bight (Fig 1). Lastly, the “Richards Bay” subsystem [RB] is located between the 30m and 40m isobaths in the northern area of the bight (Fig 1).

Fig 1

The KwaZulu-Natal Bight and location of subsystems.

DE: Durban Eddy, TM: Thukela Mouth, RB: Richards Bay.

DE: Dupan class="Chemical">rban Eddy, TM: Thukela Mouth, class="Chemical">pan class="Chemical">RB: Richards Bay.

Modelling approach

The spatial heterogeneity of ecosystem function within the KZN bight was investigated through multitrophic stoichiometric networks depicting biomass and class="Chemical">carbon (C), class="Chemical">pan class="Chemical">nitrogen (N) and phosphorus (P) exchanges within the ecosystem and across the system boundary. Networks were constructed for each subsystem representing the southern (DE), middle (TM) and northern (RB) bight. For each of the three areas, a network representing the summer season and winter season were constructed, resulting in a total of 18 networks (three areas, two seasons, three nutrients).

The class="Chemical">carbon networks were class="Chemical">partially based on class="Chemical">previously established networks reclass="Chemical">presenting the KZN Bight [44] and Thukela area [45]. For this study, networks were constructed for all three areas of the bight, and sclass="Chemical">pecies reclass="Chemical">presenting the three areas were aggregated into the same functional grouclass="Chemical">ps for all networks (S1 Table) [40,46] to account for the influence of node aggregation on ENA metrics [37,47,48]. Detritus was class="Chemical">partitioned into susclass="Chemical">pended class="Chemical">particulate organic matter (POM), sediment POM, dissolved organic matter (DOM) and dissolved inorganic matter (DIM). All data sources for all nodes are outlined in S2, S3 and S4 Tables. Cetaceans were not included in the DE network because resident class="Chemical">pan class="Species">dolphin species, which dominate the biomass of this group, do not occur beyond the 50m isobath [49,50]. Other data used for network construction were calculated from measurements taken during the two ACEP II cruises and related demersal trawls in 2010, or from literature sources (see Data sources section).

class="Chemical">Carbon networks were constructed and class="Chemical">parameterised first using Ecoclass="Chemical">path with Ecosim software, version 6 [51]. Inclass="Chemical">put class="Chemical">parameters used to class="Chemical">pan class="Disease">construct Ecopath networks include diet composition, biomass (gC m-2), production/biomass ratio (P/B, year-1), consumption/biomass ratio (Q/B, year-1) and ecotrophic efficiency (EE, proportion) which represents the proportion of production utilised in the system. If one of either B, P/B, Q/B or EE is unknown, it is estimated using two mass-balance equations which ensure the inflows into each group in the system are equal to its outflows [52]. Nitrogen and phosphorus networks were constructed using data from carbon networks and C:N, C:P and N:P ratios from [40] or literature (see 2.3 Data Sources). From these stoichiometric, mass-balanced networks, nutrient transfer efficiencies and recycling, nutrient limitations for nodes and limiting flows were calculated (see 2.4 Network Analysis Methodology), to illustrate differences and similarities in the three subsections of the KZN Bight, and link them to influences of adjacent ecosystems.

Construction of multitrophic stoichiometric networks

Data sources

class="Chemical">Carbon biomasses of most model grouclass="Chemical">ps were measured from samclass="Chemical">ples collected during the 2010 ACEP II research cruises [40]. Other basic inclass="Chemical">put data, needed to class="Chemical">pan class="Disease">construct Ecopath networks (P/B, Q/B, EE) were collected from published and grey literature (S2 Table). Exceptions to the latter were the P/B ratios of diatoms, flagellates and of prawn and shrimp. Ratios for diatoms and flagellates were calculated from biomass and production measurements in each subsystem [53]. Due to a lack of production measurements during the August [winter] cruise, P/B ratios were based on measurements from the February [summer] cruise. Diatom P/B was 199 y-1, 242 y-1 and 490 y-1 in the DE, TM and RB subsystems respectively. Flagellate P/B was 208 y-1 in the DE subsystem, 406 y-1 in the TM subsystem and 410 y-1 in the RB subsystem. In the DE networks, at the 200m isobath, a P/B of 2.5 y-1 was used for prawn and shrimp based on deep-water prawns [54], while in the TM and RB networks, which represent sites between 30m and 40m isobaths, a P/B of 2.73 y-1 was used based on shallow-water prawns [54].

Diet compositions in the networks were derived from an isotope study of the region [55,56], and from literature (S3 Table). The fate of biomass not consumed in the systems was assigned to the various [non-living] abiotic groups (S2 Table).

Imports in the form of suspended POC, DOC and DIC were included in the TM and pan class="Chemical">RB networks due to significant river outflow to the study areas. In the TM networks, susclass="Chemical">pended POC imclass="Chemical">port was calculated using a sediment concentration for the Thukela River of 4.28g L-1, the area of the Thukela mudbanks (561 km2) [45] and flow rates based on a January flow rate of 429 m3 s-1 and a June flow rate of 25 m3 s-1 [57]. Values of DOC imclass="Chemical">port were calculated based on an average global estimate for rivers of 6 mg L-1 [58], and the flow rates and area mentioned above. DIC imclass="Chemical">ports for the TM networks were calculated using DIN values from the Thukela River estuary (see below) and a DIC:DIN ratio of 8.2 [59]. In the class="Chemical">pan class="Chemical">RB networks, a sediment concentration of 0.425g L-1 was calculated based on a sediment yield for the Mfolozi River of 6.8x105 t y-1 and a mean annual flow of 1.6x1012 L y-1 [60]. A POC concentration of 2.57x1010 g y-1 was calculated using the sediment concentration and the assumption that POC was 8.4% of TSS [59]. DOC imports were calculated based on an average global estimate of 6 mg L-1 [59], a January flow rate of 6x1010 L and a June flow rate of 6x109 L for the Mfolozi River [61]. A summary of imports can be found in Table 1.

Table 1

Detrital imports [g m-2 y-1] into the Thukela Mouth and Richards Bay networks.

		Thukela Mouth		Richards Bay
Detritus group		Summer	Winter	Summer	Winter
Suspended POM	C	8 671 358	501 964	45 900	4 582
	N	1 057 483	61 215	5 598	559
	P	433 135	25 073	2 293	229
DOM	C	144 789	8 381	7 701	770
	N	7 239	419	385	39
	P	144	8	8	1
DIM	C	62 332	1 306	3 315	120
	N	7 601	159	404	15
	P	531	17	28	2

Biomasses, imports and exports in the class="Chemical">nitrogen and class="Chemical">pan class="Chemical">phosphorus networks were calculated based on the corresponding carbon data for each network and C:N and C:P molar ratios measured from samples taken during the two cruises and demersal trawls [40] or from literature (S4 Table). Seasonal DIN and DIP imports were calculated for the TM and RB networks using the flow rates and area mentioned in the carbon data sources above and nitrate and phosphate concentrations from the Thukela River [57] (Table 1). The C, N and P biomass for all 18 networks is shown in S5 Table.

Parameterisation of carbon networks

Six class="Chemical">carbon networks reclass="Chemical">presenting each subsystem in summer 2010 and winter 2010 were constructed. Initially, none of the six class="Chemical">pan class="Chemical">carbon networks were mass-balanced. The proportion of the production utilised in the system (EE) was higher than one for many groups because outflows, in the form of predation, were larger than inflows. To balance the networks, biomass data of groups which were undersampled by the ACEP II trawling gear (other cephalopods, cuttlefish, molluscs, echinoderms, large suspension feeders) [46] were removed from one or all of the carbon networks and instead estimated using an EE of 0.95 [62] in Ecopath. The proportion of cannibalism was decreased (as suggested by [62]) in large sharks and replaced with diet imported from outside the subsystem. Diet compositions were adjusted until all groups were balanced and at EE<1.

Thereafter, C networks were extracted from Ecopath to include DIC, and build the corresponding N and P networks. Ecopath does not include the uptake of DIM by primary producers. After the networks were extracted from Ecopath, DIC was added to the pan class="Chemical">carbon networks, calculated using DIN measured in each subsystem during ACEP II cruises [53] and a C:N ratio of 8.2 [59]. Flows from DIC to class="Chemical">primary class="Chemical">producer grouclass="Chemical">ps were calculated based on demand, so that inflow into each class="Chemical">primary class="Chemical">producer grouclass="Chemical">p balanced outflow. The DIC node was balanced based on demand by adding imclass="Chemical">ports or exclass="Chemical">ports across the system boundary.

Nitrogen and phosphorus networks

N and P biomasses were calculated using the stoichiometric ratios for each group (S3 Table). The N and P flows between nodes were calculated by multiplying the pan class="Chemical">carbon flow value by the C:N or C:P ratio of the class="Chemical">pan class="Species">donor biomass. Thereafter, a layered balancing procedure was used to mass balance the network while simultaneously accounting for the large difference in magnitudes of flows and preserving stoichiometry. One layer comprises all flows with a magnitude of the same power of 10. Layers are balanced one by one while keeping other layers constant, in an iterative fashion until convergence is achieved [20].

Network analysis methodology

The eighteen networks were analysed using ecological network analysis (ENA) to determine differences in size and functioning between subsystems in terms of class="Chemical">carbon, class="Chemical">pan class="Chemical">nitrogen and phosphorus and if these differences were associated with riverine nutrient sources. Nutrient limitations, cycling and flow transfer efficiencies across trophic levels were calculated to compare differences in functioning and determine the possible role of riverine nutrient sources.

System metrics

Firstly, the size of each system was determined by its activity as the Total System Throughflow (TSTf), which denotes the sum of all nodal inflows, or nodal outflows, including flows across the system boundary: where T is the flow from node i to node j, z the boundary inflow, and e and r the export, and respiration flow respectively across the boundary.

Nutrient transport between trophic levels in each subsystem was determined by assigning each node to a trophic level according to its feeding activity [63]. The trophic level of each species was calculated as column sum of the input structure matrix [64]. The individual values of the column represent the proportional feeding activity of a node on various trophic levels. Primary producers and class="Disease">detritus nodes were assumed to have a TL of one. class="Chemical">pan class="Disease">Trophic transfer efficiencies between discrete trophic levels were calculated as the fraction of input passed on to the next level via predation [65].

Nutrient cycling in proportion to the systems size was determined using the pan class="Species">Finn Cycling Index (FCI), which measures the fraction of throughclass="Chemical">put recycled [66,67]. This was calculated as follows: where T is the total throughclass="Chemical">put through grouclass="Chemical">p i and (N -1) is the throughclass="Chemical">put through grouclass="Chemical">p i resulting from cycling.

The R package enaR [68] was used to calculate the above indices.

Nutrient limitations of nodes and flows

Firstly, nodal nutrient limitations are calculated sensu Liebig's Law of the Minimum [69] from the turnover times of the three elements (C, N, P) in each compartment. Secondly, the identity (C, N or P) of the limiting flow between network components is found. Both are calculated through a biomass inclusive ascendency that is most sensitive to the limiting element [70].

For this study, the biomass inclusive ascendency was calculated in order to take into account the storage of elements [C, N, P] in each compartment of the food web, and the distribution of biomass and flows between pairs of nodes [source, receiver]. Firstly, the law of mass action was applied to estimate the exchanges between nodes, where the sum of biomass of source nodes i (ΣB), and the biomass of sink nodes j (ΣB) reflect the total biomass in the system (B. = Σ B = Σ B). A priori, mass action estimates the flow leaving node i as [Bi/B.] and the amount entering node j as (B/B.), so that the joint probability of a flow both leaving i and entering j is (BB/B.). The a posteriori conditional probability estimation of material actually leaving node i and entering node j, i.e. the flow from i to j (T) is defined from observations as T/T … The subscript '..' describes the summation over i and j. The information gained by knowing the flows a posteriori compared to estimating the exchange by the a priori probability (BB/B.) can be estimated, according to [71], as: or

The Kullback divergence, or information (I) gained by knowing the flows is then scaled by the total system throughput (i.e., K is set equal to T..) to arrive at the biomass inclusive ascendency A: or

The biomass inclusive ascendency (A) can be extended to apply to each particular nutrient [C, N or P] of interest, using the formula: where T is the flow from node i to node j of element k, and T … the total flow of all elements through all compartments [70]. The sensitivity of the compartment to its turnover time of a particular element is calculated by the partial derivative of A with respect to the biomass of element k in a particular compartment p:

The amount of change in the ascendency value resulting from changes in turnover rate for each element (k) in each compartment (p) is interpreted as follows: a comparatively slow turnover rate results in a high sensitivity value, and this denotes the element that enters the node in the least relative proportion as defined by the consumer biomass stoichiometry (the limiting nutrient), and reiterates the results of a classical Liebig test. Negative sensitivities may arise for compartments when the turnover rate of a specific element is higher compared to the turnover rate of the entire system.

The identification of the limiting source of the limiting element is an extension of the classical Liebig notion. The sensitivities attributed to the individual flows from source r to sink p can be calculated as follows [70]:

The limiting flow that depletes its source at the fastest rate compared to its availability results in the highest sensitivity value [70].

Results

The analysis of the ecosystem networks of the three different sub-regions of the KZN Bight revealed several differences. The largest subsystem, in terms of energy throughput was at the Thukela Mouth, followed by Duclass="Chemical">rban Eddy and Richards Bay (Fig 2). As exclass="Chemical">pected, the overall amount of class="Chemical">pan class="Chemical">carbon (C) flows was larger than that of nitrogen (N) and especially phosphorus (P). The Finn Cycling Index showed almost the opposite pattern–subsystems with highest throughflow had a lower degree of cycling. It was lowest in the TM region, especially in summer when more nutrients are input into this subsystem. The fraction cycled through the system was lowest for C, highest for N and intermediate for P in all three subsystems and both seasons (Fig 2).

Fig 2

Total system Throughflow (TST) and Finn Cycling Index (FCI) of each network.

DE = Durban Eddy, TM = Thukela Mouth, RB = Richards Bay, s = summer, w = winter. C = carbon, N = nitrogen, P = phosphorus.

DE = Duclass="Chemical">rban Eddy, TM = Thukela Mouth, class="Chemical">pan class="Chemical">RB = Richards Bay, s = summer, w = winter. C = carbon, N = nitrogen, P = phosphorus.

pan class="Disease">Trophic Transfer Efficiencies (TE) generally decreased from lower to higher troclass="Chemical">phic levels (TL). At the DE in the southern Bight, TE were highest at TL II and III for all three nutrients and both seasons (Fig 3). This same class="Chemical">pattern is aclass="Chemical">pclass="Chemical">parent in the middle (TM) and northern (class="Chemical">pan class="Chemical">RB) Bight for P.

Fig 3

Trophic transfer efficiencies (as fraction) between discrete trophic levels of each nutrient (C, N, P) in each region (DE, TM, RB) and season (s, w).

Nutrient limitations

The sensitivity values describe the sensitivity of the ascendency to changes in biomass (C, N, P) of individual nodes. Nodes with higher turnover rate (T/B) compared to the system’s turnover rate have a negative sensitivity value, whereas nodes with a lower turnover rate contribute positively to ascendency. A similar pattern is apparent from all 18 networks, in that small organisms with fast turnover times show highly negative values, whereas the remainder of the nodes have positive values. Nutrient nodes are excluded from this analysis, as they do not have a nutrient ‘demand’, and thus no limitation in the sense of Liebig. In our example, nodes 1–7 and 9 have low sensitivity values throughout the bight and in both seasons (Fig 4). These nodes constitute groups of small organisms (Phytoplankton, Bacteria, pan class="Disease">Small Zooplankton, Small Macrobenthos). On the other hand, the nodes with a class="Chemical">positive contribution to the system’s ascendency are larger sized organisms [various invertebrate and vertebrate grouclass="Chemical">ps] of nodes 8, 10 and above (Fig 4, S2 Table). These nodes have in addition very similar sensitivity values amongst each other, unlike the nodes reclass="Chemical">presenting smaller sized organisms, and for the three nutrients.

Fig 4

Node sensitivities for all 18 networks to changes to carbon, nitrogen and phosphorus turnover rates in each network calculated using Eq 9.

A: Sensitivity values against number of compartments. B: sensitivity value against TL of compartments. C: Sensitivity value against TL of compartments for each region and season separately. Black: primary producers, Green: bacteria and zooplankton, Red: zoobenthos, Blue: fish and invertebrate nekton.

The node sensitivities for all three subsystems in the KZN Bight showed that P was the limiting nutrient to more than 50% of nodes (S6 Table). Around 40% of nodes were N limited, whereas C limitation occurred only in the Dupan class="Chemical">rban Eddy for Bacteria and Large Sharks. Diatoms and Flagellates were P limited in all areas and both seasons, sometimes joined by Bacteria and Heterotroclass="Chemical">phic Microclass="Chemical">plankton. Zooclass="Chemical">plankton and the benthos were mainly N limited, whereas most fish were P limited. Co-limitations, denoted by similar sensitivity values were aclass="Chemical">pclass="Chemical">parent for all larger sized organisms [Nodes 8, 10 and above]. No distinct seasonal differences were aclass="Chemical">pclass="Chemical">parent, with the exceclass="Chemical">ption of class="Chemical">primary class="Chemical">producer, small class="Chemical">pan class="Disease">zooplankton and small macrobenthos nodes.

One striking difference between the three regions of the KZN Bight were the smaller sensitivity values (i.e. higher turnover rates) for diatoms and flagellates at TM, implying that nutrients were overall less limiting. This is in accordance with the comparatively high rate of nutrient delivery to this region from the Thukela River. In addition, the throughflow rate is higher in this area compared to DE and pan class="Chemical">RB. Another difference of the TM region was the comclass="Chemical">paratively higher class="Chemical">proclass="Chemical">portion of N, and thus lower P limitation, esclass="Chemical">pecially in summer (Table 2).

Table 2

Percentage of nodes limited by carbon (C), nitrogen (N) and phosphorus (P).

DE: Durban Eddy, TM: Thukela Mouth, RB: Richard Bay, s: summer, w: winter.

	C	N	P
DEs	7.1	35.7	57.1
DEw	3.6	39.3	57.1
TMs	0	44.8	51.7
TMw	0	41.4	55.2
RBs	0	37.9	58.6
RBw	0	41.4	55.2

DE: Dupan class="Chemical">rban Eddy, TM: Thukela Mouth, class="Chemical">pan class="Chemical">RB: Richard Bay, s: summer, w: winter.

We compared the nodal nutrient limitations derived from the ascendency analysis (sensitivity values) to a more traditional approach of finding the difference between the stoichiometry of supply (consumption) and that of requirement (biomass). Both approaches yielded exactly the same results (S6 Table).

In contrast to the nutrient limitation of nodes, all flows in all areas and both seasons were limited through P, except for one occasion of an N limitation. This means that flow pathways for P deplete their sources the fastest of all three nutrients, and the demands of recipient nodes are higher compared to availability at the respective source node. Among the top five source nodes of the limiting flows for all regions and seasons were Diatoms, Flagellates, Small and Medium Copepods, for between 1 and 4 flows per network. Among the top five recipient nodes of limiting flows were Cetaceans in the TM and class="Chemical">RB networks, whereas Small and Medium Coclass="Chemical">peclass="Chemical">pods as well as Small Macrobenthos were either the source or reciclass="Chemical">pients of limiting flows. In DE, limiting flows were sourced by a higher variety of nodes than in TM or class="Chemical">pan class="Chemical">RB. Neither the number of limiting flows nor the sensitivity values showed a significant relationship with the effective trophic level of either the source or the recipient node. However, higher trophic level nodes had overall lower sensitivity values, and thus deplete their source slower compared to lower trophic level nodes (Fig 5).

Fig 5

Sensitivity values of flows (calculated using Eq 10) are plotted against the effective trophic level of both the source and recipient compartment.

Discussion

The analysis of multitrophic, stoichiometric networks of the KZN Bight expanded on previous studies in the region that investigated how subsystems are influenced by riverine input and oceanographic features such as upwelling [72-74]. A major influence of riverine input [including that from the Thukela River] was previously detected through nutrient supplements reaching the nearshore [41,56,75], and energy requirements of its food web [44]. Increased macrobenthos abundance and a unique demersal fish community diversity was found in the middle region of the bight at TM [46,76]. Increased bacteria standing stock but lower primary production were associated with allochthonous inputs from rivers [53,77], and pan class="Disease">high zooplankton biomass was associated with inshore regions of the bight and the St Lucia uclass="Chemical">pwelling cell [78].

Our focus was on the functional attributes of the three subsystems, and on how stoichiometric constraints on a node level structure flows through ecosystems. The ecosystem attributes reported on here emanated from how the system deals with differences in stoichiometric demands and supplies of the various nodes. In addition to flows around individual nodes and between pairs of nodes, an important factor within the supply and demand chain is the speed with which nodes are re-supplied by recycling flows which themselves are subject to stoichiometric constraints by the nodes involved. These aspects influence the stoichiometry of supplies to nodes, and therefore play a role in the structure (e.g. predator-prey pairs, facultative feeding) and dynamics (e.g. nutrient turnover times) of ecosystems [79]. Firstly, the throughflow amount of C, N and P differed between the three subregions and nutrients, and was higher in summer compared to winter in the middle (TM) and northern bight (class="Chemical">RB). These subregions exclass="Chemical">perience inflow from large rivers, increased during the rainy season in sclass="Chemical">pring and summer from the Thukela and Mfolozi Rivers resclass="Chemical">pectively. class="Chemical">pan class="Chemical">RB was in addition influenced by the St Lucia upwelling, as also apparent from e.g. the high zooplankton biomass in the region [78]. The residence time of nutrients can be relatively high in the nearshore region of the middle to northern Bight, due to inshore currents flowing counter to the Agulhas current, and a swirl present at times in the mid-bight [73,80]. A longer residence time provides conditions to better utilise imported nutrients for primary production and by bacteria, and pass them on to higher trophic levels. The TEs were lower in the TM region, as more of the macronutrients were available compared to their demand. Although ecosystem productivity was higher in the TM region (see the TSTf, Fig 1), the primary productivity can be hindered by sediment loads from rivers which decrease light availability [53,81]. In all sub-regions, Transfer Efficiencies for P were higher compared to N, whereas the extent of cycling was higher for N. This might simply be a reflection of a more efficient recycling of N compared to P, and relatively more P being locked into biomass for a longer duration.

To our knowledge, only three other ecosystems have been presented as multitrophic networks of three macronutrients simultaneously. These are the Sylt Rømø Bight in the German/Danish Wadden Sea [82], Chesapeake Bay on the east coast of the USA [39], and an oligotrophic mangrove system offshore Belize [20]. The latter preserves the stoichiometry of nodes and flows throughout network construction and mass balancing, whereas the former two studies have constructed networks for the three nutrients in parallel. Similar to Chesapeake Bay, most nodes were P limited in our systems, whereas limitation results are not reported for the Sylt Rømø Bight. The latter exhibits comparable patterns to the KZN-Bight subsystems in transfer efficiencies (highest for P, lowest for C) and total material flow (lowest for P, highest for C). There was a difference in the amount of cycling, which was highest for P in the Sylt Rømø Bight, but highest for N in the KZN Bight. Our study area therefore experienced a different supply to demand ratio and a more efficient recycling of N relative to P.

The nutrient limitations of individual nodes were calculated from the stoichiometry of their biomass, and that of in- and outgoing links. The sensitivity values translate into differentiated nutrient limitations, as nodes with slower turnover rates for a particular nutrient have longer turnover times for that nutrient. The nutrient being available in the least relative proportion to the node can thus create a discrepancy between the stoichiometry of the biomass (C:N:P) and the turnover rates (T/B) of the elements within that stoichiometric constraint.

The limitations of individual nodes are indicated by sensitivity values that essentially describe the turnover time of nutrients within a node in comparison to that of the entire system. Slower turnover times compared to that of the system result in higher sensitivity values, which thus increased from smaller to larger organisms (Fig 4A), and from lower to higher trophic levels (Fig 4B) for all three nutrients (C, N, P). More than 50% of nodes were limited by P, which is in accordance with the higher transfer efficiencies for P compared to N and C. C limitation was apparent only from the southern bight (DE) to a small proportion.

It appears that the turnover rate of the three subsystems is held in balance by both smaller sized organisms with faster turnover rate than the system, and larger organisms with slower turnover rates. The difference in turnover rates between smaller and larger sized organisms is of course well documented in the literature for r and k strategists [83] and thus reinforces that the developmental state of ecosystems can be manifested in contrasting roles of their parts. Indeed, the contribution by both types of organisms in a network setting has been documented for their contribution to information [84,85]. Moreover, flow efficiency and redundant pathways as antagonistic agents existing in balance to increase the robustness of trophic networks have been documented previously as intrinsic features of such networks [86].

Higher trophic level (TL) organisms were less sensitive to different nutrient environments, in terms of limitations, compared to lower TL organism, expressed by more similar sensitivity values of the former. Lower TL organisms thus take a more prominent part in regulating different levels of nutrient supply for higher trophic levels. If higher TL organisms are less flexible to cope with variable stoichiometries [e.g. 87], lower TL organisms might afford their existence in a more stable stoichiometric environment and so facilitate ecosystem succession.

Network structure, data, and considerations for future stoichiometry studies

Networks for all three KZN Bight subsystems were aggregated into the same number of groups, except where the absence or presence of species demanded otherwise. This choice was made to compare functional attributes (cycling, transfer efficiency, nutrient limitations) across different ecosystems, the calculation of which is influenced by the number of nodes and flows in a network [37,47,88]. The various networks in our study are as realistic a reflection of the subsystems as possible, and the equivalent aggregation facilitated a focus on the nutrient limitations by eliminating other, confounding factors attributable to aggregation.

The calculation of system-level attributes from the information embedded in nodes [species/species groups] and their interlinkages provided an opportunity to ascertain the influence of one level of hierarchy on another. Such an approach can be implemented in the Biodiversity-Ecosystem-Function (BEF) debate [89], as the sensitivity of ecological network analysis (ENA) metrics describing functional ecosystem properties to differences in number of nodes and flow magnitudes is well known, and used especially for linking these two hierarchies [node, ecosystem] [90,91][88,89]. However, a more detailed focus on biodiversity requires a higher detail of species-specific data.

Previous studies on stoichiometry and food webs have mainly been conducted for freshpan class="Chemical">water, rather than marine ecosystems. Exclass="Chemical">periments class="Chemical">provide the bulk of data to freshclass="Chemical">pan class="Chemical">water studies [92], whereas in our data-based approach we constructed food webs based on biomass stoichiometry, and preserved this stoichiometry for resource-consumer flows, which were deducted from measurements or literature data. The former approach provides more exact data on the transformation of the stoichiometry during feeding interactions, whereas the latter approach facilitates a system’s view (e.g. system-wide recycling and transfer efficiencies) and allows to calculate the value of emerging properties of the food webs (e.g. system-wide nutrient limitations). Taking into account elemental content, nutrient ratios and fluxes has previously been mentioned as an important step in arriving at a comprehensive ecosystem view [93]. Merging both experimental and ecosystem analysis approaches could refine stoichiometric food webs by providing better constraints on recourse-consumer flow exchanges between living nodes, and with nutrient pools in multitrophic networks.

Conclusions

The three subsystems of the oligo- to mesotrophic KwaZulu-Natal Bight were influenced by terrestrial and oceanic influences, and the former provided a distinct signal apparent in functional ecosystem attributes. By tracking the C:N:P stoichiometry, we found that all three subsystems were mainly limited by pan class="Chemical">phosphorus, although in the middle bight the N limitation was comclass="Chemical">paratively higher at the exclass="Chemical">pense of P limitation. Transfer efficiencies between troclass="Chemical">phic levels were lower in the TM region, indicating a higher availability of the three nutrients comclass="Chemical">pared to their demand over all troclass="Chemical">phic levels, with an exceclass="Chemical">ption for N and P in winter. An interesting class="Chemical">pattern was aclass="Chemical">pclass="Chemical">parent for the range of sensitivity values of lower and higher troclass="Chemical">phic levels, indicating a divide between two grouclass="Chemical">ps of organisms: those with a faster turnover time than the system, and those with a slower turnover time. Such antagonistic behaviour has class="Chemical">previously been detected for information related indices and might indicate a general feature of ecosystems. Lastly, to investigate the stoichiometry of ecosystems more comclass="Chemical">prehensively and better constrain nutrient flow data in ecosystem models, a merger of various established methodologies currently investigating ecological stoichiometry (exclass="Chemical">periments, metabolic constraints, ecosystem networks) can be useful.

Demersal functional groups and representative species of each sub-region where known.

(DOCX)

Click here for additional data file.

Input data for Ecopath Carbon networks of the DE, TM and RB susbsystems.

P/B = Production/ Biomass, Q/B = Consumption/Biomass, EE = Ecotrophic Efficiency. EE and flows to detritus are proportions.

(DOCX)

Click here for additional data file.

Diet compositions used in initial unbalanced Ecopath carbon networks.

Group numbers refer to those in Table 4.1. Rows represent prey and columns represent predators. Groups 1 and 2 refer to primary producers and therefore do not require a predator column. I = imports. The first value in the column for group 28 was used in the DE networks and the second value was used in the TM and pan class="Chemical">RB networks.

(DOCX)

Click here for additional data file.

C:N and C:P ratios used to construct nitrogen and phosphorus networks.

DE = Dupan class="Chemical">rban Eddy, TM = Thukela Mouth, class="Chemical">pan class="Chemical">RB = Richards Bay.

(DOCX)

Click here for additional data file.

Carbon, nitrogen and phosphorus biomasses (g m-2) of groups in each network.

Bold values represent biomasses estimated by Ecopath.

(DOCX)

Click here for additional data file.

Nutrient limitations calculated from ascendency analysis, and turnover times in a node.

DE: Duclass="Chemical">rban Eddy, TM: Thukela Mouth, class="Chemical">pan class="Chemical">RB: Richard Bay. S: summer, w: winter. LN: Limiting Nutrient. C: carbon, N: nitrogen, P: phosphorus.

(DOCX)

Click here for additional data file.