Soil depth as a driver of microbial and carbon dynamics in a planted forest

. 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 organic C (SOC) dynamics and stabilisation is important to maximise the capacity and longevity of C sequestration. Compared to surface soil layers, little is known about the SOC 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 SOC 15 reserves to climate change. This study aimed to dive deep into the subsoils of Puruki Experimental Forest (New Zealand) and characterise the changes in SOC dynamics and the soil microbiome down to 1 metre 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 (δ 13 C) and radioactive ( 14 C) C analyses were performed to assess depth driven changes in SOC stability and age. Our research identified large declines in microbial diversity and abundance with soil depth, alongside significant structural shifts in 20 community membership. Importantly, we conservatively estimate more than 35% of total C stocks are present in subsoil layers below 30 cm. Although C age steadily increased with depth, reaching a mean radiocarbon age of 1571 yBP (years before present) in the deepest soil layers, the stability of SOC 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 comprehensively characterised the abiotic and biotic properties of a subsoil environment - a 25 frequently understudied but significant component of forest ecosystems

C means they form a vital role in mitigating the effects of climate change. Understanding the processes that regulate forest soil organic C (SOC) dynamics and stabilisation is important to maximise the capacity and longevity of C sequestration. Compared to surface soil layers, little is known about the SOC 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 SOC 15 reserves to climate change. This study aimed to dive deep into the subsoils of Puruki Experimental Forest (New Zealand) and characterise the changes in SOC dynamics and the soil microbiome down to 1 metre 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 (δ 13 C) and radioactive ( 14 C) C analyses were performed to assess depth driven changes in SOC stability and age. Our research identified large declines in microbial diversity and abundance with soil depth, alongside significant structural shifts in 20 community membership. Importantly, we conservatively estimate more than 35% of total C stocks are present in subsoil layers below 30 cm. Although C age steadily increased with depth, reaching a mean radiocarbon age of 1571 yBP (years before present) in the deepest soil layers, the stability of SOC 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 comprehensively characterised the abiotic and biotic properties of a subsoil environment -a 25 frequently understudied but significant component of forest ecosystems.

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 the total aboveground terrestrial C stocks (Batjes, 1996;Jobbágy and Jackson, 2000;Six et al., 2002). Developing forest management strategies which maximise soil C sequestration will, therefore, serve as a powerful 30 tool for offsetting fossil fuel emissions and mitigating against climate change (Jandl et al., 2007;Mukul et al., 2020). measurements of 14 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.

Study site and soil collection 70
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 ( Figure A1). Puruki Experimental Forest (henceforth referred to as Puruki Forest), is a 35-hectare Pinus radiata forest which was established in 1973, after previously being converted from podocarp/hardwood native forest to rye grass 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 75 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 all vitally important for helping us understand drivers of subsoil C dynamics (Garret et al., 2021).
Puruki Forest sits at an elevation range of 500 to 700m a.s.l., 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 80 Kelliher, 1989). Puruki Forest has an average soil pH of 5.2 at 0 to 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 to 2 m depth (Beets et al., 2004).
Soil parent materials originated from Taupo volcanic centre (1850 ± 100 BP) and older ash showers from Taupo and Okataina volcanic centres (Rijkse and Bell, 1974). The area's Orthic Pumice soils are highly permeable, being composed of loamy sand, silty sand, and gravel (Taupo lapilli) (New Zealand Soil Classification System;Hewitt, 2010). 85 A permanent sample plot (PSP) within the Rua sub-catchment of Puruki Forest was targeted for sampling. As Puruki 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 metre (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 used for 90 DNA and chemistry analysis. Following extraction, core pairs were visually compared to check for any considerable differences in soil colour and texture by depth ( Figure A2). 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 (Supplementary Methods). 95

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 manufacturer's instructions. Soils required for measurements of total C & 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 100 inorganic P analysis 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 a constant weight. Subsamples of air-dried soils were fine ground using a Spex800D ball mill for measurements of total C (%), total N (%), δ 13 C [‰], and radiocarbon 14 C analysis [‰].

Stable isotope (δ 13 C) analysis
Fine ground soil samples were sent to GNS Stable Isotope Laboratory (Wellington, New Zealand) for δ 13 C analysis by isotope ratio mass spectrometry. Soil samples were pre-treated by adding 10% HCl to 1 g of each sample and left overnight at room 115 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 δ 13 C results are reported with respect to VPDB and N-Air, normalized to GNS internal standards which were cane sugar (-10.3‰), beet sugar (-24.6‰), and EDTA (-31.1‰). The stable carbon isotope composition ( 13 C/ 12 C) is reported in δ 13 C values presented as per thousand [per mille, ‰], with the analytical precision for these measurements being 0.2‰. 120

