#Supplement 1 #JAGS code for Bayesian implementation of FvCB C3 photosynthesis model used in the paper Heberling, J.M. and J.D Fridley. Invaders do not require high resource levels to maintain physiological advantages in a temperate deciduous forest. Ecology #The model was parameterized using JAGS (Plummer, 2003) implementation of the Markov chain Monte Carlo mehtods, run from R 3.1.1 (R Core Team, 2014) using the R2jags package (Su & Yajima, 2015). model <- "model { #Classic FvCB C3 photo model (1980), excluding TPU limitation #Priors Vcmax.int ~ dnorm(25,0.001) #Peltier&Ibanez 2015 (Table 1)but bounded by 0: dnorm(25,0.001) I(0, ); values currently allowed to be negative Jmax.int ~ dnorm(55,0.001) #Peltier&Ibanzez2015 (Table 1)but bounded by zero: dnorm(55,0.001) I(0, ) Rd.int ~ dnorm(0,0.001) #Peltier&Ibanzez2015 (Table 1) dnorm(0,1) I(0, ) GammaStar.int ~ dnorm(4.275,10) #informative, from Patrick et al 2009; or treat as constant (comment out) or dnorm(3.74,10) or dnorm(3.86,10) depending upon treatment of M-M kinetic coefficients alpha.int ~ dnorm(0.24,100) #strong prior, low variation across species - Feng & Dietze 2013) #fixed effect flat priors (non-informative) beta.N ~ dnorm(0,0.0001) #leaf N effect on Vcmax beta.SLA ~ dnorm(0,0.0001) #SLA effect on alpha #beta.chl ~ dnorm(0,0.0001) #chl effect on alpha tau <- sigma^-2 #convert SD to precision (1/variance) sigma ~ dunif(0, 100) #uniform prior for normal SD # ind.tau.Gstar <- ind.sigma.Gstar^-2 #if including individual RE for Gstar # ind.sigma.Gstar ~ dunif(0, 100) spp.tau.Gstar <- spp.sigma.Gstar^-2 spp.sigma.Gstar ~ dunif(0, 100) # ind.tau.Rd <- ind.sigma.Rd^-2 #if including REs for Rd # ind.sigma.Rd ~ dunif(0, 100) # spp.tau.Rd <- spp.sigma.Rd^-2 # spp.sigma.Rd ~ dunif(0, 100) ind.tau.Vcmax <- ind.sigma.Vcmax^-2 ind.sigma.Vcmax ~ dunif(0, 100) spp.tau.Vcmax <- spp.sigma.Vcmax^-2 spp.sigma.Vcmax ~ dunif(0, 100) ind.tau.Jmax <- ind.sigma.Jmax^-2 ind.sigma.Jmax ~ dunif(0, 100) spp.tau.Jmax <- spp.sigma.Jmax^-2 spp.sigma.Jmax ~ dunif(0, 100) # ind.tau.alpha <- ind.sigma.alpha^-2 #if including REs for Rd # ind.sigma.alpha ~ dunif(0, 100) # spp.tau.alpha <- spp.sigma.alpha^-2 # spp.sigma.alpha ~ dunif(0, 100) for(i in 1:N) { #loop through observations (A/Ci and A/q data) Anet[i] ~ dnorm(An[i],tau) #Anet is the observed response, which takes on the predicted value An plus normal error An[i] <- min(Av[i],Aj[i]) - Rd[i] #An = min of two functions Rd[i]<- Rd.int#+ b.spp.Rd[spnumber[i]]+ b.ind.Rd[indiv[i]]#could add in random effects terms here Av[i] <- ((Vcmax[i]*(Ci_Pa[i]-GammaStar[i]))/(Ci_Pa[i]+(Kc[i]*(1+(O/Ko[i]))))) #Rubisco limited photosynthesis (early part of curve) Vcmax[i] <- Vcmax.int+b0.spp.Vcmax1[spnumber[i]]+b0.ind.Vcmax1[indiv[i]]+beta.N*(leafNarea[plot_sp[i]]) #Electron transport limited (latter part of curve) Aj[i] <- ((J[i]*(Ci_Pa[i]-GammaStar[i]))/((4*Ci_Pa[i])+(8*GammaStar[i]))) #J light dependency according to Tenhunen et al 1976 J[i]<-(alpha[i]*q[i]/(sqrt(1+(alpha[i]*alpha[i]*q[i]*q[i])/(Jmax[i]*Jmax[i])))) #could include Jmax-trait submodel; Feng&Dietze (2013) do alpha-trait submodel Jmax[i]<-Jmax.int+b0.ind.Jmax[indiv[i]]+b0.spp.Jmax[spnumber[i]]#+beta.SLA*(SLA[i]) alpha[i]<-alpha.int+beta.SLA*(SLA[i])#+b0.ind.alpha[indiv[i]]+b0.spp.alpha[spnumber[i]]#+beta.chl*(chl[i])# ##+beta.Thick*(Thick[i]) #could include other leaf structural traits or leaf chlorophyll GammaStar[i]<-GammaStar.int+b.spp.Gstar[spnumber[i]] #with species RE #get posteriors for Amax, defined at ambient Ci, saturating light Amax[i] <- min(((Jmax[i]*(40-GammaStar[i]))/((4*40)+(8*GammaStar[i]))),Vcmax[i]*(40-GammaStar[i])/(40+(40.49*(1+(O/27.84)))))#at 40 Ci_Pa (ambient) } #end loop # #random intercept for individual effect on Rd # for(i in 1:N.indiv) { # b.ind.Rd[i] ~ dnorm(0,ind.tau.Rd) # } # #random intercept for sp effect on Rd # for(i in 1:N.spp) { # b.spp.Rd[i] ~ dnorm(0,spp.tau.Rd) # } # #random intercept for individual effect on Gstar # for(i in 1:N.indiv) { # b.ind.Gstar[i] ~ dnorm(0,ind.tau.Gstar) # } #random intercept for sp effect on Gstar for(i in 1:N.spp) { b.spp.Gstar[i] ~ dnorm(0,spp.tau.Gstar) } #random intercept for plot (individual effect on Vcmax) for(i in 1:N.indiv) { b0.ind.Vcmax1[i] ~ dnorm(0,ind.tau.Vcmax) } #random intercept (spp effect on Vcmax) for(i in 1:N.spp) { b0.spp.Vcmax1[i] ~ dnorm(0,spp.tau.Vcmax) } #random intercept for plot (individual effect on Jmax) for(i in 1:N.indiv) { b0.ind.Jmax[i] ~ dnorm(0,ind.tau.Jmax) } #random intercept (spp effect on Jmax) for(i in 1:N.spp) { b0.spp.Jmax[i] ~ dnorm(0,spp.tau.Jmax) } #random intercept for plot (individual effect on alpha) # for(i in 1:N.indiv) { # b0.ind.alpha[i] ~ dnorm(0,ind.tau.alpha) # } # # #random intercept (spp effect on alpha) # for(i in 1:N.spp) { # b0.spp.alpha[i] ~ dnorm(0,spp.tau.alpha) # } #output will include overall Amax posterior based on fitted intercepts at ambient Ci Amax.int <-min(((Jmax.int*(40-4.275))/((4*40)+(8*4.275))),Vcmax.int*(40-4.275)/(40+(40.49*(1+(O/27.84))))) }" #end model