Belowground impacts of perennial grass cultivation for sustainable biofuel feedstock production in the tropics

Perennial grasses can sequester soil organic carbon (SOC) in sustainably managed biofuel systems, directly mitigating atmospheric CO2 concentrations while simultaneously generating biomass for renewable energy. The objective of this study was to quantify SOC accumulation and identify the primary drivers of belowground C dynamics in a zero‐tillage production system of tropical perennial C4 grasses grown for biofuel feedstock in Hawaii. Specifically, the quantity, quality, and fate of soil C inputs were determined for eight grass accessions – four varieties each of napier grass and guinea grass. Carbon fluxes (soil CO2 efflux, aboveground net primary productivity, litterfall, total belowground carbon flux, root decay constant), C pools (SOC pool and root biomass), and C quality (root chemistry, C and nitrogen concentrations, and ratios) were measured through three harvest cycles following conversion of a fallow field to cultivated perennial grasses. A wide range of SOC accumulation occurred, with both significant species and accession effects. Aboveground biomass yield was greater, and root lignin concentration was lower for napier grass than guinea grass. Structural equation modeling revealed that root lignin concentration was the most important driver of SOC pool: varieties with low root lignin concentration, which was significantly related to rapid root decomposition, accumulated the greatest amount of SOC. Roots with low lignin concentration decomposed rapidly, but the residue and associated microbial biomass/by‐products accumulated as SOC. In general, napier grass was better suited for promoting soil C sequestration in this system. Further, high‐yielding varieties with low root lignin concentration provided the greatest climate change mitigation potential in a ratoon system. Understanding the factors affecting SOC accumulation and the net greenhouse gas trade‐offs within a biofuel production system will aid in crop selection to meet multiple goals toward environmental and economic sustainability.


