Saturating effects of N fertilization on photosynthetic capacity
Evan Perkowski
2023-12-01
This script uses data from a nitrogen-by-CO2-by-inoculation growth chamber experiment to assess nonlinear saturating effects of increasing soil nitrogen fertilization on photosynthetic capacity in uninoculated Glycine max L. (Merr) seedlings. A manuscript detailing main results from this experiment has been uploaded as a pre-print to bioRxiv (DOI: 10.1101/2023.11.30.567584).
Data for the experiment are publicly available on GitHub and a version of the repository has been archived on Zenodo (DOI: 10.5281/zenodo.10177575).
This script assumes that the NxCO2xI_ms_data repository
is the root working directory.
Load libraries and import data
# Libraries
library(tidyverse)
library(nlme)
library(emmeans)
library(ggpubr)
# Load data set. Ensure that N fertilization is numeric and filter data such that all uninoculated plants where
# nodule biomass:root biomass is greater than 0.05 g/g are removed
df <- read.csv("data/NxCO2xI_data.csv") %>%
mutate(n.trt = as.numeric(n.trt)) %>%
filter(inoc == "inoc" | (inoc == "no.inoc" & nod.root.ratio < 0.05))Vcmax25 nonlinear saturating regression model
# Create regression to explain nonlinear saturating Vcmax response to fertilization.
# Note that model uses data collected from only uninoculated individuals grown under elevated CO2
vcmax.nls.elv <- nls(formula = vcmax25 ~ a + ((b * n.trt) / (c + n.trt)),
start = list(a = 24, b = 114, c = 454),
data = subset(df, inoc == "no.inoc" & co2 == "elv"))
# Inspect nonlinear model fit for plants grown under elevated CO2
vcmax.nls.elv## Nonlinear regression model
## model: vcmax25 ~ a + ((b * n.trt)/(c + n.trt))
## data: subset(df, inoc == "no.inoc" & co2 == "elv")
## a b c
## 24.07 114.85 454.90
## residual sum-of-squares: 3897
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 4.549e-06# Create regression to explain nonlinear saturating Vcmax response to fertilization.
# Note that model uses data collected from only uninoculated individuals grown under ambient CO2
vcmax.nls.amb <- nls(formula = vcmax25 ~ a + ((b * n.trt) / (c + n.trt)),
start = list(a = 3, b = 127, c = 114),
data = subset(df, inoc == "no.inoc" & co2 == "amb"))
# Inspect nonlinear model fit for plants grown under ambient CO2
vcmax.nls.amb## Nonlinear regression model
## model: vcmax25 ~ a + ((b * n.trt)/(c + n.trt))
## data: subset(df, inoc == "no.inoc" & co2 == "amb")
## a b c
## 3.294 127.467 114.013
## residual sum-of-squares: 5809
##
## Number of iterations to convergence: 2
## Achieved convergence tolerance: 6.143e-06Create data frame to draw Vcmax25 nonlinear model predictions across range in N fertilization levels
# Elevated CO2 predictions
vcmax.nls.elv.pred <- data.frame(
emmeans(vcmax.nls.elv, ~1, "n.trt",
at = list(n.trt = seq(0, 630, 1)),
data = subset(df, inoc == "no.inoc" & co2 == "elv")))
head(vcmax.nls.elv.pred)# Ambient CO2 predictions
vcmax.nls.amb.pred <- data.frame(
emmeans(vcmax.nls.amb, ~1, "n.trt",
at = list(n.trt = seq(0, 630, 1)),
data = subset(df, inoc == "no.inoc" & co2 == "amb")))
head(vcmax.nls.elv.pred)Create Vcmax25 plot
vcmax25.plot <- ggplot(data = subset(df, inoc == "no.inoc"),
aes(x = n.trt, y = vcmax25)) +
geom_point(aes(fill = co2),
alpha = 0.75, size = 4, shape = 21) +
geom_smooth(data = vcmax.nls.elv.pred,
aes(y = emmean), color = "#b2182b",
linewidth = 2, se = FALSE) +
# geom_ribbon(data = vcmax.nls.elv.pred,
# aes(y = emmean, ymax = asymp.UCL, ymin = asymp.LCL),
# fill = "#b2182b", alpha = 0.25) +
geom_smooth(data = vcmax.nls.amb.pred,
aes(y = emmean), color = "#2166ac",
linewidth = 2, se = FALSE) +
# geom_ribbon(data = vcmax.nls.amb.pred,
# aes(y = emmean, ymax = asymp.UCL, ymin = asymp.LCL),
# fill = "#2166ac", alpha = 0.25) +
geom_text(aes(500, 5, label = (paste(expression(" y = "*frac("114.85 * x", "454.90 + x")*" + 24.07")))),
parse = TRUE, size = 3.5, color = "#b2182b") +
geom_text(aes(500, 30, label = (paste(expression("y = "*frac("127.47 * x", "114.01 + x")*" + 3.29")))),
parse = TRUE, size = 3.5, color = "#2166ac") +
scale_fill_manual(values = c("#2166ac", "#b2182b"),
labels = c(expression("Ambient CO"["2"]),
expression("Elevated CO"["2"]))) +
scale_y_continuous(limits = c(0, 160), breaks = seq(0, 160, 40)) +
labs(x = "Soil N fertilization (ppm)",
y = expression(bold(italic("V")["cmax25"]*" ("*mu*"mol m"^"-2"*" s"^"-1"*")")),
fill = expression(bold("CO"["2"]*" treatment"))) +
theme_bw(base_size = 18) +
theme(legend.title = element_text(face = "bold"),
axis.title = element_text(face = "bold"),
panel.border = element_rect(linewidth = 1.25),
legend.text.align = 0)
vcmax25.plot
Jmax25 nonlinear saturating regression model
# Create regression to explain nonlinear saturating Jmax response to fertilization
# Note that model uses data collected from only uninoculated individuals grown under elevated CO2
jmax.nls.elv <- nls(formula = jmax25 ~ a + ((b * n.trt) / (c + n.trt)),
start = list(a = 46, b = 185, c = 347),
data = subset(df, inoc == "no.inoc" & co2 == "elv"))
# Inspect model fit
jmax.nls.elv## Nonlinear regression model
## model: jmax25 ~ a + ((b * n.trt)/(c + n.trt))
## data: subset(df, inoc == "no.inoc" & co2 == "elv")
## a b c
## 46.28 185.98 347.68
## residual sum-of-squares: 12020
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 4.808e-06# Create regression to explain nonlinear saturating Jmax response to fertilization
# Note that model uses data collected from only uninoculated individuals grown under ambient CO2
jmax.nls.amb <- nls(formula = jmax25 ~ a + ((b * n.trt) / (c + n.trt)),
start = list(a = 9, b = 207, c = 96),
data = subset(df, inoc == "no.inoc" & co2 == "amb"))
# Inspect model fit
jmax.nls.amb## Nonlinear regression model
## model: jmax25 ~ a + ((b * n.trt)/(c + n.trt))
## data: subset(df, inoc == "no.inoc" & co2 == "amb")
## a b c
## 9.811 207.271 96.471
## residual sum-of-squares: 20018
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 7.626e-06Create data frame to draw Jmax25 nonlinear model predictions across range in N fertilization values
# Elevated CO2 predictions
jmax.nls.elv.pred <- data.frame(
emmeans(jmax.nls.elv, ~1, "n.trt",
at = list(n.trt = seq(0, 630, 1)),
data = subset(df, inoc == "no.inoc" & co2 == "elv")))
head(jmax.nls.elv.pred)# Ambient CO2 predictions
jmax.nls.amb.pred <- data.frame(
emmeans(jmax.nls.amb, ~1, "n.trt",
at = list(n.trt = seq(0, 630, 1)),
data = subset(df, inoc == "no.inoc" & co2 == "amb")))
head(jmax.nls.amb.pred)Create Jmax25 plot
jmax25.plot <- ggplot(data = subset(df, inoc == "no.inoc"),
aes(x = n.trt, y = jmax25)) +
geom_point(aes(fill = co2),
alpha = 0.75, size = 4, shape = 21) +
geom_smooth(data = jmax.nls.elv.pred,
aes(y = emmean), color = "#b2182b",
linewidth = 2, se = FALSE) +
# geom_ribbon(data = jmax.nls.elv.pred,
# aes(y = emmean, ymax = asymp.UCL, ymin = asymp.LCL),
# fill = "#b2182b", alpha = 0.25) +
geom_smooth(data = jmax.nls.amb.pred,
aes(y = emmean), color = "#2166ac",
linewidth = 2, se = FALSE) +
# geom_ribbon(data = jmax.nls.amb.pred,
# aes(y = emmean, ymax = asymp.UCL, ymin = asymp.LCL),
# fill = "#2166ac", alpha = 0.25) +
geom_text(aes(500, 7.5, label=(paste(expression(" y = "*frac("185.98 * x", "347.68 + x")*" + 46.28")))),
parse = TRUE, size = 3.5, color = "#b2182b") +
geom_text(aes(500, 45, label=(paste(expression("y = "*frac("207.27 * x", "96.47 + x")*" + 9.81")))),
parse = TRUE, size = 3.5, color = "#2166ac") +
scale_fill_manual(values = c("#2166ac", "#b2182b"),
labels = c(expression("Ambient CO"["2"]),
expression("Elevated CO"["2"]))) +
scale_y_continuous(limits = c(0, 240), breaks = seq(0, 240, 60)) +
labs(x = "Soil N fertilization (ppm)",
y = expression(bold(italic("J")["max25"]*" ("*mu*"mol m"^"-2"*" s"^"-1"*")")),
fill = expression(bold("CO"["2"]*" treatment"))) +
theme_bw(base_size = 18) +
theme(legend.title = element_text(face = "bold"),
axis.title = element_text(face = "bold"),
panel.border = element_rect(linewidth = 1.25),
legend.text.align = 0)
jmax25.plot
Merge Vcmax25 and Jmax25 plots into single plot
ggarrange(vcmax25.plot, jmax25.plot, common.legend = TRUE,
legend = "right", labels = c("a", "b"),
font.label = list(size = 18),
label.x = 0.2, label.y = 0.975)
Figure XX Saturating effects of increasing fertilization on (a) the maximum rate of Rubisco carboxylation and (b) the maximum rate of electron transport for RuBP regeneration in uninoculated Glycine max L. (Merr) seedlings grown under one of two CO2 concentrations (420 and 1000 ppm) across nine nitrogen fertilization levels. Blue points and nonlinear trendlines indicate individuals grown under ambient CO2, while red points and nonlinear trendlines indicate individuals grown under elevated CO2. Nonlinear regression equations are included in each panel, with colored text corresponding to each CO2 treatment. Trendlines were drawn using the ‘emmeans’ function in the ‘emmeans’ R package (Lenth 2019). Figure has been redrawn using data from analyses reported in Perkowski et al. (2023).
References
Lenth R. 2019. emmeans: estimated marginal means, aka least-squares means. https://CRAN.R-project.org/package=emmeans
Perkowski EA, E Ezekannagha, NG Smith. 2023. Nitrogen demand, supply, and acquisition strategy control plant responses to elevated CO2 at different scales. bioRxiv 2023.11.30.567584 https://doi.org/10.1101/2023.11.30.567584
Perkowski EA, E Ezekannagha, NG Smith. 2023. Dataset and analysis code for: “Nitrogen demand, supply, and acquisition strategy control plant responses to elevated CO2 at different scales” (v1.1) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.10177575