Soil depth as a driver of microbial and carbon dynamics in a planted forest (Pinus radiata) pumice soil

Forest soils are fundamental in regulating the global carbon (C) cycle; their capacity to accumulate large stores of C means they form a vital role in mitigating the effects of climate change. Understanding the processes that regulate forest soil C dynamics and stabilisation is important to maximise the capacity and longevity of C sequestration. Compared with surface soil layers, little is known about soil C dynamics in subsoil layers, sensu those below 30 cm depth. This knowledge gap creates large uncertainties when estimating the distribution of global soil C stocks and assessing the vulnerability of soil C reserves to climate change. This study aimed to dive deep into the subsoils of Puruki Experimental Forest (New Zealand) and characterise the changes in soil C dynamics and the soil microbiome down to 1 m soil depth. ITS and 16S rRNA sequencing and quantitative real-time PCR were used to measure changes in soil microbial diversity, composition, and abundance. Stable (δ¹³C) and radioactive (¹⁴C) C analyses were performed to assess depth-driven changes in the stability and age of soil C. Our research identified large declines in microbial diversity and abundance with soil depth, alongside significant structural shifts in community membership. Importantly, we conservatively estimate that more than 35 % of soil C stocks are present in subsoil layers below 30 cm. Although the age of soil C steadily increased with depth, reaching a mean radiocarbon age of 1571 yr BP (years before present) in the deepest soil layers, the stability of soil C varied between different subsoil depth increments. These research findings highlight the importance of quantifying subsoil C stocks for accurate C accounting. By performing a broad range of analytical measures, this research has comprehensively characterised the abiotic and biotic properties of a subsoil environment – a frequently understudied but significant component of forest ecosystems.

1 Introduction

Forest soils are influential in regulating the global C cycle, accounting for over 70 % of the world's soil organic C (SOC) reserves and supporting over 80 % of aboveground terrestrial C stocks (Batjes, 1996; Jobbágy and Jackson, 2000; Six et al., 2002). Developing forest management strategies that maximise soil C sequestration may serve as a powerful tool for offsetting fossil fuel emissions and mitigating the effects of climate change (Jandl et al., 2007; Mukul et al., 2020).

Compared with surface soil layers, the deeper subsoil layers of forest ecosystems are comparatively under-explored (Yost and Hartemink, 2020). Although the concentration of SOC declines with depth (Schmidt et al., 2011), studies have reported that, accumulatively, up to 50 % of total SOC may be contained within soil layers below 30 cm (Balesdent et al., 2018; Jobbágy and Jackson, 2000; Gonzalez et al., 2018; Ross et al., 2020).

Deep SOC is more resistant to decomposition by soil microorganisms due to the biological, physical, and chemical conditions of the subsoil environment (Schmidt et al., 2011). Thus, when undisturbed, deep SOC is considered to represent the “stable” fraction of the SOC pool. However, the stability of sequestered SOC is complex and is related to a range of variables including soil clay content, mineralogy, structure and texture, landscape position, association with microaggregates, climate, vegetation, and soil microorganisms (Jobbágy and Jackson, 2000; Lal, 2004; Kuzyakov, 2010; Song et al., 2020). Historical biases against sampling deeper soil layers limit our understanding of how the properties of subsoil environments influence SOC sequestration and storage. This introduces uncertainty when modelling the stability of deep SOC pools to climate change (Gross and Harrison, 2019).

Soil microorganisms are responsible for the transformation and stabilisation of SOC (Lladó et al., 2017; Spohn et al., 2016) and are influential in determining patterns of C sequestration (Jastrow et al., 2007). Studying the changes that occur to soil microbial communities with depth is essential to understand the nature and extent of their role in regulating long-term SOC dynamics. Previous studies have reported that soil microbial diversity, abundance, respiration, and C turnover decline with depth (Blume et al., 2002; Eilers et al., 2012; Fang et al., 2005; Fierer et al., 2003; Spohn et al., 2016). More recent research has explored the taxonomic and functional changes that occur in microbial communities with soil depth (Brewer et al., 2019; Frey et al., 2021; Hao et al., 2021). However, particularly for coniferous forest ecosystems, our knowledge of the response of soil microbial communities to depth and the environmental properties that drive these changes remains poorly understood. Furthermore, the inherently heterogeneous nature of soil environments, both globally and locally (Curd et al., 2018; Kuzyakov and Blagodatskaya, 2015; Štursová et al., 2016), means context-specific research is required to understand the SOC dynamics and associated microbiota of forest subsoil environments. To maximise the sink capacity of forest soils and to retain the SOC already stored at depth, we need to better understand the processes that regulate C sequestration and long-term C storage within the deep soil environment.

Using Puruki Experimental Forest (central North Island, New Zealand) as a case study, this research aimed to improve our understanding of subsoil C dynamics within a production forest. Using a wide range of analytical methods, we examined the variability in subsoil C in 10 cm increments down a 1 m soil profile. Changes in the diversity, composition, and abundance of bacterial and fungal communities with depth were quantified using 16S rRNA/ITS amplicon sequencing and quantitative real-time PCR (qPCR). To investigate changes in C cycling dynamics down each soil core, measurements of soil isotopic composition (¹³C $/$ ¹²C ratio) were performed, with the ¹³C $/$ ¹²C ratio of soils previously having been used as a proxy for soil C turnover (Garten et al., 2008; Paul et al., 2020; Poage and Feng, 2004; Wang et al., 2018). Radiocarbon ageing based on measurements of ¹⁴C abundance was further used to identify older stabilised C stocks (Chabbi et al., 2009). The results of this research contribute to our understanding of the biological and chemical properties of productive forest subsoil environments, a critical research gap that is fundamental for assessing the vulnerability of forest SOC stocks to climate change.

