Dryland DIRT

Summary

In 2023, surface soils beneath sagebrush had lower total organic carbon under litter removal than under control or doubled litter (Fig. 1). Treatment responses through time also varied with depth (Fig. 2). Beneath sagebrush, particulate organic carbon was approximately 36% lower under litter removal than control, while doubled litter and control were similar. MAOC showed no clear treatment differences, although estimates were modestly higher under doubled litter (Fig. 3).

These patterns are consistent with litter inputs helping maintain surface particulate carbon, without a proportional increase when litter inputs are doubled. MicroResp can test whether litter history changes microbial respiration and substrate responses even where carbon concentrations remain similar. It cannot, by itself, determine where the additional litter carbon went.

Organic Carbon

Do litter treatments affect carbon differently across years, microsites, and depths?

These data have continuous response variable for Organic Carbon and the following 4 predictors:

Treatment: 3 levels

Year: 2 levels (included only 2013 and 2023 in this analysis)

Microsite: 2 levels

Depth: 2 levels (included only 0-5 cm & 5-10 cm in this analysis)

With Block and Plot added as random intercepts.

1| Data

I grabbed data from 2013 and 2023 only. I subset the 0-5 cm and 5-10 cookies.

dryland <- read.csv("drylanddetritus.csv")
#keep 2013 and 2023 at 0–5 and 5–10 cm

dryland <- subset(dryland, Year %in% c(2013, 2023) & Depth %in% c("0-5", "5-10"))

#convert carbon from percent to mgC/gsoil
dryland$mgC_g <- dryland$OC*10

#turn plot and year into factors
dryland$Plot <- factor(dryland$Plot)
dryland$Year <- factor(dryland$Year)
dryland$Block <- factor(dryland$Block)

#check out the response variable
ggplot(dryland, aes(mgC_g, color=mgC_g))+
  geom_density(color="plum", weight=2)+
  labs(x="Organic Carbon (mg C/g Soil)", y="Density") #values have a right skew. I will either log transform them OR gamma.

2| Explore

ggplot(dryland, aes(Year, mgC_g, fill = Treatment)) +
  geom_boxplot(position = position_dodge(width = 0.8), width = 0.7, na.rm=T)+
  facet_grid(Depth ~ Veg) +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.15))) +
  theme_classic() +
  labs(x = "Year", y = "Organic carbon (mg C/g soil)", fill = "Treatment")+
  scale_fill_manual(values=c("#d8aee5ff", "#adc3b7", "#a69e00ff"))

3| Model

Soil organic carbon concentration was analyzed using a generalized linear mixed model with a Gamma distribution and log link, fitted using the glmmTMB package in R. The Gamma distribution accommodated positive continuous concentrations with variance increasing with the expected mean. Treatment, year, microsite, depth, and their interactions were included to evaluate whether treatment responses through time differed among microsites and depths. Random intercepts for block and plot accounted for spatial grouping and repeated measurements within plots. Model adequacy was assessed using Pearson residual plots, random-intercept Q–Q plots, and convergence diagnostics.

OC_model <- glmmTMB(mgC_g ~ Treatment * Veg * Year * Depth + (1 | Block) + (1 | Plot), data = dryland, family = Gamma(link = "log"))

4| Check fit

OC_check <- model.frame(OC_model)
OC_check$fitted <- fitted(OC_model)
OC_check$residual <- residuals(OC_model, type = "pearson")

plot(residual ~ fitted, data = OC_check,
     xlab = "Fitted carbon", ylab = "Pearson residuals")
abline(h = 0, lty = 2)

qqnorm(ranef(OC_model)$cond$Plot[, 1])
qqline(ranef(OC_model)$cond$Plot[, 1])

qqnorm(ranef(OC_model)$cond$Block[, 1])
qqline(ranef(OC_model)$cond$Block[, 1])

5| Compare