Introduction
Global interest in the development of renewable energy systems stems largely from combined concerns about climate change, energy security, and environmental sustainability (Tilman et al., 2009). Biofuel production from plant feedstocks is a promising component of diverse, renewable energy plans (Liska & Perrin, 2009;Werling et al., 2014). However, biofuel production can result in net CO 2 emissions through the use of fossil fuels during land-use conversion and the production and processing of biomass (Gibbs et al., 2008;Cherubini et al., 2009). The use of conservation management practices during biofuel cultivation could sequester soil organic carbon (SOC) (Lal, 2013), providing an offset to some or all of the CO 2 emissions associated with biofuel feedstock production (Adler et al., 2007;Davis et al., 2013). However, considerable uncertainty exists regarding the magnitude and direction of changes in SOC under various biofuel feedstock production systems (Stockmann et al., 2013).
The SOC pool is a function of the dynamic balance between inputs and outputs of C belowground, with SOC increasing when inputs exceed outputs (Six & Jastrow, 2002;Cotrufo et al., 2015). Inputs of C to belowground include aboveground and belowground litter and total belowground carbon flux (TBCF), which is the sum of autotrophic C input to belowground to support root production and respiration, exudates, and microbial symbioses (Giardina & Ryan, 2002). The quantity of belowground C input is influenced by overall stand productivity and the partitioning of photosynthetically fixed C to aboveground vs. belowground (Litton et al., 2007). Primary C outputs from belowground include soil CO 2 efflux ('soil respiration') and, in some systems, leaching. Soil CO 2 efflux, the primary belowground output of C in most systems, can be decreased through reduced tillage or increased chemical recalcitrance of inputs (Davidson & Janssens, 2006).
The use of perennial grasses for biofuel production has the potential to mitigate SOC loss and greenhouse gas (GHG) emissions within the production system due to high yield and low requirements for fertilizer, pesticides, and irrigation (Liebig et al., 2008;DeLucia, 2016). In tropical and subtropical regions, perennial C4 grasses such as Saccharum officinarum (sugarcane), Saccharum spontaneum (energy cane), Pennisetum purpureum (napier grass), Megathyrsus maximus (guinea grass), and P. purpureum 9 Pennisetum glaucum (sterile napier grass hybrids) have been identified as promising candidates for biofuel feedstocks, due primarily to their potential for high biomass yields (Tran et al., 2011;Hashimoto et al., 2012;Meki et al., 2014;Mochizuki et al., 2014). Most assessments of tropical perennial C4 grasses to date have focused on aboveground yields and agronomic requirements (Kinoshita et al., 1995;Keffer et al., 2009). In turn, very little attention has been given to belowground C dynamics associated with these grasses even though SOC accrual could potentially provide a strong C sink to offset nonrenewable GHG emissions associated with their production (Gelfand et al., 2013;DeLucia, 2016).
In addition to the potential for high yield and low agricultural inputs, perennial grasses under zero-tillage cultivation have the potential to store a large amount of SOC compared to annual crops due to their extensive belowground root system (Powlson et al., 2011;Anderson-Teixeira et al., 2012). Napier grass and guinea grass species produce large stocks of root biomass, both as fine roots in the surface root zone, and at depths of up to 4.5 m (Khanal et al., 2010). Perennial grasses can be harvested by ratooning, a form of zero-tillage harvest that leaves the lower part of the plant and living roots and soil undisturbed, thereby potentially contributing to the sustainability of the production system via SOC accumulation, erosion control, and improved soil fertility . Compared to annual crops such as corn and soybeans, perennial species such as napier grass and guinea grass can be ratooned for as many as 4 years without reduction in production (Samson et al., 2005). In temperate ecosystems, the capacity for perennial grass feedstocks to mitigate climate change is fairly well documented (Anderson-Teixeira et al., 2012), but whether this will be the case in tropical regions is unknown.
Both the quantity and quality of plant inputs affect the balance of SOC. A recent meta-analysis revealed that differences in SOC pools between tillage and zero-tillage soils were due primarily to differences in the quantity of C inputs (Virto et al., 2012). However, in other cases, input quantity itself has no impact on the SOC pool (Al-Kaisi & Grote, 2007;Sanderson, 2008), suggesting a predominant influence by other factors such as input quality (i.e., decomposition dynamics). Traditionally, inputs with lower overall chemical quality and therefore lower decomposability were thought to contribute to the accumulation of SOC more than inputs with higher chemical quality due to the complex structure of low-quality molecules (Jastrow & Miller, 1996). Grandy & Neff (2008) reported that the majority of recent, plant-derived SOC was comprised of ligninrelated compounds. Recent concepts of litter transformation to SOC include rapid and slow modes of transport of different plant compounds into the mineral soil (Cotrufo et al., 2013(Cotrufo et al., , 2015, and multiple mechanisms of stabilization involving the protection of inputs from microbial degradation (Schmidt et al., 2011;Lehmann & Kleber, 2015). Ample evidence suggests that plant tissue decomposition rate, determined by the decay constant (k), is related to tissue chemical characteristics such as C : N and lignin : N (Melillo et al., 1982;Silver & Miya, 2001;Johnson et al., 2007). And, in at least one study on switchgrass, SOC positively correlated to root C : N, a ratio often used as a predictor of root decomposition (Ma et al., 2000).
The overarching objectives of this research were to quantify change in SOC over time and to identify the primary drivers of belowground C dynamics in a ratoon production system of eight perennial grass accessions under consideration as biofuel feedstocks in Hawaii. It was expected that SOC would increase following cultivation as a result of high belowground organic matter inputs. A series of hypotheses were developed to determine the primary drivers of increased SOC. First, it was hypothesized that the selected perennial grass systems would have a wide range of ecosystem C pools, fluxes, and input quality indicators that contribute collectively to the soil C balance (H1). Measured ecosystem C pools and fluxes included aboveground biomass yield, soil and root C, soil CO 2 efflux, total belowground C flux, and changes in pools over time. Input quality indicators included concentrations of cellulose, hemicellulose, lignin, C, N, and C/lignin to N ratios in roots. A negative relationship between decomposition rate constant (k) and litter lignin concentration (Johnson et al., 2007;Zhang et al., 2008) and lignin : N (Melillo et al., 1982) has been widely reported. Therefore, it was hypothesized that grasses with higher root lignin, C : N ratio, and lignin : N ratio would result in lower decomposition rates (H2). Finally, it was hypothesized that high aboveground biomass yield and high root lignin concentration would be the primary, direct drivers of SOC accumulation in this zero-tillage, tropical perennial grass system. Further, root biomass, root decay constant, and soil CO 2 efflux were expected to mediate indirect controls of these factors on SOC accumulation (H3). H3 was assessed using structural equations modeling (SEM), a multivariate approach that develops a conceptual model of causal relationships among the quantity and quality of C pools and fluxes that ultimately drive SOC accumulation.

Field site and experimental design
The study was conducted at the University of Hawaii's Experimental Research Station at Waimanalo, on the Island of Oahu, Hawaii (21°20 0 15″N, 157°43 0 30″W). The plots were located on alluvial fans at 30 m elevation with a mean annual temperature of 24.6°C and mean annual precipitation of 938 mm, 72% of which occurs between November and April (Giambelluca et al., 2013). The soil is classified as silty clay with smectitic and halloysitic mineralogy (Waialua series, very fine, mixed, superactive, isohyperthermic Pachic Haplustolls) (Soil Survey staff, accessed 11/1/14). For 24 years prior to the start of this experiment, the field site was a fallow, mechanically maintained grassy field. Information before 24 years is limited, but the presence of a dense plow layer suggests past extended periods of heavy cultivation.
The field plot design was established in October 2009 and consisted of a randomized complete block design with four replicates of eight grass accessions. The accessions studied were three napier grasses, one pearl millet 9 dwarf napier grass (PM9D) cross hybrid, and four guinea grass accessions (Table S1). These accessions were selected to encompass a wide range of aboveground yield potential based on preliminary site data. Each grass plot consisted of four rows planted in a 2 9 3 m area, with a nonplanted buffer of 0.6 m between plots. Inter-rows and buffers were covered with nylon mats to suppress weed growth. Plants were produced from stem cuttings and irrigated 3 days a week from 9 : 00 to 14 : 00 with drip irrigation for the duration of the study at~3.6 L m À2 day À1 . Fertilization with granular triple-16 formulation at the rate of 53 kg N, 23 kg phosphorus, and 44 kg potassium ha À1 occurred when the grasses were first planted in November 2009, and after the second ratooning in November 2010. The grasses were ratoon harvested in March and November 2010, and July 2011. At each harvest, grasses were cut at 10 cm above the soil surface using a brush-cutter, and all aboveground biomass was removed from the plot. Material harvested from the inner two rows, excluding the outer 0.5 m of those rows on each end, was considered the experimental unit for analysis.

Ecosystem carbon pools and fluxes
Total wet weight of the harvested material was recorded, and a subset of the material was dried at 105°C to determine moisture content, from which dry weights of the total harvested material (i.e., aboveground biomass yield) were calculated. Subsamples of leaf tissue were oven-dried at 75°C, homogenized with a ball mill (Retsch MM200 mixer mill; Retsch GmbH, Haan, Germany) to pass through a 250-lm sieve, and the C and N concentrations of leaf samples were determined by oxidative combustion on an elemental analyzer (Costech ECS 4010 CHNSO Analyzer; Costech Analytical Technologies Inc., Valencia, CA, USA). Aboveground biomass yields from the three ratoons were combined, multiplied by C concentrations of leaf materials, and averaged across years for an annualized biomass C value (g C m À2 yr À1 ).
Two perforated aluminum pans (16 9 26 9 2 cm) were placed in the inter-rows in each plot to measure aboveground litterfall (F A ). Litterfall measurements were made monthly from February to July 2011. It was assumed thereafter that, because yields over time were steady, the litterfall rate during the study period remained constant. Thus, the average monthly litterfall rate from the five measurements was multiplied by 12 to determine annual cumulative F A (g C m À2 yr À1 ). Collected litter was oven-dried at 75°C to a constant mass and weighed. Elemental C and N composition of litter was determined as described above.
Soil CO 2 efflux (F S ) was measured monthly for 1 year from August 2010 to July 2011 between 9 : 00 and 18 : 00 from each plot using a LI-6400XT portable photosynthesis system with a soil respiration chamber (LI-COR Inc., Lincoln, NE, USA). In each plot, F S was measured on five 10-cm-diameter polyvinyl chloride (PVC) collars that were inserted 2 cm into the soil 1 month prior to initial measurements. This was reduced to three collars plot À1 after November 2010 based on a lack of within plot variability in F S . Living vegetation inside the collars was clipped and removed 1 day prior to measurements. Soil temperature at 10 cm depth was measured adjacent to each collar at the time of F S measurements with a temperature probe. Volumetric soil moisture was also measured at 5 cm depths adjacent to each collar using an impedance probe calibrated to the study site soil (Hydra Soil Moisture Probe; Stevens Water Monitoring Systems Inc., Beaverton, OR, USA).
Annual cumulative F S (g C m À2 yr À1 ) was calculated using linear interpolation between monthly measurements (Litton et al., 2008). Mathematically, the interpolation can be expressed as: where F i is monthly soil CO 2 efflux rates (g C m À2 day À1 ) and D i is date of the efflux measurement, starting at 1 in July 2010 and ending as i in August 2011. Daily F i values were not corrected for diel variability because flux measurements taken every 2 h from 11 : 00 hours to 3 : 00 hours in one replicate of NG1 revealed no diel variability in F S (Tukey multiple comparison with 95% confidence) and no relationship between F S and soil temperature (r = À0.030; P = 0.704; n = 160). Previous studies have also documented lack of diel patterns in F S in tropical dry forest (Litton et al., 2008) and tropical wet forest in Hawaii (Giardina & Ryan, 2002;Litton et al., 2011). Gomez-Casanovas et al. (2013) found that linear interpolation between survey measurements was the second best method overall for estimating cumulative F S , with only linear interpolation supplemented by site-specific soil temperature dependence of F S outperforming this approach. In our study, there was no temperature dependence of F S (P = 0.704), making linear interpolation between monthly measurements the most robust approach to annual estimates of F S . Litton et al. (2008) did document an increase in F S following a large pulse precipitation event (95 mm in a system receiving~750 mm annually) that, if ignored, resulted in an underestimate of F S when using linear interpolation between monthly measurements. However, this underestimate was very small on an annual basis (0.5-7%). In addition, the focus herein is on comparison between treatments, and not with other studies, and any bias in F S resulting from the use of linear interpolations between monthly measurements would be the same across treatments. Samples for SOC (C S ) and root biomass C (C R ) pools were collected during two sampling periods: (i) after the first ratooning in April 2010 and (ii) after the 3rd ratooning in August 2011. In April 2010, two cored (5-cm-diameter) samples for C R and two augered (6-cm-diameter) samples for C S were collected at two depths (0-15 and 15-30 cm) from each plot. For August 2011, four cored (5-cm-diameter) samples (two for C R and two for C S ) were collected at the same depth increments. Soil samples were air-dried at 25°C for a week until no mass change, sieved to 2 mm, and weighed, and elemental C and N concentration was determined on a subsample as described above. Measurements were limited to 30 cm depths due to the presence of a dense, clay-rich plow layer.
The SOC pool was quantified using the equivalent mass of soil method (Ellert & Bettany, 1995;Gifford & Roderick, 2003), which negates the issue of core compaction during sampling and the effect of land management on bulk density. Equivalent mass of soil (g C m À2 ) at target depth (t) was calculated for each plot as: where M S (Z a ) and M S (Z b ) are the masses of soil for the first and second increments at depth Z a and Z b , C S (Z a ) and C S (Z b ) are the masses of C of both depth increments, and M S (t) is the target mass of soil that all samples are compared against. This linear interpolation allows comparison of SOC pools without the need of an accurate core volume. The reference soil mass selected to make comparisons within this study was 300 kg m À2 because it was approximately the mean mass of soil samples collected to 30 cm depth across all cores. As such, the SOC pools were expressed as g C m À2 in 300 kg of soil m À2 . For C R , cored soil samples were split into four parts, inserted into 250 ml Nalgene polypropylene bottles and shaken for 16 h with 100 ml of 10% sodium hexametaphosphate. Dispersed soils were wet-sieved through a 0.5 mm sieve, and collected roots were washed with deionized water, ground to pass through a 250-lm sieve using a UDY cyclone mill (Tecator Inc., Boulder, CO, USA), and root C concentrations were determined as above. No separation of live and dead roots was conducted. The root biomass in soil cores was expressed in terms of C by multiplying root C concentration by mass of root biomass extracted from wet sieving (g C m À2 to 30 cm depth). Average values across cores of C S and C R were calculated for each plot.
The belowground autotrophic input of C was estimated in each plot using the TBCF mass balance approach for nonsteady state conditions (Giardina & Ryan, 2002). The mass balance approach quantifies the sum of all pertinent C inputs that plants send belowground including root production and respiration, root exudation, and C flow to symbionts, all of which are difficult to measure independently. Use of this method for quantifying autotrophic flux of C to belowground outside forest ecosystems has been limited (Adair et al., 2009;Ford et al., 2012), but as long as all important C fluxes (inputs and outputs) and changes in C pools are quantified accurately, the mass balance approach for estimating TBCF should be applicable in any environment and is calculated as: where TBCF is the total flux of C that plants send belowground, F S is soil CO 2 efflux (g C m À2 yr À1 ), F A is aboveground litterfall (g C m À2 yr À1 ), and C S , C R, and C L are soil, root, and litter layer C pools (g C m À2 ), respectively. This approach assumes there are negligible losses of C via erosion and leaching. Soil erosion in the study site was nonexistent due to a flat topography and controlled water input. Leaching loss of C was also assumed to be negligible with <1 m precipitation annually. In addition, leaching of C has been shown to account for <2% of F S in a variety of ecosystems including grasslands and croplands (Kindler et al., 2011). It was also assumed that changes in soil and root biomass C pools below 30 cm were negligible during the measurement period. Trenches adjacent to the study plots demonstrated that minimal root growth extended below 30 cm depth in this system. Annualized changes in C S (ΔC S ) and C R (ΔC R ) were calculated as the differences between the 2010 and 2011 annualized pools. Change in litter layer C was assumed to be zero from 2010 to 2011 because all aboveground biomass was removed regularly for feedstock conversion.

Root quality and decay
Grass stalks with roots attached were excavated from soil pits (30 cm 9 30 cm, 15 cm in depth) in each plot, washed with deionized water, and air-dried for 1 week. Approximately 20 g of subsampled root tissue was analyzed for neutral detergent fiber (NDF), acid detergent fiber (ADF), cellulose, and acid unhydrolyzable compounds (lignin and other recalcitrant compounds hereafter referred to as lignin) using standard reagents (Van Soest, 1963) on a Fiber Analyser (Ankom, Macedon, NY, USA) at the Agricultural Diagnostic Service Center of the University of Hawaii Manoa. The quantities of nonfiber carbohydrates (NFC) such as organic acids, sugars, and starch were determined as 100 -NDF, and hemicellulose was calculated as NDF -ADF using results of the fiber analysis (Hall, 2003). The lignin value estimated using the sequential digestion method is a coarse estimate of lignin as polyphenolic and other unsaturated substances such as tannins and suberin may be included in the value (Van Soest & Wine, 1968). Carbon and N concentrations were determined on subsamples as described above. A common area root decomposition experiment was conducted from December 2010 to August 2011 to determine the root decomposition constant (k) using the litterbag method (Ostertag, 2001). For each replicate of each accession (four replicates 9 eight accessions = 32 experimental units), a set of five bags was made with 0.5 g subsamples of air-dried root materials 2-3 cm in length. Each bag was constructed of 5 cm 9 5 cm 0.132 mm nylon mesh. Bags were buried in a 1.5 9 7 m area adjacent to each accession plot, in a line (to facilitate future destructive sampling events) in random order such that 5-cm borders separated 32 lines of bag sets. Individual bags were buried at a 45°angle to the surface and inserted 3.5 cm into the soil. Bags from each replicate set were collected randomly at 1, 2, 3, 5, and 8 months after deployment. The site was covered with a weed mat and irrigated three times a week to simulate conditions similar to the accession plots. Collected bags were rinsed with deionized water and dried to a constant mass at 75°C, and residue C and N were determined following the procedure described above. For each experimental unit, a firstorder single-pool exponential function was applied for the root C decline over time to characterize decomposition (Wider & Lang, 1982): where L t represents the proportion of original mass at time t (year), and k represents the decay constant (t À1 ). Both firstorder single-and double-pool exponential functions were considered, but a single-pool exponential function was best suited for the data based on r 2 values (mean r 2 was 0.90 with a range of 0.69-0.99) and visual analysis of residuals.

Statistical analysis
Comparison of C pools and fluxes, root composition, and decay constant. In the statistical comparisons of C pools and fluxes (H1) and root elemental/chemical composition and decay constant (H2), three levels of comparisons were made. First, an accession effect was determined by treating each of the eight accessions independently. Analysis of variance (ANOVA) was performed to compare the effect of accessions on ecosystem C pools, fluxes, and input quality indicators using the software R 2.15 (R Development Core Team, 2012). Prior to analysis, homogeneity of variances for block and accession effects was tested with Levene's test on all response variables. Graphical assessment of normality and outliers revealed no severe nonnormality or outliers in any of the variables. A wrong accession was planted in one of the replicates of GG2 and therefore was omitted from the data analysis, giving GG2 only three replicates (Table 1). The unorthogonality created by the missing observation on sum of squares of accession effect was negligible, and therefore, type II sum of squares was used to calculate F values. When the effect of accession on response variables was significant at P ≤ 0.05, Tukey post hoc comparisons were performed to determine significant differences between accessions.
Next, two preplanned custom contrasts comparing (i) guinea grass vs. napier grass + PM9D and (ii) napier grass vs. PM9D were tested using a multcomp package in R (Hothorn et al., 2008) to determine species effect (guinea grass vs. napier grass hybrid + PM9D). Further, because PM9D is a napier grass hybrid with pearl millet, a variety effect was also determined (napier grass vs. PM9D). The two custom contrasts are independent of each other, negating the risk of inflating type I experiment-wise error rate (Seltman, 2012).
A correlation analysis was used to identify factors controlling the decay constant (k) (H2). Pearson correlation coefficients (r) were calculated (n = 32) for compositional factors (C, N, NFC, hemicellulose, cellulose, lignin, C : N ratio, and lignin : N ratio) and the root decay constant (k) calculated from the litterbag experiment.
Functional belowground carbon dynamics. Structural equation modeling (SEM), a multivariate statistical method, was used to disentangle the effect of the multiple explanatory variables into hypothesized causal pathways (Grace, 2006;Tabachnick & Fidell, 2007) (H3). Several recent ecological studies have used SEM to investigate the causal relationships among multicollinear predicting variables and their effect on soil CO 2 efflux (F S ) (Geng et al., 2012;Matias et al., 2012), the SOC pool (Jonsson & Wardle, 2010;Brahim et al., 2011), and soil microbial community composition (Eisenhauer et al., 2012). The overall purpose of SEM in this study was to delineate patterns of direct and indirect effects of explanatory variables by formulating a causal model based on both results of this study and prior knowledge of SOC accumulation. Structural equation modeling was performed using the 'sem' package in R. Variables are connected with one-way arrows indicating the flow of causal relationships. Regression coefficients are then parameterized simultaneously using a maximum-likelihood method for each arrow (Grace, 2006). The variance and covariance matrices from the parameterized coefficients were tested against the matrix from the data to determine overall fit of the hypothesized models to the data. Unstandardized and standardized path coefficients and associated standard error, Z-score, and P-value are reported for each parameter estimate.
All plots were treated as independent, and thus, all 32 plots were included in this analysis. Given the small sample size (n = 32), only observed variables and their relationships were considered; no latent (i.e., unobserved, underlying) variables were included in the model. Also, variables were limited to six to minimize estimated parameters (Grace, 2006). The six variables chosen were root lignin concentration, aboveground biomass yield, k, F S , C R 2010, and C S 2011. Because SEM relies on variance and covariance matrices of the variables, assessment of univariate and multivariate distribution of variables for outliers, linearity, and normality is crucial for subsequent inferences (Tabachnick & Fidell, 2007). Bivariate scatterplots and residual plots were used to assess linearity and normality and to check for outliers using Minitab 16 (Minitab Inc., State College, PA, USA).
The acceptability of the final model was first determined by chi-square (v 2 ) tests (P > 0.05). Nonsignificant chi-square test indicates that the variance and covariance matrices of the ) and fluxes (g C m À2 yr À1 ) of soil carbon for grass accessions (n = 4 except for GG2 where n = 3). C s and C R were expressed in terms of equivalent mass of 300 g m À2 , which was approximately equal to 0-30 cm depth. ANOVA  yr À1 ); ΔC R , difference in root C pools from 2010 to 2011 sampling dates (g C m À2 yr À1 ); C S , soil C pool (g C m À2 ); C R , root C pool (g C m À2 ); TBCF, total belowground carbon flux (g C m À2 yr À1  (Hu & Bentler, 1999). The modification indices >5, which suggests that removing a variable would improve the model fit substantially, were considered when the paths had theoretical support and estimated parameters could be interpreted (J€ oreskog & Sorbom, 1989;Grace, 2006). Additionally, a test of indirect effects through intermediate variables, called mediation, was carried out using the delta method (Sobel, 1982).

Ecosystem C pools and fluxes
In 2011, SOC pools (mean of 5434 AE 51 g C m À2 ) were, on average, 5% higher than in 2010 (mean of 5160 AE 70 g C m À2 ) (Table 1), resulting in positive ΔC S for all grass accessions (Fig. 1a). Mean SOC ranged from 4908 to 5617 g C m À2 in 2010; PM9D (5394 AE 257 g C m À2 ) was significantly greater than napier grass (5021 AE 73 g C m À2 ); but neither napier grass nor PM9D differed significantly from guinea grass (5191 AE 114 g C m À2 ). No significant differences were observed across accessions, species, or variety within napier grass in 2011. The increase in SOC pool from 2010 to 2011 (ΔC S ) was greater for napier grass (249 AE 45 g C m À2 ) compared to guinea grass (160 AE 53 g C m À2 ) and for napier grass (without PM9D) (285 AE 43 g C m À2 ) compared to PM9D (67 AE 69 g C m À2 ) ( Table 1). In 2010, root C pools (C R ) ranged from 21 to 50 g C m À2 with an average of 42 AE 2.5 g C m À2 across accessions and species, and no significant differences among accessions (Table 1). In 2011, average C R increased to 123 AE 0.8 g C m À2 across accessions, again with no significant differences among accessions. Accrual of root biomass (ΔC R ) from 2010 to 2011 averaged 61 AE 5.5 g C m À2 yr À1 across accessions and was highly variable across replicates within accessions (coefficient of variations as high as 120%), which resulted in no significant differences across accessions (Table 1).
Average aboveground biomass C ranged from 931 to 1805 g C m À2 yr À1 across accessions (Table 1) and was significantly greater for napier grasses and PM9D (1419 AE 136 g C m À2 yr À1 ) than guinea grasses (1047 AE 87 g C m À2 yr À1 ) (Fig. 1b). Aboveground production values (931-1805 g C m À2 yr À1 ) were consistent with other reported values and yields (1020-1347 g C m À2 yr À1 ) in the literature for perennial grasses in temperate and tropical regions assuming that 48% of reported yields was C (Bremer et al., 1998;Tufekcioglu et al., 2001;Dornbush & Raich, 2006). Yields of irrigated and rainfed napier grass in Hawaii were estimated to be 2352 and 1902 g C m À2 yr À1 , respectively (Kinoshita et al., 1995). Although no trials have been carried out on guinea grass in Hawai'i prior to this study, yield of guinea grass reported in Mexico reaches up to 1152 g C m À2 yr À1 (Reynoso et al., 2009).
Annualized litterfall (F A ) ranged from 110 to 229 g C m À2 yr À1 (Table 1) and was also greater for napier grasses and PM9D (128 AE 6 g C m À2 yr À1 ) than for guinea grass accessions (159 AE 20 g C m À2 yr À1 ). Within the napier grasses and PM9D, napier grass F A was greater than PM9D. Litterfall C values were 1.6-3.3 times greater than root C values, but nearly 10 times less than the aboveground C that was removed as biomass for fuel production. In addition, although some dissolved C from litter may have passed, weed mats prevented the incorporation of aboveground litter into the soil profile. Litterfall therefore was unlikely to be a substantial contributor to SOC pools.
Annualized soil CO 2 efflux (F S ) ranged from 1325 to 1788 g C m À2 yr À1 and averaged 1570 AE 58 g C m À2 White vertical bars indicate napier grass accessions, gray striped bar indicates PM9D, and gray vertical bars indicate guinea accessions. Unconnected horizontal bars indicate significant differences at P ≤ 0.05; species (guinea grass vs. napier grass) and/or variety (Napier vs. PM9D within napier grass) differences are reported as described in the Methods section. yr À1 with no significant differences across accessions (Table 1). Significant differences in ΔC S and F A canceled each other out, resulting in no significant differences in TBCF among accessions. The mean TBCF was 1614 AE 62 g C m À2 yr À1 , which was positively correlated with aboveground biomass C across accessions (r = 0.433, P = 0.015; n = 31). The range of TBCF was between 1510 and 1917 g C m À2 yr À1 , which was nearly twice as much as values reported in a nutrient poor sandy soil grassland in Minnesota (Adair et al., 2009) and in line with values for a highly productive Eucalyptus saligna plantation in Hawai'i (1400-1900 g C m À2 yr À1 ; Ryan et al., 2004).

Root quality and decay
Multiple aspects of root chemical characteristics varied significantly among accessions and by species ( Table 2). The mean concentration of root NFC, which contains soluble sugars and amino acids, varied from 21.3% to 26.8% across accessions, and root NFC concentrations were significantly greater for NG1 and NG3 than GG3 (P = 0.008). Averaged over the different accessions, napier grasses and PM9D (25.9 AE 0.8%) had significantly greater root NFC concentrations than guinea grass (22.7 AE 0.7%) (P ≤ 0.001). Mean root cellulose and hemicellulose concentrations for all accession were 33.4 AE 0.2% and 22.6 AE 0.3%, respectively, with no accession or species effects.
Root lignin concentration ranged from 17.9% to 20.9% across accessions, with GG2, GG3, and GG4 significantly greater than NG1, NG2, and PM9D ( Table 2). Comparisons of species revealed a significant and distinct pattern of root lignin concentrations in the order of PM9D (17.9 AE 0.7%) < napier grass (19.0 AE 0.4%) < guinea grass (20.8 AE 0.4%). Root C concentrations also differed significantly across accessions (Table 2), guinea grass GG1 was greater than any napier grass accessions and PM9D, and guinea grass accessions on average had significantly greater root C (44.6 AE 0.3%) compared to the average of napier grass accessions and PM9D (42.6 AE 0.2%). Root N concentrations did not vary across accessions (Table 2). However, napier grass and PM9D (0.72 AE 0.04%) had significantly greater root N concentrations than guinea grass (0.64 AE 0.04).
Root C : N ranged widely from 57 to 82 with significant differences across accessions (Table 2). On average, napier grass accessions and PM9D together had significantly lower C : N (62 AE 4) compared to guinea grass accessions (74 AE 6). Root lignin : N ratios were less variable at 24-39, but still varied significantly across accessions. GG1 had significantly greater lignin : N compared to PM9D. Similar to C : N, lignin : N was significantly lower for napier and PM9D compared to guinea grass accessions. Root chemical characteristics such as root N and lignin concentrations, and lignin : N ratio observed in this study were similar to values reported for temperate C4 grasses where % N ranged from 0.76 to 1.19, % lignin ranged from 9.3 to 29.3, and lignin : N ratio was between 11.0 and 25.8 (Vivanco & Austin, 2006).
The decay constant, k, ranged from 0.95 to 1.69, similar to the range observed for temperate grass root decomposition (0.51-1.82 yr À1 ; Vivanco & Austin, 2006). On average, roots of napier grass accessions and PxD decomposed significantly faster (1.59 AE 0.09) compared to guinea grass accessions (1.19 AE 0.09), while napier grasses and PxD showed a similar decay constant. Pearson correlation analysis determined that k was significantly correlated with all root components except cellulose and hemicellulose (Table 3).
High negative and significant correlations with k were observed in root C and lignin concentration, C : N, and lignin : N. The k value was also significantly and positively correlated with root N and NFC concentrations. Root lignin concentration was the single best predictor of k, and adding N as a second predictor did not significantly improve the fit, which is consistent with a recent litter decomposition review that manipulating a single chemical characteristic would not contribute to drastic differences in decomposition rates (Prescott, 2010).

Functional belowground C dynamics
In the initial model structure developed around our hypotheses regarding the primary drivers of SOC accumulation under cultivation, aboveground biomass yield and root lignin were assumed independent of other variables and to have a direct effect on C S 2011 and C R 2010. Fs and k served as both explanatory and response variables in the initial model, and C S 2011 did not affect any other variables (Fig. 2a). Ad hoc model modification was not considered because the value of the highest modification index was <5, suggesting that the hypothesized initial model fit the data well enough to not require major changes (J€ oreskog & Sorbom, 1989). In the final model, nonsignificant paths from lignin and yield to C S 2011 were removed because the low probability associated with these paths suggested limited chance that such an effect was present in the study (Table 4). Hypothesized paths from k to C R 2010, C R 2010 to F S , and C R 2010 to C S 2011 were retained in the model because such effects are likely to be real. As a result of these adjustments to the conceptual model, the final model fit the data well with v 2 = 3.15, df = 6, P = 0.79, RMSEA = 0, TLI = 1.18, and SRMR = 0.05 (Fig. 2b). ).
In the final model, root lignin concentration had a highly significant, negative effect on k, as expected from correlation analysis of root chemical compositions ( Table 4). The standardized path coefficient of À0.680 in the lignin to k path implies that as lignin decreases by one standard deviation, it is predicted that k increases by 68% of its standard deviation (Table 4; Fig. 2b). Expressed in terms of an unstandardized coefficient, a one unit decrease in lignin concentration would result in an increase of 0.159 in k. Expressed in terms of relevant ranges, as lignin concentrations increase its effective range from 15.88 to 22.53, k was predicted to decrease 52.8% of its effective range from 0.565 to 2.569 ( Fig. 2b). Almost half (46%) of the variation in k was accounted for by root lignin concentrations.
Aboveground yield had significant positive effects on C R 2010 and F S (Table 4), explaining 24% and 23% of variation in each variable, respectively (Fig. 2b). Aboveground yield did not have a direct effect on C S 2011, and therefore whatever effect the yield had on C S 2011 was mediated through F S and C R 2010. Soil CO 2 efflux, a pathway for loss of SOC, had a marginal negative effect on C S 2011 (Fig. 2b).
Contrary to expectation, C R 2010 was a poor predictor of both F S and C S 2011. Due to the weak effects of C R 2010 on F S or C S 2011, the indirect effect of yield on C S 2011 through F S and C R 2010 also were not significant (Table 4). Root k was the only variable in the model with a significant direct effect on C S 2011 (Fig. 2b). Together with F S , k explained 40% of the variation in C S 2011. However, k had a positive effect on C S 2011, opposite what was originally hypothesized. Moreover, a direct effect of lignin on C S 2011 was not observed in the model, while the indirect negative effect of lignin on C S 2011 through k was significant (Table 4). Therefore, the effect of lignin on C S 2011 was negatively mediated through the intervening variable k indicating that greater root lignin concentration would result in less SOC, contrary to the initial expectation. Furthermore, root k had a marginally negative effect on F S (Table 4). Contrary to initial expectations, root lignin indirectly affected F S positively at a marginal level (Table 4). The Table 3 Pearson correlation coefficients (r) and significance (P) for each chemical composition factor and the decay constant (k) from litterbag experiment (n = 31) Root composition factor r P Bolded factors indicate the largest r for positive and negative correlations. final model indicated that the indirect effect of belowground input quality, specifically, low lignin concentration, was a more important determinant of the SOC pool size than the indirect effect of aboveground biomass yield.