2 Materials and methods

2.1 Study site and soil collection

Soils were sampled from Puruki Experimental Forest (38^∘26^′ S, 176^∘13^′ E), located in the Paeroa Range of New Zealand's central North Island (Fig. S1 in the Supplement). Puruki Experimental Forest (henceforth referred to as Puruki Forest) is a 35 ha Pinus radiata forest that was established in 1973, after previously being converted from podocarp/hardwood native forest to ryegrass and clover pasture in 1957 (Beets and Brownlie, 1987; Garrett et al., 2021). Puruki Forest is now near the end of its second rotation, having previously been harvested in 1997 (Oliver et al., 2004). Puruki Forest was selected for investigation due to its high scientific impact and extensive contribution to production forestry research (Garrett et al., 2021). Furthermore, Puruki Forest has an extensive body of historical data available related to soil health, productivity, and land management – factors which are important for helping us understand the drivers of subsoil C dynamics (Garret et al., 2021).

Puruki Forest sits at an elevation range of 500–700 m a.s.l. and varies from gently (less than 12^∘) to steeply sloping (up to 30^∘) land, with a mean annual rainfall of 1500 mm and mean annual temperature of 10 ^∘C (Beets and Brownlie, 1987; Brownlie and Kelliher, 1989). Puruki Forest has an average soil pH of 5.2 at 0–10 cm depth (Beets and Brownlie, 1987). The soil pH at Puruki Forest has been reported to increase with depth, having an average pH of 5.62 at 1–2 m depth (Beets et al., 2004).

Soil parent materials at Puruki originated from the Taupo volcanic centre (1850 ± 100 BP), with older ash showers from Taupo and Okataina volcanic centres (Rijkse and Bell, 1974). The Orthic Pumice soils found at Puruki are highly permeable, being composed of loamy sand, silty sand, and gravel (Taupo lapilli) (New Zealand Soil Classification System; Hewitt, 2010).

A permanent sample plot (PSP) within the Rua sub-catchment of Puruki Forest was targeted for sampling. As Puruki Forest has a varied topology, our sampling plot was located along a 12^∘ slope. Soils were collected just outside of the plot to avoid disturbance for longer-term monitoring at the site. For sampling, 10 points spaced 2 m apart were set out along an 18 m transect. At each sampling point, two paired 1 m soil cores were extracted using a motorised percussion soil sampler capable of collecting intact cores with an internal plastic sleeve; one core was taken for bulk density analysis, and the other core was used for DNA and chemical analysis. Following extraction, core pairs were visually compared to check for any considerable differences in soil colour and texture by depth (Fig. S2). Intact soil cores were transported back to the lab in cool, dark conditions. Care was taken to minimise disruption during transport so that cores remained intact from field to lab. Following transport, soil cores were divided into 10 cm increment samples. Incidences of compaction which occurred during sample collection were adjusted for during the division of cores into depth increments (Supplement Methods).

2.2 Sample preparation

To prepare soils for DNA analysis, moist soil samples were sieved to less than 2 mm and stored at −20 ^∘C until ready for DNA extractions. DNA was extracted from 0.25 g of each soil sample using a DNeasy Powersoil Pro Kit (Qiagen) according to the manufacturer's instructions. Soils required for measurements of total C and nitrogen (N), total phosphorus (P), Bray P, pH, Mehlich 3 extractable elements (Ca, K, Mg, and Na), cation exchange capacity (CEC), exchangeable cations, organic P, and inorganic P were air dried and sieved to < 2 mm. For measurements of bulk density, air-dried fine (< 2 mm) and coarse (> 2 mm) soil fractions were oven dried separately at 104 ^∘C until constant weight. Subsamples of air-dried soils were fine ground using a Spex800D ball mill for measurements of total C (%), total N (%), δ¹³C [‰], and radiocarbon ¹⁴C analysis [‰].

2.3 Soil chemistry

Total C and N were determined on fine (< 2 mm) and coarse (> 2 mm) fraction fine-ground soils by dry combustion using a LECO CNS-2000 carbon nitrogen analyser (Rayment and Lyons, 2011). Soil pH was determined in water at a soil : water ratio of 1:2.5 (Blakemore et al., 1977). Total phosphorus (P) was measured by flow injection analyser (FIA) colorimetry after sulphuric acid digestion (Blakemore et al., 1977; Taylor, 2000). Bray P was measured by FIA colorimetry after sequential Bray 2 (NH₄F $/$ HCl) extraction at a ratio of 1:10 soil : extractant (Blakemore et al., 1977; Bray and Kurtz, 1945). Exchangeable cations (Ca, K, Mg, and Na) were measured by inductively coupled plasma mass spectrometry (ICP-MS) after 1:50 (macro) NH₄CH₃COO leaching. Cation exchange capacity (CEC) was measured by FIA colorimetry after 1 M NH₄COOH leaching, followed by 1M NaCl (Blakemore et al., 1977; Sparks et al., 1996).

2.3.1 Stable isotope (δ¹³C) analysis