OC_groups <- emmeans(OC_model, ~ Treatment | Year * Veg * Depth)
pairs(OC_groups, adjust = "tukey", type = "response")
## Year = 2013, Veg = IS, Depth = 0-5:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 0.995 0.134 Inf    1  -0.040  0.9991
##  CTL / NL 1.176 0.163 Inf    1   1.173  0.4695
##  DL / NL  1.183 0.164 Inf    1   1.212  0.4462
## 
## Year = 2023, Veg = IS, Depth = 0-5:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 0.827 0.112 Inf    1  -1.408  0.3367
##  CTL / NL 1.014 0.137 Inf    1   0.100  0.9945
##  DL / NL  1.226 0.166 Inf    1   1.509  0.2867
## 
## Year = 2013, Veg = SB, Depth = 0-5:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 1.079 0.146 Inf    1   0.560  0.8414
##  CTL / NL 1.178 0.159 Inf    1   1.215  0.4442
##  DL / NL  1.092 0.148 Inf    1   0.655  0.7895
## 
## Year = 2023, Veg = SB, Depth = 0-5:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 0.869 0.117 Inf    1  -1.043  0.5495
##  CTL / NL 1.568 0.212 Inf    1   3.331  0.0025
##  DL / NL  1.806 0.244 Inf    1   4.374 <0.0001
## 
## Year = 2013, Veg = IS, Depth = 5-10:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 1.047 0.141 Inf    1   0.337  0.9394
##  CTL / NL 1.125 0.152 Inf    1   0.873  0.6573
##  DL / NL  1.075 0.145 Inf    1   0.536  0.8535
## 
## Year = 2023, Veg = IS, Depth = 5-10:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 0.842 0.114 Inf    1  -1.272  0.4108
##  CTL / NL 0.996 0.134 Inf    1  -0.033  0.9994
##  DL / NL  1.182 0.160 Inf    1   1.238  0.4305
## 
## Year = 2013, Veg = SB, Depth = 5-10:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 1.136 0.153 Inf    1   0.944  0.6122
##  CTL / NL 1.267 0.171 Inf    1   1.754  0.1854
##  DL / NL  1.116 0.151 Inf    1   0.809  0.6973
## 
## Year = 2023, Veg = SB, Depth = 5-10:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 0.778 0.105 Inf    1  -1.862  0.1499
##  CTL / NL 0.947 0.128 Inf    1  -0.402  0.9147
##  DL / NL  1.218 0.165 Inf    1   1.458  0.3114
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

6| Results

# Generate letters from the treatment comparisons
letters_df <- as.data.frame(multcomp::cld(OC_groups, adjust = "tukey", Letters = letters))
letters_df$.group <- trimws(letters_df$.group)

# Position letters above the boxes
letter_heights <- aggregate(mgC_g ~ Year + Veg + Depth, data = dryland, FUN = max)
letters_df <- merge(letters_df, letter_heights, by = c("Year", "Veg", "Depth"))
letters_df$height <- letters_df$mgC_g + 2
**Figure 1.** Soil organic carbon concentrations in control (CTL), doubled-litter (DL), and no-litter (NL) treatments in 2013 and 2023, within interspace (IS) and sagebrush (SB) microsites at 0–5 and 5–10 cm depths. Different letters indicate treatment differences within the same year, microsite, and depth (Tukey-adjusted p < 0.05), based on a Gamma mixed model with a log link and random intercepts for Block and Plot. Letters do not compare years or panels.

Figure 1. Soil organic carbon concentrations in control (CTL), doubled-litter (DL), and no-litter (NL) treatments in 2013 and 2023, within interspace (IS) and sagebrush (SB) microsites at 0–5 and 5–10 cm depths. Different letters indicate treatment differences within the same year, microsite, and depth (Tukey-adjusted p < 0.05), based on a Gamma mixed model with a log link and random intercepts for Block and Plot. Letters do not compare years or panels.

