Retaining eucalyptus harvest residues promotes different pathways for particulate and mineral‐associated organic matter

Eucalyptus plantations have replaced other (agro)ecosystems over 5.6 Mha in Brazil. While these plantations rapidly accumulate carbon (C) in their biomass, the C storage in living forest biomass is transient, and thus, longer-term sustainability relies on sustaining soil organic matter (SOM) stocks. A significant amount of harvest residues (HR) is generated every rotation and can yield SOM if retained in the field. Yet, there is little information on how managing eucalyptus HR changes SOM dynamics. We used isotopic and molecular approaches in a 3-yr field decomposition experiment where a native grassland has been replaced by eucalyptus plantations to assess how HR management practices influence content and chemistry of two distinct SOM fractions [particulate (POM) and mineral-associated organic matter (MAOM)] at two soil depths (0–1 and 1–5 cm). The management practices investigated were HR removal (−R), only bark removal (−B), and retention of all HR (including bark, +B), combined with two levels of nitrogen (N) fertilization [0 (−N) and 200 (+N) kg/ha]. N fertilization inhibited HR decomposition (P = 0.0409), while bark retention had little effect (P = 0.1164). Retaining HR, especially with bark, increased POM-C and MAOM-C content (2.1and 1.2-fold, respectively), decreased POM-δC (1.2-fold), and increased inorganic N retention (1.7-fold) compared with plots where HR had been removed. Inorganic N applications, however, diminished the positive impacts of bark retention. Although the influence of HR management was most pronounced in POM, retaining HR reduced potential soil C mineralization by up to 20%. POM and MAOM chemistry shifted over time and revealed distinct influence of HR on the formation of these fractions. We demonstrate that HR management alters SOM dynamics and that retaining HR, particularly including bark, enhances SOM retention. With continuing conversion of native grassland ecosystems to eucalyptus, long-term sustainability will require careful HR and fertilizer management to balance total biomass harvest with sustaining belowground SOM concentrations.