Fine-ground soil samples were sent to GNS Stable Isotope Laboratory (Wellington, New Zealand) for δ¹³C analysis by isotope ratio mass spectrometry. Soil samples were pre-treated by adding 10 % hydrochloric acid (HCl) to 1 g of each sample and left overnight at room temperature. Samples were then centrifuged and rinsed to neutral pH. Following pre-treatment, samples were freeze dried and measured on an Isoprime mass spectrometer using a Eurovector elemental analyser. All δ¹³C results are reported with respect to VPDB and N-Air, normalised to GNS internal standards which were cane sugar (−10.3 ‰), beet sugar (−24.6 ‰), and EDTA (−31.1 ‰). The stable carbon isotope composition (¹³C $/$ ¹²C) is reported in δ¹³C values presented as per thousand [per mille, ‰], with the analytical precision for these measurements being 0.2 ‰.

2.3.2 Radiocarbon ¹⁴C analysis

Fine ground soil samples were graphitised at the Houghton Carbon, Water, and Soils Lab (USDA Forest Service Northern Research Station, Michigan, US), and radiocarbon measurements were conducted at the DirectAMS facility (Bothell, Washington, US). Soil samples were combusted at 900 ^∘C for 6 h in evacuated, sealed quartz tubes with cupric oxide (CuO) and silver (Ag) wire to form CO₂ gas. After combustion, CO₂ was reduced to graphite through heating at 570 ^∘C with hydrogen (H₂) gas and an iron (Fe) catalyst (Vogel et al., 1987). Graphite targets were analysed for radiocarbon abundance using accelerator mass spectrometry (AMS) and corrected for mass-dependent fractionation using measured δ¹³C values according to Stuiver and Polach (1977). Conventional Radiocarbon Age (CRA) and ¹⁴C are reported as defined by Stuiver and Polach (1977), with CRA presented as years before present (yr BP) and ¹⁴C presented as per thousand [per mille, ‰].

2.4 Microbial amplicon sequencing

Following the Earth Microbiome Project (EMP) protocols, barcoded forward and reverse primers 515F-806R were used to amplify the 16S rRNA V4 region (Caporaso et al., 2012). Nested PCRs were performed using the non-barcoded primers 27F-1492R on DNA samples that initially failed to amplify. The EMP's standard ITS amplicon protocol was followed to amplify the fungal ITS gene region using the barcoded primers ITS1f-ITS2 (Bokulich and Mills, 2013; Hoggard et al., 2018). DNA samples that initially failed to amplify were targeted for nested PCR using primers NSA3-ITS4. Full details of the PCR reaction protocols are provided in the Supplement Methods. Successful PCR amplicons were sent for Illumina 500 MiSeq 250 bp (16S rRNA) and Illumina 600 MiSeq 300 bp (ITS) paired-end sequencing at the Australian Genome Research Facility (AGRF; Melbourne, Australia). For the 16S rRNA amplicon libraries, eight no-template control (NTC) samples, which showed no visible bands during agarose gel electrophoresis, were sent for sequencing to check for contaminants.

2.5 Quantitative PCR

The absolute concentration of bacterial and fungal DNA in each DNA sample was determined from standard curves using AriaMx SYBR Green qPCR (Agilent Technologies). Broad-range forward and reverse primers 338F-518R and ITS1F-5.8s were selected to target the 16S rRNA (V3–V4) and ITS gene region for quantification (Shahsavari et al., 2016). Full details of the qPCR reaction protocol and standard curve preparation are provided in the Supplement Methods. Following quantification, 16S rRNA and ITS copy numbers were adjusted to copy numbers per nanogram of DNA.

2.6 Statistical analysis

Before analysis, slope-corrected C stocks (Mg ha⁻¹) were calculated based on soil C concentration (Total C %) and soil bulk density measurements. Slope-corrected stocks were calculated using the following formula (Jones et al., 2008):

$$(\mathrm{BD}\times {\mathrm{SOC}}_{\mathrm{conc}}\times D)\times \mathrm{SF}$$

where BD is bulk density (g cm⁻³), SOC_conc is total C (%), D is the thickness of the soil layer (cm), and SF is the slope factor. The values for total C, total N, and C stocks are reported using the sum of both the < 2 and > 2 mm soil fractions. Values for total P, Bray P (sequential extraction 1), and exchangeable cations (K, Ca, Mg, Na) were slope corrected and converted to kg ha⁻¹. The changes in soil pH, total C (Mg ha⁻¹), total N (kg ha⁻¹), C : N ratio, δ¹³C [‰], total P (kg ha⁻¹), Bray P (kg ha⁻¹), total bulk density (g cm⁻³), and exchangeable cations K, Ca, Mg, and Na (kg ha⁻¹) with soil depth were tested for significance using non-parametric Friedman tests with transect sampling position as a blocking variable. Post hoc tests were performed using pairwise Wilcoxon rank sum tests with Bonferroni correction.

To calculate differences in the allocation of C stocks down each soil profile, the proportion of C (Mg ha⁻¹) allocated to each depth increment was calculated relative to the sum of total C (Mg ha⁻¹) for each sample core. Differences in 16S rRNA and ITS copy number with depth were tested for significance using Kruskal–Wallis chi-squared tests, with post hoc testing performed using Dunn's test with Bonferroni's adjustment. To estimate the allocation patterns of microbial DNA abundance down the soil profile, the proportion (%) of 16S rRNA and ITS copy numbers allocated to each depth increment was calculated relative to the total sum of the soil core.

2.6.1 Microbial diversity analysis