**Figure 2.** Model-estimated differences in soil organic carbon change from 2013 to 2023 between doubled-litter (DL) or no-litter (NL) treatments and the control, shown by microsite (IS, interspace; SB, sagebrush) and depth. Points represent estimates and bars indicate 95% confidence intervals from a Gamma mixed model with a log link and random intercepts for Block and Plot. Estimates are log ratios of treatment versus control proportional changes. Zero indicates equal proportional change; positive values indicate greater proportional increases under treatment, and negative values indicate smaller increases or greater decreases relative to the control.

Figure 2. Model-estimated differences in soil organic carbon change from 2013 to 2023 between doubled-litter (DL) or no-litter (NL) treatments and the control, shown by microsite (IS, interspace; SB, sagebrush) and depth. Points represent estimates and bars indicate 95% confidence intervals from a Gamma mixed model with a log link and random intercepts for Block and Plot. Estimates are log ratios of treatment versus control proportional changes. Zero indicates equal proportional change; positive values indicate greater proportional increases under treatment, and negative values indicate smaller increases or greater decreases relative to the control.

Fractions

Do litter treatments differ in particulate or mineral-associated carbon in 2023 surface soils?

1| Data

fraction <- read.csv("fraction.csv")
fraction$Treatment <- factor(fraction$litter_treatment, levels = c("CTL", "DL", "NL"))
fraction$Veg <- factor(fraction$vegtype)
fraction$Block <- factor(fraction$block)

pom <- subset(fraction, fraction == "L_POM")
maom <- subset(fraction, fraction == "H_MAOM")

pom$POC <- pom$mg_POC_per_g_bulksoil
maom$MAOC <- maom$mg_MAOC_per_g_bulksoil

fraction_plot <- bind_rows(
  transmute(pom, Treatment, Veg, Fraction = "POC", Carbon = POC),
  transmute(maom, Treatment, Veg, Fraction = "MAOC", Carbon = MAOC))

2| Explore

ggplot(fraction_plot, aes(Fraction, Carbon, fill = Treatment)) +
  geom_boxplot(na.rm = TRUE) +
  facet_grid(~ Veg, scales = "free_y") +
  theme_classic() +
  labs(y = "Carbon (mg C/g bulk soil)") +
  scale_fill_manual(values = c("#d8aee5ff", "#adc3b7", "#a69e00ff"))

3| Model

POC_model <- glmmTMB(POC ~ Treatment * Veg + (1 | Block), data = pom, family = Gamma(link = "log"))
MAOC_model <- glmmTMB(MAOC ~ Treatment * Veg + (1 | Block), data = maom, family = Gamma(link = "log"))

summary(POC_model)
##  Family: Gamma  ( log )
## Formula:          POC ~ Treatment * Veg + (1 | Block)
## Data: pom
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     168.7     181.1     -76.3     152.7        27 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Block  (Intercept) 0.02786  0.1669  
## Number of obs: 35, groups:  Block, 6
## 
## Dispersion estimate for Gamma family (sigma^2): 0.049 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.79494    0.11333  15.839  < 2e-16 ***
## TreatmentDL        0.18597    0.13523   1.375  0.16904    
## TreatmentNL        0.07098    0.12797   0.555  0.57912    
## VegSB              0.78852    0.12822   6.150 7.75e-10 ***
## TreatmentDL:VegSB -0.17931    0.18795  -0.954  0.34007    
## TreatmentNL:VegSB -0.51112    0.18159  -2.815  0.00488 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(MAOC_model)
##  Family: Gamma  ( log )
## Formula:          MAOC ~ Treatment * Veg + (1 | Block)
## Data: maom
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     116.0     128.6     -50.0     100.0        28 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev. 
##  Block  (Intercept) 2.6e-11  5.099e-06
## Number of obs: 36, groups:  Block, 6
## 
## Dispersion estimate for Gamma family (sigma^2): 0.0137 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        2.101723   0.047761   44.00   <2e-16 ***
## TreatmentDL        0.093501   0.067544    1.38    0.166    
## TreatmentNL       -0.061279   0.067544   -0.91    0.364    
## VegSB             -0.008543   0.067544   -0.13    0.899    
## TreatmentDL:VegSB  0.022154   0.095522    0.23    0.817    
## TreatmentNL:VegSB  0.046264   0.095522    0.48    0.628    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4| Check fit