INTRODUCTION
Brazil is among the largest pulp and sawn wood producers in the world. Eucalyptus plantations occupy~73% of the Brazilian planted forest area of 7.8 Mha and have the highest average productivity of tree species globally (Colodette et al. 2014, IBÁ 2019. In addition to reducing pressure on native forests for timber, these plantations store~435 Mt of carbon (C) in their biomass (Sanquetta et al. 2018) and may, therefore, play a positive role in climate change mitigation. However, eucalyptus forests are managed primarily as short rotations (7-yr rotations, on average), which means C storage in the standing biomass is short-lived. The active management of eucalyptus plantations might increase or decrease soil organic matter (SOM) stocks (Epron et al. 2015, Rocha et al. 2018) and thus play a critical role in longer-term C sequestration. Plants are the primary C source for building SOM, and in short-rotation eucalyptus plantations, most inputs come from the turnover of litter, roots, and harvest residues (HR). The fate of these different C sources and their potential contributions to SOM remain uncertain but are key to understanding the potential for a more sustainable management of these forests.
Whole-tree, stem-only, or stemwood-only harvesting ( Fig. 1) are common harvesting practices that change the amount and quality of HR entering soil (Achat et al. 2015a, b). Eucalyptus HR, i.e., plant tissues with little to no commercial value, include leaves, branches, nonmarketable stem, roots, and bark, and are generated at an average rate of~30 t/ha per eucalyptus rotation, which accounts for roughly 35% of stand biomass (Gatto et al. 2010). Retaining HR in the field can reduce C and nutrient exportation and benefit soil physical (Silva et al. 2007, Jesus et al. 2015 and chemical properties , Laclau et al. 2010b, Kumaraswamy et al. 2014, and consequently tree growth (Laclau et al. 2010a, Versini et al. 2013, Rocha et al. 2016. The fate of eucalyptus HR in the soil, though, is unclear, their effect on SOM stocks is still debated (Kumaraswamy et al. 2014, Epron et al. 2015, Rocha et al. 2016, Souza et al. 2020, and possible effects of HR management on SOM chemistry and stability have rarely been investigated (but see Mathers et al. 2003).
High-productivity eucalyptus forests are commonly nitrogen (N) fertilized at planting. This manipulation of N availability might also affect the C dynamics of the retained material from previous rotations (Berg andMatzner 1997, Vivanco andAustin 2011), especially given several HR components are lignin-and tannin-rich and have a high C:N ratio (Pereira 1988, Kraus et al. 2003, Vázquez et al. 2008, Ferreira et al. 2018b. Fog (1988) proposed that chemical reactions between N and lignin polyphenols during decomposition may produce chemically complex, even toxic, compounds that hamper lignin decomposition. However, there has been no documentation for such reactions in the field. Rather, a wide range of effects of N addition on microbial community structure, physiology, and enzyme activities with feedbacks to lignin decomposition have been demonstrated (Carreiro et al. 2000, Frey et al. 2004, Allison et al. 2007, Rinkes et al. 2016, Zhou et al. 2017). On the other hand, considering the wide C:N ratio of eucalyptus HR (on average >100; Ferreira et al. 2016), N addition may alleviate N limitation and accelerate decomposition. Regardless of the mechanism, N fertilization could strongly influence the fate of litter and soil C in eucalyptus systems.
Fractionating SOM into different functional pools helps us understand the transfer of plant litter to SOM. Combining these techniques with isotope tracers provides for direct measurements of litter and SOM decomposition, transformation, and residence times. Physical fractionation is based on the premise that SOM stability and turnover are controlled by its degree of association with soil mineral particles (silt and clay;von Lützow et al. 2007, Dungait et al. 2012. Particulate organic matter (POM) is characterized by its fast turnover and low stability (Grandy et al. 2007, Gunina andKuzyakov 2014), while the mineralassociated organic matter (MAOM) pool typically has a longer residence time . It is assumed that POM is primarily formed of partially decomposed plant materials, while microbial products predominate in MAOM (Grandy et al. 2007). However, lignin-rich litter may follow progressive decomposition and be preferentially sorbed to soil minerals (Magnússon et al. 2016, Almeida et al. 2018 or soluble components may flow directly from plant litter to the MAOM fraction (Heckman et al. 2013, Sanderman et al. 2014, Cotrufo et al. 2015. Tracing the decomposition and flow of plant litter into SOM pools in forests under contrasting management regimes can provide information leading to more sustainable production forestry. Our objective was to determine whether increases in HR quantity and/or complexity (e.g., whole-tree vs stem-only vs. stemwood-only harvesting) in eucalyptus plantation soil would influence SOM content and stability. To accomplish this, we set up a decomposition experiment in a previous grassland that was recently converted to eucalyptus and used natural δ 13 C variation to quantify HR contribution to POM and MAOM pools. We also examined whether there were changes in the broad chemical patterns of distinct SOM fractions over time to provide more mechanistic insights into shifts in pool size. Finally, we investigated how increasing N availability would change the HR decomposition dynamics and their contributions to SOM pools.

Site description
The experiment was carried out in a commercial eucalyptus plantation located at São Gabriel, Rio Grande do Sul State, in southern Brazil (30°26' S; 54°31' W). The site is representative of the region where there has been a rapid conversion of Pampa native grasslands to eucalyptus plantations (Oliveira et al. 2017). The Pampa biome is a grass-dominated ecosystem, with many herb, shrub, and treelet species co-existing within the grass matrix (Overbeck et al. 2007), that spans~700,000 km 2 in South America within Brazil, Argentina, Uruguay, and Paraguay. We selected a site where the grassland (C 4 -dominated) was replaced by eucalyptus (C 3 ) two months before experimental setup to maximize the differences in 13 C natural isotope abundance between fresh input (C 3 -eucalyptus HR; δ 13 C = −28.5‰) and former SOM (δ 13 C = −13‰). Grassland vegetation was eliminated using a nonselective herbicide (glyphosate). Then, soil was prepared by subsoiling planting rows to 40 cm depth and ridge tilling, and Eucalyptus dunnii seedlings were planted at 2.2 × 3.3 m spacing.
The site is gently sloping at~150 m altitude in a subtropical zone with a mean annual temperature of 18°C and an annual rainfall of 1350 mm. Soil is an Oxyaquic Hapludalf (Soil Survey Staff 1999), formed from sedimentary materials, mainly arenite and siltite (Lemos et al. 1973, Fig. 1. Description of most common forest harvesting systems. v www.esajournals.org Roesch et al. 2009), with the following characteristics (0-10 cm soil layer): 28% clay, 55% silt, and 17% sand; bulk density, 1.15 g/cm 3 ; total porosity, 0.53 m 3 /m 3 ; pH, 4.71; and sum of bases, 7.71 cmol c /dm 3 . The characterization of SOM fractions is presented in Table 1. Detailed site and land-use change descriptions can be found in Ferreira et al. (2018a).

Experimental design
HR used in this experiment came from a two-year-old commercial E. urograndis hybrid stand. HR were separated into leaves, branches, noncommercial stem, bark, and roots. Afterward, they were dried at 45°C in a forced draft oven and chopped into 4-8 cm pieces, except leaves which were left uncut. We combined HR components based on field observations of high-productivity eucalyptus commercial plantations in Brazil (Gatto et al. 2010). The proportions used and HR chemical characterization are shown in Supplemental Table 1. More details about HR chemical composition can be obtained elsewhere , Ferreira et al. 2016.
We used micro-plots of 15 × 15 cm (diameter × height) constructed from PVC tubes as experimental units to study HR decomposition and C transfer to the soil (Appendix S1: Fig. S1). The PVC collars were buried to a depth of 10 cm between planting rows. Six lateral holes (20 mm diameter) were drilled in each tube wall and leveled with the soil surface to allow the free movement of microorganisms, fauna, and water. Any in situ litter was removed from the soil surface inside the tubes. HR were then placed on the soil, and tubes were covered with 1-mm mesh to exclude external litter input. This type of microplot minimizes some problems of litter bags, such as poor contact between litter and soil, fauna exclusion, stolen or lost bags, or the need for greater litter fragmentation to accommodate bag size, while facilitating tracking of decomposed plant residues by concentrating decomposition products (Powers et al. 2009, Shorohova andKapitsa 2014).
The effects of HR composition and N addition on SOM pools were studied in a 3 × 2 factorial experiment with sequential time samplings over three years in four fully replicated random blocks. We tested the effect of removing all residues (−R), removal of bark only (−B), and the retention of all residues (+B). HR dry mass was equivalent to 21.7 and 31.7 t/ha, for −B and +B treatments, respectively (Appendix S1: Table S1). Therefore, both HR quantity and composition were manipulated to represent different harvesting scenarios (Fig. 1). HR dynamics was studied under two levels of N: 0 kg N/ha (−N) and 200 kg N/ha (+N), which was applied as a solution of NH 4 NO 3 and ( 15 NH 4 ) 2 SO 4 enriched at 10% with 15 N tracer. N was applied only once at the start of the experiment. Experimental units were sampled as intact soil cores five times during decomposition, that is, 0, 3, 6, 12, and 36 months, sealed, and transported to the laboratory.

Harvest residues decomposition
In the laboratory, HR were carefully handpicked from the soil surface, brushed free of any adhered mineral particles, and dried in a forced air oven at 45°C. We estimated decomposition rates in −B and +B treatments by fitting a single exponential decay model (Olson 1963) to HR remaining at each time as follows: Notes: Characterization was performed in reference plots (no nitrogen and no residues) at time 0. SOM fractionation following Cambardella and Elliott (1992). POM, particulate organic matter; MAOM, mineral-associated organic matter. C, N, and δ 13 C were determined using an elemental isotope ratio mass spectrometer (EA-IRMS).
† Chemical composition determined with pyrolysis gas chromatography-mass spectrometry (Py-GC/MS) analysis and presented as relative abundance. v www.esajournals.org where X 0 = initial mass (%), X = remaining mass (%) at decomposition time t, and k = decomposition rate per day. We then estimated HR half-life (hl), i.e., the required time to decompose 50% of initial mass, in each treatment through the following equation:

Soil analysis
Micro-plot soils were carefully separated into 0-1 and 1-5 cm layers, sieved to 2 mm, and airdried for further analysis. We selected the uppermost soil layers because they tend to be more responsive to short-term management and to concentrate the first path of aboveground litter incorporation into the soil (Cotrufo et al. 2015, Mitchell et al. 2018). We measured SOM stability, C and N content, and δ 13 C and δ 15 N of SOM fractions on soils sampled at 3 yr, whereas soils sampled at 1 and 3 yr were used to investigate SOM chemistry dynamics.
SOM stability.-We assessed SOM stability by measuring respiration over a~180-d laboratory incubation of soils from 0-1 and 1-5 cm layers. 20 g of soil from each layer, treatment, and field replicate were rewetted to 60% of water-holding capacity (WHC) and placed in 100-mL jars. Jars were kept in a dark environmental chamber held at 25°C and were capped with parafilm to prevent moisture loss while allowing air flux and avoiding CO 2 saturation. Soil moisture was monitored gravimetrically and maintained at 60% of WHC by sprinkling deionized water as needed.
To measure soil respiration rates, jars were decapped, aerated, and re-capped with rubber septa, and a gas sample was immediately collected with a 10-mL syringe and injected into an infrared gas analyzer (LI-820; LI-COR, Lincoln, Nebraska, USA). Jars were then incubated in the environmental chamber for 2-12 h until a second gas sample was collected and CO 2 concentration was determined. During the initial days, jars were capped for 2-3 h to avoid excess CO 2 and respiration was measured daily. As respiration rates decreased, jars were capped for a longer time (ultimately 12 h) and the interval between samplings increased (21 d by the end of incubation). A total of 30 measurements were taken throughout the study. Soil respiration rates were calculated as the difference in CO 2 concentration between pairs of sampling points divided by the time between these pairs. SOM fractionation.-Air-dried mineral soils from 0-1 and 1-5 cm depths, sieved to 2 mm and free of visible plant material, were physically separated into POM and MAOM fractions following a protocol modified from Cambardella and Elliott (1992). Briefly, 5 g of soil was shaken in sodium hexametaphosphate (15 mL; 5 g/L) and glass beads for 16 h at 200 rpm to disperse the soil. Fractions were then separated by gently rinsing the dispersed soil onto a 53-µm sieve with deionized water until clear. The fraction passing through the sieve (<53 µm) was collected as MAOM, while the fraction remaining on the sieve was recovered as POM. Fractions were dried in an oven at 60°C until constant mass. Afterward, POM and MAOM were weighted and finely ground with a ball mill for C, 13 C, N, 15 N, and chemical characterization.
C, δ 13 C, N, and δ 15 N.-C, 13 C, total N, and 15 N were determined in fractionated samples using an elemental isotope ratio mass spectrometer (EA-IRMS GSL 20-20; Sercon, Crewe, UK). Reference gas was calibrated with Pee Dee Belemnite (PDB) and atmospheric air N 2 certified standards for isotope signature calculations. The abundance of 13 C and 15 N in samples was calculated as follows: where R = 13 C/ 12 C or 15 N/ 14 N ratios.
Based on δ 13 C and δ 15 N, we calculated the proportion of C and N in each fraction derived from HR (HR-C) and from the fertilizer (fN m ), respectively, as follows: where δ 13 C f is the δ 13 C from each treatment, replicate, and fraction, δ 13 C Ref is the reference value, herein the δ 13 C of POM or MAOM of the treatments without residues (−R) of the 3-yr. sampling, which varied upon the treatment compared (−N or +N), and δ 13 C HR is the average δ 13 C of the HR used; and v www.esajournals.org fN m ¼ ðδ 15 N f À 0:3663Þ=ðδ 15 N m À 0:3663Þ (5) where δ 15 N f is the δ 15 N from each treatment, replicate, and fraction of the N-amended treatments, 0.3663 is the N 2 air reference value, and δ 15 N m is the enrichment of the N fertilizer solution (10%). The contribution of HR-C to each SOM fraction at each soil depth, hereafter referred to as SOM yield, was calculated as follows: where C is the carbon content determined at each SOM fraction and soil depth, and C HRad is the HR decomposed (expressed in C) in each treatment at 3 yr. Pyrolysis gas chromatography/mass spectrometry. -To assess how treatments influenced POM and MAOM chemistry over time, we performed pyrolysis gas chromatography-mass spectrometry (Py-GC/MS) analysis on soil samples from 0, 1, and 3 yr following previously described protocols (Grandy et al. 2007, Wickings et al. 2012, Rinkes et al. 2016. In brief, samples were pyrolyzed with a CDS Pyroprobe 5150 pyrolyzer (CDS Analytical, Oxford, Pennsylvania, USA) at 600°C for 20 s and then transferred to a Thermo Trace GC Ultra Gas chromatograph (Thermo Fisher Scientific, Austin, Texas, USA) where they were separated on a fused silica column over the course of approximately 60 min with a starting temperature of 40°C and final temperature of 310°C. Compounds were then ionized via ion trap on a Polaris Q mass spectrometer (Thermo Fisher Scientific).
The Automated Mass Spectral Deconvolution and Identification System (AMDIS, V 2.65) and the National Institute of Standards and Technology (NIST) compound library were used to analyze and identify peaks. Peaks of each sample were expressed as the relative abundance of total sample peak area. Py-GC/MS products were grouped into seven classes: aromatics, lignin, lipids, N-bearing compounds, polysaccharides, protein, and unknown origin.

Data analysis
Data were tested for normality with the Shapiro-Wilk test and ln-or Box-Cox-transformed when necessary to meet normality assumptions.
Afterward, a two-way ANOVA was performed to assess HR decomposition with N and Residues as main factors. For SOM properties and soil respiration, a three-way ANOVA was performed separately for each SOM fraction with Depth, Residue, and N as main factors. ANOVA for SOM stability was performed on the cumulative soil respiration at the end of incubation. We performed post hoc pairwise comparisons on significant effects and interactions using Fisher's LSD test.
A multivariate description of Py-GC/MS products was performed using principal component analysis (PCA) to obtain chemical fingerprints of POM and MAOM at each treatment, depth, and time. We calculated the Euclidean distance matrixes on the seven classes of Py-GC/MS products; these were grouped to create a PCA plot. Permutational multivariate analysis of variances using distance matrices (PERMANOVA-ADO-NIS) was used to identify statistical differences in chemical composition between Fractions, Depth, Time, Residue, and N. All analysis was performed in R 3.03 (R Development Core Team 2019).

Decomposition of HR
Bark presence did not significantly alter HR field decomposition rates, although −B treatments had an average hl of~85 d longer than +B (P = 0.1164). +N treatments had an average hl of~117 d longer than −N (P = 0.0409; Fig. 2). The N × Res interaction was not significant (P = 0.3805), but we found that +N had greater effect on −B treatment alone (P = 0.0442), increasing HR-hl by 162 d, than on +B treatment, where it increased HR-hl by 71 d (P = 0.3282).
C and N content, and δ 13 C, and SOM yield of particulate and mineral-associated organic matter fractions All variables differed by soil depth, except for the C:N ratios of POM and MAOM. Both HR and N addition influenced SOM properties after 3 yr (Table 2). Changes in C content and δ 13 C were most pronounced in the 0-1 cm layer and the POM fraction (Fig. 3). Overall, HR presence increased C content and decreased δ 13 C values, v www.esajournals.org while N addition reduced δ 13 C in −R treatment and lessened +B effect on soil C accretion.
For POM-C, HR removal reduced C content and increased δ 13 C (less negative, less C 3 material) in both soil layers (Fig. 3A,B). In the 0-1 cm soil layer, −N + B resulted in the highest POM-C content (16.71 g/kg), while −N-B, +N −B, and +N+B had similar POM-C content (~11.71 g/kg). These numbers are 2.1 and 1.5 times higher, respectively, than the POM-C content in −R treatments at 0-1 cm. All treatments with HR had similar POM-δ 13 C values, which were 1.2 times  Notes: C, carbon; C:N, C:N ratio; fN m , fraction derived from mineral-applied N calculated based on 15 N content; HR, harvest residue; MAOM, mineral-associated organic matter; N, nitrogen; POM, particulate organic matter; δ 13 C, 13 C abundance. ". . ." represents P > 0.05; n/a, not applicable; ↑, ↓ indicate increase or decrease in property value, respectively, with N addition or HR retention (main effects). more negative than in −R treatment. N addition alone did not influence POM-C content and δ 13 C; however, it increased δ 13 C in −R treatments (less negative) and reduced +B effect on POM-C content in both soil layers. In the 1-5 cm depth, POM-C in −R and −B did not differ and was 2.2 times lower than in +B (for δ 13 C, −B and +B were statistically equal) if no N was applied. With N addition, −B and +B did not differ, but they had 40% more POM-C than −R. For MAOM-C, results were of lesser extent but in line with the changes observed for POM-C with HR management (Fig. 3C,D).
HR-C proportion in POM and MAOM ranged from 12 to 32 and 0% to 5%, respectively (Fig. 4). Bark retention increased HR-C in POM in the 0-1 cm depth by 30% on average. In MAOM, the greatest HR-C proportion (4.7%) was observed in the 0-1 cm layer with N addition when bark was retained. Nitrogen addition increased SOM yield of −B treatments across both SOM fractions and soil depths. HR effects on soil N were only observed in -N treatments, in both depths for POM, and in the 1-5 cm depth for MAOM (Fig. 5). +N seemed to counterbalance the positive effect of HR on POM-N in the 0-1 cm soil layer (Fig. 5A). On average, the proportion of N derived from the fertilizer (fN m ) was 1.7 and 1.8 times higher in the uppermost soil layer, for POM and MAOM fractions, respectively, compared with fN m in the 1-5 cm soil depth. The fN m was~40% greater in the POM than in the MAOM. The fN m was higher with HR maintenance in both fractions and soil depths. Both −B and +B had 1.6 times greater fN m than −R in the POM fraction in 0-1 cm depth, while in the MAOM −B and +B had 3% and 20% greater fN m than −R, respectively. In the 1-5 cm depth, −B had 1.3 and 1.4 times greater fN m than −R for POM and MAOM, respectively, while +B showed 1.7 times greater fN m than −R in both fractions.

Py-GC/MS of SOM fractions
The relative abundance of Py-GC/MS products is summarized in seven classes and is presented in Appendix S1: Table S2 (0-1 cm depth) and Appendix S1: Table S3 (1-5 cm depth). Because depth strongly influenced SOM fraction chemistry (P < 0.0001), we present data separately for each layer (Fig. 6A,B). The PCA revealed a pronounced difference in chemistry between the two SOM fractions, which also shifted over time. Lignin represented a higher proportion of POM-C in treatments with retained HR, particularly at 0-1 cm depth, and when including bark. However, distinct patterns were observed for −N and +N treatments. In −N treatments, lignin increased from 16% to 23% and 29% to 32%, in −B and +B, respectively, while in +N lignin decreased from 21% to 16% and 27% to 21%, in −B and +B treatments, respectively (Appendix S1: Table S2). When HR were removed, lignin decreased throughout the experiment, irrespective of N addition, in both soil depths and SOM fractions (Appendix S1: Tables S2, S3). Lignin played a crucial role in distinguishing 1-and 3-yr sampling times for the MAOM fraction in the 1-5 cm depth (Fig. 6B) despite its low relative abundance overall (0-2%), likely because lignin was almost absent in the 3-yr samples (Appendix S1: Table S3).
Lipid contribution to POM was higher in the 1-5 cm than in the 0-1 cm depth, while it contributed similarly to the MAOM fraction at both depths (Appendix S1: Tables S2 and S3). Over time, aromatic compounds increased in POM and decreased in MAOM at both depths. Polysaccharides and proteins varied similarly over time; they decreased in the POM fraction while increasing in MAOM, regardless of treatment.
MAOM from N-fertilized treatments at the 0-1 cm had similar SOM chemistry profiles (Fig. 6A). Treatment differentiation was more evident in the POM fraction, suggesting a constant HR flows to this fraction (Fig. 6A). Lipids were more associated with MAOM in the 0-1 cm, but more correlated with POM in the 1-5 cm layer (Fig. 6A,B). Polysaccharides, proteins, and N-bearing compounds were closely linked with the MAOM fraction in both depths, while compounds of unknown origin were more associated with POM.

SOM stability
Respiration rates and potentially mineralizable C did not differ by HR or N treatment when compared on a total soil mass basis (Fig. 7A,B), but soils from the 0-1 cm layer had higher cumulative respiration than soils from the 1-5 cm layer across treatments (P = 0.0065). When respiration was expressed on a soil C basis (inset Fig. 7A,B), no difference was observed between the two soil depths (P = 0.0943), but HR and N treatments did influence the cumulative respiration. In the 0-1 cm depth with no N addition, +B had the lowest relative respiration rate, while −B and −R did not differ (Fig. 7A, inset). In the 1-5 cm depth, +B also had a lower relative respiration rate than −B, but −B and +B did not differ from −R (Fig. 7B, inset). With N addition, in turn, −R yielded an average cumulative CO 2 across both depths 1.2-fold greater than both −B and +B, which did not differ from each other.

HR decomposition
Although litter chemistry is known to influence decomposition (Hernández et al. 2009, Vivanco and Austin 2011, Ferreira et al. 2016, we observed that HR decomposition rates were not significantly altered by the presence or absence of bark (Fig. 2). In contrast, N addition slowed decomposition. These results are the opposite of those found in the same site and treatments during initial decomposition stages, i.e., first year (Oliveira et al. 2021), when bark presence accelerated decomposition while N had no influence. The previously observed accelerated decomposition of HR with bark could have been due to faster decomposition of the bark itself , or the other components when bark was present (e.g., a priming effect of bark on these other residues), or both. We cannot separate these mechanisms, but our results suggest that bark effects were most pronounced during initial decomposition stages and became less important as decomposition progressed. On the other hand, our observation that N addition inhibited litter decomposition parallels other studies of low-quality (or ligninrich, high C:N ratio) materials (Knorr et al. 2005) or in late decomposition stages (Fog 1988). This inhibitory effect of N should be particularly  strong in bark, since it is a complex material with high contents of tannins, polyphenols, and lignin (Kögel-Knabner 2002, Pinto et al. 2013, García et al. 2014). Yet, we observed a significant effect of N addition only in the absence of bark (−B), which thus questions the hypothesis of chemical reactions between litter lignin or tannins and N as the main mechanism slowing down decomposition (Fog 1988).
Bark often has the slowest decomposition rates observed among eucalyptus HR components (Shammas et al. 2003, Epron et al. 2006, Hernández et al. 2009). However, recent studies have shown that bark retention might accelerate eucalyptus HR decomposition in nutrient-poor tropical sites because (1) it creates a physical protection that results in optimized conditions for decomposition, for example, higher humidity and protection against disturbances; and (2) its high nutrient content, especially calcium, offsets any possible inhibitory effect of chemical complexity , Ferreira et al. 2016). In our study, bark had a negligible effect on HR decomposition. We posit that this was possibly because our experimental site was located in a subtropical area with a less weathered and more nutrient-rich soil than most eucalyptus-dominated tropical soils, and thus, microorganism activity may not have been as limited by nutrient availability.
C content, δ 13 C, SOM yield , N content, and fN m in SOM fractions and stability of SOM Removing HR resulted in the lowest observed soil C concentrations (Fig. 3). Our results demonstrate that POM-C was up to~130% lower in -R plots in both 0-1 and 1-5 cm depths than in plots where HR were retained. MAOM-C was also lower in −R plots, but proportional differences were lower (up to 20%; Fig. 3). These results suggest a progressive POM-C and MAOM-C impoverishment when HR are removed, since early results from this experiment (1 yr of decomposition) indicated that POM-C was~35% lower in −R compared with more conservative HR management practices, while no differences were observed for MAOM-C (Oliveira et al. 2021). Epron et al. (2015) found a POM-C loss of 40% in the 0-5 cm soil depth with HR removal in a Congolese eucalyptus plantation, while Rocha et al. (2018) reported soil C losses of 50%, which were intensified over consecutive eucalyptus rotations in Brazil. Our numbers are slightly higher than those referred studies but could be explained by the long-term use of PVC collars, which prevented other inputs for a prolonged time. However, paralleling these studies, our results suggest an active role of HR to sustaining soil C concentrations in eucalyptus plantations.
Interestingly, δ 13 C in −R plots declined through time (Table 1; Fig. 3), indicating an unexpected C 3 contribution, which could possibly be due to fauna and microorganisms introducing C 3 -C from the surrounding eucalyptus plantation since direct root and litter input into the experimental units were excluded (Appendix S1: Fig. S1), or C kinetic isotope fractionation and preferential utilization by soil microbes (Feng 2002, Werth andKuzyakov 2010). In both soil fractions, while similar C content was observed in −R treatments (either −N or +N), a more negative δ 13 C was observed in −N−R than in +N−R (Fig. 3). The less negative δ 13 C in +N−R may be explained by +N stimulating changes in microbial community structure with consequences for either the incorporation of eucalyptus C from surrounding areas into the plots Neff 2008, Santana et al. 2015) or C consumption (Baumann et al. 2009, Oliveira et al. 2021.
Isolating the HR-C in the soil from other C sources of aboveground litter and root turnover originating during a multiyear perennial crop rotation is a technical research challenge. By setting up the experiment at the beginning of a eucalyptus rotation over a C 4 -grassland, we were able to quantify HR transference to SOM fractions through the expressive differences in 13 C natural isotope abundance (Δ 13 C~15‰). Soil conditions at the end of a eucalyptus rotation, i.e., when HR would be deposited on the soil, might differ from those here evaluated, which, in turn, would hamper the use of isotopes to quantify HR-C contribution to SOM. With this caveat in mind but considering that direct inputs from external sources were avoided, we believe our results are representative of HR impact per se on SOM pools, except for roots. To our knowledge, there has been no direct quantification of rootderived C contribution to soil C across multiple eucalyptus rotations to date. However, in line with a growing body of evidence (Rasse et al. 2005, Kramer et al. 2010, Schmidt et al. 2011, Tefs and Gleixner 2012, Berhongaray et al. 2019, it is expected eucalyptus root-derived C to be a major source of soil C across rotations, particularly to deeper soil layers; thus, root contribution to SOM was probably underestimated with our experimental design.
HR retention increased C and N content in both SOM fractions (Figs. 3, 5). Effects of HR on soil C were more pronounced in the POM than MAOM fraction, particularly in the 0-1 cm depth (Fig. 4). While this thin layer may be more susceptible to external impacts and not representative of bulk soil C stocks, it helps understanding the dynamics of C incorporation into the soil (Cotrufo et al. 2015, Mitchell et al. 2018. POM-C and MAOM-C concentration increases were higher with more conservative harvesting, eg., when all HR including bark were retained (Fig. 3). N addition, however, diminished the positive impacts of bark retention though it did not affect bark breakdown (Fig. 2), which suggests an interaction between N and bark compounds-possibly tannins (Kraus et al. 2003)-moderating SOM accrual but not decomposition. HR-C was primarily transferred to POM, accounting for 20-32% and 6-20% of POM-C in 0-1 and 1-5 cm layers, respectively, while HR-C proportion in MAOM was <5% (Fig. 4). Comparing with early decomposition stages (first year, 50% of HR decomposed; see Oliveira et al. 2021), contributions of HR-C to POM-C of 0-1 cm layer were similar, whereas POM-C concentrations and HR-C contributions to MAOM-C were increased after 3 yr, suggesting a constant and active flow of HR-C to SOM fractions. In addition, potential mineralization, i.e. cumulative respiration on a soil C basis, was lower in the treatments with HR, particularly those with bark (Fig. 7A, B), indicating that HR-C was relatively more stable, even though HR-C was predominantly found in POM-C. Our results parallel other studies showing that POM of topsoil layers is more responsive to HR management (Epron et al. 2015, Souza et al. 2020), but also confirm the potential for HR to contribute to more stable SOM fractions (Hicks Pries et al. 2017, Cusack et al. 2018, Souza et al. 2020, supporting the idea that retaining HR in the field is a valuable management practice to build and increase relative SOM stability in eucalyptus plantations. Most HR components are N-poor due to efficient internal N cycling during tree growth (Laclau et al. 2010b), and thus, they are not expected to contribute substantially to soil N pools. Surprisingly, POM-N content increased with bark retention, although N fertilizer offset these gains because it increased N content across both SOM fractions in +N treatments (Fig. 5). A substantial amount of applied N was not recovered in either HR (Ferreira et al. 2016) or SOM fractions, which suggests that much of the applied N might have been lost (leached or gaseous). Yet small (<4%), a higher proportion of the added N remained in the soil if HR were retained (Fig. 5). Since N fertilization at planting is a common practice in eucalyptus plantations in Brazil, retaining HR in the field might increase N retention and help to build soil N stocks and sustain eucalyptus tree growth through the rotation (Pulito et al. 2015).

Chemical fingerprint of SOM fractions
We characterized different SOM fractions in a decomposition chronosequence (1-and 3-yr) to understand the chemical transformations that occur when HR is incorporated into soil C. A recent incubation study with eucalyptus extracts indicated that lignin-rich litter may follow progressive decomposition and be preferentially sorbed to soil minerals (Almeida et al. 2018). We hypothesized that if this decompositionstabilization pathway also applies to in situ conditions, MAOM chemistry at 3 yr would resemble POM chemistry at 1 yr, or at least share many of the dominant compounds. This was not supported by our Py-GC/MS data; rather, POM and MAOM had distinct chemical compositions that remained distinct over time (Fig. 6). Since POM contained proportionally more lignin, we infer a progressive decomposition and direct transference from HR to this fraction. MAOM, in turn, was richer in polysaccharides, protein, and lipids indicative of microbial product dominance Neff 2008, Kallenbach et al. 2016). MAOM chemistry, however, varied by HR treatment, suggesting that HR had an indirect yet crucial role in C stabilization by providing energy for microbial processes and C incorporation into microbial biomass and mineral-stabilized v www.esajournals.org products (Cotrufo et al. 2013, Bradford et al. 2013. Indeed, eucalyptus HR retention might increase soil microbial biomass-particularly of fungi if bark is retained (Oliveira et al. 2021), and bacteria tend to forage near fresh POM (Schlüter et al. 2019), which together might result in greater microbe-derived MAOM in the long term (Kallenbach et al. 2016).
Adding N fertilizer reduced bark contributions to POM-lignin, POM-C, and MAOM-C concentrations in the 0-1 cm soil depth (Appendix S1: Table S2; Fig. 3A, C), while N effect on the decomposition of HR with bark was negligible (Fig. 2). Collectively, these results suggest active lignin breakdown but reduced lignin accumulation in the soil. Therefore, we hypothesize that N-induced shifts in microbial structure and lignin-degrading enzyme activity are more important to lignin fate in the soil than abiotic chemical reactions between N and lignin hindering its processing (Fog 1988, Rinkes et al. 2016. Given that treatments with higher POM-lignin (Appendix S1: Tables S2, S3) had reduced cumulative respiration per unit of soil C (Fig. 7), retaining HR, and thus increasing lignin content in the soil, may enhance SOM short-term stability. On the other hand, the almost complete absence of detectable lignin in 3-yr samples of 1-5 cm depth in both SOM fractions (Appendix S1: Table S3) suggests that aboveground lignin transfer and accrual are limited to the contact area with the soil (e.g., uppermost layer), and corroborates that lignin is not stable in soils in the long term (Lobe et al. 2002, Gleixner et al. 2002, Rasse et al. 2006, Thevenot et al. 2010, Klotzbücher et al. 2016).

CONCLUSIONS
The high productivity of eucalyptus plantations in Brazil places the country in a prominent position worldwide for wood supply. The sustainability of long-term production in nutrient-poor soils in the tropics relies on sustaining SOM stocks. Retaining HR in the field reduces C and nutrient exportation and can provide multiple benefits to forest soils. We demonstrate that actively managing eucalyptus HR alters SOM dynamics with potential implications for long-term C storage. Retaining HR, particularly when including bark-a more conservative harvesting-helps to build SOM of uppermost soil layers. Inorganic N applications, however, diminished the positive impacts of bark retention. Although the influence of HR management was most pronounced in POM, retaining HR reduced potential soil C mineralization. The differences between POM and MAOM chemistry over time revealed a distinct influence of HR on the formation of these fractions and provide mechanistic insights into how lignin-rich crop residues promote SOM. With continuing conversion of native grassland ecosystems to eucalyptus, long-term sustainability will require careful HR and fertilizer management to balance total biomass harvest with sustaining belowground SOM concentrations. While our conclusions remain to be tested across broader scales in longer-term studies (e.g., across multiple rotations in different regions), we provide quantitative evidence for the potential of simple management strategies to promote more sustainable forest plantations.