Introduction

The soil reservoir holds more carbon (C) than all terrestrial biomass and the atmosphere combined but is increasingly vulnerable to warming-induced losses in carbon stocks, which would be a positive feedback to climate change (Field and Raupach 2004; Davidson and Janssens 2006). The IPCC “business-as-usual-model” model (RCP 8.5) predicts global average air and soil temperatures will increase 4 °C by 2100. While more than 80% of the world's soil organic carbon (SOC) exists below 20 cm (Jobbágy and Jackson 2000), depth will not protect it from warming—CMIP5 projections indicate the same amount of warming throughout the top meter of soil (Pachauri et al. 2015; Soong et al. 2020). In response, Torn et al. (2015) recently called for an international network of deep soil manipulation experiments (iSEN) to understand the response of soil carbon throughout the soil profile to climate change. One of the first experiments warming the whole soil profile indicated deep soils have the potential to become a large positive feedback to climate change as soil respiration throughout the soil profile increased by 30–34% and subsoil carbon storage decreased by 33% in response to 4 °C of warming (Hicks Pries et al. 2017; Soong et al. 2021). Evidence from temperate and boreal soils show increasing soil temperature escalates soil microbial activity and, consequently, soil respiration rates throughout the soil profile, causing declining soil C stocks where decomposition losses outweigh input gains (Hanson et al. 2017, 2020; Soong et al. 2021).

Tropical soils hold one third of the world's SOC (Jobbágy and Jackson 2000), but few experiments warm tropical soils (Nottingham et al. 2019, 2020; Kimball et al. 2018). Historically, soil warming experiments focused on temperate and boreal systems (Plaza et al. 2019; Schädel et al. 2018). Lower latitude soils were hypothesized to be less vulnerable to rising temperatures due to the smaller absolute magnitude of warming they will experience and the increased mineral protection of organic matter in highly weathered tropical soils (Koven et al. 2017; Haaf et al. 2021). Initial predictions also hypothesized that tropical systems experiencing warm mean annual temperatures (MAT) under isothermic conditions would have little microbial response to warming as microbes might be approaching their biochemical thermal optima (Wood et al. 2012). However, model projections indicate unprecedented increases in tropical surface temperatures (Diffenbaugh and Scherer 2011), and soil incubations have found no evidence of biochemical thermal optima occurring at these projected temperatures (Wood et al. 2012). Cavaleri et al. (2015) thus highlighted the urgent need for tropical deep soil warming experiments, given the high uncertainty of their SOC response to rising soil temperatures. Limited experimental evidence now shows that tropical soils do behave more like their temperate and boreal counterparts (Hicks Pries et al. 2017; Hanson et al. 2020), respiring significantly more CO2 in response to soil warming (Zimmermann et al. 2010; Nottingham et al. 2020). The temperature sensitivity of soil respiration could lead to decreased soil carbon storage in tropical soils as observed along elevational gradients wherein carbon storage is greatest at higher, cooler elevations (Nottingham et al. 2015). In fact, one study suggests that tropical forests may be more sensitive to increases in temperature than higher latitudes, becoming large C sources (Mau et al. 2018).

Tropical soils store a larger proportion of their carbon in mineral-associations than temperate soils (Heckman et al. 2022). Recent meta-analyses (Heckman et al. 2022; Rocci et al. 2021) and experiments (Soong et al. 2021) suggest that mineral-associated carbon is less vulnerable to warming-induced losses than particulate carbon because sorption to soil minerals may offer some protection from microbial enzymatic attack. One critical end member soil lacking within soil warming experiments are volcanic-ash derived Andisols (WRB classification: Andosol), which are known for their strong mineral protection of soil organic matter (Torn et al. 1997). Andisols hold a disproportionate amount of soil carbon, storing 5% of the world's soil carbon despite making up less than 1% of the world's soils by area (Matus et al. 2014). Organo-mineral complexes are abundant in andic soils (Post and Kwon 2000; Torn et al. 1997). One recent study across a MAT gradient of Hawaiian Andisols suggested that Al and Fe poorly- and non-crystalline minerals (PNCM) have strong soil C-protection mechanisms leading to large soil C stocks (Giardina et al. 2014). However, whether the high degree of protection by PNCM will limit the response of soil C in Andisols to the intense warming predicted by 2100 is unknown.

