---
title: 'Predicting Absolute Holocellulose and Klason Lignin Contents for Peat Remains Challenging'
author: "Henning Teickner$^{1, *}$, Klaus-Holger Knorr$^1$"
date: ""
output: 
  bookdown::pdf_book:
    keep_tex: true
    citation_package: natbib
    toc: false
bibliography: references
biblio-style: unsrtnat
natbiboptions: super
header-includes:
   - \usepackage{float}
   - \usepackage{lineno}
   - \linenumbers
   - \usepackage{setspace}
   - \doublespacing
---

```{r setup-reply-2, include=FALSE}
knitr::opts_chunk$set(fig.align = "center", 
                      fig.pos = "H",
                      echo = FALSE,
                      cache = TRUE)
```

$^1$ ILÖK, Ecohydrology and Biogeochemistry Group, University of Münster, Heisenbergstr. 2, 48149 Münster, Germany  
$^*$ Corresponding author: Henning Teickner (henning.teickner@uni-muenster.de)  
  
The resilience of the large global peatland carbon stocks against climate change may depend on how peat chemistry controls decomposition rates [@Wang.2015; @Hodgkins.2018]. In their study, @Hodgkins.2018 argue that there is a latitudinal gradient in peat recalcitrance which may explain the stability of tropical peat in spite of warmer temperatures and which points toward a susceptibility of labile surface peat at high latitudes to degradation, but not deeper peat which has a chemistry more similar to tropical peat. This gradient can be described by differences in holocellulose contents --- representing relatively easily degradable carbohydrates --- and Klason lignin contents --- representing aromatics --- which @Hodgkins.2018 predicted from mid infrared spectra [@Hodgkins.2018].  
We acknowledge this conceptual model as an important step in understanding the stability of peat deposits and their susceptibility to degradation. However, we challenge the accuracy of the models used to predict peat holocellulose and Klason lignin contents. This inaccuracy has two major implications:  

1. Peat Klason lignin contents and relative accumulation during decomposition probably are underestimated by @Hodgkins.2018 which may have important consequences for decomposition of high-latitude peat deposits due to climate change.  

2. The inaccuracies imply that the predicted values are only very rough estimates for holocellulose and Klason lignin contents and at least partly biased. An interpretation of the values given in @Hodgkins.2018 in terms of quantitative concepts for organic matter chemistry (i.e. holocellulose and Klason lignin) therefore is not possible. 

Instead, we assume that if the study is reinterpreted in terms of qualitative concepts for carbohydrates and aromatics as used in the qualitative interpretation of mid infrared spectra for peat, there still is evidence for general latitudinal patterns in peat recalcitrance. However, interpretation and generalization of these findings are less straightforward. Moreover, we suggest improvements for uncertainty propagation for the models for latitudinal and temperature gradients. Remaining open questions to this approach should be addressed in future studies.

## Patterns in peat chemistry should not be interpreted in terms of holocellulose and Klason lignin contents {.unnumbered}

In a recent evaluation of the models to predict peat holocellulose and Klason lignin contents presented in @Hodgkins.2018, we found that the training data used for the models (leaves, paper products) are not representative for peat [---todo: ref]. The peat samples are more decomposed than the training data samples. Decomposition typically causes distinct features in mid infrared spectra which have complex effects on the applicability of spectral prediction models [---todo:ref]. Since these features are not present in the training data, predictions of the original model are extrapolations. Following our re-evaluation, this causes large uncertainties in model predictions for both holocellulose and Klason lignin, and a comparatively large bias towards smaller Klason lignin contents, especially for more decomposed peat with Klason lignin contents >~30 mass-% [---todo:ref]. Because of these uncertainties and biases, the analyses of @Hodgkins.2018 cannot be interpreted in terms of absolute holocellulose and Klason lignin contents.  
Instead, we suggest to interpret differences in peat chemistry presented by @Hodgkins.2018 not in terms of specific organic matter fractions (i.e. holocellulose and Klason lignin), but rather in terms of more qualitative concepts typically used in studies interpreting mid infrared spectra of peat (e.g. @Artz.2008; @Broder.2012; @Tfaily.2014; @Biester.2014), i.e. as "carbohydrates" and "aromatics". We assume that under such an interpretation, the conclusions of @Hodgkins.2018 probably are valid.  
The proposed reinterpretation is possible if values given in @Hodgkins.2018 are only interpreted in terms of the normalized peak heights used in the original prediction models, and not as absolute holocellulose and Klason lignin contents. Each model is a linear model with one independent variable (the normalized height of one peak); therefore a direct reinterpretation of the predicted absolute contents in terms of normalized peak heights is possible if the figure axes are relabeled as "Relative content of carbohydrates/aromatics" instead of "Estimated % carbohydrates/aromatics" (implying an absolute content of holocellulose and Klason lignin, respectively) [@Hodgkins.2018].   
An important open question to this approach of interpreting large-scale gradients in peat chemistry based on mid infrared spectra should be addressed in future studies: Are differences in specific mid infrared spectra peaks reliable indicators for differences in a specific concept for peat chemistry? For example, are differences in the peak around 1040 cm$^{-1}$ indicative for differences in carbohydrate contents (also across larger scales)?  
A similar problem arose in the discussion on the use of the C/N ratio as a good indicator for peat mass loss. This may only hold true if there have been no or negligible shifts in the vegetation composition which causes differences in the C/N ratio already in initial plant litter (e.g. @Biester.2014). For mid infrared spectra, such knowledge still has to be acquired in future studies, especially if large environmental gradients are considered as in @Hodgkins.2018.