plot(fitted(POC_model), residuals(POC_model, type = "pearson"),
     xlab = "Fitted POC", ylab = "Pearson residuals")
abline(h = 0, lty = 2)

plot(fitted(MAOC_model), residuals(MAOC_model, type = "pearson"),
     xlab = "Fitted MAOC", ylab = "Pearson residuals")
abline(h = 0, lty = 2)

5| Compare

POC_groups <- emmeans(POC_model, ~ Treatment | Veg)
MAOC_groups <- emmeans(MAOC_model, ~ Treatment | Veg)

pairs(POC_groups, adjust = "tukey", type = "response")
## Veg = IS:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 0.830 0.112 Inf    1  -1.375  0.3539
##  CTL / NL 0.931 0.119 Inf    1  -0.555  0.8441
##  DL / NL  1.122 0.152 Inf    1   0.850  0.6717
## 
## Veg = SB:
##  contrast ratio    SE  df null z.ratio p.value
##  CTL / DL 0.993 0.128 Inf    1  -0.052  0.9985
##  CTL / NL 1.553 0.200 Inf    1   3.421  0.0018
##  DL / NL  1.563 0.201 Inf    1   3.470  0.0015
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale
pairs(MAOC_groups, adjust = "tukey", type = "response")
## Veg = IS:
##  contrast ratio     SE  df null z.ratio p.value
##  CTL / DL 0.911 0.0615 Inf    1  -1.384  0.3491
##  CTL / NL 1.063 0.0718 Inf    1   0.907  0.6357
##  DL / NL  1.167 0.0789 Inf    1   2.292  0.0570
## 
## Veg = SB:
##  contrast ratio     SE  df null z.ratio p.value
##  CTL / DL 0.891 0.0602 Inf    1  -1.712  0.2005
##  CTL / NL 1.015 0.0686 Inf    1   0.222  0.9731
##  DL / NL  1.140 0.0770 Inf    1   1.935  0.1290
## 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## Tests are performed on the log scale

6| Results

**Figure 3.** Mineral-associated organic carbon (MAOC) and particulate organic carbon (POC) in 2023 surface soils (0–5 cm) from interspace (IS) and sagebrush (SB) microsites. Treatments are control (CTL), doubled litter (DL), and no litter (NL). Boxes represent interquartile ranges, horizontal lines indicate medians, and points beyond the whiskers indicate outlying observations. Treatments without a shared letter differ within the same fraction and microsite (Tukey-adjusted p < 0.05), based on separate Gamma mixed models with log links and random intercepts for Block. Letters do not compare panels. Vertical scales differ between fractions.

Figure 3. Mineral-associated organic carbon (MAOC) and particulate organic carbon (POC) in 2023 surface soils (0–5 cm) from interspace (IS) and sagebrush (SB) microsites. Treatments are control (CTL), doubled litter (DL), and no litter (NL). Boxes represent interquartile ranges, horizontal lines indicate medians, and points beyond the whiskers indicate outlying observations. Treatments without a shared letter differ within the same fraction and microsite (Tukey-adjusted p < 0.05), based on separate Gamma mixed models with log links and random intercepts for Block. Letters do not compare panels. Vertical scales differ between fractions.

Figures

**Figure 1.** Soil organic carbon concentrations in control (CTL), doubled-litter (DL), and no-litter (NL) treatments in 2013 and 2023, within interspace (IS) and sagebrush (SB) microsites at 0–5 and 5–10 cm depths. Different letters indicate treatment differences within the same year, microsite, and depth (Tukey-adjusted p < 0.05), based on a Gamma mixed model with a log link and random intercepts for Block and Plot. Letters do not compare years or panels.

