Statistical analyses
We used linear mixed effects models to test for effects of plant diversity and heterotroph removals on carbon fluxes using thenlme package (Pinheiro, Bates, DebRoy, Sarkar & R Core Team 2019) in R software (R Core Team 2019). We performed these analyses with both total fluxes and mass-specific fluxes obtained by dividing total fluxes by total plant biomass. Total and mass-specific fluxes were log transformed. Plant diversity was considered a factor variable with three distinct levels (1, 4 and 16 sp). In mixed model specification, plant diversity and heterotroph removal treatments were included as fixed effects. To account for the nested split plot design of the experiment, we included subplots (heterotroph removal treatments) nested within plots (plant diversity) as random intercepts in our mixed models. Full model equations are provided with the model results (Tables 1-6). All figures were generated using the ‘ggplot2’ package (Wickham 2009). To examine percent nitrogen content and relative abundances of different functional types we constructed generalized linear mixed effects models with a logit link using ‘glmer’ function in lme4 package (Bates et al. 2015). We performed path analysis with the most important variables from bivariate analyses to determine the relative importance of different pathways via which plant diversity and heterotroph removals influenced carbon fluxes. To do so, we constructed a structural equation model using the ‘piecewiseSEM’ package (Lefcheck et al. 2016).