## Klason lignin content and residual enrichment in (sub)tropical peatlands probably are larger {.unnumbered}

We suggest that @Hodgkins.2018 have underestimated peat Klason lignin contents, especially in (sub)tropical peat. As described above, the original model for Klason lignin underestimates peat Klason lignin contents. A modified model has less bias [---todo:ref]. In spite of several remaining limitations (see below and the supplementary information), we argue that this modified model can serve as rough indicator for the underestimation of peat Klason lignin contents by the original model. We therefore present values to provide tentative estimates for the underestimation of Klason lignin contents, but we do not claim that these estimates are optimal.  
With this modified model, near-surface ($\le50$ cm) peat Klason lignin contents are on average ~10 to 15 mass-% larger across the latitudinal gradient (figure \@ref(fig:reply-res-p-sites-latitude)). Specifically, average (sub)tropical near-surface peat ($<45 ^\circ$N, $\le50$ cm) Klason lignin contents are `r d_peat_pred_best_diff_mean %>% dplyr::filter(variable_y == "Klason lignin" & !high_latitude & surface_peat) %>% dplyr::pull(y_hat_diff) %>% median() %>% magrittr::multiply_by(100) %>% round(0)` [`r d_peat_pred_best_diff_mean %>% dplyr::filter(variable_y == "Klason lignin" & !high_latitude & surface_peat) %>% dplyr::pull(y_hat_diff) %>% quantile(probs = c(0.025, 0.975)) %>% magrittr::multiply_by(100) %>% round(0) %>% paste(collapse = ",")`] mass-% (median, lower and upper 95% prediction interval limit) larger than with the original model. There is also tentative indication for on average larger differences between deep (sub)tropical ($<45 ^\circ$N) and high-latitude ($\ge45 ^\circ$N) peat (figure \@ref(fig:reply-res-p-depth-profiles-latitude)).    
We also suggest that the residual enrichment of Klason lignin during decomposition probably has been underestimated. A rough estimate for the residual enrichment of Klason lignin during decomposition is the difference in Klason lignin content between near-surface peat ($\le50$ cm) and deeper peat. With the modified model, this difference is on average for (sub)tropical peatlands `r d_peat_pred_best_diff_mean1 %>% dplyr::mutate(res = y_hat_diff - y_hat_or_diff) %>% dplyr::filter(variable_y == "Klason lignin" & !high_latitude) %>% dplyr::pull(res) %>% median() %>% magrittr::multiply_by(100) %>% round(0)` [`r d_peat_pred_best_diff_mean1 %>% dplyr::mutate(res = y_hat_diff - y_hat_or_diff) %>% dplyr::filter(variable_y == "Klason lignin" & !high_latitude) %>% dplyr::pull(res) %>% quantile(probs = c(0.025, 0.975)) %>% magrittr::multiply_by(100) %>% round(0) %>% paste(collapse = ",")`] mass-% larger (figure \@ref(fig:reply-res-p-depth-profiles-latitude)). For high latitude peatlands the difference is smaller on average (`r d_peat_pred_best_diff_mean1 %>% dplyr::mutate(res = y_hat_diff - y_hat_or_diff) %>% dplyr::filter(variable_y == "Klason lignin" & high_latitude) %>% dplyr::pull(res) %>% median() %>% magrittr::multiply_by(100) %>% round(0)` [`r d_peat_pred_best_diff_mean1 %>% dplyr::mutate(res = y_hat_diff - y_hat_or_diff) %>% dplyr::filter(variable_y == "Klason lignin" & high_latitude) %>% dplyr::pull(res) %>% quantile(probs = c(0.025, 0.975)) %>% magrittr::multiply_by(100) %>% round(0) %>% paste(collapse = ",")`] mass-%).  
Consequently, in general --- and especially in (sub)tropical peatlands --- the average residual enrichment of Klason lignin due to decomposition probably has been underestimated by the original model. A consequence of this in light of climate change is that high latitude peat deposits may be more prone to shifts in decomposition rates than suggested by the original results.  