Base calling and demultiplexing of sequencing libraries to per-sample fastQ files was performed using Illumina's bcl2fastq software (version 2.20.0.422). The DADA2 version 1.18 workflow (Callahan et al., 2016) was used to process paired-end fastQ files into amplicon sequence variants (ASVs). Briefly, forward and reverse reads were quality filtered, trimmed, and de-noised before being merged into ASVs. Chimeric ASVs were removed, and taxonomies were assigned to ASVs using the Ribosomal Database Project (RDP) Classifier (Wang et al., 2007) and the UNITE (Abarenkov et al., 2021) databases. After DADA2 processing, 74.92 ± 5.60 % of 16S rRNA input reads and 44.46 ± 14.35 % of ITS input reads were retained. Following DADA2 processing, ASV tables were filtered to remove unidentified phyla, unwanted phyla (i.e. Cyanobacteria and/or Chloroplasts), and singletons. For the 16S rRNA dataset, ASVs which had a higher average count in NTC samples than eDNA samples were removed. Rarefaction curves showing ASV richness by depth increment and PCR type can be seen in Figs. S3 and S4. Following the removal of samples with low ASV counts, ASV tables were rarefied to even sampling depth before analysis (16S rRNA: 11613 ASVs per sample, ITS: 10101 ASVs per sample).

Microbial alpha and beta diversity analyses were calculated on rarefied ASV tables using the R packages phyloseq (McMurdie and Holmes, 2013) and vegan (Oksanen et al., 2020). Changes in alpha diversity with depth were tested for statistical significance using Kruskal–Wallis chi-squared tests as described previously. Pairwise Spearman's correlation tests with Bonferroni adjustment were used to correlate alpha diversity values with soil chemical parameters using psych R (Revelle, 2021). Bray–Curtis distance matrices were calculated to measure community dissimilarity and ordinated using non-metric multidimensional scaling (NMDS). Differences in community dissimilarity by depth were tested for significance using PERMANOVA, and pairwise differences were tested using pairwiseAdonis R with Bonferroni adjustment (Arbizu, 2017). For both alpha and beta diversity analyses, the impacts of transect sampling position and PCR type (nested or single) were tested for significance. Mantel tests were used to correlate Euclidean distance matrices of soil chemical properties to Bray–Curtis matrices of microbial community composition. Venn diagrams were produced to show the proportion (%) shared versus unique ASVs between depth increments using MicEco R (Russel, 2021).

Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC; Lin and Peddada, 2020) was used to identify microbial taxa which had a significant log-fold change between topsoil (0–30 cm) and subsoil layers (30–100 cm). Hierarchical clustering was performed on Euclidean distances of log-adjusted abundances, using the average linkage method (hclust R). Dendrograms were constructed using ggtree (Yu et al., 2017) to visualise the hierarchical clustering patterns of microbial orders based on their log-adjusted abundances. The log-adjusted abundances of microbial taxa obtained from ANCOM-BC analysis were correlated to soil chemical properties using pairwise Spearman's correlation tests with Bonferroni adjustment using psych R (Revelle, 2021). To assess changes in the diversity of high-level 16S rRNA functional genes with soil depth, PICRUSt2 was used to predict MetaCyc pathway abundances (Douglas et al., 2020; Caspi et al., 2018). MetaCyc pathway abundance predictions were analysed in phyloseq R using the same alpha and beta diversity methods as already outlined.

2.6.2 Microbial network analysis

ASV tables were split into samples from topsoil (0–30 cm) and subsoil (30 –100 cm) layers to obtain microbial communities specific to each soil layer. Topsoil and subsoil ASV tables were then filtered to retain only “core” communities (i.e. ASVs with a total relative abundance of at least 0.05 % in a minimum of three samples) and then summarised to genus level.

SPIEC-EASI networks were constructed using the Meinshausen–Buhlmann (MB) neighbourhood algorithm in SpiecEasi R (Kurtz et al., 2015). Networks were visualised using igraph R (Csárdi and Nepusz, 2006) and network statistics were calculated using Cytoscape 3.8.2 (Shannon et al., 2003). Significant differences (p < 0.05) in network statistics between topsoil and subsoil communities were tested using Wilcoxon signed rank tests (Mundra et al., 2021). Descriptions of these network statistics can be found in the Supplement Methods.

3 Results

3.1 Soil chemistry, C stocks, and microbial DNA abundance

Total C (Mg ha⁻¹), total N (kg ha⁻¹), total P (kg ha⁻¹), C : N ratio, and exchangeable cations (K, Ca, Mg, and Na; kg ha⁻¹) significantly declined with soil depth (p < 0.05). Soil pH and bulk density (g cm⁻³) significantly increased with soil depth (p < 0.05). The results of soil chemical analyses and the corresponding statistical tests are presented in Tables S1 and S2.

The average soil C stock of the 1 m soil cores extracted from Puruki Forest was 99.71 ± 28.78 Mg ha⁻¹ C. The average soil C stock allocated to the coarse (> 2 mm) soil fraction was 1.88 ± 1.92 Mg ha⁻¹ of C, values which did not vary by soil depth (ANOVA: F = 1.58; p > 0.05).

On average, 34.58 ± 10.25 % of the soil C stock was allocated to the 30–100 cm (subsoil) layer (Fig. 1; Table S3 in the Supplement). The proportion of soil C stock allocated to the subsoil varied across the soil cores, ranging from 46.60 % to 16.21 %. On average, 33.45 ± 10.92 % of fungal DNA and 49.11 ± 4.82 % of bacterial DNA was allocated to the subsoil (Fig. 1; Tables S3 and S4). Thus, the allocation patterns of soil C with depth corresponded well to the distribution of bacterial and fungal DNA with depth.