Ecosystem C pools and fluxes
Cultivation of fallow soil typically results in rapid decreases in SOC pools (Davidson & Ackerman, 1993), particularly if long-term, intensive tillage and planting are practiced (Dalal & Mayer, 1986). In contrast, when conservation practices such as zero-tillage or ratoon harvest are used with crops known to have extensive belowground root systems such as perennial grasses, C loss can be minimized and even accumulation of SOC is possible (Lal & Kimble, 1997;Franzluebbers, 2012;Anderson-Teixeira et al., 2009;Powlson et al., 2011;Anderson-Teixeira et al., 2013). In this tropical study, an accumulation of SOC occurred in the initial years of cultivation of a fallow field for perennial C4 grasses for biofuel production. Napier grass varieties, in particular, accumulated greater SOC during a 1-year period than guinea grass varieties, suggesting that cultivation of napier grass instead of guinea grass may have a greater benefit to the greenhouse warming potential of the production system in the form of C sequestration in soil.
In 2010, the SOC pool in adjacent fallow soils was 5423 g C m À2 (Sumiyoshi, 2011) compared to 5174 and 5442 g C m À2 for the accession plots for 2010 and 2011, respectively. Thus, it appears that any initial decrease in SOC after preparation of the field recovered under the cultivation of ratoon-harvested perennial grasses within 2 years and, in the interim, this dataset captured a dynamic window of postcultivation soil carbon dynamics. Root biomass builds rapidly after the first planting and then stabilizes following the first ratoon harvest as it comes into equilibrium with root mortality during Unstd. coeff., unstandardized path coefficient; Std. coeff., standardized path coefficient; Std relevant, standardized relevant range; Lignin, root lignin concentration (%); k, root decay constant (yr À1 ); F S , annualized soil CO 2 efflux (g C m À2 yr À1 ); C S , soil C pool (g C m À2 ); C R , root C pool (g C m À2 ). Significant effects (P < 0.05) are given in bold. Standard errors (SE) of indirect effects were estimated using the delta method (see Methods). subsequent harvests. This suggests that tilling and replanting a field that has been in zero-tillage cultivation may not lose accumulated C, as is often assumed in simulation models and life cycle analyses. Soil CO 2 efflux and TBCF were substantial, and annual C fluxes were greater than aboveground production, supporting the hypothesis that a wide range of ecosystem C pools and fluxes would contribute collectively to the soil C balance of our study systems. Consistent with the observed pattern between TBCF and aboveground net primary productivity in forest ecosystems (Litton et al., 2007), TBCF was positively related to aboveground yield in this study. This is the first study, to our knowledge, to use this mass balance approach to estimate the autotrophic flux of C to belowground in a tropical agricultural system, and suggests that coupling of aboveground and belowground C fluxes occurs universally across terrestrial ecosystems.
That total SOC pool size did not differ among the grasses was not surprising given the large residual stock of C in relation to the comparatively small change in SOC pool (<5% of total SOC) over 1 year of cultivation. Differences in the component pools and fluxes among the grasses, however, suggest that over time SOC pools may change substantially in response to cultivation of different grasses. For example, aboveground biomass yield and the annual increase in soil C differed among the accessions and more generally by species. If highyielding grasses (such as two of the tested napier grass accessions and the napier grass hybrid) also have high root biomass, and these inputs result in a net gain of soil C over time, then selecting grasses with high aboveground biomass yield may have a net greenhouse warming benefit to the production system.

