Counteracting effects of soil biota on emergence and growth of herbaceous plants
Antonio J. Perea & Álvaro López-García
1/3/2022
1.- LOADING PACKAGES AND FUNCTIONS
## Packager
library(Matrix); library(coda); library(readxl); library(lme4); library(ape); library(MCMCglmm); library(patchwork); library(carData); library(car); library(emmeans); library(lattice); library(effects); library(mvtnorm); library(MASS); library(survival); library(TH.data); library(multcomp); library(MuMIn); library(glmmTMB); library(stats4); library(bbmle); library(ggplot2); library(boot); library(sjstats); library(DHARMa); library(plyr); library(lmerTest); library(performance); library(ggeffects); library(ggbeeswarm); library(Hmisc); library(ggpubr); library(tidyverse); library(rstatix); library(broom);library(tinytex);library(latexpdf)2.- DATA
2.1- Preliminary analyses
General
## Emergence
plot_emergence<-ggplot(emer, aes(wsp, pro_emer)) +coord_cartesian(ylim = c(0,1))+
stat_summary(aes(wsp, pro_emer, col=asp),fun = mean, geom = 'point', size = 4, position = position_dodge(0.5)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.5))+
geom_quasirandom(aes(wsp, pro_emer, col=asp), alpha = 2, dodge = 0.5, size = 1) +
scale_x_discrete("", labels = c("CRA_MON" = "C. monogyna","DAP_GNI" = "D. gnidium","GEN_CIN" = "G. cinerea", "JUN_OXY" = "J. oxycedrus","JUN_PHO" = "J. phoenicia","PHI_LAT" = "P. latifolia", "ROS_OFF" = "R. officinalis","RUS_ACU" = "R. aculeatus","ULE_PAR" = "U. parviflorus"))+
scale_colour_discrete(name ="Annual plant species",
breaks=c("ana_arv", "bra_ret", "fes_rub", "lot_cor", "pla_lan", "san_min", "sil_vul"),
labels=c("A. arvensis", "B. retusum", "F. rubra", "L. corniculatus", "P. lanceolata", "S. minor", "S. vulgaris")) +
labs(title="", x = " ", y = "Emergence probability") +
theme(axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12, angle = 45,vjust = 0.5, face="italic"),
axis.title.y = element_text(size=15,vjust = 2 ),
axis.title.x = element_text(size=15, vjust = -1),
axis.title.y.left = element_text(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face= "bold", size = 14),
panel.border = element_rect(colour = "black",fill=NA),
panel.background = element_blank(),
legend.text = element_text (size = 12, face="italic"))+
stat_summary(fun = mean, geom = 'point', size = 5, position = position_dodge(0.4)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.4))
####### general mortality
plot_mortality<-ggplot(emer, aes(wsp, pro_mortality)) +coord_cartesian(ylim = c(0,1))+
stat_summary(aes(wsp, pro_mortality, col=asp),fun = mean, geom = 'point', size = 4, position = position_dodge(0.5)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.5))+
geom_quasirandom(aes(wsp, pro_mortality, col=asp), alpha = 2, dodge = 0.5, size = 1) +
scale_x_discrete("", labels = c("CRA_MON" = "C. monogyna","DAP_GNI" = "D. gnidium","GEN_CIN" = "G. cinerea", "JUN_OXY" = "J. oxycedrus","JUN_PHO" = "J. phoenicia","PHI_LAT" = "P. latifolia", "ROS_OFF" = "R. officinalis","RUS_ACU" = "R. aculeatus","ULE_PAR" = "U. parviflorus"))+
scale_colour_discrete(name ="Annual plant species",
breaks=c("ana_arv", "bra_ret", "fes_rub", "lot_cor", "pla_lan", "san_min", "sil_vul"),
labels=c("A. arvensis", "B. retusum", "F. rubra", "L. corniculatus", "P. lanceolata", "S. minor", "S. vulgaris")) +
labs(title="", x = " ", y = "Mortality") +
theme(axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12, angle = 45,vjust = 0.5, face="italic"),
axis.title.y = element_text(size=15,vjust = 2 ),
axis.title.x = element_text(size=15, vjust = -1),
axis.title.y.left = element_text(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face= "bold", size = 14),
panel.border = element_rect(colour = "black",fill=NA),
panel.background = element_blank(),
legend.text = element_text (size = 12, face="italic"))+
stat_summary(fun = mean, geom = 'point', size = 5, position = position_dodge(0.4)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.4))
###### Survival after emergence
plot_surv<-ggplot(emer, aes(wsp, (1-mort_postemer))) +coord_cartesian(ylim = c(0,1))+
stat_summary(aes(wsp, (1-mort_postemer), col=asp),fun = mean, geom = 'point', size = 4, position = position_dodge(0.5)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.5))+
geom_quasirandom(aes(wsp, (1-mort_postemer), col=asp), alpha = 2, dodge = 0.5, size = 1) +
scale_x_discrete("", labels = c("CRA_MON" = "C. monogyna","DAP_GNI" = "D. gnidium","GEN_CIN" = "G. cinerea", "JUN_OXY" = "J. oxycedrus","JUN_PHO" = "J. phoenicia","PHI_LAT" = "P. latifolia", "ROS_OFF" = "R. officinalis","RUS_ACU" = "R. aculeatus","ULE_PAR" = "U. parviflorus"))+
scale_colour_discrete(name ="Annual plant species",
breaks=c("ana_arv", "bra_ret", "fes_rub", "lot_cor", "pla_lan", "san_min", "sil_vul"),
labels=c("A. arvensis", "B. retusum", "F. rubra", "L. corniculatus", "P. lanceolata", "S. minor", "S. vulgaris")) +
labs(title="", x = " ", y = "Survival after emergence") +
theme(axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12, angle = 45,vjust = 0.5, face="italic"),
axis.title.y = element_text(size=15,vjust = 2 ),
axis.title.x = element_text(size=15, vjust = -1),
axis.title.y.left = element_text(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face= "bold", size = 14),
panel.border = element_rect(colour = "black",fill=NA),
panel.background = element_blank(),
legend.text = element_text (size = 12, face="italic"))+
stat_summary(fun = mean, geom = 'point', size = 5, position = position_dodge(0.4)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.4))
##### Aboveground Biomass
plot_AG<-ggplot(emer, aes(wsp, mean_shoot)) +coord_cartesian(ylim = c(0,0.5))+
geom_quasirandom(aes(wsp, mean_shoot, col=asp), alpha = 2, dodge = 0.5, size = 1) +
stat_summary(aes(wsp, mean_shoot, col=asp),fun = mean, geom = 'point', size = 4, position = position_dodge(0.5)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.5))+
scale_x_discrete("", labels = c("CRA_MON" = "C. monogyna","DAP_GNI" = "D. gnidium","GEN_CIN" = "G. cinerea", "JUN_OXY" = "J. oxycedrus","JUN_PHO" = "J. phoenicia","PHI_LAT" = "P. latifolia", "ROS_OFF" = "R. officinalis","RUS_ACU" = "R. aculeatus","ULE_PAR" = "U. parviflorus"))+
scale_colour_discrete(name ="Annual plant species",
breaks=c("ana_arv", "bra_ret", "fes_rub", "lot_cor", "pla_lan", "san_min", "sil_vul"),
labels=c("A. arvensis", "B. retusum", "F. rubra", "L. corniculatus", "P. lanceolata", "S. minor", "S. vulgaris")) +
labs(title="", x = " ", y = "Aboveground biomass(g)") +
theme(axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12, angle = 45,vjust = 0.5, face="italic"),
axis.title.y = element_text(size=15,vjust = 2 ),
axis.title.x = element_text(size=15, vjust = -1),
axis.title.y.left = element_text(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face= "bold", size = 14),
panel.border = element_rect(colour = "black",fill=NA),
panel.background = element_blank(),
legend.text = element_text (size = 12, face="italic"))+
stat_summary(fun = mean, geom = 'point', size = 5, position = position_dodge(0.4)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.4))
##### Belowground biomass
plot_BG<-ggplot(emer, aes(wsp, mean_root)) +coord_cartesian(ylim = c(0,0.5))+
geom_quasirandom(aes(wsp, mean_root, col=asp), alpha = 2, dodge = 0.5, size = 1) +
stat_summary(aes(wsp, mean_root, col=asp),fun = mean, geom = 'point', size = 4, position = position_dodge(0.5)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.5))+
scale_x_discrete("", labels = c("CRA_MON" = "C. monogyna","DAP_GNI" = "D. gnidium","GEN_CIN" = "G. cinerea", "JUN_OXY" = "J. oxycedrus","JUN_PHO" = "J. phoenicia","PHI_LAT" = "P. latifolia", "ROS_OFF" = "R. officinalis","RUS_ACU" = "R. aculeatus","ULE_PAR" = "U. parviflorus"))+
scale_colour_discrete(name ="Annual plant species",
breaks=c("ana_arv", "bra_ret", "fes_rub", "lot_cor", "pla_lan", "san_min", "sil_vul"),
labels=c("A. arvensis", "B. retusum", "F. rubra", "L. corniculatus", "P. lanceolata", "S. minor", "S. vulgaris")) +
labs(title="", x = " ", y = "Belowground biomass (g)") +
theme(axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12, angle = 45,vjust = 0.5, face="italic"),
axis.title.y = element_text(size=15,vjust = 2 ),
axis.title.x = element_text(size=15, vjust = -1),
axis.title.y.left = element_text(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face= "bold", size = 14),
panel.border = element_rect(colour = "black",fill=NA),
panel.background = element_blank(),
legend.text = element_text (size = 12, face="italic"))+
stat_summary(fun = mean, geom = 'point', size = 5, position = position_dodge(0.4)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.4))
##### Ratio A/B biomass
plot_ABG<-ggplot(emer, aes(wsp, ratio_mean_shoot)) +coord_cartesian(ylim = c(0,2.5))+
geom_quasirandom(aes(wsp, ratio_mean_shoot, col=asp), alpha = 2, dodge = 0.5, size = 1) +
stat_summary(aes(wsp, ratio_mean_shoot, col=asp),fun = mean, geom = 'point', size = 4, position = position_dodge(0.5)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.5))+
scale_x_discrete("", labels = c("CRA_MON" = "C. monogyna","DAP_GNI" = "D. gnidium","GEN_CIN" = "G. cinerea", "JUN_OXY" = "J. oxycedrus","JUN_PHO" = "J. phoenicia","PHI_LAT" = "P. latifolia", "ROS_OFF" = "R. officinalis","RUS_ACU" = "R. aculeatus","ULE_PAR" = "U. parviflorus"))+
scale_colour_discrete(name ="Annual plant species",
breaks=c("ana_arv", "bra_ret", "fes_rub", "lot_cor", "pla_lan", "san_min", "sil_vul"),
labels=c("A. arvensis", "B. retusum", "F. rubra", "L. corniculatus", "P. lanceolata", "S. minor", "S. vulgaris")) +
labs( x = " ", y = "Ratio A/B biomass") +
theme(axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12, angle = 45,vjust = 0.5, face="italic"),
axis.title.y = element_text(size=15,vjust = 2 ),
axis.title.x = element_text(size=15, vjust = -1),
axis.title.y.left = element_text(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face= "bold", size = 14),
panel.border = element_rect(colour = "black",fill=NA),
panel.background = element_blank(),
legend.text = element_text (size = 12, face="italic"))+
stat_summary(fun = mean, geom = 'point', size = 5, position = position_dodge(0.4)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.4))
### Total biomass ###
plot_TB<- ggplot(emer, aes(wsp, sum_shoot_root)) +coord_cartesian(ylim = c(0,1))+
geom_quasirandom(aes(wsp, sum_shoot_root, col=asp), alpha = 2, dodge = 0.5, size = 1) +
stat_summary(aes(wsp, sum_shoot_root, col=asp),fun = mean, geom = 'point', size = 4, position = position_dodge(0.5)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.5))+
scale_x_discrete("", labels = c("CRA_MON" = "C. monogyna","DAP_GNI" = "D. gnidium","GEN_CIN" = "G. cinerea", "JUN_OXY" = "J. oxycedrus","JUN_PHO" = "J. phoenicia","PHI_LAT" = "P. latifolia", "ROS_OFF" = "R. officinalis","RUS_ACU" = "R. aculeatus","ULE_PAR" = "U. parviflorus"))+
scale_colour_discrete(name ="Annual plant species",
breaks=c("ana_arv", "bra_ret", "fes_rub", "lot_cor", "pla_lan", "san_min", "sil_vul"),
labels=c("A. arvensis", "B. retusum", "F. rubra", "L. corniculatus", "P. lanceolata", "S. minor", "S. vulgaris")) +
labs(x = " ", y = "Total Biomass (g)") +
theme(axis.text.y = element_text(size=12),
axis.text.x = element_text(size=12, angle = 45,vjust = 0.5, face="italic"),
axis.title.y = element_text(size=15,vjust = 2 ),
axis.title.x = element_text(size=15, vjust = -1),
axis.title.y.left = element_text(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5, face= "bold", size = 14),
panel.border = element_rect(colour = "black",fill=NA),
panel.background = element_blank(),
legend.text = element_text (size = 12, face="italic"))+
stat_summary(fun = mean, geom = 'point', size = 5, position = position_dodge(0.4)) + # apply mean function (fun = mean) (median or other functions work too)
stat_summary(fun.data = mean_cl_normal, geom = 'errorbar', width = 0, size = 1, position = position_dodge(0.4))
plot_emergence
figure2<-ggarrange(plot_emergence, plot_mortality, plot_surv, nrow = 1, ncol = 3, align = "h", common.legend = TRUE, legend = "top", hjust = 0.5)
annotate_figure(figure2, bottom = text_grob("Soils of canopy species",color = "black", face = "bold", size = 16))
Success parameters (SI)
3.- DATA ANALYSES
3.1.- AMF Presence/absence
Model will be adjusted for random tray, asp and wsp.
3.1.1 EMERGENCE
res<- glmmTMB(pro_emer ~ mo*amf + (1|wsp) + (1|asp) + (1|tray), na.action = na.omit,
weights = n_ini, data = emer,family = betabinomial (link = "logit"))
testZeroInflation(res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.1177, p-value = 0.792
## alternative hypothesis: two.sidedsimulationOutput<-simulateResiduals(fittedModel=res)
plot(simulationOutput, add=TRUE)
summary(res) ## Family: betabinomial ( logit )
## Formula: pro_emer ~ mo * amf + (1 | wsp) + (1 | asp) + (1 | tray)
## Data: emer
## Weights: n_ini
##
## AIC BIC logLik deviance df.resid
## 2703.2 2742.5 -1343.6 2687.2 1000
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.02159 0.1469
## asp (Intercept) 2.41098 1.5527
## tray (Intercept) 0.14587 0.3819
## Number of obs: 1008, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion parameter for betabinomial family (): 6.68
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.30933 0.60152 -0.514 0.6071
## mo1 0.36894 0.14499 2.545 0.0109 *
## amf1 0.07913 0.13295 0.595 0.5517
## mo1:amf1 -0.07459 0.19411 -0.384 0.7008
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.2 PERFORMANCE PARAMETERS
Aboveground biomass
res12 <- glmmTMB(mean_shoot ~ mo*amf + (1|wsp) + (1|asp) + (1|tray) ,na.action = na.omit,
data = emer, weights = n_max_corr,family = gaussian(link= "log"))
testZeroInflation(res12)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sidedcheck_distribution(res12)## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 59%
## beta 16%
## gamma 16%
##
## Predicted Distribution of Response
##
## Distribution Probability
## beta 47%
## tweedie 22%
## gamma 12%simulationOutput<-simulateResiduals(fittedModel=res12)
plot(simulationOutput, add=TRUE)
summary(res12)## Family: gaussian ( log )
## Formula: mean_shoot ~ mo * amf + (1 | wsp) + (1 | asp) + (1 | tray)
## Data: emer
## Weights: n_max_corr
##
## AIC BIC logLik deviance df.resid
## -9975.6 -9939.0 4995.8 -9991.6 704
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.069744 0.26409
## asp (Intercept) 0.211911 0.46034
## tray (Intercept) 0.013561 0.11645
## Residual 0.001016 0.03187
## Number of obs: 712, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion estimate for gaussian family (sigma^2): 0.00102
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.98807 0.19930 -14.993 < 2e-16 ***
## mo1 -0.08181 0.02907 -2.814 0.00489 **
## amf1 0.11711 0.02648 4.422 9.78e-06 ***
## mo1:amf1 -0.01036 0.03811 -0.272 0.78565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Belowground biomass
res21 <- glmmTMB(sqrt(mean_root) ~ mo*amf + (1|wsp) + (1|asp) + (1|tray) ,
data = emer, weights = n_max_corr, family = gaussian(link= "identity"))
simulationOutput<-simulateResiduals(fittedModel=res21)
check_distribution(res21)## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 59%
## beta 12%
## weibull 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## exponential 47%
## tweedie 34%
## beta 16%plot(simulationOutput, add=TRUE)
summary(res21)## Family: gaussian ( identity )
## Formula:
## sqrt(mean_root) ~ mo * amf + (1 | wsp) + (1 | asp) + (1 | tray)
## Data: emer
## Weights: n_max_corr
##
## AIC BIC logLik deviance df.resid
## -4900.6 -4864.1 2458.3 -4916.6 694
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.0030525 0.05525
## asp (Intercept) 0.0081280 0.09016
## tray (Intercept) 0.0005628 0.02372
## Residual 0.0076539 0.08749
## Number of obs: 702, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion estimate for gaussian family (sigma^2): 0.00765
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.264255 0.039318 6.721 1.81e-11 ***
## mo1 -0.010685 0.006011 -1.778 0.0755 .
## amf1 0.008784 0.005668 1.550 0.1212
## mo1:amf1 0.005656 0.008090 0.699 0.4845
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova (res21, type = c("III"), test.statistic = c("Chisq"), vcov. = vcov(res21, complete=F))## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: sqrt(mean_root)
## Chisq Df Pr(>Chisq)
## (Intercept) 45.1709 1 1.806e-11 ***
## mo 3.1595 1 0.07549 .
## amf 2.4012 1 0.12125
## mo:amf 0.4888 1 0.48449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ratio Biomass
res31<-glmmTMB(sqrt(ratio_mean_shoot) ~ mo*amf + (1|wsp) + (1|asp) + (1|tray), na.action = na.omit,
data = emer, weights = n_max_corr,family = gaussian(link = "log"))
simulationOutput<-simulateResiduals(fittedModel=res31)
plot(simulationOutput, add=TRUE)
summary(res31)## Family: gaussian ( log )
## Formula:
## sqrt(ratio_mean_shoot) ~ mo * amf + (1 | wsp) + (1 | asp) + (1 | tray)
## Data: emer
## Weights: n_max_corr
##
## AIC BIC logLik deviance df.resid
## -1975.5 -1939.0 995.7 -1991.5 693
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.001361 0.03689
## asp (Intercept) 0.008110 0.09006
## tray (Intercept) 0.001637 0.04046
## Residual 0.025278 0.15899
## Number of obs: 701, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion estimate for gaussian family (sigma^2): 0.0253
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.164965 0.038437 -4.292 1.77e-05 ***
## mo1 -0.020025 0.012152 -1.648 0.099380 .
## amf1 0.041505 0.011328 3.664 0.000248 ***
## mo1:amf1 -0.006981 0.016250 -0.430 0.667488
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.3 SURVIVAL PARAMETERS (Exploratory analyses, not included in the manuscript)
Global Survival (nlev compare to n_ini)
res6<- glmmTMB(pro_mortality ~ mo*amf + (1|wsp) + (1|asp) + (1|tray) + (1|idpot) ,na.action = na.omit,
weights = n_ini, ziformula = ~tray, data = emer,family = betabinomial(link="logit"))
simulationOutput<-simulateResiduals(fittedModel=res6)
plot(simulationOutput, add=TRUE)
summary(res6) ## Family: betabinomial ( logit )
## Formula:
## pro_mortality ~ mo * amf + (1 | wsp) + (1 | asp) + (1 | tray) +
## (1 | idpot)
## Zero inflation: ~tray
## Data: emer
## Weights: n_ini
##
## AIC BIC logLik deviance df.resid
## 2634.7 2796.9 -1284.3 2568.7 975
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.01328 0.1152
## asp (Intercept) 2.85047 1.6883
## tray (Intercept) 0.09725 0.3119
## idpot (Intercept) 0.08070 0.2841
## Number of obs: 1008, groups: wsp, 9; asp, 7; tray, 24; idpot, 1008
##
## Dispersion parameter for betabinomial family (): 12.3
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.728120 0.649301 1.121 0.26212
## mo1 -0.381649 0.137403 -2.778 0.00548 **
## amf1 -0.010606 0.127795 -0.083 0.93386
## mo1:amf1 -0.008924 0.185883 -0.048 0.96171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.8702 7952.8341 -0.003 0.998
## tray2 -1.6963 20198.1757 0.000 1.000
## tray3 -1.4888 18450.4458 0.000 1.000
## tray4 -2.1937 21984.3380 0.000 1.000
## tray5 2.8199 9165.6073 0.000 1.000
## tray6 -2.4649 24485.2144 0.000 1.000
## tray7 -2.9222 25769.8345 0.000 1.000
## tray8 -1.5747 19993.5376 0.000 1.000
## tray9 -2.2425 22739.2074 0.000 1.000
## tray10 -3.0799 25752.9497 0.000 1.000
## tray11 -1.7714 21006.3900 0.000 1.000
## tray12 -1.1745 18640.5542 0.000 1.000
## tray13 -2.3990 23312.4991 0.000 1.000
## tray14 -1.5365 19188.0001 0.000 1.000
## tray15 -1.4182 18139.5088 0.000 1.000
## tray16 -1.8207 20895.8877 0.000 1.000
## tray17 -2.2013 23106.9724 0.000 1.000
## tray18 0.3260 14937.8291 0.000 1.000
## tray19 -2.0317 21397.0560 0.000 1.000
## tray20 -1.7978 20467.0377 0.000 1.000
## tray21 -0.5382 15563.6519 0.000 1.000
## tray22 -2.2195 22450.6500 0.000 1.000
## tray23 17.8259 7952.8341 0.002 0.998
## tray24 -2.6106 24930.0742 0.000 1.000Survival after emergency
res7<- glmmTMB((1-mort_postemer) ~ mo*amf + (1|wsp) + (1|asp) + (1|tray) + (1|idpot) ,na.action = na.omit,
weights = n_ini, ziformula = ~tray, data =emer, family = betabinomial(link="logit"))
simulationOutput<-simulateResiduals(fittedModel=res7)
plot(simulationOutput, add=TRUE)
summary(res7) ## Family: betabinomial ( logit )
## Formula:
## (1 - mort_postemer) ~ mo * amf + (1 | wsp) + (1 | asp) + (1 |
## tray) + (1 | idpot)
## Zero inflation: ~tray
## Data: emer
## Weights: n_ini
##
## AIC BIC logLik deviance df.resid
## 1431.3 1584.4 -682.7 1365.3 731
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 3.721e-08 0.0001929
## asp (Intercept) 4.117e-01 0.6416637
## tray (Intercept) 4.609e-08 0.0002147
## idpot (Intercept) 3.507e-02 0.1872692
## Number of obs: 764, groups: wsp, 9; asp, 7; tray, 24; idpot, 764
##
## Dispersion parameter for betabinomial family (): 5.66
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.06126 0.33071 9.257 <2e-16 ***
## mo1 -0.07791 0.23292 -0.334 0.738
## amf1 -0.15100 0.23028 -0.656 0.512
## mo1:amf1 0.02874 0.32305 0.089 0.929
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.890e+00 4.810e-01 -3.929 8.54e-05 ***
## tray2 -2.274e+01 3.875e+04 -0.001 0.9995
## tray3 -1.529e+00 1.141e+00 -1.341 0.1800
## tray4 1.195e+00 5.890e-01 2.030 0.0424 *
## tray5 -3.101e-01 7.769e-01 -0.399 0.6898
## tray6 -2.275e+01 3.850e+04 -0.001 0.9995
## tray7 -6.790e-01 8.796e-01 -0.772 0.4401
## tray8 -7.885e-01 8.775e-01 -0.899 0.3689
## tray9 -8.888e-01 8.767e-01 -1.014 0.3107
## tray10 -1.383e+00 1.140e+00 -1.213 0.2251
## tray11 -2.268e+01 3.945e+04 -0.001 0.9995
## tray12 -2.263e+01 3.984e+04 -0.001 0.9995
## tray13 6.241e-02 6.817e-01 0.092 0.9271
## tray14 -7.548e-01 8.790e-01 -0.859 0.3905
## tray15 -1.484e+00 1.130e+00 -1.312 0.1894
## tray16 -8.579e-01 8.778e-01 -0.977 0.3284
## tray17 -1.585e+00 1.132e+00 -1.400 0.1614
## tray18 -1.485e+00 1.132e+00 -1.312 0.1895
## tray19 -8.249e-01 8.783e-01 -0.939 0.3476
## tray20 -8.259e-01 8.789e-01 -0.940 0.3473
## tray21 -1.493e+00 1.139e+00 -1.311 0.1899
## tray22 -2.274e+01 3.875e+04 -0.001 0.9995
## tray23 -8.581e-01 8.779e-01 -0.977 0.3283
## tray24 -1.455e+00 1.137e+00 -1.280 0.2006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2 AMF RICHNESS
3.2.1 EMERGENCE
res_0<- glmmTMB(pro_emer ~ mo*rich_amf + (1|wsp) + (1|asp) + (1|tray), na.action = na.omit,
weights = n_ini, ziformula = ~tray, data = emer,family = betabinomial(link = "logit"))
simulationOutput<-simulateResiduals(fittedModel=res_0)
plot(simulationOutput, add=TRUE)
summary(res_0) ## Family: betabinomial ( logit )
## Formula: pro_emer ~ mo * rich_amf + (1 | wsp) + (1 | asp) + (1 | tray)
## Zero inflation: ~tray
## Data: emer
## Weights: n_ini
##
## AIC BIC logLik deviance df.resid
## 2739.8 2897.1 -1337.9 2675.8 976
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.01958 0.1399
## asp (Intercept) 2.47495 1.5732
## tray (Intercept) 0.15078 0.3883
## Number of obs: 1008, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion parameter for betabinomial family (): 7.42
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.226038 0.608593 -0.371 0.7103
## mo1 0.340083 0.140880 2.414 0.0158 *
## rich_amf -0.003209 0.014891 -0.216 0.8294
## mo1:rich_amf -0.006906 0.021623 -0.319 0.7494
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.203 2528.329 -0.008 0.994
## tray2 -4.658 29083.057 0.000 1.000
## tray3 -4.440 28151.249 0.000 1.000
## tray4 -4.101 23029.279 0.000 1.000
## tray5 -4.685 29351.534 0.000 1.000
## tray6 -4.561 28458.235 0.000 1.000
## tray7 16.094 2528.329 0.006 0.995
## tray8 -4.476 28168.167 0.000 1.000
## tray9 15.756 2528.329 0.006 0.995
## tray10 15.738 2528.329 0.006 0.995
## tray11 -4.264 27903.066 0.000 1.000
## tray12 -3.689 24341.296 0.000 1.000
## tray13 -4.708 29047.802 0.000 1.000
## tray14 15.888 2528.329 0.006 0.995
## tray15 -4.377 28581.406 0.000 1.000
## tray16 -4.435 28443.048 0.000 1.000
## tray17 -4.621 30574.076 0.000 1.000
## tray18 16.268 2528.329 0.006 0.995
## tray19 -4.501 29746.746 0.000 1.000
## tray20 -4.426 27914.634 0.000 1.000
## tray21 -4.331 28609.184 0.000 1.000
## tray22 -4.560 28635.467 0.000 1.000
## tray23 -4.637 29409.851 0.000 1.000
## tray24 -4.142 27218.914 0.000 1.0003.2.2 PERFORMANCE PARAMETERS
Note that, at using AMF richness model adjustment get worse than in previous cases (see DHARMAs), so in this case results must be interpreted carefully. Several models were tested in order to obtain the best fit, so here we show these ones with the best fit. Likely, because we introduced AMF richness as a gradient (i.,e. numerical), we do not have lots of intermediate values between the highest and the lowest, which promote sightly deviations from the data expected by the model and some outliers. Next approaches should consider it.
Aboveground biomass
res12_1 <- glmmTMB(mean_shoot ~ rich_amf*mo + (1|wsp) + (1|asp) + (1|tray) ,na.action = na.omit,
data = emer, weights = n_max_corr, family = gaussian(link = "identity"))
simulationOutput<-simulateResiduals(fittedModel=res12_1)
plot(simulationOutput, add=TRUE)
summary(res12_1)## Family: gaussian ( identity )
## Formula:
## mean_shoot ~ rich_amf * mo + (1 | wsp) + (1 | asp) + (1 | tray)
## Data: emer
## Weights: n_max_corr
##
## AIC BIC logLik deviance df.resid
## -10007.4 -9970.8 5011.7 -10023.4 704
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 4.342e-04 0.02084
## asp (Intercept) 9.576e-04 0.03094
## tray (Intercept) 6.675e-05 0.00817
## Residual 1.003e-03 0.03167
## Number of obs: 712, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion estimate for gaussian family (sigma^2): 0.001
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.744e-02 1.380e-02 4.162 3.16e-05 ***
## rich_amf 1.070e-03 2.301e-04 4.650 3.32e-06 ***
## mo1 -7.471e-03 2.097e-03 -3.563 0.000367 ***
## rich_amf:mo1 -1.717e-07 3.230e-04 -0.001 0.999576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Belowground biomass
res21_1 <- glmmTMB(sqrt(mean_root) ~ rich_amf*mo + (1|wsp) + (1|asp) + (1|tray) ,na.action = na.omit,
data = emer, weights = n_max_corr,family = gaussian(link = "log"))
simulationOutput<-simulateResiduals(fittedModel=res21_1)
plot(simulationOutput, add=TRUE)
summary(res21_1)## Family: gaussian ( log )
## Formula:
## sqrt(mean_root) ~ rich_amf * mo + (1 | wsp) + (1 | asp) + (1 | tray)
## Data: emer
## Weights: n_max_corr
##
## AIC BIC logLik deviance df.resid
## -4887.6 -4851.1 2451.8 -4903.6 694
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.025959 0.16112
## asp (Intercept) 0.095890 0.30966
## tray (Intercept) 0.007054 0.08399
## Residual 0.007692 0.08771
## Number of obs: 702, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion estimate for gaussian family (sigma^2): 0.00769
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4079308 0.1312969 -10.723 <2e-16 ***
## rich_amf 0.0037332 0.0020697 1.804 0.0713 .
## mo1 0.0049551 0.0188372 0.263 0.7925
## rich_amf:mo1 -0.0006561 0.0028344 -0.231 0.8170
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ratio A/B biomass
res31_1 <- glmmTMB(ratio_mean_shoot ~ rich_amf*mo + (1|wsp) + (1|asp) + (1|tray),
data = emer,weights = n_max_corr, family = gaussian(link = "identity"))
simulationOutput<-simulateResiduals(fittedModel=res31_1)
plot(simulationOutput, add=TRUE)
summary(res31_1)## Family: gaussian ( identity )
## Formula:
## ratio_mean_shoot ~ rich_amf * mo + (1 | wsp) + (1 | asp) + (1 | tray)
## Data: emer
## Weights: n_max_corr
##
## AIC BIC logLik deviance df.resid
## 1099.9 1136.3 -541.9 1083.9 693
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.002803 0.05294
## asp (Intercept) 0.016660 0.12907
## tray (Intercept) 0.003613 0.06011
## Residual 0.087831 0.29636
## Number of obs: 701, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion estimate for gaussian family (sigma^2): 0.0878
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.775209 0.055468 13.976 < 2e-16 ***
## rich_amf 0.009494 0.002141 4.434 9.27e-06 ***
## mo1 -0.027697 0.019327 -1.433 0.152
## rich_amf:mo1 -0.004114 0.002990 -1.376 0.169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2.3 SURVIVAL PARAMETERS (Data exploration, not included in the manuscript)
Global Survival (Model adjustment: Good ; nlev compare to n_ini)
res6_1 <- glmmTMB(pro_mortality ~ rich_amf*mo + (1|wsp) + (1|asp) + (1|tray) ,na.action = na.omit,
data = emer,weights = n_ini,family = betabinomial(link = "logit"))
simulationOutput<-simulateResiduals(fittedModel=res6_1 )
plot(simulationOutput, add=TRUE)
summary(res6_1 )## Family: betabinomial ( logit )
## Formula: pro_mortality ~ rich_amf * mo + (1 | wsp) + (1 | asp) + (1 |
## tray)
## Data: emer
## Weights: n_ini
##
## AIC BIC logLik deviance df.resid
## 2589.4 2628.8 -1286.7 2573.4 1000
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 0.01797 0.1340
## asp (Intercept) 2.73668 1.6543
## tray (Intercept) 0.09467 0.3077
## Number of obs: 1008, groups: wsp, 9; asp, 7; tray, 24
##
## Dispersion parameter for betabinomial family (): 9.79
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.671994 0.636241 1.056 0.29088
## rich_amf 0.007121 0.014000 0.509 0.61100
## mo1 -0.379518 0.131203 -2.893 0.00382 **
## rich_amf:mo1 0.001019 0.020264 0.050 0.95989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Similar to the previous analyses, microbiota decrease significantly the probability of mortality. Something possible due to enhancing germination, which may reduce the action of seed pathogens.
Survival after emergency
res7_1 <- glmmTMB((1-mort_postemer) ~ rich_amf*mo + (1|wsp) + (1|asp) + (1|tray) ,na.action = na.omit,
data = emer[which(emer$amf== 1),],weights = n_max_corr,family = betabinomial(link = "logit"))
simulationOutput<-simulateResiduals(fittedModel=res7_1 )
plot(simulationOutput, add=TRUE)
summary(res7_1 )## Family: betabinomial ( logit )
## Formula:
## (1 - mort_postemer) ~ rich_amf * mo + (1 | wsp) + (1 | asp) + (1 | tray)
## Data: emer[which(emer$amf == 1), ]
## Weights: n_max_corr
##
## AIC BIC logLik deviance df.resid
## 565.3 596.7 -274.7 549.3 368
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## wsp (Intercept) 1.921e-09 4.383e-05
## asp (Intercept) 1.119e+00 1.058e+00
## tray (Intercept) 3.055e-01 5.527e-01
## Number of obs: 376, groups: wsp, 9; asp, 7; tray, 21
##
## Dispersion parameter for betabinomial family (): 40.3
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.63578 0.70786 2.311 0.0208 *
## rich_amf 0.01345 0.06566 0.205 0.8377
## mo1 1.20086 0.79679 1.507 0.1318
## rich_amf:mo1 -0.10981 0.09082 -1.209 0.2266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14. PLOTS
4.1 PLOTS results AMF 1-0
EMMEANS: Estimated marginal means (Least-squares means). Compute estimated marginal means (EMMs) for specified factors or factor combinations in a linear model; and optionally, comparisons or contrasts among them. EMMs are also known as least-squares means.
4.1.1 Emergence
4.1.2 Aboveground Biomass
4.1.3 Resource allocation
4.2 PLOTS results AMF-RICHNESS
4.2.1 Aboveground biomass
4.2.2 Resource allocation