Here we present one of the first deep soil warming experiments in a tropical soil and the first in an Andisol. In contrast to most temperate and boreal systems, we hypothesized that the response of soil respiration to warming would saturate, indicating a threshold warming response controlled by poorly and non-crystalline minerals (PNCM), especially at depth. We tested this hypothesis in a wet montane tropical forest cleared of most vegetation using regression experimental design in which 50 locations across a hillslope were randomly subjected to warming along a gradient from ambient to + 4 °C and soil respiration was quantified weekly within five depth intervals (0–20, 20–40, 40–60, 60–80, and 80–100 cm).

Materials and methods

Field site description

The field site was located at the Lyon Arboretum, Honolulu, HI (21.3330° N, 157.8015° W), which has a mean annual temperature of 25.6 °C within an isothermic soil temperature regime and 4200 mm of precipitation annually (Fig. S1). The arboretum is a forest ecosystem focusing on plants native to the Pacific Islands/Asian region, but over a century ago, the land was used for cattle grazing. We selected the field site for its soil type (Medial over pumiceous or cindery, ferrihydritic, isothermic Typic Hapludands) and access to electricity. In June 2017, we manually cleared the site of invasive Hedychium flavescens understory vegetation for ease of access. Some shallow-rooted grasses (Eremochloa ophiuroides, Stenotaphrum secundatum; max rooting depth of 15 cm) were kept to mitigate soil erosion, and the site was cleared of overstory trees. We realize that removing vegetation caused a disturbance. Thus, we waited a year before beginning the warming experiment (Fig. S3) and mapped soil respiration across the hillslope to avoid any potential decomposition hotspots created by the disturbance (Fig. S1). We maintained the site by cutting any Hedychium flavescens at the base until regrowth ceased and thinning the shallow-rooted grasses to prevent them from overtaking the field site.

Experimental design

We measured surface soil respiration rates across the hillslope for 1 year before initiating the warming treatment to inform the regression experimental design. Sixty-four static chambers in a 2 × 2 m grid were sampled biweekly and interpolated into heatmaps using MATLAB to identify “hotspots” of microbial activity. Five blocks of equal surface area and initial surface respiration rates were placed across the hillslope to create a randomized block design. Areas of existing “hotspots” were not included in the blocking scheme to ensure a relatively homogenous initial surface soil respiration flux before starting the heating treatment (Fig. S1). Each block contained a single 2.5 m deep heating rod designed to create a gradient of heating across each block. Two sampling locations for each sampling depth interval (0–20, 20–40, 40–60, 60–80, and 80–100 cm) were randomly assigned to each block for measuring soil temperature, soil moisture, and soil respiration. Because this gradient design did not evenly heat across space, each sampling point experienced diffusion of CO2 from its warmer side that was likely balanced by diffusion of CO2 towards its cooler side. At each gas well set, soil temperature was recorded by the novel soil temperature sensors at the depth of interest (bottom of depth interval). The temperature sensors consisted of a single high accuracy digital sensor mounted on a board (TSYS01 TE Connectivity, Switzerland, custom made by UH Manoa MESH lab) and enclosed in a 3 ml polystyrene test tubes filled with epoxy that reported soil temperature every 30 s to an online server. Soil moisture was measured using a Sentek Diviner 2000 1.6 m rod (Sentek, Australia) at 30 static locations across the hillslope due to cost limitations. We assigned each gas well to its closest Diviner rod.