Figure 1 The relative proportion (%) of (a) total carbon stocks (C Mg ha⁻¹), (b) bacterial DNA abundance, and (c) fungal DNA abundance allocated to each depth increment. The values for soil C stocks are reported using the sum of < 2 and > 2 mm soil fractions. The transect position (m) shows where each soil core was extracted from along the sampling transect. The red line shows the typical cutoff point of the topsoil subsoil boundary, where measures below 30 cm are rarely taken.

3.2 δ¹³C and ¹⁴C natural abundance measurements

Overall, δ¹³C [‰] significantly increased by 1.58 ‰ (p < 0.05; Table S1) between depth increments of 0–10 and 90–100 cm. However, between soil depths of 70–90 cm the enrichment of ¹³C was lower than the rest of the subsoil (Fig. 2). Values for ¹⁴C [‰] and the CRA [yBP] of soil C declined with depth (Fig. 2; Table S1). The age of soil C steadily increased with depth and reached a mean age of 1570.83 ± 316.28 yr BP in the deepest soil layers of 90–100 cm.

Figure 2 The mean ± SD changes in total C stock (Mg ha⁻¹), δ¹³C [‰], Δ¹⁴C [‰], and conventional radiocarbon age (CRA) [yr BP] with soil depth.

3.3 Changes in microbial diversity, abundance, and composition with depth

The richness, diversity, and evenness of fungal and bacterial communities significantly declined with depth (p < 0.05; Fig. 3; Tables S6 and S7). The same trends were evident for the predicted abundances of bacterial functional genes with depth (Table S7). Transect position had no relationship to bacterial alpha diversity; however, fungal evenness did vary across sampling locations (Table S7). The DNA abundance of bacteria and fungi significantly declined with soil depth (Fig. 3; Table S7) but was unaffected by transect position. Most of the significant differences in alpha diversity and DNA abundance occurred between the upper 0–20 cm depth increment and the rest of the soil; few significant pairwise differences were evident among soil layers sampled below 30 cm (Figs. S5–S8). When tested using pairwise Spearman's rank correlations (Fig. 4), microbial alpha diversity and abundance were strongly positively correlated to total C, total N, total P, C : N ratio, and ¹⁴C (ρ > 0.77; p < 0.05). Soil pH was strongly negatively correlated to microbial alpha diversity and abundance (ρ > −0.85; p < 0.05). δ¹³C was strongly negatively correlated to fungal richness (ρ = −0.85, p < 0.05) and bulk density (g cm⁻³) was strongly negatively correlated to bacterial abundance (ρ = −0.90, p < 0.05).

Figure 3 The mean ± SD changes in microbial richness, Shannon diversity, Pielou's evenness, and 16S rRNA/ITS copy number with soil depth.

Figure 4 A pairwise Spearman's correlation matrix showing the correlation between microbial richness, diversity, DNA abundance, and the associated soil chemical properties. Blank cells indicate non-significant (p > 0.05) correlations.

Bacterial (PERMANOVA: R²=0.21, p < 0.001) and fungal (R²=0.26, p < 0.001) community composition varied significantly between depth increments (Fig. 5). Pairwise adonis tests identified that the community composition of soils between 0 and 30 cm contributed to the majority of significant differences between groups, with the differences in microbial composition between soil increments past 30 cm being non-significant (Fig. S9). The composition of bacterial and fungal communities significantly differed by sampling transect position (bacteria: R²=0.13, p < 0.001; fungi: R²=0.04, p < 0.001). However, this accounted for less variation in community composition than soil depth. The composition of bacterial functional genes was significantly different between soil depth increments (R²=0.34, p < 0.001; Fig. S10).

Figure 5 NMDS ordination plots that were calculated using Bray–Curtis dissimilarity matrices showing the differences in (a) bacterial and (b) fungal community composition between the different soil depth increments.

Mantel tests identified that ¹⁴C (R=0.38, p < 0.001), C : N ratio (R=0.18, p < 0.05), bulk density (R=0.20, p < 0.05), Total p (R=0.16, p < 0.001), δ¹³C (R=0.13, p < 0.05), total N (R=0.13, p < 0.05), and exchangeable Ca (R=0.12, p < 0.01) shared a significant but weak correlation to bacterial community composition. Variables that shared a significant but weak correlation to fungal community composition included ¹⁴C (R=0.32, p < 0.001), δ¹³C (R=0.19, p < 0.01), total P (R=0.19, p < 0.01), total N (R=0.17, p < 0.01), total C (R=0.14, p < 0.05), and exch. Ca (R=0.21, p < 0.01). Furthermore, total N (R=0.17, p < 0.01), total P (R=0.17, p < 0.01), exchangeable Ca (R=0.13, p < 0.01), bulk density (R=0.12, p < 0.05), δ¹³C (R=0.18, p < 0.01), and ¹⁴C (R=0.38, p < 0.001) had a significant but weak correlation to bacterial functional genes.

3.4 Microbial taxonomic responses to soil depth

Only 9 % of bacterial ASVs and 2 % of fungal ASVs were present in all depth increments down the entire soil profile (Fig. S11). Upper depth increments (0–20 cm) had the highest number of unique bacterial (20 %) and fungal (32 %) ASVs, suggesting that these upper layers supported a greater range of unique bacterial and fungal taxa than the lower depths. Less than 3 % of bacterial ASVs and 1 % of fungal ASVs were found to be unique within depth increments from the 40 to 100 cm soil layer, suggesting that subsoils did not harbour vastly unique or distinctive microbial communities from those of topsoil layers.

