library(tidyverse)

Initial values and assumptions

#Percentage of carbon (C) per fresh mass
cperw = 0.2 * 0.52

##1. Correct for the expected difference between CO2 evolved and O2 consumed using a respiratory quotient (RQ) of 0.95 
##Klekowski 1972 (doi: 10.1163/187529272x00665) reports a aO2 value of 1.40 nl O2 per hour. This can be converted to a aCO2 value of 2.43 using:
##aCO2 = mass.CO2 * (44 x 10^9)
a.O2 = 1.4
V.CO2 = 0.95 * a.O2

##2. Calculate the number of moles (m) of CO2
##We assume a temperature of 15 degrees Celcius.
t = 15
m = V.CO2 / (0.0821 * (273 + t) * 1E9)

##3. Convert from aO2 to aCO2 
a.CO2 = m * (44 * 1E9)

##4. Convert to daily rate, scale to relative mass of C in CO2. Divide by 1000 to convert from ng to ug.
a.CO2.daily = a.CO2 * 24 * (12/44) / 1000

Data source: For adult biomass and carbon budget calculations we use data from Nemaplex (http://nemaplex.ucdavis.edu/) and Mulder and Vonk 2011, doi:10.2307/23034835. For juveniles, we use data from Mulder and Vonk 2011.

Data import

# Nemaplex data
setwd("~/Work/ETH/Projects/Nematodes/Calculations")
Nemaplex_data_raw <- read.csv("Family Ecophysiology qParameters.csv", stringsAsFactors = FALSE, fileEncoding = "latin1") 

# Mulder 2011 data
mulder_2011 <- read.csv("Mulder2011.csv")
mulder_data <- mulder_2011 %>% 
  select(trophic.group, Length, Width, Juv_Adult)

# Biome abundances
Biome_Abundances <- read.csv("20180827_Biome_Abundances.csv")

# Lowest estimates
numbers_min <- read.csv("20180912_Biome_Abundances_Minimum.csv") %>% 
  select(-Total_Number) %>% 
  gather(key, value, -Biome) %>% 
  rename(trophic.group = key) %>% 
  rename(Total.number = value)

# Highest estimates
numbers_max <- read.csv("20180912_Biome_Abundances_Maximum.csv") %>% 
  select(-Total_Number) %>% 
  gather(key, value, -Biome) %>% 
  rename(trophic.group = key) %>% 
  rename(Total.number = value)

adult_data_1 <- Nemaplex_data_raw %>% 
  select(cp.value, feeding.code, Length.micm, Width.micm) %>% 
  mutate(feeding.code = replace(feeding.code, feeding.code == "p", "Predators")) %>% 
  mutate(feeding.code = replace(feeding.code, feeding.code == "o", "Omnivores")) %>% 
  mutate(feeding.code = replace(feeding.code, feeding.code == "b", "Bacterivores")) %>% 
  mutate(feeding.code = replace(feeding.code, feeding.code == "h", "Herbivores")) %>% 
  mutate(feeding.code = replace(feeding.code, feeding.code == "f", "Fungivores")) %>% 
  filter(feeding.code != "e") %>% 
  filter(feeding.code != "d") %>% 
  filter(feeding.code != "") %>% 
  rename(trophic.group = feeding.code) %>% 
  rename(Length = Length.micm) %>% 
  rename(Width = Width.micm)

cp.values <- adult_data_1 %>% group_by(trophic.group) %>% summarise(cp.value = mean(cp.value)) 

mulder_data$trophic.group <- substr(mulder_data$trophic.group,1,4)  
adult_data_2 <- mulder_data %>%   
  mutate(trophic.group = replace(trophic.group, trophic.group == "bact", "Bacterivores")) %>% 
  mutate(trophic.group = replace(trophic.group, trophic.group == "fung", "Fungivores")) %>% 
  mutate(trophic.group = replace(trophic.group, trophic.group == "herb", "Herbivores")) %>% 
  mutate(trophic.group = replace(trophic.group, trophic.group == "omni", "Omnivores")) %>% 
  mutate(trophic.group = replace(trophic.group, trophic.group == "pred", "Predators")) %>% 
  filter(Juv_Adult == "adult")
adult_data_2 <- adult_data_2 %>% select(-Juv_Adult)
adult_data_2 <- merge(cp.values, adult_data_2[1:3], by = "trophic.group")

adult_data <- rbind(adult_data_1, adult_data_2)

Bodymass calculations

Adults

#mass as function of length and width: W = (L ∙D^2)/(1.6 ∙ 10^6 )
adult_data$mass <- with(adult_data,(Length * Width * Width)/(1.6*1E6))

#Respiration rate in ug C per day as funciton of mass and aC02. b-value is set to 0.75
adult_data$Respiration.rate <- with(adult_data,(a.CO2.daily * ((adult_data$mass)^0.75)))

#Production rate in ug C per day as function of mass and cp value
adult_data$Production.rate <- with(adult_data, (cperw* ((adult_data$mass)/(12 * adult_data$cp.value))))

#Daily carbon budget in ug C; sum of respiration rate and production rate
adult_data$Carbon.budget <- with(adult_data, (adult_data$Respiration.rate + adult_data$Production.rate))

summary_adult <- adult_data %>% 
  group_by(trophic.group) %>% 
  summarise(n.measurements = n(), 
            mean.mass.ug = mean(mass),
            mean.daily.respiration.rate.ug = mean(Respiration.rate),
            mean.daily.production.rate.ug = mean(Production.rate),
            mean.daily.carbon.budget.ug = mean(Carbon.budget)) %>% 
  na.omit() 
summary_adult
## # A tibble: 5 x 6
##   trophic.group n.measurements mean.mass.ug mean.daily.resp…
##   <chr>                  <int>        <dbl>            <dbl>
## 1 Bacterivores            3645        0.824          0.0106 
## 2 Fungivores               390        0.361          0.00674
## 3 Herbivores              4386        1.94           0.0160 
## 4 Omnivores                507        7.27           0.0567 
## 5 Predators                569        4.29           0.0441 
## # … with 2 more variables: mean.daily.production.rate.ug <dbl>,
## #   mean.daily.carbon.budget.ug <dbl>

Juveniles

juv_data <- mulder_2011 %>% 
  select(Taxonomy, trophic.group, Length, Width, Juv_Adult)

juv_data$trophic.group <- substr(juv_data$trophic.group,1,4)  
juv_data <- juv_data %>%   
  mutate(trophic.group = replace(trophic.group, trophic.group == "bact", "Bacterivores")) %>% 
  mutate(trophic.group = replace(trophic.group, trophic.group == "fung", "Fungivores")) %>% 
  mutate(trophic.group = replace(trophic.group, trophic.group == "herb", "Herbivores")) %>% 
  mutate(trophic.group = replace(trophic.group, trophic.group == "omni", "Omnivores")) %>% 
  mutate(trophic.group = replace(trophic.group, trophic.group == "pred", "Predators")) %>% 
  filter(Juv_Adult == "juveniles")

cp.values <- adult_data %>% group_by(trophic.group) %>% summarise(cp.value = mean(cp.value)) 

juv_data<- merge(cp.values, juv_data[1:5], by = "trophic.group")

#Mass as function of length and width: W = (L x D^2)/(1.6 x 10^6 )
juv_data$Mass <- with(juv_data, (Length * Width * Width)/(1.6*1E6))

#Respiration rate in ug C per day as function of mass. aCO2 is as calculated above b-value assumed at 0.75
juv_data$Respiration.rate <- with(juv_data,(a.CO2.daily * ((juv_data$Mass)^0.75)))

#Production rate in ug C per day as function of mass and cp value
juv_data$Production.rate <- with(juv_data, (cperw* ((juv_data$Mass)/(12 * juv_data$cp.value))))

#Daily carbon budget in ug C; sum of respiration rate and production rate
juv_data$Carbon.budget <- with(juv_data, (juv_data$Respiration.rate + juv_data$Production.rate))

summary_juv <- juv_data %>% 
  group_by(trophic.group) %>% 
  summarise(n.measurements = n(), 
            mean.mass.ug = mean(Mass),
            mean.daily.respiration.rate.ug = mean(Respiration.rate),
            mean.daily.production.rate.ug = mean(Production.rate),
            mean.daily.carbon.budget.ug = mean(Carbon.budget)) %>% 
  na.omit() 
summary_juv
## # A tibble: 5 x 6
##   trophic.group n.measurements mean.mass.ug mean.daily.resp…
##   <chr>                  <int>        <dbl>            <dbl>
## 1 Bacterivores           12275       0.156           0.00365
## 2 Fungivores              1193       0.0595          0.00180
## 3 Herbivores              8477       0.122           0.00299
## 4 Omnivores                770       0.453           0.00772
## 5 Predators                516       1.22            0.0163 
## # … with 2 more variables: mean.daily.production.rate.ug <dbl>,
## #   mean.daily.carbon.budget.ug <dbl>

Number of nematodes per biome

The numbers per biome are outputted from the geospatial models.

#Import calculated totals of individuals
numbers <- Biome_Abundances %>% 
  select(-Total_Number) %>% 
  gather(key, value, -Biome) %>% 
  rename(trophic.group = key) %>% 
  rename(Total.number = value)

#Total number of individuals, per trophic group
numbers_tot <- numbers %>% 
  group_by(trophic.group) %>% 
  summarise(tot.numb.indiv = sum(Total.number))

#Total number of individuals, per trophic group, per biome
numbers_biome <- numbers%>% group_by(trophic.group, Biome) %>% summarise(tot.numb.indiv = sum(Total.number)) 

Per-trophic group sums

#Total number of individuals, per trophic group
numbers_tot_print <- numbers_tot  %>% 
  rename("Trophic group" = 1) %>% 
  rename("Number of individuals" = 2) %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))