Root quality and decay
The link between high-yielding grasses and soil C accumulation is dependent in part on the fate of the root biomass as an input to the belowground soil system. The fate (i.e., whether C is released from the system or retained), in turn, is in part dependent on the chemical composition of detrital inputs (Cotrufo et al., 2015). Here, the hypothesis that the selected perennial grasses would have a wide range of root input quality indicators was supported. The concentration of carbohydrates, lignin, C, N, and their relative abundance varied among grasses, and often by species. Specifically, many indicators of low C quality, such as high lignin concentration, C : N ratio, and lignin : N ratio, were greatest for the guinea grasses compared to napier grasses and the napier grass hybrid. A view emerges of complex interactions among the potential factors controlling soil C accumulation when considering variation within these chemical indicators of input quality in combination with variation described above in ecosystem C pools and fluxes.
Root decay was strongly related to root lignin concentration, which is consistent with a recent review by Prescott (2010) that identified lignin as the most reliable predictor of root decomposition rate. Mass loss during the litterbag study, and the resultant decay constant, was interpreted to represent losses of potential C inputs from the soil system through microbial respiration or leaching. Accordingly, we expected root lignin concentration to have a direct, positive effect on SOC pool size or an indirect effect via the pathway of high lignin concentration driving low root decay (and therefore low losses of C outputs via respiration and leaching) resulting in the accumulation of root fragments or residues as SOC.