Archaeal taxa belonging to Thermoplasmata, Desulfurococcales, Thermoprotei, Woesearchaeota incertae sedis (AR15, AR18, AR16, AR20), Euryarchaeota, Candidatus Aenigmarchaeum, and Thermoproteales exhibited significant positive log-fold changes with depth (Fig. 6). Furthermore, bacterial taxa belonging to Acidobacteria (e.g. Bryobacter, Gp20, Gp18, Gp25, Gp19, Gp6, Acidicapsa), Elusimicrobia (Blastocatella, Candidatus Endomicrobium), Planctomycetes (Phycisphaerales, Tepidisphaerales, Phycisphaerae), Firmicutes (Lactobacillales, Clostridia, Thermoanaerobacterales), Betaproteobacteria (Neisseriales, Ferritrophicales, Hydrogenophilales, Methylophilales), and Armatimonadetes (Gp2, Fimbriimonadales) exhibited a positive log-fold change in subsoils (Fig. 6). Fungal taxa which exhibited a positive log-fold change in subsoils (Fig. 6) included Calcarisporiellomycota, Exobasidiomycetes (Entylomatales), Sordariomycetes (Microascales, Trichosphaeriales, Annulatascales, Ophiostomatales), Agaricomycetes GS29, Sareomycetes (Sareales), Orbiliomycetes (Orbiliales), Malasseziomycetes (Malasseziales), and Kickxellomycetes (Kickxellales and GS29).

Figure 6 Cluster dendrograms of bacterial and fungal orders based on their bias-adjusted abundances (ANCOM-BC analysis). Hierarchical clustering was performed on Euclidean distances of bias-adjusted abundances using the average linkage method. The log-fold change (LFC) value of each taxon between topsoil and subsoil is displayed. Taxa with positive LFC values were higher in subsoils and taxa with negative LFC values were lower in subsoils. The statistical significance of each taxa's differential abundance between topsoil and subsoil is displayed, with TRUE = p < 0.05 and FALSE = p > 0.05.

In contrast, bacterial taxa which had a significantly negative log-fold change in subsoils (Fig. 6) included Chlamydiae (Chlamydiales), Bacteroidetes (Sphingobacteriales, Cytophagales), Acidobacteria (Edaphobacter, Candidatus Solibacter, Candidatus Koribacter, Gp5, Gp16, Gp1, Terriglobus), Gammaproteobacteria (Chromatiales), Firmicutes (Clostridiales, Bacillales), Gemmatimonadetes (Gemmatimonadales), Verrucomicrobia (subdivision 3, Spartobacteria), Chloroflexi, Latescibacteria, and Actinobacteria (Ktedonobacterales, Actinomycetales, Gaiellales). The fungal phyla Mucoromycota, Rozellomycota, Mortierellomycota, and Glomeromycota exhibited large negative log-fold change values in subsoils, with members of the classes Tremellomycetes (e.g. Tremellales, Filobasidiales, Trichosporonales), Umbelopsidomycetes (Umbelopsidales), Agaricomycetes (Thelephorales, Boletales), Leotiomycetes (Helotiales, Thelebolales), Sordariomycetes (Xylariales, Chaetosphaeriales, Sordariales), and Microbotryomycetes exhibiting strong negative responses to increasing soil depth (Fig. 6). More detailed outputs of the ANCOM-BC analysis are presented in Tables S8–S13.

When correlating the abundance of microbial taxa to soil chemical properties, the microbial phyla Bacteriodetes, Chlamydiae, Chloroflexi, Gemmatimonadetes, Mortierellomycota, and Mucoromycota were strongly correlated (Spearman ρ > 0.7, p < 0.05) with ¹⁴C abundance (Table S14). Furthermore, of all the soil chemical properties measured, radiocarbon ¹⁴C abundance exhibited the greatest number of strong and significant correlations with the abundances of taxonomic groups (Table S14).

3.5 Structural changes in the soil microbiome with depth

The average node degree and neighbourhood connectivity of the topsoil interkingdom network was significantly higher than that of the subsoil network. In contrast, the average shortest path length and clustering coefficient were significantly higher in the subsoil interkingdom network (Fig. S12; Table S15). Compared with subsoils, the topsoil interkingdom network had a more connected structure with a greater degree of co-occurrences between genera. For both topsoil and subsoil interkingdom networks, bacterial genera were the dominant community members compared with fungi and archaea (Fig. 7). Networks of fungal-only topsoil communities had a significantly (p < 0.05) higher node degree, neighbourhood connectivity, average shortest path length, and betweenness centrality than fungal-only subsoil networks (Fig. S12; Table S15). For prokaryote-only networks (i.e. bacteria and archaea), the average node degree and average shortest path length was significantly higher in the topsoil than in the subsoil network (p < 0.05; Table S15). The fungal-only subsoil network exhibited a sparse structure, with 80.65 % of nodes being single nodes, i.e. they shared no co-occurrence with any other node (versus 30.34 % in the fungal topsoil network). In prokaryote-only networks, single nodes compromised 7.87 % of nodes in subsoils and 1.18 % of nodes in topsoils. In all of the networks, there was a greater proportion of positive edges (average 70.41 %) than negative edges (29.59 %) connecting the nodes (Table S15).

Figure 7 SPIEC-EASI networks that were constructed on core microbial genera associated with topsoil (0–30 cm depth) and subsoil (30–100 cm depth) layers. Node size represents the normalised (centred log ratio) mean abundance of each genus. Blue nodes represent bacterial genera, red nodes represent archaeal genera, and yellow nodes represent fungal genera. Positive edge weights are shown by black lines connecting nodes and negative edge weights are shown by red lines connecting nodes.