## Warning: funs() is soft deprecated as of dplyr 0.8.0
## please use list() instead
## 
## # Before:
## funs(name = f(.)
## 
## # After: 
## list(name = ~f(.))
## This warning is displayed once per session.
numbers_tot_print
## # A tibble: 6 x 2
##   `Trophic group` `Number of individuals`
##   <chr>                             <dbl>
## 1 Bacterivores      192471000000000000000
## 2 Fungivores         63764000000000000000
## 3 Herbivores        125001000000000000000
## 4 Omnivores          39411000000000000000
## 5 Predators          19734300000000000000
## 6 Total             440381299999999983616

Relative abundances per biome

#Relative abundance per biome
abundances <- Biome_Abundances %>%
  mutate(Biome = replace(Biome, Biome == 1, "Tropical Moist Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 2, "Tropical Dry Forest")) %>%
  mutate(Biome = replace(Biome, Biome == 3, "Tropical Coniferous Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 4, "Temperate Broadleaf Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 5, "Temperate Conifer Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 6, "Boreal Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 7, "Tropical Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 8, "Temperate Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 9, "Flooded Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 10, "Montane Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 11, "Tundra")) %>%
  mutate(Biome = replace(Biome, Biome == 12, "Mediterranean Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 13, "Deserts")) %>% 
  mutate(Biome = replace(Biome, Biome == 14, "Mangroves"))

abundances_biome <- abundances %>% select(-Total_Number)
abundances_biome[,2:6] <- prop.table(as.matrix(abundances_biome[,2:6]), 1)
abundances_biome[,2:6] <- round(select_if(abundances_biome[,2:6], is.numeric), 3)
abundances_biome
##                          Biome Bacterivores Fungivores Herbivores
## 1       Tropical Moist Forests        0.443      0.120      0.249
## 2          Tropical Dry Forest        0.447      0.124      0.285
## 3  Tropical Coniferous Forests        0.423      0.135      0.235
## 4  Temperate Broadleaf Forests        0.445      0.153      0.248
## 5    Temperate Conifer Forests        0.390      0.201      0.239
## 6               Boreal Forests        0.446      0.158      0.293
## 7          Tropical Grasslands        0.371      0.140      0.354
## 8         Temperate Grasslands        0.478      0.138      0.278
## 9           Flooded Grasslands        0.410      0.141      0.308
## 10          Montane Grasslands        0.395      0.166      0.260
## 11                      Tundra        0.462      0.130      0.303
## 12       Mediterranean Forests        0.435      0.155      0.263
## 13                     Deserts        0.423      0.127      0.296
## 14                   Mangroves        0.443      0.137      0.221
##    Omnivores Predators
## 1      0.116     0.071
## 2      0.096     0.047
## 3      0.127     0.080
## 4      0.108     0.045
## 5      0.114     0.057
## 6      0.074     0.029
## 7      0.094     0.041
## 8      0.078     0.027
## 9      0.095     0.045
## 10     0.119     0.060
## 11     0.071     0.033
## 12     0.098     0.049
## 13     0.079     0.075
## 14     0.127     0.072
abundances_tot <- abundances %>% select(Biome, Total_Number) 
abundances_tot[,2] <- prop.table(as.matrix(abundances_tot[,2]), 2)
abundances_tot[,2] <- round(abundances_tot[,2], 3)
abundances_tot
##                          Biome Total_Number
## 1       Tropical Moist Forests        0.103
## 2          Tropical Dry Forest        0.016
## 3  Tropical Coniferous Forests        0.005
## 4  Temperate Broadleaf Forests        0.124
## 5    Temperate Conifer Forests        0.042
## 6               Boreal Forests        0.250
## 7          Tropical Grasslands        0.081
## 8         Temperate Grasslands        0.079
## 9           Flooded Grasslands        0.005
## 10          Montane Grasslands        0.042
## 11                      Tundra        0.138
## 12       Mediterranean Forests        0.014
## 13                     Deserts        0.098
## 14                   Mangroves        0.002
rel.abundance.trophic <- format(round(prop.table(data.frame(as.list(apply(subset(abundances, select= c(-Biome, -Total_Number)), 2, sum)))), 3), nsmall = 3)
rel.abundance.trophic
##   Bacterivores Fungivores Herbivores Omnivores Predators
## 1        0.437      0.145      0.284     0.089     0.045
abundances_final <- rbind(rel.abundance.trophic)

write.csv(abundances_final, "TableS5_rel_abundance_biome.csv")

Carbon budget

Adults

#Adults
summary_tot_adult <- merge(summary_adult, numbers_tot, by = "trophic.group")
#Total biomass in Gt: mean mass per trophic group x individuals
summary_tot_adult$total.biomass.Gt <- with(summary_tot_adult, (summary_tot_adult$mean.mass.ug * summary_tot_adult$tot.numb.indiv)/1E21)
#Total biomass in Gt C: biomass x 20% fresh mass x carbon content 52% of dry mass
summary_tot_adult$total.Cbiomass.Gt <- with(summary_tot_adult, (summary_tot_adult$total.biomass.Gt * cperw))
#Total monthly respiration in Gt C: daily rate x individuals x 31 days
summary_tot_adult$total.mon.resp.Gt <- with(summary_tot_adult, (summary_tot_adult$mean.daily.respiration.rate.ug * summary_tot_adult$tot.numb.indiv * 31)/1E21)
#Total monthly production in Gt C: daily rate x individuals x 31 days
summary_tot_adult$total.mon.prod.Gt <- with(summary_tot_adult, (summary_tot_adult$mean.daily.production.rate.ug * summary_tot_adult$tot.numb.indiv * 31)/1E21)
#Total monthly carbon budget in Gt C: sum of respiration + production
summary_tot_adult$total.C.budget.Gt <- with(summary_tot_adult, (summary_tot_adult$total.mon.prod.Gt + summary_tot_adult$total.mon.resp.Gt))

#Final output
summary_tot_adult <- summary_tot_adult %>%  
  select(trophic.group, tot.numb.indiv, total.biomass.Gt, total.Cbiomass.Gt, total.mon.resp.Gt, total.mon.prod.Gt, total.C.budget.Gt) %>% 
  rename("Trophic group" = 1) %>% 
  rename("Computed individuals" = 2) %>% 
  rename("Fresh biomass (Mt)" = 3) %>% 
  rename("Biomass (Mt C)" = 4) %>% 
  rename("Monthly respiration (Mt C)" = 5) %>% 
  rename("Monthly production (Mt C)" = 6) %>% 
  rename("Monthly carbon budget (Mt C)" = 7) %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

summary_tot_adult[,3:7] <- summary_tot_adult[,3:7] * 1000
summary_tot_adult[,-1] <- round(select_if(summary_tot_adult, is.numeric), 2)
summary_tot_adult[,2] <- formatC(summary_tot_adult[,2], format = "e", digits = 2)
summary_tot_adult
##   Trophic group Computed individuals Fresh biomass (Mt) Biomass (Mt C)
## 1  Bacterivores             1.92e+20             158.60          16.49
## 2    Fungivores             6.38e+19              23.02           2.39
## 3    Herbivores             1.25e+20             242.48          25.22
## 4     Omnivores             3.94e+19             286.69          29.82
## 5     Predators             1.97e+19              84.75           8.81
## 6         Total             4.40e+20             795.54          82.74
##   Monthly respiration (Mt C) Monthly production (Mt C)
## 1                      63.04                     31.19
## 2                      13.33                      2.18
## 3                      62.14                     20.31
## 4                      69.26                     17.56
## 5                      26.99                      6.07
## 6                     234.76                     77.32
##   Monthly carbon budget (Mt C)
## 1                        94.23
## 2                        15.51
## 3                        82.45
## 4                        86.83
## 5                        33.06
## 6                       312.08

Juveniles

#Juveniles
summary_tot_juv <- merge(summary_juv, numbers_tot, by = "trophic.group")
#Total biomass in Gt: mean mass per trophic group x individuals
summary_tot_juv$total.biomass.Gt <- with(summary_tot_juv, (summary_tot_juv$mean.mass.ug * summary_tot_juv$tot.numb.indiv)/1E21)
#Total biomass in Gt C: biomass x 20% fresh mass x carbon content 52% of dry mass
summary_tot_juv$total.Cbiomass.Gt <- with(summary_tot_juv, (summary_tot_juv$total.biomass.Gt * cperw))
#Total monthly respiration in Gt C: daily rate x individuals x 31 days
summary_tot_juv$total.mon.resp.Gt <- with(summary_tot_juv, (summary_tot_juv$mean.daily.respiration.rate.ug * summary_tot_juv$tot.numb.indiv * 31)/1E21)
#Total monthly production in Gt C: daily rate x individuals x 31 days
summary_tot_juv$total.mon.prod.Gt <- with(summary_tot_juv, (summary_tot_juv$mean.daily.production.rate.ug * summary_tot_juv$tot.numb.indiv * 31)/1E21)
# #Total monthly carbon budget in Gt C: sum of respiration + production
summary_tot_juv$total.C.budget.Gt <- with(summary_tot_juv, (summary_tot_juv$total.mon.prod.Gt + summary_tot_juv$total.mon.resp.Gt))

summary_tot_juv <- summary_tot_juv %>%  
  select(trophic.group, tot.numb.indiv, total.biomass.Gt, total.Cbiomass.Gt, total.mon.resp.Gt, total.mon.prod.Gt, total.C.budget.Gt) %>% 
  rename("Trophic group" = 1) %>% 
  rename("Computed individuals" = 2) %>% 
  rename("Fresh biomass (Mt)" = 3) %>% 
  rename("Biomass (Mt C)" = 4) %>% 
  rename("Monthly respiration (Mt C)" = 5) %>% 
  rename("Monthly production (Mt C)" = 6) %>% 
  rename("Carbon budget (Mt C)" = 7) %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

summary_tot_juv[,3:7] <- summary_tot_juv[,3:7] * 1000
summary_tot_juv[,-1] <- round(select_if(summary_tot_juv, is.numeric), 2)
summary_tot_juv[,2] <- formatC(summary_tot_juv[,2], format = "e", digits = 2)

Total carbon budget

# Final table of biomass and carbon budget

#Assuming a relative proportion of 30% adults and 70% juveniles
adult_frac <- 0.30
juv_frac <- 0.70

#Multiply adult and juvenile summaries by relative proportions
summary_tot_final <- summary_tot_juv[,1:2]
summary_tot_final[3:7] <- ((adult_frac * summary_tot_adult[3:7]) + (juv_frac * summary_tot_juv[3:7])) 
summary_tot_final[3:7] <- round(select_if(summary_tot_final[3:7], is.numeric), 2)
summary_tot_final
##   Trophic group Computed individuals Fresh biomass (Mt) Biomass (Mt C)
## 1  Bacterivores             1.92e+20              68.57           7.13
## 2    Fungivores             6.38e+19               9.56           0.99
## 3    Herbivores             1.25e+20              83.41           8.67
## 4     Omnivores             3.94e+19              98.50          10.25
## 5     Predators             1.97e+19              42.25           4.39
## 6         Total             4.40e+20             302.30          31.44
##   Monthly respiration (Mt C) Monthly production (Mt C)
## 1                      34.16                     12.22
## 2                       6.49                      0.91
## 3                      26.74                      7.01
## 4                      27.38                      6.08
## 5                      15.05                      3.00
## 6                     109.82                     29.24
##   Monthly carbon budget (Mt C)
## 1                        46.39
## 2                         7.40
## 3                        33.75
## 4                        33.46
## 5                        18.06
## 6                       139.06
#Write to file
setwd("~/Work/ETH/Projects/Nematodes/Calculations")
write.csv(summary_tot_final, "Table1_carbon_budget_2.csv")

outer bounds (i.e. ± values) in our estimates

#to calculate the spread of the data, we take the minimum and maximum predicted values across all ensembled models
numbers_tot_min <- numbers_min %>% 
  group_by(trophic.group) %>% 
  summarise(tot.numb.indiv = sum(Total.number))

numbers_tot_max <- numbers_max %>% 
  group_by(trophic.group) %>% 
  summarise(tot.numb.indiv = sum(Total.number))

numbers_tot_min_print <- numbers_tot_min %>% 
  rename("Trophic group" = 1) %>% 
  rename("Number of individuals" = 2) %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

numbers_tot_max_print <- numbers_tot_max %>% 
  rename("Trophic group" = 1) %>% 
  rename("Number of individuals" = 2) %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

variation <- numbers_tot_print[,1]
#variation$avg <- (((numbers_tot_max_print$"Number of individuals" - numbers_tot_print$"Number of individuals") + (numbers_tot_print$"Number of individuals" - numbers_tot_min_print$"Number of individuals"))/2)
variation$mean.min <- (numbers_tot_print$"Number of individuals" - numbers_tot_min_print$"Number of individuals")
#variation$max.mean <- (numbers_tot_max_print$"Number of individuals" - numbers_tot_print$"Number of individuals") 

#Adults
summary_tot_min_adult <- merge(summary_adult, numbers_tot_min, by = "trophic.group")
#Total biomass in Gt: mean mass per trophic group x individuals
summary_tot_min_adult$total.biomass.Gt <- with(summary_tot_min_adult, (summary_tot_min_adult$mean.mass.ug * summary_tot_min_adult$tot.numb.indiv)/1E21)
#Total biomass in Gt C: biomass x 20% fresh mass x carbon content 52% of dry mass
summary_tot_min_adult$total.Cbiomass.Gt <- with(summary_tot_min_adult, (summary_tot_min_adult$total.biomass.Gt * cperw))
#Total monthly respiration in Gt C: daily rate x individuals x 31 days
summary_tot_min_adult$total.mon.resp.Gt <- with(summary_tot_min_adult, (summary_tot_min_adult$mean.daily.respiration.rate.ug * summary_tot_min_adult$tot.numb.indiv * 31)/1E21)
#Total monthly production in Gt C: daily rate x individuals x 31 days
summary_tot_min_adult$total.mon.prod.Gt <- with(summary_tot_min_adult, (summary_tot_min_adult$mean.daily.production.rate.ug * summary_tot_min_adult$tot.numb.indiv * 31)/1E21)
#Total monthly carbon budget in Gt C: sum of respiration + production
summary_tot_min_adult$total.C.budget.Gt <- with(summary_tot_min_adult, (summary_tot_min_adult$total.mon.prod.Gt + summary_tot_min_adult$total.mon.resp.Gt))

#Final output
summary_tot_min_adult <- summary_tot_min_adult %>%  
  select(trophic.group, tot.numb.indiv, total.biomass.Gt, total.Cbiomass.Gt, total.mon.resp.Gt, total.mon.prod.Gt, total.C.budget.Gt) %>% 
  rename("Trophic group" = 1) %>% 
  rename("Computed individuals" = 2) %>% 
  rename("Fresh biomass (Mt)" = 3) %>% 
  rename("Biomass (Mt C)" = 4) %>% 
  rename("Monthly respiration (Mt C)" = 5) %>% 
  rename("Monthly production (Mt C)" = 6) %>% 
  rename("Monthly carbon budget (Mt C)" = 7) %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

summary_tot_min_adult[3:7] <- summary_tot_min_adult[3:7] * 1000
summary_tot_min_adult[,-1] <- round(select_if(summary_tot_min_adult, is.numeric), 3)
summary_tot_min_adult[,2] <- formatC(summary_tot_min_adult[,2], format = "e", digits = 3)

#Juveniles
summary_tot_min_juv <- merge(summary_juv, numbers_tot_min, by = "trophic.group")
#Total biomass in Gt: mean mass per trophic group x individuals
summary_tot_min_juv$total.biomass.Gt <- with(summary_tot_min_juv, (summary_tot_min_juv$mean.mass.ug * summary_tot_min_juv$tot.numb.indiv)/1E21)
#Total biomass in Gt C: biomass x 20% fresh mass x carbon content 52% of dry mass
summary_tot_min_juv$total.Cbiomass.Gt <- with(summary_tot_min_juv, (summary_tot_min_juv$total.biomass.Gt * cperw))
#Total monthly respiration in Gt C: daily rate x individuals x 31 days
summary_tot_min_juv$total.mon.resp.Gt <- with(summary_tot_min_juv, (summary_tot_min_juv$mean.daily.respiration.rate.ug * summary_tot_min_juv$tot.numb.indiv * 31)/1E21)
#Total monthly production in Gt C: daily rate x individuals x 31 days
summary_tot_min_juv$total.mon.prod.Gt <- with(summary_tot_min_juv, (summary_tot_min_juv$mean.daily.production.rate.ug * summary_tot_min_juv$tot.numb.indiv * 31)/1E21)
# #Total monthly carbon budget in Gt C: sum of respiration + production
summary_tot_min_juv$total.C.budget.Gt <- with(summary_tot_min_juv, (summary_tot_min_juv$total.mon.prod.Gt + summary_tot_min_juv$total.mon.resp.Gt))

summary_tot_min_juv <- summary_tot_min_juv %>%  
  select(trophic.group, tot.numb.indiv, total.biomass.Gt, total.Cbiomass.Gt, total.mon.resp.Gt, total.mon.prod.Gt, total.C.budget.Gt) %>% 
  rename("Trophic group" = 1) %>% 
  rename("Computed individuals" = 2) %>% 
  rename("Fresh biomass (Mt)" = 3) %>% 
  rename("Biomass (Mt C)" = 4) %>% 
  rename("Monthly respiration (Mt C)" = 5) %>% 
  rename("Monthly production (Mt C)" = 6) %>% 
  rename("Carbon budget (Mt C)" = 7) %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total")))

summary_tot_min_juv[3:7] <- summary_tot_min_juv[3:7] * 1000
summary_tot_min_juv[,-1] <- round(select_if(summary_tot_min_juv, is.numeric), 3)
summary_tot_min_juv[,2] <- formatC(summary_tot_min_juv[,2], format = "e", digits = 3)

#Multiply adult and juvenile summaries by relative proportions
summary_tot_min_final <- summary_tot_min_juv[,1:2]
summary_tot_min_final[3:7] <- ((adult_frac * summary_tot_min_adult[3:7]) + (juv_frac * summary_tot_min_juv[3:7]))

variation <- variation %>% 
  rename("Trophic group" = 1) %>% 
  rename("Number ±" = 2)
variation$'Fresh biomass ±' <- summary_tot_final$`Fresh biomass (Mt)` - summary_tot_min_final$`Fresh biomass (Mt)`
variation$"Biomass C ±" <- summary_tot_final$`Biomass (Mt C)` - summary_tot_min_final$`Biomass (Mt C)`
variation$"Respiration ±" <- summary_tot_final$`Monthly respiration (Mt C)` - summary_tot_min_final$`Monthly respiration (Mt C)`
variation$"Production ±" <- summary_tot_final$`Monthly production (Mt C)` - summary_tot_min_final$`Monthly production (Mt C)`
variation$"Budget ±" <- summary_tot_final$`Monthly carbon budget (Mt C)` - summary_tot_min_final$`Monthly carbon budget (Mt C)`
variation
## # A tibble: 6 x 7
##   `Trophic group` `Number ±` `Fresh biomass … `Biomass C ±` `Respiration ±`
##   <chr>                <dbl>            <dbl>         <dbl>           <dbl>
## 1 Bacterivores       2.08e19            7.42         0.770            3.69 
## 2 Fungivores         6.46e18            0.969        0.0966           0.661
## 3 Herbivores         1.14e19            7.59         0.785            2.43 
## 4 Omnivores          4.56e18           11.4          1.19             3.17 
## 5 Predators          3.08e18            6.59         0.682            2.35 
## 6 Total              4.63e19           34.0          3.54            12.3  
## # … with 2 more variables: `Production ±` <dbl>, `Budget ±` <dbl>
write.csv(variation, "variation.csv")

#Biomass per biome

#Total biomass per biome, per trophic group

#adult
summary_biome_adult <- merge(summary_adult, numbers_biome, by = "trophic.group")
summary_biome_adult$total.mon.resp.Gt <- with(summary_biome_adult, (summary_biome_adult$mean.daily.respiration.rate.ug * summary_biome_adult$tot.numb.indiv * 31)/1E21)
summary_biome_adult$total.biomass.Gt <- with(summary_biome_adult, (summary_biome_adult$mean.mass.ug * summary_biome_adult$tot.numb.indiv)/1E21)
summary_biome_adult$total.biomass.C.Gt <- with(summary_biome_adult, (summary_biome_adult$total.biomass.Gt * cperw * 1000))
summary_biome_adult$total.mon.prod.Gt <- with(summary_biome_adult, (summary_biome_adult$mean.daily.production.rate.ug * summary_biome_adult$tot.numb.indiv * 31)/1E21)
summary_biome_adult$total.C.budget.Gt <- with(summary_biome_adult, (summary_biome_adult$total.mon.prod.Gt + summary_biome_adult$total.mon.resp.Gt))


biomesum_adult <- summary_biome_adult %>% group_by(Biome, trophic.group) %>% summarise(sum(total.biomass.C.Gt)) %>%
  spread(key = trophic.group, value = 'sum(total.biomass.C.Gt)') 
biomesum_adult <- transform(biomesum_adult, Total.biomass=rowSums(subset(biomesum_adult, select= -Biome))) %>% 
  mutate(Biome = replace(Biome, Biome == 1, "Tropical Moist Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 2, "Tropical Dry Forest")) %>%
  mutate(Biome = replace(Biome, Biome == 3, "Tropical Coniferous Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 4, "Temperate Broadleaf Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 5, "Temperate Conifer Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 6, "Boreal Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 7, "Tropical Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 8, "Temperate Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 9, "Flooded Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 10, "Montane Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 11, "Tundra")) %>%
  mutate(Biome = replace(Biome, Biome == 12, "Mediterranean Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 13, "Deserts")) %>% 
  mutate(Biome = replace(Biome, Biome == 14, "Mangroves")) %>% 
  mutate(Biome = replace(Biome, Biome == 105, "Total"))

# juveniles
summary_biome_juv <- merge(summary_juv, numbers_biome, by = "trophic.group")
summary_biome_juv$total.mon.resp.Gt <- with(summary_biome_juv, (summary_biome_juv$mean.daily.respiration.rate.ug * summary_biome_juv$tot.numb.indiv * 31)/1E21)
summary_biome_juv$total.biomass.Gt <- with(summary_biome_juv, (summary_biome_juv$mean.mass.ug * summary_biome_juv$tot.numb.indiv)/1E21)
summary_biome_juv$total.biomass.C.Gt <- with(summary_biome_juv, (summary_biome_juv$total.biomass.Gt * cperw * 1000))
summary_biome_juv$total.mon.prod.Gt <- with(summary_biome_juv, (summary_biome_juv$mean.daily.production.rate.ug * summary_biome_juv$tot.numb.indiv * 31)/1E21)
summary_biome_juv$total.C.budget.Gt <- with(summary_biome_juv, (summary_biome_juv$total.mon.prod.Gt + summary_biome_juv$total.mon.resp.Gt))

biomesum_juv <- summary_biome_juv %>% group_by(Biome, trophic.group) %>% summarise(sum(total.biomass.C.Gt)) %>%
  spread(key = trophic.group, value = 'sum(total.biomass.C.Gt)') 
biomesum_juv <- transform(biomesum_juv, Total.biomass=rowSums(subset(biomesum_juv, select= -Biome))) %>% 
  mutate(Biome = replace(Biome, Biome == 1, "Tropical Moist Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 2, "Tropical Dry Forest")) %>%
  mutate(Biome = replace(Biome, Biome == 3, "Tropical Coniferous Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 4, "Temperate Broadleaf Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 5, "Temperate Conifer Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 6, "Boreal Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 7, "Tropical Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 8, "Temperate Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 9, "Flooded Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 10, "Montane Grasslands")) %>%
  mutate(Biome = replace(Biome, Biome == 11, "Tundra")) %>%
  mutate(Biome = replace(Biome, Biome == 12, "Mediterranean Forests")) %>%
  mutate(Biome = replace(Biome, Biome == 13, "Deserts")) %>% 
  mutate(Biome = replace(Biome, Biome == 14, "Mangroves")) %>% 
  mutate(Biome = replace(Biome, Biome == 105, "Total"))

biomesum <- biomesum_adult[,1, drop = FALSE]
biomesum[2:7] <- ((adult_frac * biomesum_adult[2:7]) + (juv_frac * biomesum_juv[2:7]))

biomesum_print <- biomesum %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total"))) %>% 
  rename("Bacterivores (Mt C)" = 2) %>% 
  rename("Fungivores (Mt C)" = 3) %>% 
  rename("Herbivores (Mt C)" = 4) %>% 
  rename("Omnivores (Mt C)" = 5) %>% 
  rename("Predators (Mt C)" = 6) %>% 
  rename("Total biomass (Mt C)" = 7) 
biomesum_print[-1] <- format(round(biomesum_print[-1], 3), nsmall = 5)
biomesum_print
##                          Biome Bacterivores (Mt C) Fungivores (Mt C)
## 1       Tropical Moist Forests             0.74500           0.08500
## 2          Tropical Dry Forest             0.11700           0.01400
## 3  Tropical Coniferous Forests             0.03600           0.00500
## 4  Temperate Broadleaf Forests             0.90400           0.13100
## 5    Temperate Conifer Forests             0.26500           0.05700
## 6               Boreal Forests             1.81500           0.27000
## 7          Tropical Grasslands             0.48900           0.07800
## 8         Temperate Grasslands             0.61900           0.07500
## 9           Flooded Grasslands             0.03500           0.00500
## 10          Montane Grasslands             0.27300           0.04800
## 11                      Tundra             1.03700           0.12300
## 12       Mediterranean Forests             0.10200           0.01500
## 13                     Deserts             0.67800           0.08600
## 14                   Mangroves             0.01500           0.00200
## 15                       Total             7.13100           0.99400
##    Herbivores (Mt C) Omnivores (Mt C) Predators (Mt C)
## 1            0.78400          1.37200          0.71900
## 2            0.13900          0.17600          0.07400
## 3            0.03700          0.07600          0.04100
## 4            0.94400          1.54100          0.55200
## 5            0.30400          0.54600          0.23200
## 6            2.23500          2.11600          0.71000
## 7            0.87400          0.86600          0.32300
## 8            0.67500          0.71000          0.21400
## 9            0.05000          0.05800          0.02400
## 10           0.33700          0.58000          0.25200
## 11           1.27700          1.12000          0.44800
## 12           0.11600          0.16200          0.06900
## 13           0.88800          0.89200          0.72400
## 14           0.01400          0.03000          0.01500
## 15           8.67500         10.24400          4.39500
##    Total biomass (Mt C)
## 1               3.70500
## 2               0.51900
## 3               0.19400
## 4               4.07200
## 5               1.40400
## 6               7.14600
## 7               2.63000
## 8               2.29200
## 9               0.17200
## 10              1.49000
## 11              4.00600
## 12              0.46500
## 13              3.26700
## 14              0.07600
## 15             31.43900
write.csv(biomesum_print, "TableS6_biomass.csv")

#Relative biomass of each trophic group

#Relative biomass of each trophic group 
rel.biomass.trophic <- format(round(prop.table(data.frame(as.list(apply(subset(biomesum, select= c(-Biome, -Total.biomass)), 2, sum)))), 3), nsmall = 3)
rel.biomass.trophic
##   Bacterivores Fungivores Herbivores Omnivores Predators
## 1        0.227      0.032      0.276     0.326     0.140

#Relative biomass of each trophic group per biome

#Relative biomass of each trophic group per biome
rel.biomass.biome <- as.data.frame(prop.table(as.matrix(biomesum[2:6]),1))
rel.biomass.biome <- cbind(biomesum$Biome, format(round(rel.biomass.biome, 3), nsmall = 3)) %>% 
  rename(Biome = 1)

rel.biomass.trophic$Biome <- "Total"
rel.biomass <- rbind(rel.biomass.biome,rel.biomass.trophic)
#relative biomass per biome
rel.biomass.sum <- as.data.frame(prop.table(as.matrix(subset(biomesum, select = Total.biomass))),1) %>% 
  bind_rows(summarise_all(., funs(if(is.numeric(.)) sum(.) else "Total"))) %>% 
  rename("Across trophic groups - fraction of total biomass" = 1)

rel.biomass.final <- cbind(rel.biomass, format(round(rel.biomass.sum, 3), nsmall = 3)) %>% 
  rename("Bacterivores - fraction of total" = 2) %>% 
  rename("Fungivores - fraction of total" = 3) %>% 
  rename("Herbivores - fraction of total" = 4) %>% 
  rename("Omnivores - fraction of total" = 5) %>% 
  rename("Predators - fraction of total" = 6)
rel.biomass.final
##                          Biome Bacterivores - fraction of total
## 1       Tropical Moist Forests                            0.201
## 2          Tropical Dry Forest                            0.225
## 3  Tropical Coniferous Forests                            0.184
## 4  Temperate Broadleaf Forests                            0.222
## 5    Temperate Conifer Forests                            0.189
## 6               Boreal Forests                            0.254
## 7          Tropical Grasslands                            0.186
## 8         Temperate Grasslands                            0.270
## 9           Flooded Grasslands                            0.206
## 10          Montane Grasslands                            0.183
## 11                      Tundra                            0.259
## 12       Mediterranean Forests                            0.220
## 13                     Deserts                            0.208
## 14                   Mangroves                            0.198
## 15                       Total                            0.227
##    Fungivores - fraction of total Herbivores - fraction of total
## 1                           0.023                          0.212
## 2                           0.026                          0.269
## 3                           0.025                          0.192
## 4                           0.032                          0.232
## 5                           0.041                          0.217
## 6                           0.038                          0.313
## 7                           0.030                          0.333
## 8                           0.033                          0.294
## 9                           0.030                          0.290
## 10                          0.032                          0.226
## 11                          0.031                          0.319
## 12                          0.033                          0.249
## 13                          0.026                          0.272
## 14                          0.026                          0.185
## 15                          0.032                          0.276
##    Omnivores - fraction of total Predators - fraction of total
## 1                          0.370                         0.194
## 2                          0.338                         0.142
## 3                          0.390                         0.209
## 4                          0.379                         0.136
## 5                          0.389                         0.165
## 6                          0.296                         0.099
## 7                          0.329                         0.123
## 8                          0.310                         0.093
## 9                          0.337                         0.137
## 10                         0.389                         0.169
## 11                         0.280                         0.112
## 12                         0.349                         0.149
## 13                         0.273                         0.222
## 14                         0.399                         0.193
## 15                         0.326                         0.140
##    Across trophic groups - fraction of total biomass
## 1                                              0.118
## 2                                              0.017
## 3                                              0.006
## 4                                              0.130
## 5                                              0.045
## 6                                              0.227
## 7                                              0.084
## 8                                              0.073
## 9                                              0.005
## 10                                             0.047
## 11                                             0.127
## 12                                             0.015
## 13                                             0.104
## 14                                             0.002
## 15                                             1.000
write.csv(rel.biomass.final, "TableS7_rel_biomass.csv")

#Comparing nematodes to humans

#total nematode biomass (freshweight) in Gt:
tot_nema_biomass <- summary_tot_final$`Fresh biomass (Gt)`[6]

#Human biomass is calculated using assumptions from Bar-On et al. 2018 PNAS 
#human population from 2015 UN data
human_pop <- 7383008.82 * 1000 #https://esa.un.org/unpd/wpp/Download/Standard/Population/

#human population versus nematode population in billions
cat("For each human on Earth,", format(round((numbers_tot_print$`Number of individuals`[6]/human_pop)/1E9, 1), nsmall = 1), "billion nematodes are present.")
## For each human on Earth, 59.6 billion nematodes are present.
#Average human weight assumed at 50 kg by Hern et al. 1999 https://link.springer.com/article/10.1023/A:1022153110536
human_weight <- 50
human_biomass <- (human_pop * human_weight)/1E12

#nematode biomass as fraction of human biomass
cat("Total soil nematode biomass represents", format(round((tot_nema_biomass/human_biomass*100), 1)),"% of total human biomass")
## Total soil nematode biomass represents  % of total human biomass