Data analyses were performed using RStudio version 3.4.0 (R Development Core Team 2017) statistical software. Linear regression analysis was used to examine changes in soil chemical and physical properties (i.e., major cations and BD) relative to increase time after pasture conversion into OP plantations. Assumption of normality was checked for all analyses with Shapiro-Wilk test and visual inspection of normality plots. If assumption of normality was not satisfied, then permutation tests were performed as in the case of Ca, sum of cations, and BS. Patterns in changes of soil C stocks over the OP chronosequence time were examined for the bulk soil and for each of the four sampled soil layers using regression models. The nonlinear least square “nls” and linear model “lm” functions in R were used to fit nonlinear and linear regression models, respectively. In addition, the “segmented” function was used to perform segmented (broken line) regression analysis. Statistical significance was declared at P < 0.05. Similarly, model fit for changes in C3-OP–derived C and C4 pasture–derived C over the chronosequence time was also examined for each of the four studied soil layers by testing the above-mentioned models. These regression analyses allowed estimating rates of C3 accumulation and C4 decomposition and decay constants (k), rates of total soil C decrease, and break points in soil C stock changes (time at which a change in the direction of change in C stocks occurred).

After testing various models (i.e., mono-exponential and bi-exponential), model performance assessment was based on Akaike information criterion (AIC) and coefficient of determination (R2) values. Models with the largest R2 and lowest AIC values were selected (see table S1). Accordingly, (i) changes in whole soil C stocks and C stocks in each soil layer were described by fitting segmented regression models (except linear regression on the deepest soil layer of 30 to 50 cm), which yielded estimated break points; (ii) C3-derived C in the 10- to 20-cm, 20- to 30-cm, and 30- to 50-cm soil layers were described by linear models, while pattern in C3-derived C in the soil surface layer (0 to 10 cm) was described by an exponential rise to equilibrium model of the formy=((k*y0A)*exp(k*t+A)/k(4)where y is the C stock, k is the annual decay constant of the pool, y0 is the C3-C stock before OP cultivation started (thus, 0), A is the C3 annual input to the C3 pool, and t is time after conversion; (iii) A first-order decay model was fitted to the obtained pasture-derived C data in the four studied soil layersy=c*exp(k*t)(5)where y is the C stock, c is the original C stock before OP cultivation, k is the decay rate constant, and t referred to time. The half-life (HL) of the original C stock in Eq. 5 and of the OP input in Eq. 4 was calculated asHL=ln(2)/k(6)

A PCA was carried out for further exploration of the relations between soil parameters that can be affected by cultivation (BD, C:N, 13C, C, N, EA, Na, pH, K, Mg, 15N, Ca, BS, and P).