Common structure in the heterogeneity of plant-matter decay

Carbon removed from the atmosphere by photosynthesis is released back by respiration. Although some organic carbon is degraded quickly, older carbon persists; consequently carbon stocks are much larger than predicted by initial decomposition rates. This disparity can be traced to a wide range of first-order decay-rate constants, but the rate distributions and the mechanisms that determine them are unknown. Here, we pose and solve an inverse problem to find the rate distributions corresponding to the decomposition of plant matter throughout North America. We find that rate distributions are lognormal, with a mean and variance that depend on climatic conditions and substrate. Changes in temperature and precipitation scale all rates similarly, whereas the initial substrate composition sets the time scale of faster rates. These findings probably result from the interplay of stochastic processes and biochemical kinetics, suggesting that the intrinsic variability of decomposers, substrate and environment results in a predictable distribution of rates. Within this framework, turnover times increase exponentially with the kinetic heterogeneity of rates, thereby providing a theoretical expression for the persistence of recalcitrant organic carbon in the natural environment.


INTRODUCTION
Greater than 90 per cent of the carbon dioxide input to the atmosphere -ocean system each year derives from the natural decay of organic carbon [1,2]. Decay is heterogeneous in space and time: organic molecules vary in lability [3,4]; micro-environmental heterogeneity such as the aggregation of minerals in soil and sediments interfere with degradation [5]; humification and repolymerization [4,6] result in polymers that are difficult to degrade; and decomposer communities are diverse and varied [4,7]. Physical and chemical changes in local environment [8] can speed up or prevent decomposition altogether. Spatial heterogeneity of soil nutrient concentrations at the meter scale [9] may also influence degradation rate heterogeneity. Combined, these diverse effects yield kinetic heterogeneity: older compounds appear to decay at slower rates than younger compounds [10,11], and carbon stores and their turnover times are larger than predicted from initial decomposition rates [12].
The sizes of the organic carbon stores and their rates of turnover are required for quantifying feedback between climate and the carbon cycle [13,14] in order to predict changes in carbon dioxide levels and climate [1,2,8,15]. Because decay time scales vary widely, from minutes to millions of years, estimates of carbon stocks and turnover times require knowledge of all decay rates, from fast to slow [16,17]. However, rate distributions and the mechanisms that determine them remain unknown. Identification of rate distributions should provide insight, not only for predictive purposes [17][18][19] but also for understanding the ecological dynamics [20] of decomposition.
Because the wide range of time scales makes it impossible to directly measure decay over all phases of decomposition, we focus on plant litter decay and early transformations to young soil organic matter. Specifically, we investigate measurements from the Long Term Intersite Decomposition Experiment Team (LIDET) study [2,12,[21][22][23]. This study monitored the decomposition of 27 different types of litter, including needles, leaves, roots, wood, grass and wheat distributed among 28 different locations across North America ranging from Alaskan tundra to Panamanian rainforests.
Litter was collected and then re-distributed in litter bags at different sites in order to investigate the effect of composition, ecosystem and climatic parameters on decomposition. Litter bags were collected and analysed at least once per year for up to 10 years. We show how to estimate kinetic heterogeneity from these observations of decay. We also analyse how these rate distributions are related to climatic conditions and litter composition.
The remainder of this paper is organized as follows. In §2, we pose and solve an inverse problem to find the rate distributions corresponding to the decomposition of plant matter from the LIDET study. We then show in §3 that the distributions are lognormal on average. Subsequently, in §4, we show how the two parameters of the lognormal distribution depend on composition and environment. In §5, we derive the relation of these parameters to the turnover time of carbon stocks. These results show that turnover times grow exponentially as the heterogeneity of rates increase, thereby highlighting the dependence of carbon stocks on their slowest rates of decay [17].

DISORDERED KINETICS
Organic matter decomposition may be viewed as the relaxation of a kinetically disordered system. In this section, we specify a model for the influence of disorder on decay. We then describe how we invert it to obtain a distribution of decay rates.

Model
We suppose that decay rate constants k derive from stochastic reactions between heterogeneous substrates and ecological communities in a random environment. We describe this scenario using a 'static' model of 'disordered kinetics' [24][25][26]. In this model, the mass g(t) is a decreasing function of time that derives from a continuous superposition of exponential decays e 2kt weighted by the probability pðkÞdk that the rate constant k is present at the onset of decay. Given these assumptions, decay proceeds as where p(k) ! 0 and Ð 1 0 pðkÞdk ¼ 1. Models similar or identical to equation (2.1) have been previously employed to describe organic matter decay [17][18][19][27][28][29][30][31][32][33]. When the distribution p(k) is discrete, the integral in equation (2.1) becomes a sum known as a 'multi-G' or 'multi-pool' model [21,34,35]. Although equation (2.1) lacks detailed mechanisms of the processes involved in decomposition, its simplicity and common application suggest that it is a reasonable first attempt at characterizing decomposition dynamics. Dispersion of the rates k in this model are probably associated with variations in the quality of plantmatter compounds [36,37], which range from highly labile simple sugars to more refractory lignin, waxes and phenolic compounds [35]; local spatial heterogeneity in soil moisture and nutrients [9]; chemical transformation of compounds [6] and decomposer and metabolic diversity [38,39]. Rather than attempting a detailed characterization of these individual mechanisms, we simply seek the distribution of rates associated with the minimal description of decomposition given by equation (2.1).
Although equation (2.1) represents a system of parallel steady decays, decomposition also involves temporal disturbances and serial processes. However, serial transformation processes can be mathematically expressed as parallel decays [17]. Regarding temporal fluctuations, we interpret the steady distribution p(k) as the probability that decay occurs at an effective first-order rate k that is averaged over seasonal and other disturbances [17]. We also note that if the difference between the time scales of two serial processes is large, the system effectively relaxes at the time scale of the slower process. For example, the degradation time scale of a particle attached to a mineral surface may be much larger than the duration of the transient period before attachment; similarly, the time scale of humification is probably short relative to the lifetime of the slowly degrading humic substance [17]. Decomposition may be approximated as proceeding initially from the mineral-associated or humic state [17]. A consequence of a parallel decay model is that resulting decays g(t) are convex (concave-up). Specifically, any completely monotone decay g(t)/g(0) can be described by a linear superposition of rates weighted by a probability density function p(k) [40].

Inverse problem
Under certain physical conditions, distributions p(k) of reaction rates can be calculated analytically [24,27,28] and evaluated by comparing g(t) to experimental data. However, given the complex nature of decomposition, purely physical models may not be appropriate. We therefore seek to identify the distribution that best fits the data without resorting to assumptions beyond those implied by equation (2.1). Once the best distribution is found, physical reasoning then allows the identification of mechanisms that can generate this distribution.
Mathematically, equation (2.1) is a Laplace transform and p(k) can be found from its inverse. However, the inverse Laplace transform is ill-posed [44], meaning that small changes in the data g(t) can result in large changes in the solution p(k). A standard method to solve such ill-posed problems is to seek solutions p(k) that are minimally 'rough' [44]. Here, we use Tikhonov regularization [44,45] to identify an optimally smooth p(k) that best fits the data ( §2). Such methods have been previously applied to problems of NMR spin relaxation [41] to probe the structure of porous media [46,47] and the properties of biological tissue [42].

RATES ARE DISTRIBUTED LOGNORMALLY
We apply this procedure to litter decomposition data from the LIDET study. An example of decay from an LIDET dataset is shown in figure 1a. The corresponding estimate of the rate distribution in logarithmic space, expressed as rðxÞ ¼ p½kðxÞdk=dx, where x ¼ ln k, is shown in figure 1b. The rate k is rescaled by the period of seasonal forcing (1 year) and is therefore non-dimensional. The good fit of r(ln k) to a Gaussian indicates that the distribution of rates is lognormal, characterized by the parameters m and s, where m is the mean of ln k and s 2 is the variance of ln k.
To investigate the extent to which the lognormal distribution applies to the remainder of the LIDET data, we identify the 234 LIDET datasets that contain at least five measurements with replicates. These datasets contain 11 different litter types distributed among 26 sites. We then employ several tests on each of these 234 datasets to check whether equation (2.1) is an appropriate description of these datasets. First, we find that seven of these datasets show insignificant mass loss between the first and last field measurement, rendering equation (2.1) irrelevant. Six of these datasets are associated with root decay, suggesting that roots can persist for long times in certain conditions. We next identify the datasets that decay faster than exponentially, counter to the assumption of first-order kinetics and a superposition of exponential decays. We use two tests to identify these datasets. First, we check the curvature of the datasets and find that nine of the datasets have negative curvature (concave-down). Such superexponential decay cannot be consistent with the Laplace transform relation (2.1) [40]. These datasets are primarily located at sites having low precipitation, indicating that decay dynamics may be limited by moisture or decomposer activity, rather than substrate availability. Second, we apply our inversion procedure to the remaining datasets and find that three datasets have a significant trend in the residual error. These datasets decay faster than exponentially, but are not concave down. All three of these datasets are associated with wood decay. In summary, our tests disqualify 7 þ 9 þ 3 ¼ 19 datasets from further consideration. Further details of the tests are given in appendix A.1.
Of the remaining 215 datasets, our inversion procedure indicates that 33 datasets are characterized by a single rate constant and decay exponentially. Guided by the result of figure 1b, we then fit a Gaussian to the 182 estimates of r(ln k) exhibiting a non-zero variance and rescale each by the fitted parameters m and s. We plot the mean of the rescaled distributions in figure 1c. Although there is scatter and skew among the individual estimates of r(ln k), figure 1c shows that the mean of the rescaled distributions of ln k is very similar to a Gaussian distribution. Because the 33 single-rate datasets correspond to a lognormal distribution with zero variance, our results indicate that the lognormal represents the average rate distribution of the 215 datasets for which the model (2.1) applies. Lognormally distributed variables arise naturally from multiplicative stochastic processes [48]. Here, lognormally distributed rates may result from the multitude of seemingly stochastic requirements for decomposition, such as the presence of water, the presence of an appropriate microbe, the lack of predation, the conditions for expression of hydrolytic enzymes, the encounter of enzymes with the organic matter, etc. [6]. More generally, the probability of completing any task that relies on the successful completion of many subtasks is lognormal [49]. In this context, the lognormal can be viewed as a null hypothesis in which decomposition rates result from the occurrence of a large number of independent decay requirements [6]. Mathematically, if we assume that the probability P of decomposing a parcel of organic matter over a time span Dt is the product of independent probabilities of satisfying various requirements for decay over that interval, then the first-order rate constant k ¼ P/Dt becomes asymptotically lognormal as the number of requirements increases. In this manner, the multiplicative stochasticity of a decay system results in the lognormal distribution. This general description suggests that attempts to precisely model the individual mechanisms that stochastically interact to form this broader pattern would be overly complex. It also agrees with the idea that decay rates are the product of many compositional and environmental effects [22].
Previously suggested forms of the rate distribution p(k) are the gamma distribution [17,27] and the loguniform distribution [28]. The log-uniform distribution, for which r(ln k) is constant and pðkÞ / 1=k between prescribed limits, approximates the lognormal when ðln k À mÞ=2s 2 ( 1 [49]. Moreover, its Laplace transform asymptotically approaches the Laplace transform of the lognormal distribution as s ! 0. The gamma distribution, however, differs significantly from the lognormal. We find that the lognormal distribution predicts both the data g(t) and describes the inferred rate distribution p(k) better than the gamma distribution for 177 out of 215 datasets (electronic supplementary material, §2). We have also compared the lognormal to exponential and multi-pool models. Because our inversion procedure indicates that only 33 of 234 datasets are described by a simple exponential decay, we find that simple exponential decays are generally under-parametrized for describing litter decay datasets, consistent with previous studies [12,35]. The best fitting type of multiple pool model varies widely among the datasets, with no single model type describing all datasets [12,50]. A universal multiple pool model, containing pools of various types (leached, labile, refractory, inert, etc.), would be overparametrized. Furthermore, the number of pools and rates associated with each pool are sensitive to noise, as different combinations of pools can represent the same decay [51,50]. This sensitivity makes understanding the constitutive relationships between pools and environmental and compositional parameters difficult [31].
An advantage of the lognormal is that it parametrizes decay by only two variables, m and s. We proceed in §4 to identify relations between the lognormal parameters m and s, and the climatic and compositional parameters associated with the LIDET study.

CONTROLS ON THE LOGNORMAL PARAMETERS
We seek an understanding of the controls on m (the mean order of magnitude of rates) and s 2 (the variance of those orders of magnitude). Before analysing all 215 estimates of these parameters, we identify values of m and s that are highly uncertain by disregarding the small fraction of datasets having anomalously long turnover times t. Assuming a soil carbon store is in steady state with a constant litter input, its turnover time t is equal to its mean residence time [52], which in the random-rate model (2.1) equals the mean time constant kk À1 l [31]; thus After evaluating the turnover times associated with all 215 datasets using equation (4.1), we find that there is a distinct group of datasets associated with excessively long turnover times greater than 1000 years (figure 5). These datasets contain a significant mass fraction that is effectively inert, having unknown decay dynamics. Extrapolating the kinetics of such slow processes therefore has considerable uncertainty. There are 24 datasets in this outlying cluster. These data are typically associated with root decay at certain locations (table 2), suggesting that the soils of certain ecosystems can enable the persistence of roots for long times. We do not consider these 24 datasets in our subsequent analysis of m and s. Further discussion of these outliers can be found in appendix A.1.
Owing to the nature of the LIDET study, many different litter types were placed at the same site and we therefore have many estimates of m and s at each value of temperature, precipitation and other climatic variables. Similarly, because each litter type was planted at many sites, there are many different estimates of m and s for each value of initial lignin concentration, nitrogen concentration, etc. In the following section, we study how the average values m and s of the lognormal parameters m and s vary with measured independent variables such as temperature, lignin, nitrogen, etc. When analysing the effects of climatic variables, m and s represent the averages over all litters at each site, and when analysing the effects of compositional variables, m and s represent the averages over all sites where the litter was deployed. Similarly, s 2 represents the average variance s 2 , and t represents the average turnover time, etc. Analagous depictions of the unaveraged data can be found in appendix A.3.

The mean m
We first investigate how climatic conditions and composition affect m. Figure 2a shows a significant positive correlation between m and temperature. From this trend, we find that the median decomposition rate e m increases by a factor Q 10 ¼ 2.0 + 0.3 (1 s.d.) with a 108C increase in temperature, in agreement with previous estimates [2]. All other measured and synthetic climatic parameters also significantly correlate with m, with the climate decomposition index [21,22] exhibiting the highest correlation (table 1).
The parameter m is also related to composition: figure 2b shows that m decreases as the initial ligninto-nitrogen ratio ('/N) increases. The observed trend indicates that increases in the lignin concentration, a refractory component of plant matter, are associated with a reduction in m, while increases in organic nitrogen, an important nutrient for microbial decomposers [36,37], are associated with an increase in m. This is consistent with the use of '/N as a measurement of litter quality [6,21,22,53]. The carbon-to-nitrogen ratio (C/N) and other nutrient measures are also correlated with m (table 1). Concentrations of lignin, N, S, P, etc., represent initial values associated with each type of litter. Figure 2b also indicates that needles have lower values of m than leaves. This effect, however, may be related to the difference in '/N between the two tissue types.

The variance s 2
We next investigate the relation of climatic conditions to the heterogeneity of decomposition rates, represented by s. Figure 2c shows that temperature has no significant effect on s. Moreover, s is uncorrelated with all climatic parameters monitored in the LIDET study (table 1); thus climatic conditions appear unrelated to s. We therefore find no evidence from the decadal LIDET data that the Q 10 of refractory components is significantly different than the Q 10 of labile components. This supports respiration models such as CENTURY [16,54], which uses the same temperature and soil moisture factor for each pool of organic matter, independent of lability. We note that if rate  dispersion reflects the variation in activation energies of decay processes [55], then Arrhenius kinetics suggest that s only slightly decreases with temperature over the 358 temperature range associated with the LIDET sites. This is consistent with the data presented in figure 2c, but the wide variation in s indicates that this trend is not significant and that kinetic heterogeneity is controlled by other variables.
Although s exhibits no relation to climate, it does vary with composition. Figure 2d indicates that s decreases as the initial lignin-to-nitrogen ratio ('/N) increases. Because '/N correlates negatively with m, decreasing the ratio of these components tends to both shift and stretch the rate distribution, increasing the rate constants k of the faster decay processes, while the rate constants of slower, more refractory processes are relatively unchanged. Nutrients such as N and S, and to a lesser extent P and K, exhibit similar relationships with m and s (table 1). Physically, these relationships indicate that nutrient limitation is present at early times as faster processes appear to depend strongly on the nutrient content of the litter. Slower, more refractory processes take place at rates probably sustained by the transport and immobilization of nutrients from the surrounding soil [53] and are not nutrient-limited. In fact, increased nitrogen content may inhibit the degradation of transformed plant compounds [6], widening the slow tail of the distribution and increasing s. Lignin, on the other hand, may reduce the rate constants k of more labile compounds by shielding them via a ligno-cellulose polymer matrix [21], suggesting that '/N measures a resistivity to initial decay. The effect of '/N on s also appears to saturate at low '/N, suggesting that these mechanisms lose control after crossing a threshold [6] of high N or low lignin content is reached. We also observe in figure 2d that roots and leaves tend to have higher s than needles, yet the effect of '/N on s appears less strong for roots and wood, both of which decompose underground. Roots and wood do however follow the trend of m versus '/N, suggesting that the effect of initial composition may persist over time in roots and wood, effecting a wider portion of their rate distribution, not just the fast rates. This behaviour may be related to components in their tissues, underground decomposition or both.
Collectively, the results of figure 2 and table 1 suggest that climate variability changes the median rate of decay, e m , whereas the variance of decay time scales, s 2 , appears to be a property of the litter sample itself and its relationship to the decomposer community inhabiting it. Table 1 identifies additional correlations between climatic, compositional and the lognormal parameters. Sulphur, another important microbial nutrient, is highly correlated with both m and s. Potassium exhibits a similar trend as well. The causality of the trends in table 1 however is not always clear. For example, ash also has a significant positive correlation with both m and s. However, this is most probably explained by the strong rank correlation (r ¼ 0.89) between ash and sulphur, as well as a strong correlation between ash and metals which also have a positive correlation with m and weak correlation with s. Ash is composed of sulphates, K, P, Ca and other metals [4]. Phosphorus surprisingly does not show as strong a signal as N or S and its large p-values suggest that trends with m and s may not be significant. It is possible that the initial phosphorus concentrations may contain errors because phosphorous, as with N, S, K and Ca, is present in lower concentrations in conifer needles than in deciduous leaves [56]; the values of phosphorus measured in the needles and leaves of the LIDET study do not follow this pattern. Metals contain some important rare nutrients for microbial decomposers; we find that they are more significantly correlated with m than s. The lack of a significant trend for organic compound types (other than lignin) is also surprising, as we would expect water soluble carbohydrates to affect faster decomposition time scales, and cellulose to also play a role in dynamics.

Further trends
Latitude, used as a proxy for the variability in seasonal temperature, does not show a correlation with s, indicating that temporal fluctuations in temperature do not contribute to the rate heterogeneity. A comparison of average monthly temperature and precipitation data with s also supports this finding. This result provides further evidence that rate heterogeneity is set by non-climatic factors, and that climate scales the time scale of both labile and refractory processes roughly equally.

SCALING UP TO THE CARBON CYCLE
The heterogeneity of decomposition rates has strong implications for the dynamic properties of carbon stocks. The derivative of equation ( The mean kkl is exponentially greater than the median e m because of the heavy tail of p(k). A similar amplification acts to exponentially increase the turnover time t to values much greater than e Àm . Using equation (4.1) and assuming p(k) is lognormal, one finds ¼ kkl À1 e s 2 : ð5:3Þ These relations show that rate heterogeneity has a profound effect: kkl À1 underestimates t by a factor that grows exponentially with the variance s 2 . As the distribution widens, fast rate-constants weigh heavily on the calculation of kkl, whereas slower rate-constants set the mean residence time kk À1 l. The upshot is that both the size of organic carbon stocks ( proportional to t in the steady state) and the time scale of the transient response to a disturbance (also related to t) grow exponentially with the heterogeneity s 2 of rates. These effects are a consequence of the heavy tail of the lognormal distribution. We calculate kkl and t for each dataset from our inversion using equations (5.1) and (4.1) and find the average of the log of their values, ln kkl and ln t, for each litter type. Focusing on the effects of composition, figure 3a shows a strong negative correlation between ln kkl and '/N, whereas figure 3b shows no significant correlation between the average order of magnitude of turnover time ln t and '/N. Physically, these relations reflect the unequal influence of composition on faster and slower rate constants k. Because kkl is also the initial decomposition rate, we conclude that the initial '/N exhibits strong control over early decomposition [21,53]. This influence of initial composition is eventually lost, not only at later times [6] but also in the steady state. Mathematically, these observed trends follow from equations (5.2) and (5.3), given that '/N correlates negatively with both m and s (figure 2b,d). We find that leaves, needles and roots on average have roughly the same turnover times: 10 years, 11 years and 14 years, respectively. The geometric mean turnover time of all 191 datasets is 11.5 years, but deviations from this characteristic value appear not to be controlled by initial composition. Recall from §4 that roots may also have uncharacterizably long residence times in certain locations and these are not analysed in figure 3, suggesting a larger departure of root turnover time from needles and leaves. Conditions resulting in extremely persistent root organic matter are unclear (see appendix A.1). Because the turnover time is unaffected by initial nitrogen concentration, we cannot claim that changes in the nitrogen content of the litter (perhaps through changes in nitrogen deposition) will affect the turnover time of plant matter or carbon storage in soils. It is possible that soil composed of the parent material (as opposed to the LIDET transplant study) may show a different relationship between nitrogen and turnover time. Changes in temperature and precipitation on the other hand affect m only and therefore do influence turnover time and soil carbon storage. The patterns observed in figure 3a -c suggest the following physical interpretation: initial litter composition tends to change the faster rates in the continuum, which affect both m and s. The slower rates associated with a long-term behaviour and turnover time are less related to initial litter chemistry and are more likely to be determined by soil and microbial community properties. Therefore, during the later stages of litter decay, continued transformation to soil organic matter and its subsequent decay are less a function of the parent material and more a function of semi-transformed compounds and its local interaction with soil [6]. Furthermore, early degradation may be nutrient-limited and depend on the nutrient content of the litter, whereas the slower paced degradation of more recalcitrant materials may be sustained by immobilization of nutrients from the surrounding soil. The departure of roots from the trend in figure 3c, specifically the relative constancy of s under a change in m, suggests that the effect of initial composition may persist during root decay or decomposition below ground, influencing the rates of slower processes as well. Figure 4 depicts our main findings: (i) decomposition rates are distributed lognormally; (ii) environmental change acts as a catalyst that scales all rates similarly, consistent with the models (such as CENTURY) that assign the same temperature and moisture sensitivity across all pools of organic matter; and (iii) faster processes are more sensitive to litter composition (e.g. '/N, tissue type) than slower processes. The first result, made possible by inverting equation (2.1), identifies the structure of the kinetic heterogeneity associated with decomposition. The second addresses an ongoing debate concerning the temperature sensitivity of decomposition at different time scales [21,55]. The third result identifies a control for the dispersion of decomposition time scales and shows why composition affects initial decay without changing the turnover time. Each conclusion is separate and independent of the others.

CONCLUSION
Ecosystem models are often coupled with global circulation models [14,[57][58][59][60][61] in order to provide an insight into the climate system. Incorporation of lognormally distributed decay rates in popular ecosystem models and use of the lognormal to precisely predict carbon turnover and storage would require careful parametrizations [6,14,16,54] between the lognormal parameters m and s and climatic, soil and compositional parameters. We have provided a first approach for quantifying these relations. However, a more detailed analysis incorporating known mechanisms [6] is required to provide a more comprehensive picture.
We note that the wide range of conditions under which lognormal rates are expected suggests that our results are general, applicable to other degradation processes in natural environments. Evidence of this generality is seen in the decay of marine sedimentary organic matter, which is well described by the quantitatively similar log-uniform distribution [28]. The ubiquity of lognormally distributed degradation rates suggests that a focus on factors that affect rate heterogeneity, rather than specific rates themselves, will lead to a greater understanding-and improved predictions [6,13]-of the ways in which the carbon cycle interacts with climate.

A.1. Data screening
Our analysis uses litter bag data from the LIDET study [12,21 -23]. An important measurement taken during the study is the fraction of original ash-free mass remaining after a given type of plant matter decomposes for a set amount of time in a particular location. Each measurement represents the average of up to four replicates at each site, litter type and removal time. Litter bags were typically collected and analysed each year for up to 10 years, except for bags at tropical and sub-tropical sites, which were more frequently collected at three to six month intervals.
We call a data point the average mass fraction remaining of all replicates of a given plant-matter type, site and duration. If a data point consists of only one litter bag (no replicates), we average that litter bag with the litter bags of the temporally nearest data point (typically 1 year before or 1 year after). We call the time series of data pointsgðtÞ ¼ gðtÞ=gð0Þ associated with a given combination of site and litter type a dataset. Before analysis, we subject the LIDET datasets to a series of four tests.
(1) First, we consider only the 234 LIDET datasets that each contained at least five data points (not including time zero) with at least one replicate each. (2) The datasets are subjected to a Mann -Kendall test of trend [62] to check for a significant trend of mass loss between the time of the first data point (not time zero) and the final data point. According to this criterion, seven datasets have no significant decay over this duration and are therefore eliminated from further consideration. (3) The remaining datasets are then checked against the assumption that mass loss is described by the superposition of exponential decays given in equation (2.1) of the main text and equation (A 1) given below. To be consistent with the superposition, any decaygðtÞ must be completely monotone, meaning that all even derivatives ofgðtÞ must be !0 and all odd derivatives must be 0 (see §XIII.4 of Feller [40]). Because, the Mann-Kendall trend test identifies significant decreasing trends, we additionally check only the average sign of the second derivative. To do so, we fit a quadratic function a þ bt þ ct 2 togðtÞ and eliminate the nine datasets for which c , 0. (4) Of the remaining 218 datasets, three more datasets are eliminated because the inversion algorithm cannot identify a solution p(k) satisfying pðkÞ . 0 and Ð 1 0 pðkÞdk ¼ 1 without having a significant trend in the residual error. The combined application of these four criteria leaves 215 datasets for analysis. As discussed in the text, we find that the average r(ln k) determined from the inversion of these 215 datasets appears lognormal.
Before identifying trends in the parameters m and s, we check their values for outliers. Outliers are identified by checking the turnover times of each dataset by substituting the inversion into using equation (4.1). While 155 out of these 215 datasets had turnover times less than 50 years, a few datasets were characterized by turnover times extremely large for litter decay, over 1000 years. As shown in figure 5, a histogram of turnover times identified two clusters of datasets: a main cluster of 191 datasets having turnover times less than 1000 years, and 24 datasets having turnover times greater than 1000 years. For the remaining analysis of the parameters m and s, we do not analyse the outlying cluster of 24 datasets with t . 1000 for the reasons discussed in the main text. Inclusion of these datasets in our analysis adds noise to our estimates of m and s but does not change general trends.
A.1.1. Trends in data not well described by a superposition of exponential decays. We present the datasets that do not appear well described by equation (2.1) in table 2 and find that there are noticeable patterns.
The datasets flagged by test 2 (negligible mass loss) and test 5 (turnover time) are associated with slow degradation and the persistence of organic carbon. Table 2 shows that these datasets are predominantly associated with the decay of roots located at the sites CPR, HFR, VCR, NIN, NWT, ARC and BSF. While it is known that roots can be highly recalcitrant [6], it is unclear from this data what controls the persistence of roots; there appear to be no common patterns among the associated sites; they vary highly in climatic parameters and biome type. Although all three types of roots (ANGE, DRGL, PIEL) were also planted at the sites AND, CDR, CWT, HBR, KNZ, LBS, LUQ, LVW, OLY, SMR and UFL, these sites do not contain a single dataset showing extremely slow root decay. There also appears to be no common patterns among these sites.
As discussed in the main text, the datasets filtered by tests 3 (negative curvature) or 4 (trend in residual error) both identify faster than exponential decays. These were consistently located at sites that had low precipitation or were wood, indicating that decay in dry conditions and decay of wood may not always be substrate-limited.

A.2. Regularized inversion
As described in the main text, we assume that a superposition of exponential decays describes a typical decay experiment, wheregðtÞ is the mass fraction remaining at time t and p(k) is the initial probability distribution of first-order decay rates k. Mathematically,gðtÞ is the Laplace transform of a given rate distribution p(k). As mentioned in the main text, distributions p(k) have been hypothesized and evaluated by inserting them into equation (A 1) and comparing the resultinggðtÞ with data [27,28]. In principle, we could calculate the inverse Laplace transform ofgðtÞ exactly to determine p(k). However, the inverse Laplace transform is ill-posed [44,63], meaning that small changes ingðtÞ can result in large changes in p(k), effectively resulting in non-unique solutions. We address this problem as follows.
First, we transform the probability density function p(k) to rðxÞ, where x ¼ ln k: The Laplace transform (A 1) then becomes rðln kÞe Àkt d ln k: ðA 3Þ To solve for r(ln k) from a given set of datagðtÞ, we discretize equation (A 3) as where g is a vector of data pointsgðtÞ of length m, A is the discretized Laplace transform operator of size m Â n and r is the solution vector of length n associated with n prescribed intervals of ln k.
Because the system is both constrained and underdetermined (n . m), we seek the r that minimizes the norm of the residual error. In other words, we find the distribution r that solves min r k g À Ar k subject to the non-negativity constraint  and the constraint that r(ln k) integrates to unity: X n j¼1 A 1j r j ¼ 1: ðA 7Þ This constrained inversion of the Laplace transform, however, remains ill-posed. Regularization techniques are commonly used to solve under-determined and ill-posed inverse problems [44,45,63]. These techniques typically work by attempting to find the 'simplest' solution that fits the data signal but not the noise. Here, we use a Tikhonov regularization technique [44,45,63], which favours smooth solutions by minimizing both the residual error and a quantitative measure of the complexity of r(ln k). Solution complexity is commonly associated with the roughness or the intensity of fluctuations present in the solution. Here, we measure roughness by the squared norm of the first derivative of the solution vector, i.e. drðln kÞ d ln k ¼ X i r iþ1 À r i ln k iþ1 À ln k i 2 ! 1=2 : ðA 8Þ The problem then becomes min r fk g À Ar k 2 þl k Rr k 2 g; ðA 9Þ where R is the bi-diagonal first derivative operator (i.e. the discretization of equation (A 8)) and l is the regularization parameter, which controls the relative weights of the solution roughness and the residual error norm. The method proceeds by identifying the value of l that sets an optimal balance between the residual error k g À Ar k and the roughness k Rr k. A common approach is to use an L-curve technique [44,63]. An L-curve is generated by parametrically varying l and solving equation (A 9) for r(ln k), obtaining values for the residual error norm k g À Ar k and the roughness norm k Rr k for each value of l. The parametric relationship of k Rr k versus k g À Ar k typically has the shape of an 'L'. Although each point on the L-curve may be considered an optimal solution for a certain value of l, the corner of this curve is associated with the regularized solution to the inverse problem. The corner represents the value of l for which increasing l strongly increases solution roughness, while providing little reduction in residual norm, and also for which decreasing l greatly sacrifices residual norm with little reduction in roughness.
Before choosing the corner, we check the solution with l ¼ 0 in case it is just a single delta function, indicating that only one rate is present during decay. If greater than 0.99 of the mass fraction is contained within one rate, we conclude that the simplest solution is an exponential decay having a single rate constant. This occurs for 33 of the datasets.

A.3. Unaveraged parameter analysis
We present here the unaveraged data shown in figures 2 and 3 of the main text. Figure 6 shows how the values of m and s for the 191 LIDET datasets vary with temperature and the lignin-to-nitrogen ratio '/N. Figure 7a,b show the unaveraged variation of kkl and t with '/N. Figure 7c shows the unaveraged plot of s 2 versus m. While there is much scatter among the data, the general trends remain the same as in figures 2 and 3 of the main text.