## Larger uncertainties due to a lack of uncertainty propagation and unrepresentative data {.unnumbered}

Besides a reinterpretation in terms of qualitative concepts for organic matter recalcitrance, we also believe that a more comprehensive quantification and propagation of uncertainties is important to further scrutinize the results of @Hodgkins.2018. In the original analysis, the following uncertainties are omitted: (1) As described above, the training dataset is not representative for peat and therefore predictions are biased. This uncertainty is not considered in the analyses of @Hodgkins.2018. (2) The second source of uncertainty which was omitted is the prediction uncertainty of the models used to predict holocellulose and Klason lignin contents. (3) A third source of uncertainty not considered yet is that in computing the regression model to analyze the latitudinal gradient based on average core near-surface contents, the uncertainty for these average values was not propagated.  
The first two sources of uncertainty can be omitted if only normalized peak heights (as qualitative indicators) are used instead of predicted holocellulose and Klason lignin contents, as described above. But the third source of uncertainty should also be considered when reinterpreting the results in terms of area normalized peak heights.  
To support our concerns, we reanalyzed the latitudinal gradient for holocellulose and Klason lignin, respectively, using the modified models described in [---todo:ref] and fully considering the last two sources of uncertainty using a Bayesian approach (described in more detail in the supplementary information). With the modified models, the estimated slopes for the latitudinal gradient are more uncertain (95% confidence intervals are [`r s1_res_holocellulose %>% dplyr::pull(b_latitude) %>% quantile(probs = c(0.025, 0.975)) %>% round(1) %>% paste0(collapse = ",")`] and [`r s1_res_klason_lignin %>% dplyr::pull(b_latitude) %>% quantile(probs = c(0.025, 0.975)) %>% round(1) %>% paste0(collapse = ",")`] for holocellulose and Klason lignin, respectively; figure \@ref(fig:reply-res-p-sites-latitude) and supporting figure 3<!---\@ref(fig:reply-res-p-latitude-poc-slope)--->). A similar reanalysis of the depth profiles between the two latitude categories shows that differences are also more uncertain (figure \@ref(fig:reply-res-p-depth-profiles-latitude)). Our analysis therefore indicates larger uncertainties both in the depth profiles and the latitudinal gradient. Considering these uncertainties is important when peat properties predicted with a model are analyzed.
Again, the values presented here are only tentative because the modified models probably are also biased.

## Future directions {.unnumbered}

Notwithstanding our critical remarks, the idea of @Hodgkins.2018 to relate mid infrared spectra to precisely quantifiable and defined organic matter fractions is important. Such fractions are relatively well understood and therefore allow us to better deduce what effects climate shifts can have on peat decomposition and recalcitrance. It is much more informative if we are able to precisely estimate a change in Klason lignin contents of e.g. 5 mass-% if average annual temperatures increase by 1$^{\circ}$C than if we are only able to make qualitative statements that the relative amount of aromatics (without describing in detail what this fraction means) increases, as is currently the case.  
When interpreting mid infrared spectra, distinguishing between precisely quantifiable concepts for organic matter chemistry, such as holocellulose and Klason lignin contents, and qualitative concepts which currently are used in interpreting peat MIRS, such as "carbohydrates and aromatics" (without specifying how this fraction could actually be measured), is crucial: Models relying on quantitative concepts are more testable and allow to actually quantify fractional changes in organic matter chemistry, whereas models built on qualitative concepts can only be refuted if there are blatant discrepancies relative to expectations (i.e. a decrease in "carbohydrates" with latitude). Ultimately, we need more precise models for peat variations in chemistry than can be developed based on qualitative concepts.  
We suggest the following steps to gain deeper insights into latitudinal gradient in peat chemistry: As outlined above, we suggest to interpret the results of @Hodgkins.2018 in terms of qualitative concepts for peat chemistry, i.e. "carbohydrates" and "aromatics". These suggestions also apply to a follow-up study using the same models [@Verbeke.2022]. Further studies are needed to check if specific peaks are good indicators for these qualitative concepts for peat chemistry, how the relation may be confounded e.g. by vegetation composition or other organic matter fractions, as could be the case especially for peaks related to aromatics [---todo:ref]. This could be done by computing models to reliably predict holocellulose and Klason lignin contents of peat from mid infrared spectra and analyzing which spectral variables are relevant for this purpose and if there are confounding variables. Such models will also help us to get more precise insights into the degree to which high latitude peat accumulation may be decreased by decomposition due to rising temperatures, because they allow replacing more costly measurements. Finally, future studies should carefully consider uncertainties and fully propagate these into model estimates. Predicting absolute peat holocellulose and Klason lignin contents from infrared spectra remains challenging.