Radiocarbon 14 C analysis
Fine ground soil samples were graphitized 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 hours in evacuated, sealed quartz tubes with cupric oxide (CuO) and silver (Ag) wire to form CO2 gas. After combustion, CO2 was reduced to graphite through heating at 570°C with 125 hydrogen (H2) 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 δ 13 C values according to Stuiver and Polach (1977). Conventional Radiocarbon Age (CRA) and 14 C are reported as defined by Stuiver and Polach (1977), with CRA presented as years before present (yBP) and 14 C presented as per thousand [per mille, ‰].

Microbial amplicon sequencing 130
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 on DNA samples which initially failed to amplify using the non-barcoded primers 27F-1492R. Following nested PCR, amplicons were amplified using barcoded 515F-806R primers. 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 which initially failed to 135 amplify were targeted for nested PCR using primers NSA3-ITS4. Nested PCR products were PCR amplified using the

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 145 the qPCR reaction protocol and standard curve preparation are provided in the Supplementary Methods. Following quantification, 16S and ITS rRNA copy numbers were adjusted to copy number per ng DNA.

Statistical analysis
Prior to analysis, slope corrected C stocks (Mg ha -1 ) were calculated based on total C (%) and soil bulk density measurements.
Slope corrected stocks were calculated using the following formula: '(BD x SOCconc x D) x SF' (Gattinger et al., 2012;Jones 150 et al., 2008). Where 'BD' is bulk density (g cm -3 ), 'SOCconc ' is total C (%), 'D' is thickness of the soil layer (cm), and 'SF' is slope factor. The values for total C, total N, and C stocks are reported using the sum of <2 mm 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 -1 . The changes in soil pH, total C (Mg ha -1 ), total N (kg ha -1 ), C:N ratio, δ 13 C [%o], total P (kg ha -1 ), Bray P (kg ha -1 ), total bulk density (g cm -3 ), and exchangeable cations K, Ca, Mg, and Na (kg ha -1 ) with soil depth were tested for 155 significance using non-parametric Friedman's tests with soil core transect 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 -1 ) allocated to each depth increment was calculated relative to the sum of total C (Mg ha -1 ) for each sample core. Differences in 16S and ITS rRNA copy number with depth were tested for significance using Kruskal Wallis chi-squared, with post hoc testing performed using 160 Dunn's test with Bonferroni's adjustment. To estimate the allocation patterns of microbial DNA abundance down the soil profile, the proportion (%) of 16S and ITS rRNA copy numbers allocated to each depth increment were calculated relative to the total sum of the soil core.

Microbial diversity analysis
Base calling and demultiplexing of sequencing libraries to per-sample fastQ files was performed using Illumina's bcl2fastq 165 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 denoised 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 rRNA input reads were retained. 170 Following DADA2 processing, ASV tables were filtered to remove unidentified phyla, unwanted phyla (i.e. Cyanobacteria/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 Figures A3 and A4. Following the removal of samples with low ASV counts, ASV tables were rarefied to even sampling depth prior to analysis (16S rRNA: 11613 ASVs per sample, ITS rRNA: 10101 ASVs per sample). 175 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 180 multidimensional scaling (NMDS). Differences in community dissimilarity by depth were tested for significance using PERMANOVA, and pairwise differences tested using pairwiseAdonis R with Bonferroni adjustment (Arbizu, 2017). For both alpha and beta diversity analyses, the impacts of transect sampling position and nested vs single PCR 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 vs unique 185 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 identify microbial taxa which had a significant log-fold change between topsoil (0 to 30 cm) and subsoil layers (30 to 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 190 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 outlined 195 above.

Microbial network analysis
ASV tables were split into samples from topsoil (0 to 30 cm) and subsoil (30 to 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. those with a total relative abundance of at least 0.05% in a minimum of 3 samples) and then summarised to genus level. 200 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 Supplementary Methods. 205