We installed five heating rods vertically into the soil to create a gradient of soil warming across the hillslope depending on distance from the rods. The rods were fabricated from 1.9 cm electric metallic tube conduit enclosing 2.5 m of heating cable (120 V,10 W/ft, Brisk Heating SLCAB Self Regulating cable) and backfilled the conduit with sand to promote conductance. The heating rods were controlled by a PID controller that maintained the level of heating at + 4 °C above ambient with a monitoring sensor 70 cm away from the rod placed 60 cm into the soil profile. Ambient temperature was calculated by averaging four temperature sensors outside the expected radius of heating (3 m away from rods and 60 cm deep). The temperature sensors consisted of a single high accuracy digital sensor mounted on a board (TSYS01, custom made by UH Manoa MESH lab) and enclosed in a 3 ml polystyrene test tubes filled with epoxy that reported soil temperature every 30 s to an online server.

We were able to maintain the targeted range of 0 to + 4 °C differences in mean soil temperatures between ambient and heated areas from November 2018 through November 2019 (Fig. S2). There were several inconsistencies in the soil warming across time because the heaters were turned off due to site vandalization in December 2018 and internet connectivity problems in February and July 2019. Even though the average ambient soil temperature varied throughout the year, soil temperature data shows we were able to maintain + 4 °C at the target radius (70 cm) and depth of heating (60 cm) for > 90% of the time preserving the target gradient of soil temperatures.

Soil CO2 production

Weekly sampling of soil CO2 concentrations, soil temperature, and soil moisture took place November 2018 through November 2019. At each sampling location, gas wells (0.635 cm stainless steel seamless tubing) were installed at a 45° angle in sets of three: at the depth of interest, the depth of interest − 20 cm, and the depth of interest + 20 cm. We sampled 5 ml of gas from each gas well into evacuated vials, whose CO2 concentrations were then analyzed on a Shimadzu GC-2014 (Shimadzu USA Manufacturing Inc, Canby, OR). A surface static chamber (10 cm diameter × 30 cm high PVC collar covered with a PVC cap installed with a septum during gas sampling) was used to measure the flux density across 0 cm (i.e., the surface respiration flux) needed to calculate CO2 production within the 0- 20 cm depth interval (see below).

Using the CO2 concentrations at each depth, we calculated CO2 production within each depth using the following series of equations based on Fick’s Law (Tang et al. 2003; Moyes and Bowling 2013; Contosta et al. 2016; Hicks Pries et al. 2017) (Eq. 1).

$$F = -D\frac{dC}{dz}$$
(1)

where F is the flux density (g m−3 h−1), D is the diffusion coefficient, and dC/dz is the change in concentration over depth. The diffusion coefficient was determined by Eq. 2:

$$D = {D}_{ao} {\left(\frac{T}{293.15}\right)}^{1.75}\left(\frac{101.3}{P}\right)\xi$$
(2)

where D is the diffusion coefficient of CO2 in soil, T is the temperature in °C at the time of sampling, P is local atmospheric pressure in kPa, and \(\xi\) is a dimensionless tortuosity factor. The tortuosity equation was an empirical relationship chosen from Moyes and Bowling (2013) (Eq. 3).

$$\xi =0.95{\varepsilon }^{1.93}$$
(3)

where ε is the air-filled porosity derived from the total soil porosity, bulk density each depth, and volumetric water content at each sampling depth during the time of measurement. To calculate CO2 production across a given depth interval, the flux densities across that interval were calculated by subtracting the flux density of the lower depth from the flux density of the upper depth and then dividing by the 20 cm depth increment. CO2 production is the soil respiration occurring within that depth interval, which is a mixture of root, root-derived, and microbial respiration.

Soil characterization

