2 MATERIALS AND METHODS
2.1 Experimental design and soil sample collections
Two of the experimental sites used for this study were located in Wisconsin (Oregon and Hancock) and three in Michigan (Lux Arbor, Lake City, and Escanaba), USA. The soils of Oregon, Lux Arbor, and Escanaba sites are Alfisols, and of Hancock and Lake City sites are Entisol and Spodosol, respectively. At each site a randomized complete block design experiment with 3 (Hancock) or 4 (the rest of the sites) replications has been established in 2013. Details on the research sites and soil descriptions have been reported by Kasmerchak and Schaetzl (2018) and Lee et al. (2023).
The two studied vegetation systems were: (1) non-fertilized monoculture switchgrass (Panicum virgatum L.; Cave-in-Rock variety); and (2) restored prairie, which consisted of 18 plant species of grasses (including switchgrass), forbs, and legumes. Soil sampling was conducted in 2020 (Oregon site) and 2019 (the other 4 sites). Two types of soil samples were collected from each replicated plot. First, three intact soil cores (5 cm in height and 5 cm Ø) were collected from 5 to 10 cm depth for X-ray computed micro-tomography (µCT) scanning. Then, the loose soil surrounding the cores was also collected for measurements of soil texture and mineralogy. All samples were stored at 4 °C until scanning and measurements.
2.2 Measurements of soil porosity, texture, and mineralogy
Soil porosity was calculated from bulk and particle densities of the collected samples (i.e., total porosity), and texture was determined using the hydrometer method (Gee and Or, 2002) in all replicated plots of all sites. Soil mineralogy composition was measured using X-ray powder diffraction (XRD) carried out at the Illinois State Geological Survey (Champaign, Illinois, USA). Because of the high costs of XRD analyses only three replicated samples from each system in each site were subjected to these measurements. Prior to XRD analyses, the samples were cleaned, dried in a vacuum oven, and ground to < 44 µm. Then, one subset of the prepared sample was powdered by the McCrone mill (MBP) (McCrone Accessories & Components, Westmont, IL, USA) for quantification of non-clay minerals (quartz, clay, K-feldspar, P-feldspar, calcite, dolomite, siderite, and pyrite/marcasite), and the other was fractionized into < 2 µm powder for clay minerals (smectite, mica, kaolinite, and chlorite). The prepared two types of powders were then spread on a glass slide and analyzed using a Siemens/Bruker D5000 X-ray Powder Diffraction instrument (Billerica, MA, USA). The JADE™ software was used to identify percentages of constituents in each powdered sample from XRD patterns.
2.3 Soil core scanning and image analysis
Pore structure assessments were performed via X-ray µCT. Soil cores were drained at -28 kPa using a 5-bar pressure plate extractor (Soilmoisture Equipment Corp., Goleta, CA) prior to μCT scanning to remove water from pores Ø >10 μm and increase the contrast between the solids and air on X-ray μCT images. Then, the cores were scanned using X-ray μCT machine (North Star Imaging, X3000, Rogers, MV, USA) at the Department of Horticulture facility, Michigan State University. The energy settings were 75 kV and 450 μA. The scanning resolution of 18 μm was achieved using the Subpix-mode of the scanner, combining four individual scans shifted half pixel in vertical and horizontal directions. Scanned images of 3014 projections were reconstructed by efX software (North Star, Rogers, MN, USA).
A schematic summary of the steps involved in the image processing for this study is outlined in Figure 1. The image pre-processing was conducted using ImageJ-Fiji software (Schindelin et al., 2012) to remove artifacts and noises. First, to exclude sampling artifacts near the soil core walls, images were centered and cropped into prisms (1500 × 1500 × 2240 pixels corresponded to 2.7 cm in length, 2.7 cm in width, and 4.1 cm height). Then, ‘Remove Background’ tool in Xlib/Beat plugin was used to remove shadowing effects from the images, followed by the removal of ring artifacts on the image’s polar domain using a stripe filter of the Xlib/Beat plugin. After that, a 3D non-local mean filter (σ = 0.1) implemented in scikit-image (Walt et al., 2014) was used to reduce the noise (Buades, Coll, & Morel, 2011; Darbon, Cunha, Chan, Osher, & Jensen, 2008). The pre-processing steps dropped the resolution of images from 18 to 36 μm.
Root residues, which we will refer to as particulate organic matter (POM), were segmented from the filtered images with Ilastik software, a machine learning-based tool (Berg et al., 2019). A random forest classifier was used on a multi-dimensional feature space of the filtered gray scale images. The classifier was trained using two cores from each combination of different vegetation and sites (20 of total 114 cores) and then applied on entire cores. The training dataset produced out-of-bag error rate estimates less than 1.8% in overall, and all segmented POM images were visually inspected to ensure the accuracy and integrity of the segmentation process. The outcome of POM segmentation was denoised by removing objects smaller than 4 voxels in diameter from the images.
Segmentation of the filtered grayscale images into pore and solid binary images was performed to identify the pores visible at the image resolution, referred further on as image-based pores. For each sample the segmentation threshold was estimated as an ensemble of six segmentation methods (i.e., Otsu, Triangle, Huang, IsoData, Li, and Moments). The global thresholds for the stack of images in each individual core, estimated using the six segmentation methods, were averaged and applied to that stack to separate the solid and air-filled voxels in the images using SimpleITK in Python (Beare, Lowekamp, & Yaniv, 2018; Lucas et al., 2022). Obtained images were used to compute pore size distributions using ‘Local Thickness’ tool, an approach based on the maximum inscribed sphere method (Hildebrand & Rüegsegger, 1997; Vogel, Weller, & Schlüter, 2010).
Biopores were identified as described by (Lucas et al., 2022; Lucas, Schlüter, Vogel, & Vetterlein, 2019b). Specifically, to employ tubular-shaped features of biopores in differentiating them from other irregularly shaped pores, we used the Tubeness plugin in ImageJ-Fiji for shape detection. As rising σ-values significantly increased the computational time, binary images were scaled down to 50% and 20% for Tubeness filtering with σ-values ranging from 1-4 and 2-30, for each scale respectively, with a step size of 1. Gaussian blurring was applied to the entire binary image with varying σ-values in order to efficiently identify biopores of various diameters. The resulting tubular channels were slightly smaller than the root channel itself due to the exclusion of rough surface on biopore walls. Thus, to better capture the actual width of biopores, 3D dilation steps were employed as a postprocessing measure. After combining all elongated objects, misclassified objects were removed (Phalempin, Lippold, Vetterlein, & Schlüter, 2021a). Proportions of biopores in the entire pore system were calculated. After that, proportions of biopores occupied by POM were computed by first calculating the volume of POM located in biopores and then dividing this volume by the entire volume of biopores.
Nine masks corresponding to interval regions nine distances away from the segmented POM (0-0.25, 0.25-0.5, 0.5-1.0, 1.0-1.5, …, 3.0-3.5, and 3.5-4.0 mm) were created using 3D distance transform in ImageJ-Fiji. Then, masks of interval regions were applied to the pore-solid segmented image and to the pore size distribution image of the entire sample to calculate the porosities and the size distributions individually for each interval region. Contributions of pores of different size classes to image-based porosity of the distance interval regions were expressed as pore fractions (%). We considered three pore size classes, namely 36-150 μm, 150-300 μm, and > 300 μm Ø. The 36 μm Ø corresponded to the smallest pore size that could be reliably detected on the studied images. Pores < 150 μm Ø are known to have especially high microbial activity and strongly contribute to the C processing (Kravchenko & Guber, 2017; Kravchenko et al., 2019; Strong et al., 2004), and pores < 300 μm Ø function as the secondary pathways for water and nutrient supplies to resident microorganisms (Franklin et al., 2021).
A Connectivity tool of BoneJ plugin in the ImageJ-Fiji was used for the pore connectivity calculations: first, Euler numbers (χ) were computed, and the numbers were divided by the total volume of corresponding regions (V ) (Odgaard & Gundersen, 1993; Vogel & Roth, 2001):
\(\chi_{V}=\ \frac{N-C+H}{V}\) (1)
Where N is the number of isolated objects, C is the number of redundant connections or loops, and H is the number of completely enclosed cavities, which are typically negligible in soil pore system (Lucas, Vetterlein, Vogel, & Schlüter, 2021; Vogel, 2002). The minimum size of the object was 2 × 2 × 2 voxel. Higher, e.g., positive, χV values calculated via Eq. (1) correspond to lower connectivity, while lower, e.g., negative, to the higher connectivity. To simplify the presentation of the connectivity data we report the results as negative values of χV , that is, the high values of -χV correspond to high connectivity while the low values to the low one.
Since pore connectivity can be affected by the volume of the soil in which it is calculated, we did not assess it at the same distance intervals as those that were used for the image-based porosity and pore size distributions described above. Instead, we only calculated it in immediate vicinity of the residue, i.e., the region of 0-0.25 mm away from the POM, and for the entire soil volume. The resultant two estimates of the connectivity were used for comparisons among the five experimental sites and plant systems.
2.4 Statistical analysis
The data were analyzed using SAS 9.4 (SAS Institute Inc., NC, USA) procedures of PROC MIXED and PROC GLIMMIX. Since we did not expect that 6-7 years of disparate vegetation influenced soil texture and mineralogy, the statistical models for texture and mineralogy characteristics included only the fixed effects of the experimental sites. For the other soil properties, the statistical models included the fixed effects of sites, plant systems, and their interaction. Statistical models for analyses of image-based porosity data at different distances from POM additionally included the same fixed effects as the soil properties and individually tested by distance intervals, as the interaction among the sites, plant systems, and distances was significant. All models included the random effects of experimental blocks nested within the sites and, when necessary, the random effects of cores nested within the blocks, plant systems, and sites. The latter were used as an error term for testing the plant system effects. The assumptions of normality and variance homogeneity were assessed using normal probability plots, plots of residuals vs. predicted values, and Levene’s tests for equal variances.
Additionally, we grouped the sites into two soil texture classes for comparing pore size distributions and connectivity between finer-textured soils and coarser-textured soils. The first group included Oregon, Lux Arbor, and Escanaba sites, the three soils with < 66 % sand content, and the second group consisted of Hancock and Lake City sites with > 82 % sand content (Table 1). Models for analysis of pore size distribution data within each distance interval and of connectivity data within 0-0.25 mm distance and entire image stack included the fixed effects of plant systems, soil groups, and their interaction, and random effects of experimental sites nested within soil groups, blocks nested within the sites, and cores nested within the blocks, plant systems, sites, and groups.
Linear relationships among soil texture variables, mineralogical and clay mineralogy variables, proportions of biopores and POM in biopores, and the distance-based porosities and connectivity were assessed using Pearson’s correlation coefficients.