4 Discussion

4.1 Changes in soil C quantity, stabilisation, and age with depth

Previous research has identified large stores of unaccounted C stocks located in subsoil layers below 30 cm (Balesdent et al., 2018; Gonzalez et al., 2018; Ross et al., 2020). Our research identified that within Puruki Forest, 35 % of soil C stocks were allocated to the subsoil layers below 30 cm. Within our sampling site, subsoil C storage patterns were highly variable and ranged from 19 % to 46 % of the total soil C stocks. Previously, Oliver et al. (2004) remarked upon the large variability in soil C storage patterns across Puruki Forest, reporting average C stocks of 143 Mg ha⁻¹ to 1 m soil depth; values which are higher than those observed in our study (99.71 Mg ha⁻¹). High variability in subsoil C storage within Puruki Forest may be a signature of the natural textural and structural variation observed in the pumice soils (Silver et al., 2000; Telles et al., 2003). Furthermore, this variability may be a signature of Puruki's previous land use history (native forest > pasture > plantation) and the harvesting disturbances associated with historical management practices (Beets et al., 2002; Garrett et al., 2021; Rijkse and Bell, 1974). It is important to consider that our study was performed across a small spatial scale, with our 10 soil cores extracted from one transect within Puruki Forest. Given the high variability in soil C values observed in our data, we recommend that future studies perform soil sampling at a higher replication level, with soil cores extracted over a larger spatial extent within the study area under investigation.

We observed an overall enrichment of soil δ¹³C with depth which aligns with previous research (Garten, 2011; Wang et al., 2018). However, this was not a clear linear enrichment and δ¹³C was highly variable between depths of 70–90 cm. This corresponded to the high structural variability in our soil cores, as there were layers of coarse pumice material (Taupo lapilli) present throughout the subsoil horizons below 50 cm (Froggatt, 1981; Rijkse and Bell, 1974). Wynn et al. (2005) identified soil texture as an important factor governing soil δ¹³C enrichment, with finer-textured soils being associated with a greater degree of δ¹³C enrichment with soil depth than coarse-textured soils. As discussed prior in reference to subsoil C variability, the variability in soil δ¹³C with depth may be a consequence of previous disturbances from land management practices which caused a mixing of topsoil and subsoil layers (Oliver et al., 2004). Down to 1 m soil depth, ¹⁴C abundance declined by 229 ‰, which is lower than the median global decline of 502 ‰ (Mathieu et al., 2015). However, the mean radiocarbon age of our deepest soil layers (∼ 1570 yr BP) corresponds well to the age of the soil parent materials. Soil parent materials were dominated by the Taupo eruption and this can be directly observed by the layers of Taupo lapilli on Waimihia, Rotoma, Waiohau, and Rotorua ash beds present in several of our subsoils. As these layers of Taupo lapilli have been dated as 1760 ± 80 yr BP before 1950 (Rijkse and Bell, 1974), we would not have expected any soil radiocarbon ages to have exceeded this value.

4.2 Soil microbial diversity and abundance changed with soil depth

Although microbial DNA abundance declined with depth, accumulatively subsoils held 49 % of bacterial and 33 % of fungal DNA abundance which corresponded to values for subsoil C stocks. Soil microorganisms are responsible for the mineralisation of SOC and are influential in determining SOC formation, composition, and persistence (Domeignoz-Horta et al., 2021; Jansson and Hofmockel, 2020). Although subsoil C is considered more resistant to microbial decomposition (Schmidt et al., 2011), an accumulation of microbial “stocks” in proximity to subsoil C stores presents concerns, particularly given our uncertainty on the effects of climate change on soil microbial activity and C substrate use (Yang et al., 2021). However, the contribution of microbial necromass C to total soil C is known to increase with soil depth (Ni et al., 2020). As our research does not provide information on the viability or activity of subsoil microbial DNA, much of the accumulated bacterial and fungal DNA in the subsoil may be non-viable microbial cells. Much like incorporating measures of soil C age and stability when quantifying subsoil C stocks, further research should measure the abundance of subsoil microbial cells which are viable and active, as well as those which are dormant but have the potential to become active under global change. Doing so will help us better investigate the vulnerability of subsoil C stocks to microbial decomposition under climate change.

Soil microbial diversity and abundance declined sharply with depth which fits in agreement with previous research findings (Blume et al., 2002; Eilers et al., 2012; Mundra et al., 2021; Rosling et al., 2003). These declines occurred alongside clear shifts in microbial community composition. The magnitude of the depth-driven changes slowed past 30 cm, which marked a clear transitional change between the topsoil–subsoil boundary. Although fungal diversity and DNA abundance were generally much lower than that of bacteria, patterns in response to depth were consistent across both microbial kingdoms. The changes in microbial diversity and abundance were highly correlated with the concentration (total C %) and radiocarbon age of soil C. These findings support research that has proposed declines in soil C quantity, quality, and availability to be strong drivers of declines in microbial diversity and abundance (Eilers et al., 2012). Due to their proximity to aboveground vegetation, topsoil layers typically receive a wider load and range of fresh C inputs from surface litter and plant roots (Fierer et al., 2003; Spohn et al., 2016). This may explain why soil depths 0–20 cm supported the greatest number of unique ASVs (i.e. those which occurred only in that soil layer), with the number of unique ASVs declining sharply with depth.

4.3 Shifts in community co-occurrence patterns with depth