Figure 1. Soil organic carbon concentrations in control (CTL), doubled-litter (DL), and no-litter (NL) treatments in 2013 and 2023, within interspace (IS) and sagebrush (SB) microsites at 0–5 and 5–10 cm depths. Different letters indicate treatment differences within the same year, microsite, and depth (Tukey-adjusted p < 0.05), based on a Gamma mixed model with a log link and random intercepts for Block and Plot. Letters do not compare years or panels.

**Figure 2.** Model-estimated differences in soil organic carbon change from 2013 to 2023 between doubled-litter (DL) or no-litter (NL) treatments and the control, shown by microsite (IS, interspace; SB, sagebrush) and depth. Points represent estimates and bars indicate 95% confidence intervals from a Gamma mixed model with a log link and random intercepts for Block and Plot. Estimates are log ratios of treatment versus control proportional changes. Zero indicates equal proportional change; positive values indicate greater proportional increases under treatment, and negative values indicate smaller increases or greater decreases relative to the control.

Figure 2. Model-estimated differences in soil organic carbon change from 2013 to 2023 between doubled-litter (DL) or no-litter (NL) treatments and the control, shown by microsite (IS, interspace; SB, sagebrush) and depth. Points represent estimates and bars indicate 95% confidence intervals from a Gamma mixed model with a log link and random intercepts for Block and Plot. Estimates are log ratios of treatment versus control proportional changes. Zero indicates equal proportional change; positive values indicate greater proportional increases under treatment, and negative values indicate smaller increases or greater decreases relative to the control.

**Figure 3.** Mineral-associated organic carbon (MAOC) and particulate organic carbon (POC) in 2023 surface soils (0–5 cm) from interspace (IS) and sagebrush (SB) microsites. Treatments are control (CTL), doubled litter (DL), and no litter (NL). Boxes represent interquartile ranges, horizontal lines indicate medians, and points beyond the whiskers indicate outlying observations. Treatments without a shared letter differ within the same fraction and microsite (Tukey-adjusted p < 0.05), based on separate Gamma mixed models with log links and random intercepts for Block. Letters do not compare panels. Vertical scales differ between fractions.

Figure 3. Mineral-associated organic carbon (MAOC) and particulate organic carbon (POC) in 2023 surface soils (0–5 cm) from interspace (IS) and sagebrush (SB) microsites. Treatments are control (CTL), doubled litter (DL), and no litter (NL). Boxes represent interquartile ranges, horizontal lines indicate medians, and points beyond the whiskers indicate outlying observations. Treatments without a shared letter differ within the same fraction and microsite (Tukey-adjusted p < 0.05), based on separate Gamma mixed models with log links and random intercepts for Block. Letters do not compare panels. Vertical scales differ between fractions.

References

Bolker, B. M., Brooks, M. E., Clark, C. J., Geange, S. W., Poulsen, J. R., Stevens, M. H. H., & White, J.-S. S. (2009). Generalized linear mixed models: A practical guide for ecology and evolution. Trends in Ecology & Evolution, 24(3), 127–135. https://doi.org/10.1016/j.tree.2008.10.008

Brooks, M. E., Kristensen, K., van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., Skaug, H. J., Mächler, M., & Bolker, B. M. (2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal, 9(2), 378–400. https://doi.org/10.32614/RJ-2017-066

Lenth R, Piaskowski J (2026). emmeans: Estimated Marginal Means, aka Least-Squares Means. doi:10.32614/CRAN.package.emmeans https://doi.org/10.32614/CRAN.package.emmeans. R package version 2.0.4, https://CRAN.R-project.org/package=emmeans.

Hothorn T, Bretz F, Westfall P (2008). “Simultaneous Inference in General Parametric Models.” Biometrical Journal, 50(3), 346-363.