Soil chemistry, C stocks, and microbial DNA abundance
Total C (Mg ha -1 ), total N (kg ha -1 ), total P (kg ha -1 ), C:N ratio, and exchangeable cations (K, Ca, Mg, and Na; kg ha -1 ) significantly declined with soil depth (p < 0.05). Soil pH and bulk density (g cm -3 ) significantly increased with soil depth (p < 0.05). Results of soil chemical analyses and the corresponding statistical tests are presented in Table A1 and A2. 210 The average soil C stock of the 1 m soil cores extracted from Puruki Forest was 99.71 ± 28.78 Mg ha -1 C. The average soil C stock allocated to the coarse (>2 mm) soil fraction was 1.88 ± 1.92 Mg ha -1 of C. The amount of soil C stock allocated to the coarse soil fraction did not vary by soil depth increment (ANOVA: F-value = 1.58; p > 0.05).
On average, 34.58 ± 10.25% of the soil C stock was allocated to the 30 to 100 cm (subsoil) layer (Figure 1; Table A3). The proportion of soil C stock allocated to the subsoil varied across the soil cores, ranging from 46.60% to 16.21%. On average, 215 33.45 ± 10.92% of fungal DNA and 49.11 ± 4.82% of bacterial DNA was allocated to the subsoil layers past 30cm depth ( Figure 1; Tables A3 and A4). Thus, the allocation patterns of soil C with depth corresponded well to the distribution of bacterial and fungal DNA with depth.

δ13C and 14 C natural abundance measurements
Overall, δ 13 C [‰] significantly increased by 1.58‰ (p < 0.05; Table A1) between depth increments of 0 to 10 and 90 to 100 220 cm. However, between soil depths of 70 to 90 cm the enrichment of 13 C was lower than the rest of the subsoil (Figure 2).
Values for 14 C [‰] and CRA [yBP] of SOC declined with depth ( Figure 2; Table A1). The age of SOC steadily increased with depth and reached a mean age of 1570.83 ± 316.28 yBP in the deepest soil layers of 90 to 100 cm.

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; Figure  225 3; Table A6 and A7). The same trends were evident for the predicted abundances of bacterial functional genes with depth (Table A7). Transect position had no relationship to bacterial alpha diversity, however fungal evenness did vary across sampling location (Table A7). The DNA abundance of bacteria and fungi significantly declined with soil depth (Figure 3; Table A7) but were unaffected by transect position. Most of the significant differences in alpha diversity and DNA abundance occurred between the upper 0 to 20 cm depth increment and the rest of the soil; few significant pairwise differences were 230 evident among soil layers sampled below 30 cm ( Figures A5 to A8). When tested using pairwise Spearman's rank correlations ( Figure 4), microbial alpha diversity and abundance were strongly positively correlated to total C, total N, total P, C:N ratio, and 14 C (rho > 0.77; p < 0.05). Soil pH was strongly negatively correlated to microbial alpha diversity and abundance (rho > -0.85; p < 0.05). δ 13 C was strongly negatively correlated to fungal richness (rho= -0.85, p < 0.05), and bulk density (g cm -3 ) was strongly negatively correlated to bacterial abundance (rho= -0.90, p < 0.05). 235 Bacterial (PERMANOVA: R 2 = 0.21, p < 0.001) and fungal (R 2 = 0.26, p < 0.001) community composition varied significantly between depth increments ( Figure 5). Pairwise adonis tests identified that the community composition of soils between 0 to 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 ( Figure A9). The composition of bacterial and fungal communities significantly differed by sampling transect position (bacteria: R 2 = 0.13, p < 0.001; fungi: R 2 = 0.04, p < 0.001). However, this 240 accounted for less variation in community composition than soil depth. The composition of bacterial functional genes was significantly different between soil depth increments (R 2 =0.34, p < 0.001; Figure A10).