Initial soil cores (0.8 in diameter, 100 cm depth) were taken in May 2018, after site clearing and 6 months before soil warming at each of the 50 sampling locations to characterize the carbon (C), nitrogen (N), poorly- and non-crystalline minerals (PNCM), organo-mineral complexes, active Fe ratio, pH, ∆pH, and the richness of the microbial communities based on 16S and ITS sequences. 16S rRNA analysis provided genus-level information on bacteria while ITS analysis resulted in fungal inventory (see below for information on DNA extraction). We present only the data from samples that corresponded with the depth interval measured at each sampling location. Accurately sampling 50 1-m deep soil profiles for bulk density was not logistically possible and would create disturbances. Therefore we estimated bulk density at each of the 50 sampling locations by dividing the dry weight of the core sample (g) by the volume of the core (cm3). Soil cores were split into 0–20, 20–40, 40–60, 60–80, and 80–100 cm depths for the analyses. C and N were analyzed on dried, ground soil with a Carlo Erba NC 2500 Elemental Combustion System (Thermo Fisher Scientific, Waltham, MA). Three different extraction methods were used to determine the Fe and Al content: hydroxylamine hydrochloride (Ross et al. 1985), citrate dithionite (Holmgren 1967), sodium pyrophosphate (McKeague 1967). Hydroxylamine hydrochloride extraction was chosen in place of ammonium oxalate because it dissolves sustainably lower amounts of crystalline oxides (like magnetite common in Hawaii's volcanic soils) (Silva et al 2015; Chao and Zhou 1983; Ross et al. 1985). The extracted analytes were then measured using Inductively Coupled Plasma analysis (ICP). PNCM was calculated as the hydroxylamine hydrochloride extractable Al + 0.5Fe molar quantities per gram soil. Fe and Al organo-mineral complexes were calculated from the sodium pyrophosphate extraction. The active Fe ratio was calculated from the hydroxylamine hydrochloride extractable Fe divided by the citrate dithionite extractable Fe. pH was determined by using a 2:1 MilliQ water to soil ratio and measured using a Thermo Scientific benchtop pH meter. ∆pH was determined by the deionized water method—KCl method (1 M KCl, 2:1 KCl to soil ratio). ∆pH is a proxy of the organo-mineral association potential with a negative ∆pH indicating adsorption is more likely. The relative charges of the mineral and OM particles are dependent on soil pH via the point of net zero charge (pzc). Soils with OM above the pzc but minerals below the pzc will have oppositely charge particles resulting in adsorption. Soil minerals above the pzc (indicated by ∆pH greater than 0 in more acidic soils) will have positively charged OM and mineral particles, making adsorption less likely (Bailey et al. 2019).

DNA was extracted from 0.25 g of each soil sample using the Qiagen DNeasy PowerSoil Kit. Bacterial 16S (V4 region) rRNA genes were sequenced using the primer pairs 515F (Parada et al. 2016) and 806R (Apprill et al. 2015) using a dual-barcoding approach. For fungi, the ITS1 region was targeted using primer pairs ITS1F (White et al. 1990) and ITS2 (Gardes and Bruns 1993) using a single barcode on the ITS2 primer (Smith and Peay 2014). Amplified samples were pooled in equimolar concentrations and sequenced together as a single run using the Illumina MiSeq PE250 at the University of California, Berkeley QB3 Genomics center. Bioinformatics processing of the sequences was performed using QIIME2 (Bolyen et al. 2019) and an established pipeline (Nguyen et al. 2016) that uses stringent filtering criteria to minimize noise. For bacterial sequences, DADA2 exact sequence variance were considered as OTUs (operational taxonomic units) and fungal exact sequence variance were further clustered at 97% to obtain the best fungal OTUs.

Data analysis

First, we performed one-way ANOVAs with depth as a categorical main effect in R (R Core Team 2021). A significant depth effect was followed by a tukey post hoc test to determine which depths differed from one another. We used a Bonferroni correction of our alpha to account for the number of tests performed (13), so a significant effect had a p < 0.0036. Our response variables included the initial soil characterization data and the mean values of soil moisture (VWC), soil temperature, and soil CO2 production at the 50 sampling locations averaged over time. For the soil temperature data, Δ soil temperature (mean soil temperature − mean ambient soil temperature) at each sampling location was used as a measure of the soil warming treatment. Soil CO2 production was log transformed to meet assumptions of the ANOVA. Next, we completed a principal component analysis (PCA) using the same response variables to inform the parameters chosen for a multiple regression (prcomp function, R Core Team 2021). The PCA highlighted potential drivers of spatial variability and the importance of depth, supporting the need for regressions that included soil depth.

In order to model the most significant factors affecting the response of CO2 production across depths, we used generalized additive modeling (GAM) since it is robust to both parametric and non-parametric data (Hastie and Tibshirani 1986). Furthermore, our hypothesis predicted a non-linear response of CO2 production to temperature, which a GAM can capture. The GAM structure can be written as Eqs. 4:

$$g(E\left(Y\right)= \propto + {s}_{1}\left({x}_{1}\right)+\dots +{s}_{p}\left({x}_{p}\right)$$
(4)

where Y is the response variable, E(Y) is the predicted value, and \({s}_{1}\left({x}_{1}\right)+\dots +{s}_{p}\left({x}_{p}\right)\) are the estimated non-parametric variables fit with a smoothing function.

First, a global model combining the data from all depths was developed utilizing the “mgcv” package in R (Wood 2006) with the temporal means of CO2 production, VWC, and soil warming as well as C, N, pH, ΔpH, PNCM, active Fe ratio, OM complexed Al, OM complexed Fe, 16S richness, and ITS richness across all 50 sampling locations, each of which is a specific depth interval. To optimize the fit of the model, we log transformed the response variable, CO2 production. We checked the concurvity of predictors in this global model and dropped one from each pair of variables that had a concurvity value > 0.8. We then reran the global model without N (related to C), active Fe ratio (related to PNCM), and OM-complexed Al (related to OM-complexed Fe).

Next, we wanted to test how the controls on CO2 production differ with depth. The global model was limited by the number of smoothers (k) and therefore was unable to accommodate modeling by each of the sampling depths. When optimizing models for each depth, the number of smoothers limited the number of terms included in the model to four. Fitting by AIC determined that three smoothers (k) per model term limited underfitting while including the maximum terms for the optimal model. Given these limitations, four terms could be included in the depth-wise model. We set the soil warming treatment, Δ soil temperature, as one of these terms and then iterated the model through all possible model combinations with three other terms (chosen from among VWC, C, pH, ΔpH, PNCM, Fe-organic complexes, 16S richness, and ITS richness), each with three smoothers, by five depth groupings (0–20, 20–40, 40–60, 60–80, and 80–100 cm) for a total of 56 model combinations. Model selection was completed by Akaike information criterion (AIC) utilizing the “MuMIn” package in R (Barton 2019; Peters et al. 2019). Within this manual model iteration by depth grouping, the top model was selected based on the corrected Akaike information criterion (AICc) threshold below used for small sample sizes (n < 100) (Takezawa 2014):

$$\Delta i = AICc_{i} - AICc_{min}$$
(5)

where AICci is the AICc for the ith model, and AICcmin is the minimum of AICc for all model candidates. The optimal model is Δi = Δmin = 0 (Burnham and Anderson 2003). Therefore, the model with the lowest AICc was selected as the “best model” for representing the system.

Results

Many of the soil characteristics significantly differed by depth including %C (p < 0.0001), %N (p < 0.0001), Δ pH (p < 0.0001), ITS richness (p = 0.00017), active Fe ratio (p = 0.0006), PNCM (p = 0.0001), and OM-complexed Fe (p = 0.003); however, the temperature treatment (Δ Temp; p = 0.20), VWC (p = 0.046), pH (p = 0.043), OM-complexed Al (p = 0.006), and the 16S richness (p = 0.0097) did not based on the bonferroni-adjusted alpha of 0.0036(df = 4, n = 50 for all; Fig. 1). Soil respiration decreased significantly with depth; the most significant differences were between the surface 0–20 cm and bottom 60–100 cm depths (p < 0.0001, df = 4, n = 50; Fig. 2a). The amount of PNCM per gram of soil C increased significantly with depth (p = 0.003, df = 4, n = 50; Fig. 2b). The PCA of the temporally averaged (November 2018 to November 2019) and initial soil parameters explained 59% of the variance in the data within the first two principal components (Fig. 3; Table S1). The first principal component, which explained 42% of the variance, shows a gradient across depth with the shallowest sampling depth, 20 cm, having a strong correlation with CO2 production and soil C and N. Sampling depths below 40 cm were strongly correlated with the concentration of PNCM, the active Fe ratio, and ΔpH, the proxy index for net mineral charge.

Fig. 1
figure 1

Depth wise changes in warming treatment (∆ soil temperature), volumetric water content (VWC), carbon (C), nitrogen (N), pH, ΔpH, 16S and ITS-based microbial community richness, active Fe ratio, poorly crystalline and non-crystalline minerals (PNCM), and organically complexed (OM Comp.) Al (k) and Fe (l). The warming treatment and VWC were averaged over November 2018 through November 2019 for each depth interval. The middle lines are the medians, the boxes are the 25th to 75th percentiles, and the whiskers are 1.5 times the interquartile range

Fig. 2
figure 2

The natural log of CO2 production averaged over November 2018 through November 2019 for each depth interval (a) and the amount of poorly and non crystalline minerals per g soil carbon within each depth interval (b). The middle lines are the medians, the boxes are the 25th to 75th percentiles, and the whiskers are 1.5 times the interquartile range. One way ANOVA, depth effect, Df = 4, n = 50, p < 0.0001 for a and p = 0.003 for b)

Fig. 3
figure 3

First and second principal components from the mean timeseries and initial soil characterization datasets for the Lyon Arboretum field site. Top three positive and negative variables from the correlation matrix are listed for each principal component with corresponding correlation values

The global model of the controls on CO2 production explained 91.7% of the model deviance (R2 = 0.87). The most significant term from the global model was depth (p < 0.0001; Table 1). Other significant terms included soil warming (Δ Temp; p = 0.0033), VWC (p = 0.024), Δ pH (p = 0.027), and %C, which was nearly significant (p = 0.066; Table 1). As expected, CO2 production increased with the amount of soil warming and carbon content of the soil (Fig. 4a, c). In contrast, CO2 production decreased with increasing volumetric water content and depth (Fig. 4b, e). ΔpH increased the amount of CO2 produced until ΔpH around 0 after which it had little effect (Fig. 4d).

Table 1 Approximate significance of the smoothed terms from the best fit model by depth grouping (0–20, 20–40, 40–60, 60–80 and 80–100 cm), which had the following formula: log(CO2) ~ s(Δ soil temperature, k = 3, by = depth) + s(vwc, k = 3, by = depth) + s(P-NCM, k = 3, by = depth) + s(S_16S, k = 3, by = depth)
Fig. 4
figure 4

Global GAM predictions for all significant predictors of log-transformed CO2 production (n = 50). ∆ soil temperature (a), volumetric water content (VWC; b), ∆ pH (d), and depth (e) are all significant predictors (p < 0.05) and carbon (c) was marginally significant (p < 0.07). The points are the data colored by depth (see e for values), the lines are the predictions of each significant predictor with the other model variables held at their mean, and the grey ribbons are the 95% confidence intervals

Given the high significance of depth in the global model and in the PCA visualization, further model optimization by depth was needed. Model optimization provided three key predictors when included with the soil warming treatment and factored by five depth groupings (0–20, 20–40, 40–60, 60–80, and 80–100 cm). These key predictors were VWC, PNCM concentration, and the bacterial OTU (16S) richness (R2 = 0.93, Deviance explained = 96.8%, n = 50; Fig. 4, Table 1). Soil CO2 production increased significantly with the warming treatment in the surface 0–20 and 20–40 cm depths, but not in the deeper soil (Fig. 5a). As in the global model, CO2 production decreased with increasing VWC; an effect that was significant in the 20–40, 40–60, and 60–80 cm depths (Fig. 5b). CO2 production significantly declined with PNCM content in the surface soil only and had idiosyncratic relationships in the deeper soil (Fig. 5c). For example, CO2 production increased until about 40 g PNCM kg−1 after which it decreased at 20–40 and 60–80 cm depths, but had no significant relationships with PNCM content at 40–60 and 80–100 cm depths. Lastly, CO2 production significantly increased with increasing bacterial OTU richness at 0–20, 60–80, and 80–100 cm depths (Fig. 5d).

Fig. 5
figure 5

Plot of each smoothed term for the top GAM model by depth (top to bottom: 0–20, 20–40, 40–60, 60–80, 80–100 cm) as predictors of log transformed CO2 production. ∆ soil temperature (a), VWC (volumetric water content; b), poorly crystalline and non-crystalline minerals (PNCM; c), and 16S-based richness (d) were significant predictors in the depth wise GAM (p < 0.05). The points are the data colored by depth, the lines are the model predictions with the other model variables held at their mean for the given depth interval, and the grey ribbons are the 95% confidence intervals of the model predictions

Discussion

The importance of soil minerals in determining the soil warming response

Warming the whole profile of a tropical Andisol increased CO2 production within the top 40 cm by 50–300% per 1 °C increase in temperature, but had no significant effect on CO2 produced below 40 cm. Mineral protection, attributed to high concentrations of poorly and non-crystalline Al and Fe at depth, was likely the main factor inhibiting the temperature-dependent response of soil CO2 production. The lack of a warming response at depth contrasts the results of the only other whole-soil-profile warming experiment to partition soil respiration by depth where a 4 °C increase in soil temperature led to a 37% increase in surface CO2 flux with 40% of that response coming from deeper soils (Hicks Pries et al. 2017; Soong et al. 2021). In another whole-soil-profile tropical soil warming experiment, SWELTR, depth-resolved fluxes were not quantified; however, 4 °C of warming caused a substantial 55% increase in surface CO2 flux (Nottingham et al. 2020). Our surface soils had an even larger response, with CO2 production doubling in response to 2 °C of warming at 0–20 cm and increasing by a factor of 6 in response to 2 °C of warming at 20–40 cm. This strong response may have been driven in part by decaying root biomass as a result of our site preparation. Despite their much warmer climate, tropical soils can be just as responsive to warming as their temperate and arctic counterparts, at least in the short term (Lu et al. 2013; Wood et al. 2019).

The unresponsiveness of deep soil respiration to experimental warming demonstrates the importance of mineral-associations in reducing the susceptibility of SOM to climate change-driven losses (Rocci et al. 2021). Shallow SOM is often more accessible to soil microbes who can degrade it readily; whereas in deep soils, a greater proportion of the carbon is associated with soil minerals that can stabilize and protect SOM from microbial decomposition (Adhikari et al. 2019; Porras et al. 2018; Jackson et al. 2017). Hawaiian Andisols are known for their strong organo-mineral associations driven by high concentrations of poorly and non-crystalline minerals that protect carbon through sorption, aggregation, or co-precipitation reactions (Cismasu et al. 2015; Lalonde et al. 2012) leading to long soil carbon residence times (Torn et al. 1997). These PNCM increased with depth in our soils and were a significant control on soil respiration at several depths such as at the surface where soil respiration declined with increasing PNCM. While we do not have direct measures of the amount of mineral-associated carbon in these soils, the average amount of PNCM per carbon increased by an order of magnitude from 0.2 g PNCM per g C at 0–20 cm to 3 g of PNCM per g C at 80–100 cm, so carbon at depth had a greater likelihood of being mineral-associated (Fig. 2b). Soils with high concentrations of PNCM such as allophane, imogolite, and primary Al-OM complexes are common in Andisols and other soils derived from volcanic ash (Nanzyo 2002). Even at low concentrations, PNCM can protect OM from destabilization and decomposition (Kramer and Chadwick 2016). High concentrations of both PNCM and soil C typically persist over time deep in the soil profile of Hawaiian Andisols (e.g., Giardina et al. 2014). While the importance of PNCM in controlling carbon abundance and stabilization is well known (Giardina et al. 2014; Torn et al. 1997), this is the first study to demonstrate its importance in protecting soil carbon from loss when faced with climate warming.

The significant relationship of soil respiration to ∆pH in the global GAM model further supports the importance of organo-mineral associations in controlling the response of soil respiration to warming. Soil respiration decreased with more negative ∆pH values below zero. ∆pH values less than zero indicates stronger mineral associations because minerals are more likely to have a net negative charge to facilitate adsorption of positively charged organic matter. At ∆pH greater than zero, where minerals and organic matter are less likely to have opposing charges weakening mineral protection of organic matter, there was an increase in soil respiration rates.

Organo-mineral associations also explain why our results contrast those of a previous whole soil warming experiment in a temperate forest. At that site, there was a very strong response of deeper soil respiration to warming driven by the large proportion of free particulate carbon at depth, which was preferentially lost relative to mineral-associated carbon over 5 years (Soong et al. 2021). That soil was developed from a granitic parent material, which had weaker mineral protection of organic carbon than nearby andesitic soils (Rasmussen et al. 2005) and had similar proportions of particulate and mineral-associated carbon throughout the soil profile (Hicks Pries et al. 2017).

Other drivers of CO2 production throughout the profile

Soil moisture reduced CO2 production. At this tropical site there are 4200 mm of precipitation a year leading to volumetric water contents of 20.8–49.4% throughout the year and greater than 35% annually. Other tropical soil respiration studies found soil moisture had a significant parabolic effect on CO2 produced (Hashimoto et al. 2004; Wood et al. 2013). Interestingly, the tipping point of that relationship from positive to negative was a volumetric water content of 35% in Wood et al. (2013), the lowest average annual water content measured in this study. High volumetric water content indicates increased water-filled pore space that limits oxygen diffusion leading to anaerobic sites within the soil profile that can reduce microbial respiration (Silver et al. 1999; Henry 2012). Reduced diffusion of CO2 out of the soil due to increased water-filled pore space, a mechanism sometimes cited for reductions in surface CO2 production at high moisture levels (e.g., Wood and Silver 2012; Wood et al. 2013), was likely less of a factor in this study wherein CO2 production was measured within individual depth increments via gas wells. In these very wet soils, soil moisture was not affected by soil warming as in other experiments (linear regression, df = 48, p = 0.51, n = 50; Xu et al. 2013).

As expected, soil respiration increased with increasing soil carbon content. Unfortunately, we were unable to quantify the amount of carbon within the volume of each depth interval since a proxy was used to estimate bulk density due to logistical reasons. Bacterial richness, which was positively related to soil carbon content (linear regression, p = 0.016, n = 50), was a significant predictor in the depth-wise model and, as with soil carbon concentrations, increased CO2 production. It is unknown whether this effect was driven by bacterial diversity itself (Delgado-Baquerizo et al. 2016; Liu et al. 2018) or another metric that scales with bacterial diversity such as microbial biomass or carbon amount (Bastida et al. 2021). However, it should be noted that the relationship between microbial diversity and biomass is negative in tropical forests (Bastida et al. 2021). Overall, bacterial richness being significant in the depth-wise model instead of soil carbon indicates that bacterial richness may be a better indicator of carbon availability at our site than carbon itself. As previously stated, soil carbon at depth is more likely to be protected from microbial decomposition by minerals, decoupling the relationship between soil carbon and respiration in this tropical Andisol. The mechanisms by which microbial communities interact with mineralogical controls and soil carbon will be tested in future research at this site.

Conclusion

Understanding the feedbacks of soil respiration to rising global air temperatures provides critical information to climate models on how climate changes may impact the terrestrial carbon cycle. Climate change experiments in underrepresented soil and climate regimes must be expanded upon in order to develop an accurate projection of future of soil organic carbon storage and sequestration. This study piloted a low-cost methodology for warming and monitoring soils while adding to the lack of empirical evidence on the effects of soil warming in tropical soils and Andisols. Andisols are critical soils to be integrated into models given the demonstrated role of their unique mineralogy on SOC stabilization, which can be used to refine the relationship between mineral-protected carbon and temperature sensitivity.

Unlike any other deep soil warming experiment to date, there was no significant CO2 response to soil warming below 40 cm of the profile. Multimodal analyses confirmed the hypothesis that high concentrations of PNCM was the primary driver in the lack of CO2 response, followed by high relative soil moisture and low bacterial richness at depth. The lack of CO2 response due to potential PNCM protection suggests a limited positive feedback to warming within Andisols, especially at depth. Therefore, Andisols could be explored as a priority for restoration or protection areas when considering management and policy decisions for climate action.