(ref:reply-res-p-sites-latitude-caption) Predicted average surface ($\le 50$ cm) peat holocellulose (a) and Klason lignin (b) contents plotted against latitude (compare with Fig. 3 in @Hodgkins.2018). Lines and shaded areas represent average predictions from regression models and 95% confidence intervals. "Modified" is the modified Bayesian hierarchical regression model which simultaneously models individual samples' contents from mid infrared spectra and the latitudinal gradient of average core near-surface peat contents. This model considers prediction uncertainty from the mid infrared spectra and from computing per-core averages. "Original" is the original linear regression model [@Hodgkins.2018] computing only the latitudinal gradient of average core near-surface contents. This model does not consider prediction uncertainty from mid infrared spectra, nor uncertainty from computing per-core averages. Points are average core near-surface contents predicted from the model ("Modified") or computed from the average predictions for individual samples ("Original") with error bars representing 95% confidence intervals. Points for "Modified" are shifted by $+0.1^\circ$.

```{r reply-res-p-sites-latitude, echo=FALSE, out.width="70%", fig.height=3.5, fig.width=6.5, fig.cap='(ref:reply-res-p-sites-latitude-caption)'}
p_poc_latitude6
```

(ref:reply-res-p-depth-profiles-latitude-caption) Predicted median holocellulose (A) and Klason lignin (B) depth profiles of peat core data classified into two latitude categories following @Hodgkins.2018 (compare with Fig. 2 in @Hodgkins.2018). Lines are averages of LOESS smoothers fitted to the predicted values by the models. Shaded regions are corresponding 95% confidence intervals, comprising prediction uncertainty in the holocellulose and Klason lignin estimates, respectively. Vertical dashed lines in each columns represent approximate surface and average mean contents as predicted using the improved models ("Modified").

```{r reply-res-p-depth-profiles-latitude, echo=FALSE, warning=FALSE, out.width="70%", fig.height=5.5, fig.width=6.5, fig.cap='(ref:reply-res-p-depth-profiles-latitude-caption)'}
p_depth_profiles_latitude
```

```{r reply-res-p-veg-peat-differnece, echo=FALSE, eval=FALSE, out.width="90%", fig.height=3.5, fig.width=7.5, fig.cap='Predicted surface ($\\le 50$ cm) peat and vegetation median holocellulose (A) and Klason lignin (B) contents for the peat cores shown in figure \\ref{fig:res-p-gaussian-beta-depth-profile-best} grouped according to \\cite{Hodgkins.2018}. Predictions of the (biased) original models are equivalent to the data shown in \\cite{Hodgkins.2018} (figure 5). Averages are median values of predictions, error bars are 95\\% prediction intervals. "Difference" is the difference between predictions for vegetation and peat.'}
p_veg_peat_difference
```

# Acknowledgements {.unnumbered}

We like to thank Jeffrey Chanton and Suzanne Hodgkins for a constructive dicussion which greatly helped to improve a first version of this manuscript. This study was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) grant no. KN 929/23-1 to Klaus-Holger Knorr and grant no. PE 1632/18-1 to Edzer Pebesma.

# Competing interests {.unnumbered} 

The authors declare no competing interests.

# Author contributions {.unnumbered}

HT performed the calculations, prepared the figures, and wrote the original text. Both HT and KHK revised and edited the text. KHK provided funding for this study.

# Data and Code availability {.unnumbered}

Data and code to reproduce our analyses are available via <https://doi.org/10.5281/zenodo.6325760.>