Functional belowground C dynamics
Linking ecosystem C pools and fluxes to accurately represent soil function remains a challenge, yet is critical for predicting the effect of land-use change and management on the net C balance of agroecosystems. In a renewable energy system, environmental sustainability is often assessed from the perspective of the net global warming potential for the collective feedstock production system, fuel conversion, and delivery in a life cycle analysis. Understanding functional belowground C dynamics and linking soil C accumulation to measureable indices can be important in guiding feedstock and management choices during the planning and implementation of a production system where climate change mitigation is a desirable outcome. Based on the results from H1 and H2, we developed a conceptual model for the factors controlling soil C accumulation in the study system. However, hypothesis H3, that aboveground biomass and root lignin concentration (both easily measured factors) provide direct positive control on SOC pool size, was not supported. Further, the roles of root C, soil CO 2 efflux, and decay constant as drivers of the indirect effects of yield and root lignin on SOC pool size were confirmed, but also were not as expected.
The marginal and negative correlation between root decay rate (k) and soil CO 2 efflux (F S ) was unexpected. However, F S is the cumulative flux of CO 2 from soils to the overlying atmosphere that is comprised of both autotrophic and heterotrophic components, and k is only one component of heterotrophic respiration contributing to F S in addition to decomposition of aboveground litter and SOM. Prior studies have shown that F S is largely driven by recent photosynthetic C input to belowground (see Kuzyakov & Gavrichkova, 2010) and that the root respiration contribution to F S is~30-50% and increases asymptotically with increasing F S (Bond-Lamberty et al., 2004). As such, k is likely to be a small component of overall F S (~10-15%) and it is entirely feasible that k could be marginally and negatively correlated with F S . In this study, dissolved organic matter transport and fragmentation contributed to k, in addition to biochemical oxidation.
In the final model describing the functional belowground C dynamics of this system, high aboveground biomass yields did result in high root biomass C, but did not directly result in high SOC pools. Instead of high root biomass C resulting in high SOC, the chemistry of roots, primarily lignin concentration, was the critical factor in controlling root decomposition and subsequent SOC accumulation. Specifically, low root lignin concentration drove high root decomposition rates, as hypothesized. High decomposition rates, however, resulted not in greater soil CO 2 losses as expected, but rather in greater SOC accumulation. Thus, the hypothesis that grasses with higher aboveground productivity and higher lignin concentration (and lower decomposition rate) would accumulate the greatest amount of SOC over the study period was not supported. Instead, biofuel grasses cultivated with zero-tillage that had high productivity and roots with low lignin (i.e., high quality) and high decomposition rates accumulated the most SOC in this study. In temperate switchgrass systems, Ontl et al. (2015) reported that crops that maximize root productivity would lead to the largest gains in protected soil C. On a longer time frame, as roots die and decompose with repeated ratooning, it seems likely that accessions with high root biomass, high proportion of root death upon ratoon, low total lignin concentration, and high proportion of readily degradable lignin monomers would accumulate the greatest amount of SOC.
Soil OC formation occurs via both a rapid dissolved organic matter-microbial path and a slower physical transfer path as litter fragments enter the soil (Cotrufo et al., 2013(Cotrufo et al., , 2015. A likely mechanism for SOC accumulation may be through incorporation of root fragments and microbial biomass C, necromass, and by-products into aggregates in the absence of tillage (Six et al., 2000;Ontl et al., 2015;Tiemann & Grandy, 2015). Huang et al. (2011) showed in a forest environment that inputs with lower lignin concentration resulted in more SOC accumulation due to greater inputs through fragments, similar to results from this study. Once within aggregates, fragmented tissue, soluble organic matter, and microbial biomass by-products are physically protected from further decomposition. Taken together, the results from this study agree with the recent understanding that microbial transformation of plant residue and physical protection within aggregates can be a more important driver of SOC stabilization than chemical recalcitrance per se (Schmidt et al., 2011;Lehmann & Kleber, 2015).
The effect of biofuel biomass cultivation on SOC may be positive or negative depending on the type of cropping system. Perennial crops and low-or zero-tillage systems with little residue removal are the most likely to result in a net increase or maintenance of SOC (Lewandowski, 2013). Much of the research on C sequestration in bioenergy crop cultivation has focused on temperate and subtropical regions (Gelfand et al., 2013;Agostini et al., 2015), but increasingly tropical perennial grasses such as sugarcane and napier grass are being utilized (Lemus & Lal, 2005;Anderson-Teixeira et al., 2009). Recent work completed on Miscanthus (switchgrass) cultivation in Europe showed that the accumulation of SOC was dependent on the genotypespecific root distribution of each variety , further highlighting the importance of understanding the underlying mechanisms driving SOC dynamics in a given system. On similar systems in the US Midwest, soil type affected the efficacy of aggregates to form and protect C derived from giant Miscanthus inputs under cultivation for biofuel (Tiemann & Grandy, 2015). Identifying characteristics of plants (such as biomass production, rooting morphology, and tissue chemistry) and of soil (such as texture, mineralogy, and organic matter dynamics) and how they will interact under cultivation is critical to predicting the capacity of a production system to sequester C in soil.
The most beneficial, sustainable biofuel systems should focus on feedstocks that complement food production, do not result in direct or indirect losses of C through land-use change, and do not cause net greenhouse gas emissions (Tilman et al., 2009;Anderson-Teixeira et al., 2012). Productivity and sustainability of systems cultivated for biofuel feedstocks often are positively associated with SOC content (Lehmann, 2012;Anderson-Teixeira et al., 2013), and increasing SOC in agricultural systems simultaneously restores soil quality through improvements in the biological, physical, and chemical functions provided by organic matter (Powlson et al., 2011;Franzluebbers, 2012;Murphy, 2014) while sequestering atmospheric CO 2 (Lal, 2013;Werling et al., 2014;DeLucia, 2016). In this study of tropical perennial C4 grasses, high-yielding varieties with low root lignin cultivated in zero-tillage systems were the most promising biomass feedstock to meet the complex needs of society for renewable fuels, food, and climate change mitigation through net carbon sequestration during feedstock production. Genetic improvements are possible through technology and breeding to maximize climate change mitigation of renewable fuels systems, and these results reveal plant traits that may be the best predictors of soil carbon accumulation in the tropics. Identifying soil management practices that ensure productivity and profitability for farmers while building SOC remains a challenge (Hill et al., 2006;Meki et al., 2014), but isolating readily measured indicators of potential SOC accumulation such as yield, root morphology, and tissue chemistry improves our ability to make system-level decisions that ensure both environmental and economic sustainability.