Whilst bacterial diversity and abundance declined sharply with depth, bacterial subsoil communities exhibited relatively minor differences in network structure compared with topsoil communities. Furthermore, subsoil fungal networks exhibited greater structural differences between topsoil and subsoil layers compared with bacterial communities. Subsoil fungal networks were almost entirely disconnected and there was no suite of co-occurring fungal genera that characterised the subsoil microbiome. The breakdown in subsoil fungal network structure may be explained simplistically by reduced community size limiting the ability for fungal co-occurrences to establish. Although still in low abundance, the relatively larger bacterial community abundance may have contributed to the lesser effects of soil depth on their co-occurrence patterns compared with fungi. Furthermore, as the properties of soil fungal communities are strongly driven by characteristics of aboveground vegetation (Likulunga et al., 2021; Urbanová et al., 2015), the breakdown in subsoil fungal co-occurrences may be due to their increasing physical distance from the surface litter and plant rhizosphere.

4.4 Differential responses of microbial taxa between topsoil and subsoils

Consistent with previous research, we observed an increased abundance of archaeal taxa in subsoils, mainly the phyla Euryarchaeota, Crenarchaeota, and Woesearchaeota (Brewer et al., 2019; Eilers et al., 2012; Frey et al., 2021; Hartmann et al., 2009). Proposed reasons for the increased abundance of archaea with depth include their adaptation to chronic energy stress (Valentine, 2007), involvement as ammonia oxidizers in driving autotrophic nitrification in deep soil layers (Eilers et al., 2012), preference as methanogens for anaerobic environments (Feng et al., 2019), and adaptation to nutrient-poor environments as slow growing oligotrophs (Turner et al., 2017).

Several representatives of the bacterial phyla Acidobacteria and Firmicutes increased in abundance with depth which may be explained by their tolerance for nutrient-limited environments (Bai et al., 2017; Brewer et al., 2019; Frey et al., 2021; Feng et al., 2019; Lladó et al., 2017). However, there were several representatives of Acidobacteria and Firmicutes that also declined with depth, such as Candidatus Solibacter, Candidatus Koribacter, and Clostridiales. Hansel et al. (2008) reported high variability in the distribution of Acidobacteria subdivisions among different soil horizons. The broad phylogenetic coverage and ecological diversity of Acidobacteria can make it difficult to infer metabolic and functional adaptations based on their taxonomic position alone (Hansel et al., 2008). Such issues highlight the need to adopt more trait-based approaches when investigating the microbial taxa that exhibit adaptive responses to soil depth. For example, the large genome size of Candidatus Solibacter (Challacombe et al., 2011) may explain its high abundance in surface soils, which have a wide range of C and nutrient inputs and greater exposure to environmental flux (Barberán et al., 2014). Theoretically, soil depth may drive a reduction in microorganisms with large genome sizes, as specialisation becomes more favourable under harsh conditions and the maintenance of a larger genome becomes too energetically costly.

Strong declines in the abundance of Bacteroidetes with depth have been previously attributed to their copiotrophic behaviour (Eilers et al., 2012; Feng et al., 2019; Mundra et al., 2021). In our study, the abundance of Bacteroidetes was positively correlated with ¹⁴C abundance. This may indicate the phylum's preference for fast-cycling, bioavailable C which typically declines with depth as more chemically resistant, slow-cycling C forms a greater proportion of the C pool (Koarashi et al., 2012; Torn et al., 1997). Previous research has observed an increased abundance of Gemmatimonadetes with depth due to their preference for low soil moisture (Bai et al., 2017; Chu et al., 2016; Frey et al., 2021; Li et al., 2020; Seuradge et al., 2016; Zhang et al., 2019; Debruyn et al., 2011). This is inconsistent with our findings and may be explained by the highly permeable soils within Puruki Forest that exhibit a greater potential for soil moisture to increase with depth, with the subsoil horizons having a high moisture capacity (Beets and Brownlie, 1987; Will and Stone, 1967; Beets and Beets, 2020). Carteron et al. (2021) observed large shifts in the ratio of saprophytic : mycorrhizal abundance with soil depth, which correlated to changes in soil chemistry. In our study, soil depth was associated with large declines in fungal taxa (i.e. Mucoromycota, Rozellomycota, Mortierellomycota, and Glomeromycota) previously reported in the literature as saprophytic and rhizospheric fungi (Bonfante and Venice, 2020; Hoffmann et al., 2011; Wagner et al., 2013; Yurkov et al., 2016).

5 Conclusions

The findings of our research address knowledge gaps related to subsoil C dynamics for production forests such Puruki Forest. A considerable proportion (at least 35 %) of soil C stocks within Puruki Forest were allocated to the subsoil layers below 30 cm, with this value reaching up to 49 % in some soil cores. In the deepest soil layers, the mean C age was 1571 yr BP, which corresponded closely to the estimated age of the soil parent materials. Although there were large declines in bacterial and fungal DNA abundance and diversity with soil depth, the magnitude of these changes slowed past 30 cm depth. This marked a topsoil–subsoil transitional boundary in the soil profile. In total, 49 % of bacterial and 33 % of fungal DNA abundance was allocated to the subsoil. Quantifying the viable activity of this microbial “biomass” is essential for predicting the vulnerability of subsoil C stocks to microbial decomposition under climate change. Numerous archaeal taxa exhibited an increased abundance in subsoils, indicating their potential ecological importance in deep soil environments. Furthermore, numerous fungal and bacterial taxa exhibited significant differential abundances between topsoil and subsoil layers. This supports the notion that subsoil environments may act as an environmental filter, with the vertical distribution of soil microorganisms reflecting the sharp changes in soil physical and chemical properties with depth.