Microbial taxonomic responses to soil depth 250
Only 9% of bacterial ASVs and 2% of fungal ASVs were present in all depth increments down the entire soil profile ( Figure   A11). Upper depth increments (0 to 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 255 layers.

Structural changes in the soil microbiome with depth
The average node degree and neighbourhood connectivity of topsoil interkingdom networks was significantly than that of subsoils. In contrast, the average shortest path length and clustering coefficient was significantly higher in subsoil interkingdom networks than topsoil interkingdom networks ( Figure A12; Table A15). Compared to subsoils, the topsoil interkingdom network had a more connected structure with a greater degree of co-occurrences between genera. For both topsoil and subsoil 285 interkingdom networks, bacterial genera were the dominant community members compared with fungi and archaea ( Figure   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 ( Figure A12; Table   A15). 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 subsoil network (p < 0.05; Table A15). The fungal-only subsoil network exhibited a 290 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 A15).

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 total C stocks. Previously, Oliver et al. (2004) remarked upon the large variability in soil C storage 300 patterns across Puruki Forest, reporting average C stocks of 143 Mg ha -1 to 1 m soil depth; values which are higher than observed in our study (99.71 Mg ha -1 ). 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 an artefact 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). 305 High variability in soil C stocks across a relatively short spatial distance reaffirms previous research that attributes soil-carbon variability to much of the uncertainty in our ability to quantify subsoil C stocks (Ross et al., 2020). Soils are a significant contributor to the total C sequestration capacity of Puruki Forest and can hold 38% of total planted forest C stocks at end of forest rotation (Garrett et al., 2021). It is important to consider that this study was performed across a small-spatial scale and soil samples were obtained from 10 soil cores extracted along one transect within Puruki Forest. Given the high variability in 310 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 δ 13 C with depth which aligns with previous research (Garten, 2011;Wang et al., 2018). However, this was not a clear linear enrichment and δ 13 C was highly variable between depths of 70 to 90cm. This corresponded to the high structural variability in our soil cores, as there were layers of coarse pumice material (Taupo lapilli) 315 present throughout the subsoil horizons below 50cm (Froggatt, 1981;Rijkse and Bell, 1974). Wynn et al. (2005) identified soil texture as an important factor governing soil δ 13 C enrichment, with finer textured soils being associated with a greater degree of δ 13 C enrichment with soil depth than coarse textured soils. As discussed, Puruki Forest was once formerly covered by native forest prior to being converted to pasture and then planted forest (Beets and Brownlie, 1987). Thus, variability in soil δ 13 C with depth may be a consequence of previous disturbances from land management practices which caused a mixing 320 of topsoil and subsoil layers (Oliver et al., 2004). Down to 1 m soil depth, 14 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 yBP) correspond 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 yBP before 1950 (Rijkse and 325 Bell, 1974), we would not have expected any soil radiocarbon ages to have exceeded this value.

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 the formation, composition and persistence of SOC (Domeignoz-Horta et al., 2021;330 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 . 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 subsoil microbial DNA, much of the accumulated bacterial and fungal 335 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 340 (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. 345 These findings support research which 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 to 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. 350

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 to topsoils. Compared to bacterial communities, subsoil fungal networks exhibited greater structural differences between the topsoil and subsoil. 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 355 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 to 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 due to their increasing physical distance from the surface litter and plant 360 rhizosphere.

Differential responses of microbial taxa to topsoil vs 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 365 (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 370 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 amongst 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 375 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 maintenance of a larger genome becomes too energetically costly. 380 Strong declines in the abundance of Bacteroidetes with depth has been previously reported and attributed to Bacteroidetes' copiotrophic behaviour (Eilers et al., 2012;Feng et al., 2019;Mundra et al., 2021). In our study, the abundance of Bacteroidetes positively correlated with 14 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). Another bacterial phylum which declined in abundance with depth, Chlamydiae, was positively 385 correlated to younger SOC ages. Chlamydiae has also been positively associated with lower pH, organic matter content and C:N ratio (Zhalnina et al., 2015;Ma et al., 2021), properties which commonly change with soil depth. The declines in Gemmatimonadetes and Chloroflexi with depth observed in our study is inconsistent with previous research (Bai et al., 2017;Chu et al., 2016;Frey et al., 2021;Li et al., 2020;Seuradge et al., 2016;Zhang et al., 2019). The increased abundance of Gemmatimonadetes' with depth has previously been attributed to their preference for low soil moisture conditions (Debruyn 390 et al., 2011;Frey et al., 2021). However, the highly permeable soils within Puruki Forest 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). Thus, the decline of Gemmatimonadetes with soil depth at Puruki Forest may be due to its preference for low soil moisture conditions. 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 395 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).

Conclusions
The findings of our research address knowledge gaps related to subsoil C dynamics for production forests such Puruki Forest. 400 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 yBP, corresponding closely to the estimated age of 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 topsoilsubsoil transitional boundary in the soil profile. In total, 49% of bacterial and 33% of fungal DNA abundance was allocated to 405 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 act as an environmental filter; with the vertical distribution of soil microorganisms reflecting the sharp 410 changes in soil physical and chemical properties with depth.

Code/Data availability
The R code and datasets generated during this current study are available from the corresponding author on reasonable request.     Hierarchical clustering was performed on Euclidean distances of bias-adjusted abundances using the average linkage method. The log fold change (LFC) value of each taxa 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 450 between topsoil and subsoil is displayed, with TRUE = p < 0.05 and FALSE = p > 0.05. Figure 7. SPIEC-EASI networks constructed on core microbial genera associated with topsoil (0 to 30cm depth) and subsoil (30 to 100cm depth) layers. Node size represents each genus' normalized (centred log ratio) mean abundance. Blue nodes represent bacterial genera, red nodes represent archaeal genera, and yellow nodes represent fungal genera. Positive edge weights are shown 455 by black lines connecting nodes, negative edge weights are shown by red lines connecting nodes.