Supplementary Material
Energy flux: The link between multitrophic biodiversity and ecosystem functioning
Andrew D. Barnes1,2,3*, Malte Jochum4, Jonathan S. Lefcheck5, Nico Eisenhauer1,2, Christoph Scherber3, Mary I. O‘Connor6, Peter de Ruiter7,8, & Ulrich Brose1,9
1German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Deutscher Platz 5e, 04103 Leipzig, Germany
2Institute of Biology, Leipzig University, Johannisallee 21, 04103 Leipzig, Germany
3Institute of Landscape Ecology, University of Muenster, Heisenbergstrasse 2, 48149 Muenster, Germany
4Institute of Plant Sciences, University of Bern, Altenbergrain 21, 3013 Bern, Switzerland
5Bigelow Laboratory for Ocean Sciences, East Boothbay, ME, 04544 USA
6Dept of Zoology and Biodiversity Research Centre, Univ. of British Columbia, Vancouver, BC, V6T 1Z4, Canada
7Institute for Biodiversity and Ecosystem Dynamics, University of Amsterdam, Science Park 904, 1098XH Amsterdam, The Netherlands
8Biometris, Wageningen University, Droevendaalsesteeg 1, 6708 PB Wageningen, The Netherlands
9Institute of Ecology, Friedrich Schiller University Jena, Dornburger-Str. 159, 07743 Jena, Germany
Corresponding author: Barnes, A.D. ()
Introduction
Here, we will provide a step-by-step example of how to get from ecological community data to the calculation of whole-community energy flux.
Data come from [S1]. The experimental design nested species combinations based on their functional traits within 4 levels of species richness: 1, 3, 6, or 9 species. We will ignore the traits component for this exercise and focus only on the levels of richness.
The experimental communities included draws from a pool of estuarine consumers, including: 5 herbivore species (the snail Bittiolum varium, two amphipods, Gammarus mucronatus and Ampithoe spp., an isopod Erichsonella attenuata, and the shrimp Hippolyte pleuracanthus) and five predators (the blue crab Callinectes sapidus, the grass shrimp Palaemonetes sp., the pipefish Syngnathus sp. and the mummichog Fundulus heteroclitus). All communities received a 30 g wet-weight of the red alga Gracilaria sp..
The experiment was a combined substitutive/additive design with biomass of herbivores equalized and initial biomass estimated for all predators and controlled statistically.
Experimental communities were assembled in filtered flow-through mesocosms and allowed to interact for 3 weeks, after which time the setup was broken down. All animals were processed in the lab: for larger predators and algae, individuals were weighed and combusted to obtain an ash-free dry weight (in mg), while ash-free dry weight of smaller herbivores was estimated through size fractionation and the equations in [S2]. Final algal biomass summed contributions from recruiting red and green algae, mostly Polysiphonia sp. and Cladophora sp..
Steps 1-6 following Box 1, Figure I
General setup
First, load the required packages, set the working directory, and load the data file.
rm(list=ls())
library(splitstackshape) # used to convert data to individual-level data frame
library(igraph) # used to visualize the network structure
# Set local folder
setwd("...")
# Read the mesocosm data
mesocosmdata <-read.csv("Lefcheck & Duffy 2015.csv", sep =",")
The new dataframe mesocosmdata provides community data with each row containing information for a certain size class (sieve) within each species and sample. The columns abundance, biomass(mg dry weight), trophic.level, and mlO2.mg.h (metabolic demand in ml O2 mg-1 h-1) contain further details for this level of resolution. Column richness defines the richness treatment based on the respective sample id.
For our purpose and in order to provide a step-by-step example of the whole workflow, we will now transform the initial data structure in mesocosmdata to yield a single row in the data frame for each individual. The new individual-level data frame is called idat. This can then be used to calculate a body mass and metabolic rate for each individual.
idat <-expandRows(mesocosmdata, "abundance", drop =FALSE) # 23341 individuals
Note that, in this new dataframe, information on abundance and biomass is information on the level of the size class per species and sample. It helps us to calculate individual-level body mass, but is no longer relevant on the individual level.
Step 1: Measure or assign individual body masses
As our example dataset provided biomass and abundance per size-class, species, and sample id, we simply divide biomass by abundance to get the individual body mass in mg dry weight.
idat$dw <-idat$biomass/idat$abundance # dry weight
Subsequently, we convert these dry body masses to fresh body masses (mg) based on a general conversion factor (fresh weight = dry weight * 4) from [S3]. Dry-fresh mass scaling relationships are also available from various literature sources, but we use this general conversion factor here for simplicity.
idat$fw <-idat$dw *4# fresh weight
Step 2: Calculate individual metabolic rates
Here, we will calculate individual metabolic rates for herbivores and predators based on measured metabolic demand and metabolic-rate regressions from the literature, respectively. Wherever possible, we use metabolic rate information measured for the same species, individual body mass, experimental temperature (23 °C) and salinity (20 ppt) as in our dataset. For most of the herbivores, the example dataset already contains measured metabolic rates for the species in the experiment. Therefore, we only need to apply these metabolic rates to each individual and convert the rates to the desired units (in this case joules per hour, or J h-1). To do this, we multiply the given metabolic rates in column idat$mlO2.mg.h (in ml O2 mg-1 h-1) by body mass and multiply the product by 20.1, according to the conversion factor for ml O2 to J provided by [S3].
idat$mr <-idat$mlO2.mg.h *idat$fw *20.1# metabolic rate in J h-1
There are five remaining species in the dataset that do not have measured metabolic rates. For these, we must rely on literature data.
Fundulus heteroclitus
For this mesopredator fish, we use parameters from a regression of metabolic rate on body mass and temperature from [S4]. The regression equation is
where M is individual fresh body mass and T is temperature in °C. Metabolic rate Q is in ml O2 g-1 h-1, so we multiply it by body mass and then again by the conversion factor 20.1 to get metabolic rate in J h-1. Note that, where necessary, further unit conversions of e.g. mg to g are made in the below calculations without specifically mentioning it in the text every time.
idat$mr[idat$species=="Fundulus heteroclitus"] <-
exp(-0.919-0.204 *log(idat$fw[idat$species=="Fundulus heteroclitus"]/1000)
+(0.028 *23)) *(idat$fw[idat$species=="Fundulus heteroclitus"]/1000) *20.1
Syngnathus spp.
For the pipefish, we use data from [S5], who provide parameters for a regression of metabolic rate R (µmol O2 h-1) on fresh body mass M in g.
We use a conversion of 1 µmol O2 = 0.022391 ml O2 from
idat$mr[idat$species=="Syngnathus spp."] <-
exp(0.8 *log(idat$fw[idat$species=="Syngnathus spp."]/1000)-5.43) *22.391 *20.1
Callinectes sapidus
The Blue Crab individuals in our example dataset are juveniles, for which we found suitable data on metabolic rates from [S6]. Metabolic rates at 20 °C are found to be 0.0741 ml O2 g-1 h-1. Again, we need to use the 20.1 conversion factor to get from ml O2 g-1 h-1 to J h-1.
idat$mr[idat$species=="Callinectes sapidus"] <-
0.0741 *(idat$fw[idat$species=="Callinectes sapidus"]/1000) *20.1
Palaemonetes spp.
For the shrimp Palaemonetes spp., we found data from [S7]. They report metabolic rates of 46.5 µmol O2 g-1 h-1 at 25 °C and 17 ppt salinity, which is the closest we could get to the experimental conditions of our example dataset. We then use the above mentioned conversion factor from µmol O2 to ml O2 and the conversion from ml O2 to J.
idat$mr[idat$species=="Palaemonetes spp."] <-
46.5 *(idat$fw[idat$species=="Palaemonetes spp."]/1000) *0.022391 *20.1
Hippolyte pleuracanthus
For this second shrimp species, we did not manage to find species-specific metabolic rate information. We thus decided to use the same scaling as for the other shrimp species, Palaemonetes spp..
idat$mr[idat$species=="Hippolyte pleuracanthus"] <-
46.5 *(idat$fw[idat$species=="Hippolyte pleuracanthus"]/1000) *0.022391 *20.1
Step 3: Assign network topology and feeding preferences
As described above, there are three trophic levels in this example dataset: predators, herbivores and primary producers. For the sake of simplicity, we use these three trophic levels as a three-node food chain (for an illustration, see Figure S1). In this example there is only a single node per trophic level, and each trophic node only has one potential resource. Therefore, we do not need to specify feeding preferences.
Step 4: Group individuals into network nodes and calculate node metabolism for each experimental community
Along with the node metabolism for predators and herbivores (sum of all individual metabolic rates for each node), the sample id, richness, predator, herbivore and primary producer abundance, and fresh biomass per sample now get stored in a joint data frame sampledata.
ids <-unique(idat$sample) # 83 sample ID's
ids <-sort(ids) # sort the ID's
richness <-as.numeric(tapply(idat$richness,idat$sample,mean))
abundance <-as.numeric(tapply(idat$sample,idat$sample,length))
pred.abund <-as.numeric(tapply(idat$sample[idat$trophic.level=="Predator"],
idat$sample[idat$trophic.level=="Predator"],length))
pred.abund[is.na(pred.abund)] <-0# replace Na's by 0's
herb.abund <-as.numeric(tapply(idat$sample[idat$trophic.level=="Herbivore"],
idat$sample[idat$trophic.level=="Herbivore"],length))
herb.abund[is.na(herb.abund)] <-0# replace Na's by 0's
prim.abund <-as.numeric(tapply(idat$sample[idat$trophic.level=="Primary producer"],
idat$sample[idat$trophic.level=="Primary producer"],length))
biomass.fresh <-as.numeric(tapply(idat$fw,idat$sample,sum))
pred.biomass <-as.numeric(tapply(idat$fw[idat$trophic.level=="Predator"],
idat$sample[idat$trophic.level=="Predator"],sum, na.rm=TRUE))
pred.biomass[is.na(pred.biomass)] <-0# replace Na's by 0's
herb.biomass <-as.numeric(tapply(idat$fw[idat$trophic.level=="Herbivore"],
idat$sample[idat$trophic.level=="Herbivore"],sum, na.rm=TRUE))
herb.biomass[is.na(herb.biomass)] <-0# replace Na's by 0's
prim.biomass <-as.numeric(tapply(idat$fw[idat$trophic.level=="Primary producer"],
idat$sample[idat$trophic.level=="Primary producer"],sum, na.rm=TRUE))
prim.biomass[is.na(prim.biomass)] <-0# replace Na's by 0's
# Calculate summed metabolic demand of predators per sample
pred.met <-as.numeric(tapply(idat$mr[idat$trophic.level=="Predator"],
idat$sample[idat$trophic.level=="Predator"],sum, na.rm=TRUE))
# Calculate summed metabolic demand of predators per sample
herb.met <-as.numeric(tapply(idat$mr[idat$trophic.level=="Herbivore"],
idat$sample[idat$trophic.level=="Herbivore"],sum, na.rm=TRUE))
# Joint data frame
sampledata <-data.frame(ids=ids, richness=richness, abundance=abundance,
pred.abund=pred.abund, herb.abund=herb.abund, prim.abund=prim.abund,
biomass.fresh=biomass.fresh, pred.biomass=pred.biomass,
herb.biomass=herb.biomass, prim.biomass=prim.biomass, pred.met=pred.met,
herb.met=herb.met)
Step 5: Assign assimilation efficiencies to each interaction
Here, for simplicity, we define assimilation efficiencies at the consumer level. Average values for predators and herbivores are taken from [S8].
# Predator assimilation efficiency:
ePre <-0.906
# Herbivore assimilation efficiency:
eHer <-0.545
It should be noted that the level of resolution used to apply parameters such as assimilation efficiency or metabolic rate will always be a trade-off between the availability of data at the desired resolution (e.g., species-specific assimilation efficiencies will not be available for many species) and the research question one is trying to answer. For example, if the research focus is the effect of temperature on energy flux in food webs, then scaling assimilation efficiencies and metabolic rates with temperature would be very important, given the established dependence of these parameters on temperature [S8,S9].
The importance of resolution is true for various other parameters used in the calculation of energy fluxes. Depending on the research question, parameters such as body masses, metabolic rates, feeding preferences, and assimilation efficiencies should ideally be applied at a suitable level of resolution (i.e., individual, species, or functional group level) and, where appropriate and possible, scaled with environmental variables such as temperature or other relevant parameters (e.g., resource quality [S10]).
Step 6: Calculate energy fluxes including losses to predation
To calculate energy flux in a food web, we must start by calculating flux to the highest trophic level and then sequentially calculate the next lowest trophic level toward the bottom of the food web. Thus, in our case, we first calculate flux to predators, and then to herbivores. This way, we account for losses to consumption by higher trophic levels. Before we actually calculate the fluxes, we will first remove some experimental replicates from our data frame. This is because, for the purposes of this example, we only want to compare replicates where the full tri-trophic community is present (i.e., replicates with all trophic levels present). Therefore, we remove all replicates that are missing either predators, herbivores, or both.
# Retain only samples where all three trophic levels are present
sampledata <-subset(sampledata, pred.abund!=0herb.abund!=0)
Now, energy flux will be calculated based on the formula , where e is assimilation efficiency, X is node metabolic demand and L is loss to consumers (i.e., energy flux to higher-level nodes).
First, we define empty vectors to store the fluxes in:
F_Ps <-c() # flux from herbivores to predators per sample
F_Hs <-c() # flux from primary producers to herbivores per sample
samplefluxs <-c() # community energy flux per sample (sum of all fluxes per sample)
Then, we run a for-loop over all sample id’s (experimental replicates) that loads the respective predator and herbivore metabolic demand (X), checks which trophic levels are present (if-conditions), calculates the fluxes based on the above equation and assigned parameters, and stores the calculated fluxes in the prepared vectors.
for(i in 1:length(sampledata$ids)){
temp.id <-sampledata$ids[i] # current sample
X_P <-sampledata$pred.met[i] # metabolic demand of predators in sample i
X_H <-sampledata$herb.met[i] # metabolic demand of herbivores in sample i
if(sampledata$pred.abund[i]>0sampledata$herb.abund[i]>0){# pred. & herb. present
F_P <-(1/ePre) *X_P # predators have no loss to consumption
F_H <-(1/eHer) *(X_H+F_P) # herbivore energetic loss due to predation by predators
}
if(sampledata$pred.abund[i]==0sampledata$herb.abund[i]>0){# no pred., but herb.
F_P <-0# no flux to predators, because there are none in these samples
F_H <-(1/eHer)*X_H # no loss to predation
}
F_Ps <-c(F_Ps, F_P) # write calculated fluxes in vectors: flux to predators
F_Hs <-c(F_Hs, F_H) # write calculated fluxes in vectors: flux to herbivores
sampleflux <-sum(F_P,F_H) # calculate sum of all fluxes
samplefluxs <-c(samplefluxs, sampleflux) # write fluxes in vectors: total flux
}
Finally, we add the calculated whole-community flux, flux to herbivores and to predators to the sample data frame.
sampledata <-cbind(sampledata, flux.to.pred=F_Ps, flux.to.herb=F_Hs,
energyflux=samplefluxs)
Visualize results
Finally, we visualize the results in a figure where panel a) shows the trophic network (in this case, a tri-trophic food chain) and panel b) shows the relationship between the treatment levels of species richness and energy flux, with multitrophic flux (of the whole, tri-trophic network) and individual fluxes to predators (predation) and herbivores (herbivory) displayed seperately.
pdf("Energyflux_results_lefcheck_example.pdf", width =10, height =5)
par(mfrow=c(1,2), oma=c(0,0,2,1), mar=c(4,4,0,0))
# Define variables for models and plotting
richness <-sampledata$richness
energyflux <-sampledata$energyflux
predflux <-sampledata$flux.to.pred
herbflux <-sampledata$flux.to.herb
##########################################
# a) plot igraph food chain
chain <-graph( c("Primary producers","Herbivores","Herbivores","Predators"))
layout_chain <-matrix(c(1,1,1,1,2,3,1,1,1), nrow=3, ncol=3)
plot(chain,
layout=layout_chain,
vertex.color=c( "grey","#fbbf17","#ef3729"),
vertex.frame.color=FALSE,
vertex.size=c(log10(mean(prim.biomass)) *20,
log10(mean(herb.biomass)) *20,log10(mean(pred.biomass)) *20),
edge.color=c("#fbbf17","#ef3729"),
edge.arrow.size=0,
edge.width=c(mean(herbflux),mean(predflux)),
vertex.label.cex=1.2,
vertex.label.dist=c(14,11,10),
vertex.label.degree=0)
mtext("a)", side=3, line=1.5, cex=1.2, adj=0, font=2)
##########################################
# b) plot energy flux and richness
modfull <-lm(energyflux~richness)
summary(modfull)
modpred <-lm(predflux~richness)
summary(modpred)
modherb <-lm(herbflux~richness)
summary(modherb)
# Plot energyflux vs richness
plot(richness, energyflux, type="n", xlab="Richness treatment", ylab="Energyflux [J/h]")
points(richness, energyflux, col="black", pch=19)
points(richness+0.05, predflux, col="#ef3729", pch=19) # 0.05 offset (overlapping points)
points(richness+0.1, herbflux, col="#fbbf17", pch=19) # 0.1 offset (overlapping points)
legend("topleft", legend =c("Multitrophic flux","Herbivory","Predation"), bty ="n",
col =c("black","#fbbf17","#ef3729"), lty=c(1,2,2), lwd=2.5, pch=NULL, cex=1)
x <-seq(min(richness),max(richness),length.out=length(richness))
yfull <-predict(modfull, newdata=list(richness=x))
lines(x,yfull, col="black", lwd=2.5)
x <-seq(min(richness),max(richness),length.out=length(richness))
ypred <-predict(modpred, newdata=list(richness=x))
lines(x,ypred, col="#ef3729", lwd=2.5, lty=2)
x <-seq(min(richness),max(richness),length.out=length(richness))
yherb <-predict(modherb, newdata=list(richness=x))
lines(x,yherb, col="#fbbf17", lwd=2.5, lty=2)
mtext("b)", side=3, line=1.5, cex=1.2, adj=0, font=2)
dev.off()
Figure S1.Trophic structure (a) and relationship between energy flux and the experimental richness treatment (stocked species richness) (b) for example data from [S1]. Edge width and node size in panel (a) are scaled by the grand mean of all fluxes to consumers and their biomasses, respectively. Energy flows from the bottom to the top of the network. X-axis values in panel (b) are slightly offset to better visualize overlapping points due to fixed richness levels.
References
S1Lefcheck, J.S. and Duffy, J.E. (2015) Multitrophic functional diversity predicts ecosystem functioning in experimental assemblages of estuarine consumers. Ecology 96, 2973–2983
S2Edgar, G.J. (1990) The use of the size structure of benthic macrofaunal communities to estimate faunal biomass and secondary production. Journal of Experimental Marine Biology and Ecology 137, 195–214