*Article* **Soil Autotrophic Bacterial Community Structure and Carbon Utilization Are Regulated by Soil Disturbance—The Case of a 19-Year Field Study**

**Chang Liu 1,2, Junhong Xie 1,2, Zhuzhu Luo 1,3, Liqun Cai 1,3 and Lingling Li 1,2,\***


**Abstract:** The roles of bacterial communities in the health of soil microenvironments can be more adequately defined through longer-term soil management options. Carbon dioxide (CO2) fixation by autotrophic bacteria is a principal factor in soil carbon cycles. However, the information is limited to how conservation tillage practices alter soil physiochemical properties, autotrophic bacterial communities, and microbial catabolic diversity. In this study, we determined the changes in autotrophic bacterial communities and carbon substrate utilization in response to different soil management practices. A replicated field study was established in 2001, with the following soil treatments arranged in a randomized complete block: conventional tillage with crop residue removed (T), conventional tillage with residue incorporated into the soil (TS), no tillage with crop residue removed (NT), and no tillage with residue remaining on the soil surface (NTS). Soils were sampled in 2019 and microbial DNA was analyzed using high-throughput sequencing. After the 19-year (2001–2019) treatments, the soils with conservation tillage (NTS and NT) increased the soil's microbial biomass carbon by 13%, organic carbon by 5%, and total nitrogen by 16% compared to conventional tillage (T and TS). The NTS treatment increased the abundance of the *cbbL* gene by 53% in the soil compared with the other soil treatments. The *cbbL*-carrying bacterial community was mainly affiliated with the phylum *Proteobacteria* and *Actinobacteria*, accounting for 56–85% of the community. Retaining crop residue in the field (NTS and TS) enhanced community-level physiological profiles by 31% and carbon substrate utilization by 32% compared to those without residue retention (T and NT). The 19 years of soil management lead to the conclusion that minimal soil disturbance, coupled with crop residue retention, shaped autotrophic bacterial phylogenetics, modified soil physicochemical properties, and created a microenvironment that favored CO2-fixing activity and increased soil productivity.

**Keywords:** carbon fixation; *cbbL* gene abundance; *cbbL*-harboring bacteria; carbon substrate utilization; soil management

#### **1. Introduction**

Soil microbial communities serve as the reservoir of biological processes and play a critical role in converting organic carbon into plant nutrients [1,2]. It is estimated that approximately 80 to 90% of the soil biogeochemical processes are mediated by microbial activities, including organic carbon degradation [3]. An essential biosynthetic process in soil is carbon dioxide (CO2) fixation by autotrophic bacteria, a key process of carbon cycling [4] which is responsible for the sequestration of atmospheric CO2 into the soil [5]. In agroecosystems, autotrophic CO2-fixing bacteria are present in various types, including distinct physiological groups, such as cyanobacteria and purple photosynthetic and chemoautotrophic bacteria [6]. The groups of bacteria can incorporate 14C into their microbial biomass through a series of biochemical reactions [7,8]. Those reactions are related

**Citation:** Liu, C.; Xie, J.; Luo, Z.; Cai, L.; Li, L. Soil Autotrophic Bacterial Community Structure and Carbon Utilization Are Regulated by Soil Disturbance—The Case of a 19-Year Field Study. *Agriculture* **2022**, *12*, 1415. https://doi.org/10.3390/ agriculture12091415

Academic Editor: Luciano Beneduce

Received: 21 July 2022 Accepted: 20 August 2022 Published: 8 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to the community composition of autotrophic bacteria [9], influencing the intensity of CO2 fixation.

The Calvin–Benson–Bassham cycle (i.e., the CBB cycle) is the most important autotrophic CO2 fixation pathway [10,11] which influences the primary carbon source in an ecosystem [12]. Ribulose-1,5-bisphosphate carboxylase/oxygenase (known as RubisCO), a key enzyme in the CBB cycle, regulates the process of CO2 reduction and the oxygenolysis of ribulose-1,5-bisphosphate [10,11]. RubisCO exists in four forms (referred as I, II, III, and IV) with different structures, catalytic activities, and O2 sensitivities, of which form I's RubisCO mainly occurs in plant and photo- and chemo-autotrophic bacteria [13,14]. In the plant–soil–microbiome environment, autotrophic microbial diversity can be reflected by the functional biomarkers (i.e., the *cbbL* genes). These genes encode the large subunit of RubisCO form I that has a large scale of sequences for phylogenetic analysis [15,16].

Many anthropogenic activities, such as tillage and crop rotation, affect the composition and diversity of the microbial community in soil [17,18]. Soil tillage affects soil water conservation [19], fertility enrichment [20], and carbon sequestration [21,22]. Reduced or zero tillage (i.e., minimal to no disturbance to the soil profile) coupled with crop residue retention can increase the nutrient supply potential [23] and enhance carbon cycling [24] while decreasing the carbon footprint [25,26], which helps alleviate the challenge of global climate change [27]. Diversified crop rotation increases soil carbon [28,29], improves microbial diversity [30], and enhances the system's productivity and resiliency. However, little is known about the long-term impact of multiple years of soil tillage and cereal– legume alternate rotation on the function of autotrophic CO2-fixing bacteria. It is unclear how the nexus of autotrophic bacterial community composition—CO2 fixation—tillage practices impact soil microenvironments.

To elucidate those effects, we established the field experiment in 2001 on the western Loess Plateau of China. Different tillage treatments have since been applied to soil that has a cereal–legume crop rotation system. After the 19-year treatments (in 2019), soil was sampled from each plot in three replicates and analyzed for its biophysiochemical properties. The abundance and community composition of autotrophic CO2-fixing bacterial communities were analyzed using real-time quantitative PCR and high-throughput sequencing. Soil carbon substrate utilization and its relation to autotrophic bacterial community were also determined. In a previous paper [31], we reported the diversity of bacterial communities and the characterization of the phylogenic composition in relation to soil management. In the present paper, we present the findings on (1) the changes in soil physicochemical parameters, autotrophic bacterial community composition, and soil catabolic diversity in response to the 19-year tillage treatments, and (2) the possible mechanisms responsible for shifts in the soil's autotrophic bacterial communities in relation to tillage. These findings from the long-term field study bring us to a new level of understanding of the interaction among autotrophic bacterial community, carbon fixation, and soil management practices.

#### **2. Materials and Methods**

#### *2.1. Site Description*

A field experiment was started in 2001 at the Rainfed Agricultural Education Center (35◦28 N, 104◦44 E; 1971 m a.s.l.) of Gansu Agricultural University in China. The site is on the western Loess Plateau and has a temperate, semiarid, continental monsoon climate [32]. Its solar radiation is 5.68 KJ m−<sup>2</sup> and sunshine duration is 2476 h annually [33]. The mean temperature is 6.5 ◦C and the 10 ◦C-based accumulative temperature is 2339 ◦C. The average precipitation was 392 mm per year, with 54% occurring from June to September in most years. The experimental plots were on level terrain with dark loessal Calcaric Cambisol soil (FAO/UNESCO soil classification system, 1990). Prior to the start of the experiment in 2001, the site had a long history of the continuous cropping of flax (*Linum usitatissimum* L.) under conventional tillage practices [34].

#### *2.2. Experimental Design*

A randomized complete block design was used to accommodate the four soil management treatments in each of the three replicates. The description of the treatments and the implementation of the tillage practices and crop residue retention are summarized in Table 1. More details on the tillage and field operations of the treatments can be viewed in a previous publication [34]. In brief, all the treatments were implemented in both the wheat (*Triticum aestivum* L.) and field pea (*Pisum arvense* L.) phases of the rotation each year. For the TS and T treatments, the soil was plowed to a depth of 0.20 m in the fall and harrowed in the following spring's seeding time. The experiment included four tillage and crop residue management treatments per spring. In each of the three replicates, the following treatments were arranged in plot sizes of 80 m2 (4 m wide × 20 m long). The crops were managed according to the recommended agronomic practices. The crops were planted in late March to early April, varying slightly from one year to the other due to weather conditions. The seeding rate was 187.5 kg ha−<sup>1</sup> for wheat and 180 kg ha−<sup>1</sup> for pea, on average, varying slightly each year due to seed size and percent germination. The seeding rates were aimed to achieve a plant population of 250 plants per m−<sup>2</sup> for wheat and 80 plants per m−<sup>2</sup> for field pea. The wheat was fertilized with urea nitrogen (46% N) at 105 kg N ha−<sup>1</sup> and the pea at 20 kg N ha−1. Both crops received a phosphorus fertilizer (Ca(H2PO4)2.2H2O, with 16% P2O5) at 46 kg P2O5 ha−1. The fertilizers were broadcast at sowing. No potassium (K) was applied as the soil contained exchangeable K greater than 300 mg kg<sup>−</sup>1. No herbicides were applied to any crop because of the dry weather and little weed pressure. No irrigation was applied. At maturity (usually in mid-July to early August), the entire plot was harvested using a plot harvester and the biomass and yields were recorded.

**Table 1.** Treatment structure and the description of the tillage practices and crop residue implementation in the long-term field experiment started in 2001.


<sup>a</sup> The treatment descriptions and abbreviations will be used in the other Tables and Figures throughout the article.

#### *2.3. Soil Sampling and Physicochemical Property Measurement*

Bulk soils were sampled in June 2019 when pea plants were at the mid-flowering stage and wheat plants were at the late stage of flowering. During flowering, crop plants are most active, with roots interacting with soil microbes under semiarid environments [35,36]. Sampling at this stage can provide researchers with invaluable information about rhizospheres and their interactions with environmental factors [37]. From each plot, five soil samples were randomly taken from a 0 to 0.2 m depth using a 50 mm diameter iron soil corer, and the five samples were bulked for each plot. After passing through a 2 mm mesh sieve, two subsamples were taken from each plot: one was for DNA extraction and carbon substrate utilization assessment and the other was for soil property analysis.

The soil pH was determined with a pH meter (Mettler Toledo FE20, Shanghai, China) using a deionized suspension with 1:2.5 soil: water ratio (mass: volume). Soil organic carbon (SOC) was analyzed using the Walkley–Black wet oxidation method and total N (TN) was measured with the Kjeldahl method. Olsen phosphorus was examined using the colorimetric method with 0.5 M NaHCO3. A UV-1800 spectrophotometer (Mapada Instruments, Shanghai, China) was used to determine soil NH4 +-N and NO3 −-N. Soil microbial biomass carbon and nitrogen (SMBC and SMBN, respectively) were measured using the Chloroform method. In each plot, soil moisture was measured gravimetrically and converted into volumetric water content using soil bulk density, which was measured by sampling the intact core of the soil using known-volume metal rings, drying the soil, and weighing the dried soil.

#### *2.4. Soil Microbial Catabolic Diversity*

For each soil sample from each plot (i.e., three replicates), community-level physiological profiles (CLPP) were determined using Biolog Eco-PlatesTM (Biolog Inc., Hayward, CA, USA) which contained 31 widely used *C* substrates along with a control well. The 31 C sources included nine carboxylic acids, seven carbohydrates, three other substrates, six amino acids, and two amines/amides [38]. Each well of the Eco-plates was inoculated with 150 μL of soil suspensions before being incubated at 25 ◦C. The carbon substrate utilization based on well color development was measured with optical density (OD) values at a wavelength of 590 nm (color development plus turbidity) and 750 nm (turbidity only) every 24 h for 7 d using the Biolog method (Biolog Inc., Hayward, CA, USA). Average well-color development (AWCD) was used to assess the microbial metabolic activity, as follows [39]:

$$\text{AWCD} = \sum \frac{C - R}{31} \text{\AA}$$

where *C* is the OD value of each well and *R* is the absorbance reading of the control well. Substrate richness and evenness, defined as the Shannon–Weaver index (*H'*), were calculated as follows:

$$H' = -\sum\_{i=n} p\_i (
ln p\_i)\_{..}$$

where *pi* represents the ratio of color development of a certain well to the sum of the color development of all wells in a micro-plate.

#### *2.5. DNA Extraction and Quantification of cbbL Gene Abundance*

Soil microbial DNA was extracted from each approximately 0.5 g of fresh soil sample, in triplicate, using a DNeasy PowerSoil Kit (QIAGEN, Inc., Hilden, Germany ) following the manufacturer's protocol. The concentrations of the extracted DNA were determined with a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and the quality of the extracted DNA was determined using agarose gel electrophoresis.

The *cbbL* gene copy number was assessed using real-time quantitative PCR (qPCR) on a LightCycler 480 II real-time PCR system (Roche Diagnostics, Mannheim, Germany). The *cbbL* genes were amplified with the forward primer K2f (5 -ACCAYCAAGCCSAAGCTSGG-3 ) and the reverse primer V2r (5 -GCCTTCSAGCTTGCCSACCRC-3 ) (reference). The reaction mixture contained 1 <sup>μ</sup>L of DNA template, 5 <sup>μ</sup>L of 2 × LightCycler® 480 SYBR Green I Master, 0.2 μL of forward primer, 0.2 μL of reverse primer, and 3.6 μL of nucleasefree water, and it made up a final volume of 10 μL. The thermal protocol was 10 min of incubation at 95 ◦C, followed by 40 cycles of 10 s at 95 ◦C and 30 s at 60 ◦C. Melting curve analysis was performed to validate the specificity of the amplified PCR products. The standard curves were produced using triplicate 10-fold dilutions of plasmid DNA with inserted *cbbL* gene. PCR amplification efficiencies were 94% with an *R*<sup>2</sup> value of 0.99.

#### *2.6. High-Throughput Sequencing of the cbbL Gene*

The primer set (K2f-V2r) for *cbbL* gene real-time qPCR was used for the *cbbL* gene amplification for each sample, respectively. The PCR components contained 5 μL of Q5 reaction buffer (5×), 5 μL of Q5 High-Fidelity GC buffer (5×), 0.25 μL of Q5 High-Fidelity DNA polymerase (5 U/μL), 2 μL (2.5 mM) of dNTPs, 1 μL (10 uM) of each forward and reverse primer, 2 μL of DNA template, and 8.75 μL of ddH2O. The PCRs were performed with the following procedures: 2 min of initial denaturation at 98 ◦C, followed by 25 cycles of 15 s at 98 ◦C, 30 s for annealing at 55 ◦C, 30 s for elongation at 72 ◦C, and

a final extension at 72 ◦C for 5 min. PCR amplicons were extracted from 1% agarose gels and purified with Agencourt AMpure beads (Beckham Coulter, Indianapolis, IN, USA) and quantified using the PicoGreen dsDNA assay kit (Invitrogen, Carlsbad, CA, USA). High-throughput sequencing data were obtained through commercial laboratory services (Personal Biotechnology Co., Shanghai, China) and deposited in the NCBI database with accession PRJNA689959.

#### *2.7. Data Processing and Bioinformatics Analysis*

Raw fastq files were quality-filtered and merged using VSEARCH pipeline (version 2.13.4) on the platform of QIIME 2. The sequences with a quality score <20, having ambiguous bases, or those containing mononucleotides repeats of >8 bp were removed. The operational taxonomic units (OTUs) were clustered with a 97% sequence similarity cut-off using VSEARCH pipeline (v. 2.13.4). All representative reads in the OTUs were taxonomically classified and annotated based on the GenBank® nucleotide sequence database (http://ncbi.nlm.nih.gov; accessed on 15 March 2020).

#### *2.8. Statistical Analysis*

Data were analyzed using SPSS software (SPSS Inc., Chicago, IL, USA). Significant differences between treatments in soil properties, *cbbL* gene abundance, and catabolic diversity were determined at the probability of <0.05.

The QIIME2 package was employed to determine alpha diversity, and the community evenness and richness were expressed using the Shannon and Simpson indexes. The variation of *cbbL*-harboring bacterial community composition across the samples was determined using principal co-ordinates analysis (PCoA) based on the genus-level compositional profiles. Hierarchical clustering analysis was performed with the unweighted pair-group method of arithmetic means (UPGMA) based on Bray–Curtis distance matrices. The differences in the taxonomic and phylogenetic communities were compared by clustering them using the R "*vegan*" package (R Foundation for Statistical Computing, Vienna, Austria). The significant differences in microbial structure among treatments was determined by permutational multivariate analysis of variance (PERMANOVA) and analysis of similarities (ANOSIM).

Network analysis was performed to investigate the co-occurrence patterns among the *cbbL*-harboring bacterial communities. The OTUs with a relative abundance of higher than 0.01% were retained and the "*psych*" package of R version 4.0.2 was used to analyze the preprocessed data and calculate the spearman correlation coefficient matrix. Statistical correlations with a cut-off at an absolute *r* value of higher than 0.6 and a *p* value of lower than 0.05 were retained for further analysis. The software Gephi (version 0.9.2) was employed for network visualization and network properties measurement. To identify the keystone bacterial taxa driving community turnover under different treatments, the abundance of bacterial taxonomic OTUs was regressed against the four field treatments based on random forest classification analysis using the "*randomForest*" package of R version 4.0.2. Lists of biomarkers taxa were ranked to compare the importance of OTUs using 10-fold cross-validation, implemented using R '*rfcv*' function. The top 15 bacterial biomarkers were chosen according to the stabilized cross-validation error curve. The relative abundance of biomarker taxa was illustrated using the "*heatmap*" package of R version 4.0.2. The Spearman's correlation coefficients and the Mantel test with 999 permutations in the "*vegan*" package of R were used to assess the correlations of soil properties and the distribution of the dominant genera of *cbbL*-harboring bacteria based on the distance matrix.

#### **3. Results**

#### *3.1. Soil Physiochemical and Biological Properties*

The soil's biophysiochemical properties significantly differed between treatments (Table 2). The contents of SOC, TN, TK, SMBC, and SMBN were significantly greater for the NT and NTS treatments compared to the T and TS treatments (*p* < 0.05). Soil organic carbon was 3.2%, 13.4%, and 6.1% greater for the NT, NTS, and TS treatments, respectively, compared with the T treatment. The highest soil TN concentration (1.07 g kg−1) was obtained in the NTS soil. The TK for the NTS treatment was 10.9% (*p* < 0.05) greater than that for the NT, TS, and T treatments. The soil's bulk density was 4.6%, 9.2%, and 6.1% lower for the NT, NTS, and TS treatments, respectively, compared with the T treatment. The NT, NTS, and TS treatments increased SMBC by 12.6%, 31.8%, and 17.0%, respectively, compared to the T treatment. In addition, SMBN for the NTS treatment was greater (*p* < 0.05) than that for the T treatment (*p* < 0.05), with no difference among NT, TS, and T. There was no difference in soil pH, TP, NO3 −-N, NH4 +–N, soil temperature, and moisture among the treatments (*p* > 0.05).

**Table 2.** Biophysiochemical properties of the soils under different management practices.


<sup>a</sup> SOC, soil organic carbon; TN, total nitrogen; TP, total phosphorus; TK, total potassium; SC, saturation conductivity; SMBC, soil microbial biomass carbon; SMBN, soil microbial biomass nitrogen. <sup>b</sup> The treatment descriptions are summarized in Table 1. <sup>c</sup> Data (means ± SD, n = 3) in each row followed by different letters are significantly different at *p* < 0.05.

#### *3.2. cbbL Gene Abundance and the Relation to Soil Properties*

The soils under the NTS treatment had more copy numbers of *cbbL* genes than the soils under the other treatments, and the increment was 52% when NTS was compared with the T treatment (*p* < 0.01; Table 3). Correlation analysis across the treatments showed that the soil *cbbL* gene was positively correlated with soil moisture, saturation conductivity, SMBC, and SMBN, but it was negatively correlated with bulk density (*p* < 0.05; Table S1).

**Table 3.** The *cbbL* gene copy numbers of the soil samples under the different soil management practices.


<sup>a</sup> The treatment descriptions are summarized in Table 1. <sup>b</sup> Values in each column followed by different letters are significantly different at *p* < 0.05.

#### *3.3. Autotrophic Bacterial Community Diversity and OTUs Richness*

After quality control, a total of 468,535 filtered sequences with 32,250 to 43,646 sequences per sample remained (Table S2). The rarefaction curves were close to the saturation phase with the increase in sample size. There was sufficient sequencing depth and the OTUs

were representative of the overall bacterial community libraries. To compare soil carbon fixation and microbial community diversity, the same survey effort level of 30,637 sequences was randomly selected from each soil sample (95% of the smallest number of reads). In total, 4988 OTUs were observed across all samples, with the number of OTUs ranging from 503 to 1692, varying with the soil samples (Table S2). The four tillage and crop residue management treatments shared 451 OTUs, accounting for 9.04% of the core OTUs (Figure 1a). There were differences in the OTU-sharing between treatments: the NTS and NT treatments shared the highest number of OTUs (1107 OTUs, 22.19%), whereas the TS and T treatments shared the least number of OTUs (683 OTUs, 13.69%) and the NTS and T treatments shared 18.04% of the OTUs (900 OTUs). Moreover, the NTS treatment harbored the highest number of unique OTUs (1458 OTUs, 29.23%) and the T treatment harbored the lowest number of unique OTUs (308 OTUs, 6.17%).

**Figure 1.** OTUs richness and diversity of soil CO2-fixing bacteria communities across all the treatments. (**a**) Venn diagram containing the unique and shared and unique OTUs among the treatments. (**b**) Alpha diversity index including Chao1, Simpson, Shannon, and observed OTUs. The asterisks indicate significant differences among the treatments. (**c**) Principal coordinate analysis (PCoA) based on bacterial community composition. The treatment descriptions are summarized in Table 1. The numbers following the treatment name denote the sampling replications. For example, NT1, NT2, and NT3 mean the soil sampling, respectively, from replicate 1, 2, and 3 of the NT plots.

The alpha diversity of bacterial populations differed significantly among soil treatments, as shown by the Chao1, Simpson, Shannon, and observed OTUs indexes (*p* < 0.05; Figure 1b). The Shannon index for the NTS treatment was significantly higher compared to the T treatment and the diversity was ranked NTS > NT > TS > T. No differences in the Chao1 and observed OTUs indexes were observed among the four treatments (*p* > 0.05; Figure 1b). The PCoA results showed similarities and differences in CO2-fixing bacterial community structure across the different treatments (Figure 1c), where total variance explained by PCoA1 and PCoA2 was 35.5 and 13.6%, respectively. The bacterial communities of the soil under the NT and TS treatments were clustered closely, whereas the community compositions in the soil under the NTS and T treatments differed from those of the NT and TS treatments. Permutational multivariate analysis of variance (PERMANOVA) and analysis of similarities (ANOSIM) confirmed the significant differences in the compositions of the soil *cbbL*-carrying bacterial communities under different tillage practices and crop residue managements (*p* < 0.05; Table S3).

#### *3.4. Composition and Networks of Soil Autotrophic Bacteria Communities*

The composition of the bacterial communities was determined for soils from each of the three field replicates (e.g., T1, T2, and T3 for the T treatments in Figure 2 and Table S4). At the phylum level, the bacterial communities were dominated primarily by the phyla *Proteobacteria* (38.52–80.94%) and *Actinobacteria* (5.72–22.65%), both accounting for an average of 56–85% of the relative abundance of the bacterial communities (Figure 2a; Table S4). The phyla *Euryarchaeota* and *Cyanobacteria* had low relative abundances. The tillage treatment had a significant (*p* < 0.05) impact on the relative abundance of a bacterial community. The relative abundance of the phylum *Proteobacteria* was 34.2, 41.6, and 33.5% lower in the soil under NT, NTS, and TS, respectively, while the relative abundance of the phylum *Actinobacteria* was 158.4, 245.7, and 142.5% higher, respectively, than that measured in the soil under the T treatment.

At the class level, *alpha-Proteobacteria* and *beta-Proteobacteria* dominated the bacterial communities, containing 14.89–30.71% and 5.71–22.65% of the total *cbbL* gene sequences, respectively (Figure 2b; Table S4). The class *Rubrobacteria* was the most dominant within the phylum *Actinobacteria*, which accounted for 0.57–4.00% of the total relative abundance. At the genus level, the *cbbL*-carrying bacteria were mainly affiliated with *Bradyrhizobium*, *Variovorax*, *Gaiella*, and *Steroidobacter*, accounting for 20–35% of the identified communities (Figure 2c; Table S4). The tillage treatments had an impact on the community composition, as shown in the top 10 genera; compared to the T treatment, the NTS treatment increased the relative abundance of genus *Bradyrhizobium* by 59.3% and reduced the relative abundance of genus *Variovorax* by 60.7%. However, no significant differences in the taxonomic composition were observed among the NT, TS, and T treatments. The hierarchical cluster analysis (Figure 2c) had results similar to those obtained from PERMANOVA and ANOSIM (Table S4), where the bacterial communities formed two main clusters (T and TS vs. NT and NTS), and the largest differences were between the NTS and T treatments.

Co-occurrence network analyses were conducted to evaluate the different co-occurrence patterns of the *cbbL*-harboring bacterial communities under the different treatments (Figure 3). In total, there were 35 nodes and 71 links, with an average degree of 4.057 and a graph density of 0.119 in the autotrophic bacterial network (Table S5). The ratio of positive edges was higher than the ratio of negative edges in the *cbbL*-harboring bacterial communities, with the ratio of positive edges accounting for 85.92%. The genera *Methylocella* and *Advenella* were the most connected keystone genera in the *cbbL*-harboring bacterial communities, and both were from the phylum *Proteobacteria*, with significant correlations with other genera (Figure 3; Table S5).

**Figure 2.** Microbial community composition under the four tillage practices and crop residue management treatments. (**a**) Relative abundance (%) of the most abundant phylum across all soil samples. (**b**) Taxonomic differences of *cbbL*-harboring microbial community. The largest circles represent the phylum level and the inner circles represent class, order, family, and genus. (**c**) Hierarchical cluster analysis and relative abundance (%) of the top 10 *cbbL*-harboring bacterial genera. The treatment descriptions are summarized in Table 1. The numbers following the treatment name denote the sampling replications. For example, NT1, NT2, and NT3 mean the soil sampling, respectively, from the NT plots in each of the three replicates.

#### *3.5. Biomarkers in the Autotrophic Bacterial Communities*

To determine the key biomarkers discriminating the *cbbL*-harboring bacterial communities under different treatments, we performed a regression of the abundances of bacterial OTUs against the tillage treatments based on the random forests machine learning algorithm. The top 15 keystone OTUs were chosen as the respective bio-marker taxa because the cross-validation error curve was stabilized when using these OTUs for both bacterial communities (Figure 4). The biomarker taxa were mainly affiliated with the classes of *alpha-Proteobacteria* and *beta-Proteobacteria* in the *cbbL*-harboring bacterial communities (Table S6). The heatmaps in Figure 4 further illustrate the scaled abundances of these biomarker taxa in response to soil tillage management. It was found that the phylum *Proteobacteria* played more crucial roles than the other phyla across soil samples for the soil's *cbbL*-harboring bacterial communities.

**Figure 3.** Spearman correlation analysis showing the *cbbL*-harboring bacterial co-occurrence networks across the soil samples. Each connection line represents a statistically significant correlation (*p* < 0.05) with a Spearman correlation coefficient of *r* > 0.6. The nodes of the networks are colored at the phylum level and labeled at the genus level. The size of the node is proportional to the number of connections (degrees), and the thickness of the edge is proportional to the value of the Spearman's correlation coefficients. The blue edges stand for positive correlations between two bacterial nodes, while the red edges stand for negative correlations.

**Figure 4.** Bacterial taxonomic biomarkers in *cbbL*-carrying bacterial communities. The top 15 biomarker taxa were ranked in descending order of the importance to the accuracy of the models. The heatmap illustrates the variation in the abundance of the 15 predictive biomarker taxa (OTUs). The treatment descriptions are summarized in Table 1. The numbers following each treatment name in the heatmaps denote the sampling replications. For example, NT1, NT2, and NT3 mean the soil samples were taken from replicate 1, 2, and 3 of the NT plots, respectively.

#### *3.6. Soil Microbial Catabolic Diversity and Carbon Utilization Pattern*

The soil's microbial metabolic activity was expressed as the average well-color development (AWCD), an indicative measurement of the physiological profiles at the community level. The AWCD differed significantly between the treatments (Table S7). The NTS and TS

treatments increased the AWCD index by 21.6 and 19.1%, respectively, while the NT treatment reduced the AWCD index by 15.4% in comparison with the T treatment. However, no significant differences for the *H'* (Shannon–Wiener) index were observed among the four treatments.

The soil's microbial carbon utilization patterns based on six fundamental groups differed among the treatments (Figure 5a). Compared to the T treatment, NT and TS reduced carbon resource utilization based on amines/amides by 63.6 and 35.2%, respectively, while NTS and TS increased carboxylic acid utilization by 37.1 and 44.9%, respectively. The overall carbon utilization capabilities differed significantly among the six carbon substrate groups (Figure 5b). The carbon utilization for carbohydrates, carboxylic acids, and amino acids were significantly higher than that for the polymers, miscellaneous, and amines/amides groups. The highest carbon utilization was obtained in the carbohydrates group, with a mean value of 8.808, while the lowest carbon utilization was found in the amines/amides group, with a mean value of 0.815.

**Figure 5.** Carbon utilization patterns of soil bacterial communities under the four tillage and crop residue management treatments, with (**a**) the carbon utilization based on six carbon substrate groups, (**b**) the carbon utilization of six carbon substrate groups, and (**c**) the utilization capabilities of 31 carbon substrate resources among the treatments. The different lowercase letters in (**b**) indicate significant differences among the substrates at *p* < 0.05. The asterisks in (**c**) indicate significant differences among the tillage treatments. Data are means ± standard deviation (n = 3). The treatment descriptions are detailed in Table 1.

The carbon utilization capability of soil microbiomes, based on seven specifical substrates (L-Arginine, L-Phenylalanine, α-D-Lactose, N-Acetyl-D-Glucosamine, D-Mannitol, D-Malic Acid, and α-Cyclodextrin), differed significantly among the tillage treatments (Figure 5c). The utilization of L-Arginine substrate in the TS treatment was 73.9% higher (*p* < 0.05) than that in the T treatment. The NTS treatment increased the utilization of L-Phenylalanine substrate by 85.5%, compared with the T treatment (*p* < 0.05). The NT treatment reduced the utilization of α-D-Lactose and N-Acetyl-D-Glucosamine substrates by 28.9 and 34.4%, respectively, compared with the T treatment (*p* < 0.05). The NT treatment increased the utilization of D-Mannitol substrate by 18.67% (*p* < 0.05) compared with

the NTS, TS, and T treatments, with the latter three having a similar utilization. For the utilization of D-Malic Acid substrate, the NTS and TS treatments respectively had 38.8 and 45.6% greater (*p* < 0.05) utilization than the T treatment. The utilization of α-Cyclodextrin substrate was 68.9% greater in the NTS and 102.7% greater in the TS treatments compared to the T treatment, while the NT treatment reduced the α-Cyclodextrin substrate utilization by 64.8% compared with the T treatment (*p* < 0.05). Across the four tillage treatments, the mean carbon utilization capabilities of the soil microbiomes differed significantly among the 31 carbon substrates. D-Mannitol substrate had the highest carbon utilization, averaging 1.624, and 2-Hydroxy benzoic Acid substrate had the lowest utilization, averaging 0.023.

#### *3.7. Correlation among Physiochemical Traits, Catabolic Diversity, and the Autotrophic Bacterial Community*

There were significant (*p* < 0.05) correlations among the soil physicochemical properties and *cbbL*-harboring bacterial community compositions (Figure 6a). Spearman correlation analysis across the tillage treatments revealed that SOC was positively correlated with TN, saturation conductivity, SMBC, and SMBN, but it was negatively correlated with bulk density (Table S8). Total soil N was correlated positively with TK, saturation conductivity, and SMBN, while soil bulk density was negatively correlated with saturation conductivity, SMBC, and SMBN. Further, the Mantel test revealed that the structure of the autotrophic bacterial communities was highly correlated with SOC (*r* = 0.315, *p* < 0.05), TK (*r* = 0.328, *p* < 0.05), bulk density (*r* = 0.724, *p* < 0.05), total porosity (*r* = 0.415, *p* < 0.05), saturation conductivity (*r* = 0.349, *p* < 0.05), and SMBC (*r* =0.377, *p* < 0.05) (Figure 6a).

**Figure 6.** A complex of correlations (**a**) between the autotrophic bacterial compositions (Table S4) and environmental factors (including soil physiochemical properties) and (**b**) between autotrophic bacterial community and microbial catabolic diversity in response to soil tillage treatments. In the (**a**) analysis, the pairwise comparisons of environmental factors are shown with the color gradient denoting Spearman's correlation coefficients (range 1 to −1); Spearman's rank correlations (r) of key soil parameters; and Mantel tests for the correlation between environmental factors and the distribution of the autotrophic bacterial community compositions. The edge width in the connection lines denotes the Mantel's r statistic for the corresponding distance correlations, and the edge color represents the statistical significances based on 999 permutations. In the (**b**) analysis, the red lines represent carbon substrate utilization, the blue lines represent the bacterial genus-level taxonomy, and the different colored circles represent the soil samples from all replications (n = 3) of each treatment. The treatment descriptions are summarized in Table 1.

Redundancy analysis, showing the relationship between carbon source utilization and the structure of *cbbL*-harboring communities at the genus level, revealed that the community composition varied from 14.20 to 20.68% (Figure 6b). The carbon utilization of polymers, amino acids, carbohydrates, and carboxylic acid were clustered to the edge of the genus *Bradyrhizobium* and *Methyloversatilis*, and they were negatively correlated to the genus *Rhodopila*. Among the carbon substrates, polymers explained 16.6% of the variance (*F* = 2.0, *p* = 0.178) and amino acids explained 8.1% of the variance (*F* = 0.9, *p* = 0.414) across the soil samples, and they were the key indicators of carbon utilization capacity (with 999 permutations) (Table S9). The abundance and diversity of the genus *Gaiella* were closely related to miscellaneous. The 19 years of tillage treatments had an impact on bacterial community structure, with the genera *Bradyrhizobium* and *Gaiella* presenting highly in the soil under the NTS treatment, while the genus *Variovorax* was high in the soil with continuous T practices.

#### **4. Discussion**

#### *4.1. Diversity of Soil cbbL-Carrying Bacterial Communities Shaped by Crop Residue Retention*

The biomarker *cbbL* genes, which encode a large subunit of ribulose-1,5-bisphosphate (RUbP) carboxylase/oxygenase in the Calvin–Benson–Bassham cycle, are considered important indicators of CO2-fixing autotrophic microbes in soil. We found a significant shift in the phylogenetic diversity of *cbbL*-carrying bacteria in response to soil management, and the largest differences in diversity existed between the NTS and the T treatments. The NTS treatment reduced the relative abundance of *Proteobacteria* while increasing the abundance of *Actinobacteria* at the phylum level, and at the genus level, the NTS treatment increased the relative abundance of *Bradyrhizobium* while reducing the relative abundance of *Variovorax*. The mechanisms for the shift in the bacterial community composition due to different soil management treatments are unclear, but we suggest the following possible mechanisms: **(a)** the minimal physical disturbance to the soil, coupled with continuous crop residue input (in each growing season, year-after-year), provides cumulative feedback to the soil's microenvironment, forming a balanced C to N ratio over the years, which favors microbial activities. A high N availability in soil is found to promote *Proteobacteria* while prohibiting *Actinobacteria*. A high C availability in soil is found to stimulate *Proteobacteria* while prohibiting other microbes compared to low N or low C soils. **(b)** The 19 years of continuous input of crop residue as organic material in the NTS treatment increased C sources for energy, essential for the cell synthesis of organisms [17], leading to the increased abundance of *cbbL* genes encoding ribulose-1,5-bisphosphate carboxylase oxygenase for CO2 fixation. The decomposition of the input material increased the substrate sources, leading to the great shift in the taxonomic composition of *cbbL*-carrying bacterial communities. **(c)** The shift in the microbial community composition due to soil management led to the formation of a microenvironment favorable to the humification process of SOC, increasing the size of the organic C pools. Additionally, **(d)** the facultative CO2-fixing bacterial members in soil may grow heterotrophically with the addition of crop residue (besides autotrophic activities) [40,41]. Nevertheless, our results, in combination with the findings reported by others [42–44], bring us to a new level of understanding that long-term crop residue input with minimal disturbance can shape microbial community composition, benefiting CO2 fixation processes and soil C cycling.

#### *4.2. cbbL-Carrying Bacterial Catabolic Diversity in Relation to Soil Management*

Bacterial OTU richness and catabolic diversity were highest in the soil under no till with crop residue retention among the soil treatments evaluated, suggesting that soil management does alter the diversity of soil microbiota, provided that the soil is treated for a longer period of time (19 years in the present study). Other researchers have found that the catabolic diversity is reflected by the shift in microbial community composition accompanied by the different capacities of carbon substrate utilization [45]. However, the catabolic diversity found in bulk soils may not be detectable in rhizosphere soils [46], since the expression of genes for soil nutrient reduction varies with soil microenvironments. Further, we found that the 19 year soil treatments did not affect soil pH (Table 1). The *cbbL*-carrying taxa were mainly affiliated with the types of bacteria, regardless of the level of pH tolerance. The diversity of autotrophic bacterial communities was highest in soil under the NTS treatment and lowest in soil under the T treatment, mainly due to the continuous lack of disturbance and crop residue input that created an ideal microenvironment that favored microbial activities. Our results support the findings of some researchers that soil pH does not have a direct role in impacting bacterial community diversity [47,48], and they do not agree with some others [49,50]. We suggest that the relationship between soil pH and bacterial community diversity is likely driven by pH-related soil nutrients. For example, the heavy use of N fertilizers can reduce the diversity and abundance of *cbbL* gene-containing microorganisms [51], mainly due to soil N decreasing pH, which prohibits some beneficial microbial activities [17,52].

#### *4.3. Relationships among Physiochemical Properties, CO2 Fixation, and Bacterial Diversity*

Soil is a comprehensive ecosystem in which microbial communities have a strong correlation with soil management options, host plants, and edaphic factors [53,54], including nutrient availability, ion exchange capacity, and SOC turnover [55,56]. In our study, the 19 year soil management practices led to substantial changes in physiochemical and biological properties. No till coupled with crop residue retention significantly increased *cbbL* gene abundance, SOC, SMBC, and SMBN while decreasing bulk density compared to conventional tillage. Our results add new knowledge that in a no-till soil microenvironment, bulk density, SMBC, and SMBN impact the *cbbL* gene abundance, and that the increased *cbbL* gene abundance, in turn, increases CO2 fixation, affirming the positive role of *cbbL*-carrying bacteria in soil's C cycling. The *cbbL*-carrying facultative autotrophs played a more prominent role than the obligate autotrophs in soil CO2 fixation, as we identified both facultative and obligate bacteria through sequencing, and we found that the relative abundance of *Proteobacteria* at the phylum level was highest in all the soil treatments evaluated. *Proteobacteria* is a well-known facultative aerobic autotroph that assimilates CO2 via the Calvin–Benson cycle. However, the CO2-fixing capability of autotrophs may vary with other factors, such as plant species [57], crop fertilization [51], and soil amendments [58].

#### *4.4. Implication of Long-Term Field Studies*

Bacteria play an important role in soil C cycling, some of which serve as the decomposers of root exudates, plant litter, and other carbon sources which break down the organisms and convert organic energy into forms useful to the rest of the organisms in the food web. The conversion process releases nutrients and thus improves soil fertility. However, limited information is available regarding how the bacterial community composition and diversity are influenced by long-term soil management options, such as tillage practices and crop residue input. In the scientific literature, some aspects of the relationship between soil bacterial communities and anthropogenic activities are documented, for example, soil tillage practices shape microbial community composition [59,60], change microbial diversity [46,61], and provide feedback to soil physiochemical properties [62]. However, most published findings are primarily drawn from short-term (< 10 years) experiments, such that the conclusions are often inconclusive, inconsistent, or even controversial from one study to another. Running a large field experiment continuously, year after year, is undoubtedly costly, but the experience of conducting the 19-year field experiment has provided us with the high confidence for making concrete conclusions on the subject area.

The improved soil microenvironmental conditions due to no till and crop straw retention for 19 years has led to increased soil productivity. Pea yield ranged from 1030 to 1980 kg ha−<sup>1</sup> and aboveground plant biomass ranged from 4300 to 6400 kg ha−<sup>1</sup> (Figure S1). Spring wheat followed a similar trend of treatment effect as that obtained in field pea. The effect of soil management on crop yields in 2019 was similar to that averaged cross the 19 study years; the NTS, TS, and NT treatments increased (*p* < 0.01) yields by 32.5, 25.1, and 10.0%, on average, of wheat and pea, respectively, compared to the T treatment. The soils receiving residue retention treatments (i.e., NTS and TS) had significantly higher productivity (22.7% more) compared to the soil without crop residue retention (i.e., NT and T).

Knowledge of the complex nexuses among soil's physiochemical properties, microbiomes, and anthropogenic activities (such as tillage practices and crop residue input) is highly demanded in the development of strategies to improve the health of soil microenvironments. Our findings are well positioned to meet the needs. One of the important aspects in the assessment of the impact of long-term soil management practices is to quantify the relative weights or percent impact between the two subfactors—tillage practice and crop straw retention. There is a need to understand whether the combination of the two subfactors actually generates additional synergy between the two. For example, the NTS treatment in our study had the greatest benefits in terms of soil microbial diversity, C-cycling, and nutrient enrichment compared to the rest of the treatments; however, we were unable to differentiate the relational weight of the impact between tillage practice and straw retention. Did the benefits stem from the NTS treatment due to straw input (x% of the benefits), no tillage (y% of the benefits), or the combination of the two? Nevertheless, our long-term experiment has been running since 2001 at a site representative of the 1.56 million hectares of the China Loess Plateau. The results of the study can provide some guidelines for improving soil health and productivity in local regions, and possibly extended to areas with similar climatic conditions to the experimental site.

#### **5. Conclusions**

The results of the 19-year field experiment show that continuous disturbance to the soil profile elevated the selective pressure to the soil CO2-fixing bacteria community, disturbed microbial catabolic activity, and created hardship for the microenvironment. Conversely, no disturbance to the soil, coupled with residue-retaining on the soil surface or crop residue being incorporated into the soil, increased the soil's organic carbon and nutrients and lowered the soil's bulk density. The underlying micro-driven functional changes led to enhanced soil carbon utilization. We conclude that bacterial community structure, composition, and CO2-fixing capability are highly regulated by soil management practices, and that a minimal disturbance to the soil's microenvironment, coupled with the retention of crop residues in the soil, will improve bacteria-involved biological activities and increase nutrient cycling and soil productivity.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/agriculture12091415/s1, Table S1: Pearson correlation coefficients between *cbbL* gene abundance and soil physicochemical properties, Table S2: Detailed sequencing depth and OTUs number of the *cbbL*-harboring bacterial communities under different soil treatments, Table S3: PERMANOVA and ANOSIM results for the effects of tillage and crop residue management treatments on soil *cbbL*-harboring bacterial communities, Table S4: Relative abundance of soil *cbbL*-harboring bacterial taxonomic composition at different levels for all samples as affected by tillage and crop residue management treatments, Table S5: Co-occurrence network properties of soil *cbbL*-harboring bacterial communities under the four tillage and crop residue management treatments, Table S6: Lists of the top 15 bacterial biomarkers taxa of *cbbL*-harboring bacterial communities under the different treatments, Table S7: Soil microbial catabolic diversity response to the four tillage and crop residue management treatments, Table S8: Spearman correlation coefficients among the soil's physicochemical and biological properties, and Table S9: Redundancy analysis of the relationship between carbon source utilization and the *cbbL*-harboring bacterial community at the genus level and clustering of the soil samples. Figure S1: Annual biomass and grain yield of spring wheat and field pea grown under different soil management practices.

**Author Contributions:** Conceptualization, L.L.; validation, C.L.; resources, L.L.; writing—original draft preparation, C.L.; writing—review and editing, L.L.; supervision, L.L.; project administration, J.X., L.C. and Z.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** The project was financially supported by the National Natural Science Foundation of China (31801320 and 31761143004), the Department of Science and Technology of Gansu Province (20JR10RA545), and the Fostering Foundation for Excellent Ph.D. Dissertations of Gansu Agricultural University (YB2018002).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data and material produced in the present study are provided in this manuscript and the supplementary materials. All high-throughput sequencing data were deposited in the NCBI Sequence Read Archive (SRA) database under accession number PRJNA689959.

**Acknowledgments:** We appreciate the excellent technical assistance for the field work and laboratory tests provided by undergraduate and graduate students at the Gansu Agricultural University Rainfed Agricultural Experimental Station.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Spectroscopic Investigation on the Effects of Biochar and Soluble Phosphorus on Grass Clipping Vermicomposting**

**Etelvino Henrique Novotny 1,\*, Fabiano de Carvalho Balieiro 1, Ruben Auccaise 2, Vinícius de Melo Benites <sup>1</sup> and Heitor Luiz da Costa Coutinho 1,†**


**Abstract:** Seeking to evaluate the hypothesis that biochar optimises the composting and vermicomposting processes as well as their product quality, we carried out field and greenhouse experiments. Four grass clipping composting treatments (only grass, grass + single superphosphate (SSP), grass + biochar and grass + SSP + biochar) were evaluated. At the end of the maturation period (150 days), the composts were submitted to vermicomposting (Eisenia fetida earthworm) for an additional 90 days. Ordinary fine charcoal was selected due to its low cost (a by-product of charcoal production) and great availability; this is important since the obtained product presents low commercial value. A greater maturity of the organic matter (humification) was observed in the vermicompost treatments compared with the compost-only treatments. The addition of phosphate significantly reduced the pH (from 6.7 to 4.8), doubled the electrical conductivity and inhibited biological activity, resulting in less than 2% of the number of earthworms found in the treatment without phosphate. The addition of soluble phosphate inhibited the humification process, resulting in a less-stable compound with the preservation of labile structures, primarily cellulose. The P species found corroborate these findings because the pyrophosphate conversion from SSP in the absence of biochar may explain the strong acidification and increased electric conductivity. Biochar appears to prevent this conversion, thus mitigating the deleterious effects of SSP and favouring the formation of organic P species from SSP (78.5% of P in organic form with biochar compared to only 12.8% in the treatments without biochar). In short, biochar decreases pyrophosphate formation from SSP, avoiding acidification and salinity; therefore, biochar improves the whole composting and vermicomposting process and product quality. Vermicompost with SSP and biochar should be tested as a soil conditioner on account of its greater proportion of stabilized C and organic P.

**Keywords:** 13C nuclear magnetic resonance; 31P nuclear magnetic resonance; charcoal; *Eisenia fetida*; pyrogenic carbon

#### **1. Introduction**

Biochar is a C-rich product distinct from charcoal because it is produced for soil application purposes in order to improve its quality (organic matter stability, cation exchange capacity (CEC), soil fertility, water retention capacity, porosity and biological activity are commonly increased), prevent nutrient leaching, improve C storage or purify the soil from pollutants [1–3]. The uses of and studies regarding this pyrolyzed biomass for agriculture have increased in recent decades, in part due to recognition of the role of carbonised biomass in the high fertility of *terra preta de indios* soils, which are anthropic soils that had a high input of carbonised biomass during the pre-Colombian period [4–6]. Although a number of studies have investigated the effects of different pyrolyzed residues on plant

**Citation:** Novotny, E.H.; Balieiro, F.d.C.; Auccaise, R.; Benites, V.d.M.; Coutinho, H.L.d.C. Spectroscopic Investigation on the Effects of Biochar and Soluble Phosphorus on Grass Clipping Vermicomposting. *Agriculture* **2022**, *12*, 1011. https:// doi.org/10.3390/agriculture12071011

Academic Editor: Borbála Biró

Received: 25 April 2022 Accepted: 7 July 2022 Published: 13 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

productivity and chemical and physical soil properties [7–11], little is known about *biochar* interactions with soil microbiota and fauna or with the native microorganisms from the raw materials used in composting [8,12–14].

Composting is a technique in which microbial activity during organic material decay is optimised and the material is transformed into either a fertiliser or soil conditioner. Vermicomposting, on the other hand, is a process by which earthworms aid in recycling and increase the velocity of organic matter stabilisation through residue homogenisation. The degradation of cellulose and inoculation of other organisms are stimulated by these macrofauna representatives, and humification tends to increase in the presence of these organisms [15–17].

Although earthworms occur in most environments, these organisms are sensitive to soil management practices [18,19], trace metals (mainly Cd) [20,21] and carbon polyaromaticcontaining substances [21,22]; thus, they can be used as quality indicators for soil or different organic substrates [21,23,24]. Fertilisers can affect the abundance and density of earthworms in the soil, although limited studies have investigated the effect of soluble sources, including phosphorus (P), on these soil organisms. Sarathchandra et al. [25] performed field studies and did not observe changes in the populations of *Lumbricus rubellus* or *Aporrectodea caliginosa* in either the adult or young forms after three years of single superphosphate (SSP) and rock phosphate application in North Carolina. However, laboratory studies have shown that direct contact with SSP can be lethal to the Californian red earthworm (*Eisenia fetida*), which has led to studies recommending the monitoring of these animal populations after applications of inorganic and soluble sources of P to the field [26].

In order to find a way to minimise the external input (fertilisers and soil conditioner) needed to maintain the extensive grasslands at the Rio de Janeiro International Airport in Brazil, which produces approximately 2 Mg of grass clippings daily, and also to reduce the disposal cost of this low-density raw material and hence reduce the environmental impact of the airport, the potential use of grass clippings for compost was proposed [27,28]. Therefore, this study aimed to evaluate the impact of using soluble phosphate and biochar on the establishment and development of the Californian red earthworm *E. fetida* in compost and vermicompost that included grass clippings from the Rio de Janeiro International Airport as well as evaluate the final products using 13C and 31P solid-state nuclear magnetic resonance spectroscopy. We hypothesized that the introduction of P in the composting process would increase the nutritional status of the compost/vermicompost, since P is a limiting element to plant growth in tropical oxidized soils, and that biochar would improve the composting and vermicomposting processes as well as the soil conditions for grass growth in the airport field area (not tested).

#### **2. Materials and Methods**

#### *2.1. Composting Site and Origin of Raw Material*

The experiment was conducted in a covered hangar of the Rio de Janeiro/Galeão International Airport (−22.804114◦; −43.229841◦), which is located in Rio de Janeiro, Brazil. Composting was performed using grass clippings obtained by mowing the airport's islands, i.e., the unused grassy areas between taxiways, between runways or between a taxiway and a runway. The clippings were ground to homogenise the particle size. The chemical composition of the residue (mean ± standard deviation, n=20) was as follows [28]: C = 415 ± 2 g kg<sup>−</sup>1; N = 11 ± 2 g kg−1; K = 10.7 ± 0.7 g kg−1; Ca = 4.2 ± 0.7 g kg−1; Mg = 1.3 ± 0.1 g kg<sup>−</sup>1; P = 0.5 ± 0.1 g kg<sup>−</sup>1; and C/N ratio = 47 ± 7.5. At the beginning of the experiment, the average moisture and density of the grass clippings were 70–80% and ~80 kg m<sup>−</sup>3, respectively.

#### *2.2. Composting of Grass Clippings*

Four treatments arranged in a full factorial design (22: two levels of phosphate and of *biochar* in grass clipping composting) were tested: grass (G); grass + biochar (GB); grass + phosphate (GP); and grass + biochar + phosphate (GBP). The phosphate used was

single superphosphate (SSP) with 8.7% P, mainly in the form of monocalcium phosphate (). Previous composting was performed in 2.5 m3 piles. The charcoal used in the piles was the residue of vegetable charcoal production (*Eucalyptus* sp., conventional carbonisation at 450 ◦C for six hours). This residue was a fine-grained charcoal, with a grain size usually smaller than 5 mm. The used sample exhibited non-uniform and small grain sizes, with only 10% of the mass retained on a 4-mm sieve. The quantities of each input used to prepare the compost piles are presented in Table 1.


**Table 1.** Quantities (kg) of each material/input used in the composting piles (2.5 m3).

At the beginning of the process, the piles were manually turned every two days. After 30 days, the piles were turned every 10 days. After 120 days of composting, all of the piles were ground to 5 cm, and the piles were turned twice a week for 30 days.

After 150 days of composting, six sub-samples of each pile were collected and transferred to plastic boxes for vermicomposting.

#### *2.3. Vermicomposting*

Vermicomposting was performed for 90 days using a completely randomised experimental design with six replicates. Twenty adult *E. fetida* were added to 10 L of the different composts and placed in 20 L plastic boxes. *E. fetida* earthworms were chosen because they exhibit high prolificacy, precociousness, survival rates and adaptability to captivity [4]. The experimental design was a 2 × 4 full factorial design (with and without earthworms × 4 former treatments, see Section 2.2) with 6 replicates and a total of 48 worm boxes.

Compost moisture was kept constant with frequent wettings. After 90 days, the earthworms from each box were counted.

#### *2.4. Physical, Chemical and Spectroscopy Characterisation*

Sub-samples from each box were dried at 45 ◦C for 48 h and ground, and the C and N concentrations were determined using a PerkinElmer 2400 CHNS (Waltham, MA, USA) elemental analyser. The pH and electrical conductivity (EC) of the different composts were determined as follows: approximately 15 g of compost was weighed in lidded plastic containers and then 85 mL of distilled water was added. The samples were stirred using a horizontal agitator for 30 min at 40 rpm and left to stand for 20 min, after which the EC and pH were measured.

For spectroscopic characterisation, sub-samples of each treatment were freeze-dried and ground in liquid N2. Solid-state 13C and 31P nuclear magnetic resonance (NMR) spectra were obtained with a Varian Inova (11.74 T) spectrometer at 13C, 31P and 1H frequencies of 125.7, 202.5 and 500.0 MHz, respectively. For this, a T3NB HXY with a 4-mm probe was utilised to detect the 13C and 31P nuclei and the rotors were spun using dry air at 15 kHz. All experiments were carried out at room temperature.

Two NMR pulse procedures were applied: variable amplitude cross-polarisation (CP) for 13C and direct polarisation (DP) for 31P.

In the 13C CP-MAS experiment an optimised recycle delay (d1) of 500 ms was used; the 1H 90-pulse was set to 3 μs, the contact time value to 1 ms and the acquisition time to 15 ms. High-power two-pulse phase modulation (TPPM) 1H decoupling at 70 kHz was employed. The cross-polarisation time was selected after variable contact time experiments and the recycle delays were selected to be five times longer than the longest spin–lattice relaxation time (T1), as determined by inversion recovery experiments.

In the 31P DP-MAS experiment, the recycle delay was 10 min (after inversion-recovery experiments), the 31P 90-pulse was set to 3 μs and the acquisition time to 15 ms.

#### *2.5. Data Analysis*

The data from the vermicomposting experiment (pH, EC and number of earthworms) were subjected to a multivariate analysis of variance (MANOVA), while the data from composting and further vermicomposting were analysed using repeated measures MANOVA. The normality and homoscedasticity of the residuals were tested and, when necessary, appropriate data transformation (log for C and C/N ratio and square root for the number of earthworms) was applied. When a statistically significant effect was detected, the means were compared using Duncan's test at *p* < 0.05. The software Statistica 7.1 [29] was used for all statistical analyses.

The multivariate analysis of the spectral data was performed using the software Unscrambler 10.4 (CAMO, Norway). Principal component analysis (PCA) was carried out using the full 13C NMR spectra obtained, after area normalisation and mean-centring of the data. To aid in the analyses of the 31P NMR results, the multivariate curve resolution (MCR) procedure was carried out. The basic goals of MCR are: the determination of the number of components co-existing in the chemical system; the extraction of the pure spectra of the components (qualitative analysis); and the extraction of the concentration profiles of these components (quantitative analysis). This analysis was preceded by PCA to estimate the number of components in the mixture. After this, the rotation of the PC was calculated without orthogonality constraints (in this way, it would have infinite solutions). To solve this, new constraints were adopted (e.g., non-negative concentrations and/or non-negative spectra). In this way, when the goals of MCR are achieved, it is possible to unravel the "true" underlying sources of data variation, after which the results with physical meaning are easily interpretable [18].

#### **3. Results**

#### *3.1. Earthworm Reproduction with Different Substrates*

The application of only phosphate drastically inhibited earthworm reproduction (Figure 1a). However, the interaction between the main factors of phosphate and biochar was statistically significant (*<sup>p</sup>* = 2.6 × <sup>10</sup>−6), with this detrimental effect mitigated by biochar, resulting in significantly larger populations than the GP treatment but still lower populations than the treatments without phosphate (G and GB).

#### *3.2. Chemical and Spectroscopy Analysis of the Composts and Vermicomposts*

The addition of phosphate fertiliser increased the medium's acidity and the EC of the vermicompost; however, biochar decreased these effects significantly (Figure 1b,c). Biochar prevented the pronounced acidification caused by phosphate; however, the pH was still lower than that of the treatments without phosphate (significant interaction, *<sup>p</sup>* = 2.0 × <sup>10</sup><sup>−</sup>6). Meanwhile, for EC just the main factors (phosphate and biochar) were statistically significant (*<sup>p</sup>* = 5.5 × <sup>10</sup>−<sup>6</sup> and 9.9 × <sup>10</sup>−4, respectively), without significant interaction; therefore, the biochar decreased the EC, fully mitigating the deleterious effect of phosphate application on the vermicompost EC.

**Figure 1.** Square root of earthworm number after inoculation for 90 days in different composted substrates (**a**); pH (**b**) and electrical conductivity (EC) (**c**) of the vermicompost. G: grass clippings; G+B: grass + biochar; G + P: grass + phosphate; and G + B + P: grass + biochar + phosphate. Columns with the same lower-case letters (for earthworms, pH and biochar factor for EC) and upper-case letters (for P factor for EC) indicate no statistical difference at *p* < 0.05 using Duncan test. Vertical bars denote standard errors.

The introduction of biochar increased C concentrations in all substrates (Figure 2a), and no significant effect was observed for phosphate or vermicomposting. Concerning N concentration and C/N ratio, the interaction of phosphate x biochar was significant (*<sup>p</sup>* = 2.8 × <sup>10</sup>−<sup>4</sup> and 2.5 × <sup>10</sup><sup>−</sup>4, respectively), being that phosphate and biochar decreased N and increased the C/N ratio (Figure 2b), with the strongest effect being from biochar. In addition, further vermicomposting increased the N content and decreased the C/N ratio (Figure 2c), indicating further humification.

**Figure 2.** Biochar effect on C concentration (**a**), phosphate and biochar interaction effects on N concentration and on C/N ratio (**b**), and effects of vermicomposting on these variables (**c**). G: grass clippings; G + B: grass + biochar; G + P: grass + phosphate; and G + B + P: grass + biochar + phosphate; Comp: composting; and Vermic: further vermicomposting. Columns with the same lower-case letters (C and N) and symbols with the same upper-case letters (C/N ratio) indicate no statistical difference at *p* < 0.05 using Duncan test. Vertical bars denote standard errors.

The change in humification was confirmed by the 13C and 31P NMR spectra, which are shown below. The C/N ratio increased in all treatments with biochar (Figure 2b), which was also explained by the inclusion of material rich in recalcitrant C and poor in N. Meanwhile, a significant effect of vermicomposting we observed was that it decreased the C/N ratio. The addition of phosphate to the substrates in the absence of biochar resulted in higher C/N ratios in the compost and vermicompost (no interaction), confirming the inhibition of humification by phosphate.

The organic chemical structures detected using 13C NMR (Figure 3) indicate that, following composting, the material still exhibited significant amounts of labile structures, mainly cellulose (O-alkyl, with signals at ~65 and 73 ppm and di-O-alkyl, with signal at ~105 ppm). The relative contribution of these regions varied between substrates, indicating selective degradation of cellulose (higher biological activity and higher humification) in certain substrates. For the substrates with biochar, a pronounced increase of the C-aryl signal was observed (~125 ppm), which is a typical signal of condensed aromatic rings. These differences were summarised and highlighted using PCA and are discussed in detail below.

**Figure 3.** Solid state 13C NMR spectra of the different compost and vermicompost obtained. Without biochar (**A**) and with biochar (**B**).

#### *3.3. Progress of Humification According to Principal Component Analysis (PCA)*

Using PCA, we were able to satisfactorily model the data (97% of variance captured) with two principal components (PCs). The first PC (92% of total variance) identified the different substrates exhibiting bipolar loadings (Figure 4A), with positive loadings for the sp2 C signals of biochar (aryl) and negative for the sp<sup>3</sup> hybridised C of grass cellulose (Oalkyl and di-O-alkyl). The second PC (5% of total variance) was characterised by negative loadings for the signals typical of labile structures (Figure 4A), primarily cellulose, that were partially oxidised to uronic acids (O-alkyl, di-O-alkyl and carboxyl with a signal at ~173 ppm). Therefore, this PC served as an indicator of the progress of humification, since less labile structures indicate the progression of the humification process.

The samples with biochar (GB, VGB, GBP, VGBP) presented the highest scores for PC1 (Figure 4B), confirming the contribution of polycondensed aryl structures towards their 13C NMR spectra.

The samples containing only grass, with and without earthworms, exhibited the highest PC2 scores, indicating a lower proportion of labile structures and more advanced humification (Figure 4B). Treatments with phosphate and without biochar (compost—GP and vermicompost—VGP) exhibited the lowest scores for PC2 (Figure 4B), indicating that phosphate inhibited the advance of humification by suppressing the biological activity of micro- and macro-organisms (after the earthworm population decrease, Figure 1a), resulting in a material rich in labile structures, i.e., partially oxidised cellulose. This finding confirmed the results for the N concentrations and C/N ratios as well as the detrimental effects on earthworms. The samples with biochar (GB, VGB, GBP, VGBP) exhibited intermediate scores (Figure 4B) for this component (uronic acids), indicating that the maturation of the compost and vermicompost were similar for the samples with biochar, regardless of phosphate presence.

**Figure 4.** (**A**) PCA loadings obtained from the 13C NMR spectra. (**B**) PCA scores. G: grass compost; VG: grass vermicompost; GP: grass + phosphate compost; VGP: grass + phosphate vermicompost; GBP: grass + biochar + phosphate compost; VGBP: grass + biochar + phosphate vermicompost.

The 31P NMR and its multivariate curve resolution analysis (Figure 5) indicate that, in the absence of biochar, the addition of phosphate to the compost significantly changed the distribution of P species, with a likely predominance of pyrophosphate (with a chemical shift of approximately −1.3 ppm). Moreover, in samples without phosphate or with phosphate and biochar, the predominant form was likely organic (mono and di-ester phosphate, with chemical shifts of approximately 0.96 ppm and 2 ppm, respectively) or mixed with non-hydrolysed monobasic phosphate (chemical shift of 0.96 ppm).

**Figure 5.** Multivariate curve resolution results obtained from the 31P NMR spectra showing the two component solutions. (**A**): estimated component spectra; discontinued lines show the 31P NMR spectra of 31P standards. (**B**): estimated concentrations. G: grass compost; VG: grass vermicompost; GP: grass + phosphate compost; VGP: grass + phosphate vermicompost; GBP: grass + biochar + phosphate compost; and VGBP: grass + biochar + phosphate vermicompost.

#### **4. Discussion**

The grass clipping compost obtained from the International Airport of Rio de Janeiro did not exhibit limitations for composting, which was consistent with many other works [27,28]. We found that grass clippings could also be used for vermicomposting because a statistically significant increase of 400% of the earthworm population was observed during the studied period (Figure 1a). Curry and Schmidt [30] cited, amongst a wide range of organic materials, the presence of grass in the digestive tract of these soil animals, which indicates that grass clipping material can be used to feed earthworms and then in vermicomposting.

However, when SSP was added to the composting substrate, a detrimental effect was observed. The acid reaction from the sulphate present in the phosphate inhibited microbial activity, with the resulting compost exhibiting lower N concentrations, higher C/N ratios (Figure 2b) and higher levels of labile structures (uronic acids—Figure 4). Moreover, the phosphate severely inhibited earthworm reproduction and led to a significant decrease in the earthworm population (i.e., death or escape of the earthworms, Figure 1a). The phosphate dosage (1 kg m−3) used in the present experiment resulted in a significant increase in compost acidity (pH < 5) and in the ionic strength of the substrates (Figure 1b,c), which are factors that negatively affect the development and reproduction of earthworms of this genus [31]. This acidification may result from the conversion of applied phosphate into pyrophosphate in the absence of biochar, as detected by 31P NMR (Figure 5). On the other hand, in the absence of phosphate or in the presence of biochar, the predominant forms of P were organic or non-hydrolysed monobasic phosphate, indicating the highest biological activity and incorporation of P into biological tissues.

In addition, the presence of certain metals (Cd, Cu, Ni, As, Cr and Pb) and radionuclides (226Ra, 228Ra and 210Pb) has been reported in phosphate fertilisers commercialised in Brazil [32], albeit in reduced concentrations. However, their solubility may have increased because of the pronounced pH decrease, thus contributing to the detrimental effect on the earthworms. Nevertheless, additional specific studies are required to confirm this hypothesis. It should be noted that in a review on the uptake and accumulation of heavy metals by earthworms, the metal concentration in the soil was found to be a poor predictor of metal concentration in the earthworms' tissues, and pH was considered to be the main factor associated with metal uptake by earthworms [33].

The presence of biochar mitigated these detrimental effects of phosphate. The pyrolysed biomass acquired different properties depending on the pyrolysis temperature and nature of the feedstock [13,34]. Biochar application has been found to promote changes in nutrient availability (usually pH increase, N immobilisation etc.) [19]; interference with organism and plant signalling compounds; xenobiotic detoxification; and refuge availability for pathogens or growth-promoting organisms, etc. These changes are examples of multiple beneficial or detrimental effects (direct and indirect) of biochar on soil quality and plant production [24,35–37].

The litter earthworm is sensitive to high acidity, high salinity and certain toxic elements, and the treatments with biochar exhibited increased pH and decreased EC, similar or close to the values from the treatments without phosphate. These characteristics may have contributed to the mitigation of detrimental effects of phosphate on the earthworm population and also towards microbiological activity, as indicated by the humification of the biomass (Figure 4).

The 31P NMR spectra and their MCR analysis (Figure 5) show that, in the absence of biochar, the addition of phosphate to the compost significantly changed the distribution of P species, with a likely predominance of pyrophosphate. Moreover, in samples without phosphate or with phosphate and biochar, the predominant form was likely organic (mono and di-ester phosphate) or mixed with monobasic phosphate. This conversion of applied phosphate into pyrophosphate in the absence of biochar may explain the strong acidification of the medium as well as the increase in EC. Biochar was able to prevent this conversion and favoured the formation of P organic species, which may have been caused by the highest microbial activity in the presence of biochar as indicated by the highest compost humification (Figure 3).

Vermicomposting increased the N concentration of the tested substrates (Figure 2c), probably due to the labile organic matter evolving to CO2 and resulting in a relative enrichment in N.

Biochar addition to the substrates led to lower N concentrations in the compost and vermicompost (Figure 2b), which may have been caused by a dilution effect related to the inclusion of material with high recalcitrant C concentrations and low N concentrations.

A decrease in the C/N atomic ratios was observed after vermicomposting (Figure 2c), indicating a further humification process promoted by the earthworms. The relative increase of N concentrations likely occurred because the raw materials were N-poor and cellulose degradation (CO2 emission) was favoured by the presence of earthworms; we confirmed this using the 13C NMR spectra (Figures 3 and 4).

#### **5. Conclusions**

A detrimental effect of SSP on macro- and micro-biota was observed with a drastic decrease in earthworm population and lower compost and vermicompost humification.

This detrimental effect of SSP was likely caused by acidic conditions resulting from pyrophosphate formation and high salinity, inhibiting humification. This was indicated by the presence of materials with increased labile structures, low N concentrations and high C/N ratios. However, biochar mitigated this negative effect during composting by maintaining microbial activity, which was indicated by higher humification, and during vermicomposting, by mitigating the detrimental effect of SSP on the earthworm population. Biochar probably reduces the deleterious acidity and salinity induced by single superphosphate by decreasing pyrophosphate formation, avoiding acidification and salinity.

SSP alone is not recommended for the composting and vermicomposting of grass clippings, but in combination with biochar it should be tested as a soil conditioner, on account of its greater proportion of stabilized C (from grass clippings and pyrogen) and organic P.

**Author Contributions:** Conceptualisation (E.H.N.; F.d.C.B.; V.d.M.B. and H.L.d.C.C.); methodology (E.H.N.; F.d.C.B.; R.A. and V.d.M.B.); investigation (E.H.N.; F.d.C.B.; R.A. and V.d.M.B.); formal analysis (E.H.N.); visualisation (E.H.N.); writing—original draft, review and editing (E.H.N.; R.A.; F.d.C.B. and H.L.d.C.C.); supervision and funding acquisition (E.H.N.; F.d.C.B.; V.d.M.B. and H.L.d.C.C.); project administration (F.d.C.B.). All authors have read and agreed to the published version of the manuscript.

**Funding:** Funding was provided by: the Funding Authority for Studies and Projects (Financiadora de Estudos e Projetos—FINEP) grant project "Production of organic fertilisers and substrates from residues of grass maintenance in urban areas"; the National Counsel of Technological and Scientific Development (Conselho Nacional de Desenvolvimento Científico e Tecnológico—CNPq) for Novotny's and Balieiro's research fellowships (309391/2020-2 and 307434/2020-6, respectively); the Brazilian Federal Agency for Support and Evaluation of Graduate Education (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—CAPES) for Auccaise's post-doctoral research fellowships (PNPD Nº 03020/09-6); and Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) for Novotny's research fellowship (E-26/202.874/2018).

**Institutional Review Board Statement:** Not applicable since the work does not involve humans or animals.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** The authors wish to thank the Nuclear Magnetic Resonance Laboratory of the Brazilian Centre of Physics Research (Laboratório de Ressonânica Magnética Nuclear do Centro Brasileiro de Pesquisas Físicas). The authors wish to thank the Embrapa Soils trainees (scholarship granted - PIBIC/Embrapa Solos): João Paulo Moura Barata, Carolina Araújo de Queiroz Costa and Natália Evaristo Belo do Santos, as well as the airport workers who helped conduct the experiments at the airport and at Embrapa Soils. I thank Adriana Aquino from Embrapa Agrobiology for the donation of earthworms and recommendations for the management of the experiment.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article Trichoderma* **Bio-Fertilizer Decreased C Mineralization in Aggregates on the Southern North China Plain**

**Lixia Zhu, Mengmeng Cao, Chengchen Sang, Tingxuan Li, Yanjun Zhang, Yunxia Chang and Lili Li \***

College of Life Science and Agronomy, Zhoukou Normal University, Zhoukou 466001, China; 20181006@zknu.edu.cn (L.Z.); 202007030031@zknu.cn (M.C.); 202007030055@zknu.cn (C.S.); 202007030014@zknu.edu.cn (T.L.); 202007030047@zknu.cn (Y.Z.); 20061015@zknu.edu.cn (Y.C.) **\*** Correspondence: 19961011@zknu.edu.cn

**Abstract:** *Trichoderma* bio-fertilizer is widely used to improve soil fertility and carbon (C) sequestration, but the mechanism for increasing C accumulation remains unclear. In this study, effects of *Trichoderma* bio-fertilizer on the mineralization of aggregate-associated organic C were investigated in a field experiment with five treatments (bio-fertilizer substitute 0 (CF), 10% (BF10), 20% (BF20), 30% (BF30) and 50% (BF50) chemical fertilizer nitrogen (N)). Aggregate fractions collected by the dry sieving method were used to determine mineralization dynamics of aggregate-associated organic C. The microbial community across aggregate fractions was detected by the phospholipid fatty acid (PLFA) method. The results indicated that *Trichoderma* bio-fertilizer increased organic C stock across aggregate fractions and bulk soil compared with CF. Cumulative mineralization of aggregate-associated organic C increased with the increasing bio-fertilizer application rate. However, the proportion of organic mineralized C was lower in the BF20 treatment except for <0.053 mm aggregate. Moreover, the PLFAs and fungal PLFA/bacterial PLFA first increased and then decreased with increasing bio-fertilizer application rates. Compared with CF, the increases of bacteria PLFA in >2 mm aggregate were 79.7%, 130.0%, 141.0% and 148.5% in BF10, BF20, BF30 and BF50, respectively. Similarly, the PLFAs in 0.25–2, 0.053–0.25 and <0.053 mm aggregates showed a similar trend to that in >2 mm aggregate. Bio-fertilizer increased the value of fungi PLFA/bacteria PLFA but decreased G+ PLFA/G− PLFA, and BF20 shared the greatest changes. Therefore, appropriate *Trichoderma* bio-fertilizer application was beneficial to improving soil micro-environment and minimizing risks of soil degradation.

**Keywords:** *Trichoderma* bio-fertilizer; soil aggregate stability; mineralization; microbial community

#### **1. Introduction**

Increasing carbon (C) sequestration in cropland has been recognized as an effective way to reduce CO2 emissions, improve soil structure and promote soil microbial diversity [1,2]. Though protected from microbial decomposition [3], organic C in each aggregate fraction could be affected by agronomic practices (i.e., bio-fertilizer application). Furthermore, the inherent heterogeneity of microbes existing in different aggregate fractions was also influenced by organic material application [2]. Generally, soil aggregate fractions would affect the soil microbial diversity, while soil C contributed to the microbial community structure [4,5]. However, sequestration of soil organic carbon (SOC) was closely related to the C input and output, which were affected by the microbial mineralization of organic materials. Clay soil with high SOC content was more vulnerable to decomposition, showing that dynamics of SOC might also be affected by soil texture which was influenced by aggregate fraction [6]. Therefore, evaluating the sequestration and mineralization of organic C in aggregate fractions and its related microbial community is vital for a deeper understanding of soil C stability.

Studies have demonstrated that long-term fertilization significantly affected SOC stability in different aggregate fractions [7–9], and mineralization of organic C would vary

**Citation:** Zhu, L.; Cao, M.; Sang, C.; Li, T.; Zhang, Y.; Chang, Y.; Li, L. *Trichoderma* Bio-Fertilizer Decreased C Mineralization in Aggregates on the Southern North China Plain. *Agriculture* **2022**, *12*, 1001. https://doi.org/10.3390/ agriculture12071001

Academic Editors: Yinglong Chen, Masanori Saito and Etelvino Henrique Novotny

Received: 9 June 2022 Accepted: 8 July 2022 Published: 11 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

with soil aggregates [6,10]. The mineralization of C in small aggregates was vulnerable to soil conditions in an experiment in which soil types concluded Alfisols, Inceptisols, Mollisols and Ultisols [11]. Reeves et al. [12] showed that the mineralization of C in microaggregates was higher in Vertisol, while Wang et al. [13] found an opposite trend in the soil of Udic Ferralsols. No difference in organic C mineralization between different aggregate fractions was also observed in Dermosols [14]. Such diverse results clearly indicated that mineralization of aggregate-associated organic C varied with soil profiles. Soil microbes, the driving force of organic C decomposition, were affected by soil textures and aggregate class [13,15]. Therefore, it is important to investigate changes in soil microbial community across aggregate fractions to deeply understand soil C stability. Generally, organic materials clearly increased fungi decomposing biosolids and decreased the abundance of CO2 emission-related fungi [16]. However, little is known about variations of soil microbial community in aggregate fractions and how different classes of aggregates protect and sequester C in soil cultivated with *Trichoderma* bio-fertilizer.

The winter wheat-summer maize cropping system is the major planting pattern in the North China Plain (NCP), which occupies around 16 M ha [17]. The Shajiang Calci-Aquic Vertosol, an important soil resource, is widely distributed in the southern NCP. However, research information is lacking about the mineralization of aggregate-associated organic C in this area. Thus, this study was conducted to investigate (1) the aggregate distribution and aggregate-associated organic C under the dry sieving method, (2) the mineralization dynamics of C in aggregates in soil treated with bio-fertilizer and (3) the soil microbial community within aggregate fractions in Shajiang Calci-Aquic Vertosol. We hope this research will provide fundamental evidence to strengthen the understanding of microbial regulation of C dynamics across aggregate fractions and guide appropriate bio-fertilizer application in the southern NCP.

#### **2. Materials and Methods**

#### *2.1. Study Area and Experimental Design*

The study was conducted in Xiping county in 2018 in the southern NCP, China (N 113◦12 , E 33◦27 , with an altitude of 64 m above sea level), with a mean annual precipitation of 786 mm. The field had been cultivated for more than 50 years before this study started. The mean annual temperature is 14.7 ◦C, and the sunshine duration has been 2181 h over the past 50 years. Based on the Chinese Soil Taxonomy, the soil in this area is classified as a Shajiang Calci-Aquic Vertosol with a pH value of 6.9. The soil bulk density, SOC, total nitrogen (TN), mineral nitrogen, available phosphorus and available potassium of the 0–20 cm depth are 1.31 g cm−3, 9.4 g kg−1, 0.88 g kg−1, 14.38 mg kg−1, 6.06 mg kg−<sup>1</sup> and 179.23 mg kg<sup>−</sup>1, respectively.

Five treatments were included in this experiment: CF (100% chemical fertilizer N), BF10 (chemical fertilizer N supplemented with 10% *Trichoderma* bio-fertilizer N), BF20 (chemical fertilizer N supplemented with 20% *Trichoderma* bio-fertilizer N), BF30 (chemical fertilizer N supplemented with 30% *Trichoderma* bio-fertilizer N) and BF50 (chemical fertilizer N supplemented with 50% *Trichoderma* bio-fertilizer N). Each treatment received the same amounts of N, P2O5 and K2O (220, 90 and 90 kg ha<sup>−</sup>1, respectively) during the maize growing season. The chemical fertilizers were urea (N 46%), calcium super phosphate (P2O5 12%) and potassium sulphate (K2O 51%). The nitrogen, phosphorus and potassium contents of *Trichoderma* bio-fertilizer were deducted and supplemented with corresponding chemical fertilizer in each treatment. The details of the *Trichoderma* bio-fertilizer were listed in Zhu et al. [9]. Fermented with wheat straw, the bio-fertilizer had approximately <sup>2</sup> × <sup>10</sup><sup>8</sup> colony forming units per gram of *Trichoderma asperellum* ACCC30536, with 297.2 g kg−<sup>1</sup> organic matter and 120.1 g kg−<sup>1</sup> total nitrogen. Both the *Trichoderma* bio-fertilizer and chemical fertilizer were spread evenly on the surface of each plot and thoroughly mixed with the top 0–20 cm layer by a rotary cultivator a week before sowing the maize. Summer maize in each plot was planted at a density of 65,000 plants ha−<sup>1</sup> in early June and harvested in early October each year.

#### *2.2. Soil Sampling and Analysis*

Soil cores were sampled after the maize harvest in October 2020. Five random soil cores from the top 20 cm layer of each plot were taken and mixed thoroughly as one composite sample. Visible materials were removed and stored at 4 ◦C before analysis. Subsamples of 300 g were shaken on a motorized vibratory sieve-shaker (8411; Zhejiang Chenxin Machine Equipment Co., Ltd., Shaoxing, China) for 3 min with a mesh size of 2, 0.25 and 0.053 mm to obtain four size fractions: >2, 0.25–2, 0.053–0.25 and <0.053 mm.

The incubation experiment was prepared as follows: 30 g samples of the four aggregate fractions were incubated in 500 mL buckets at 25 ◦C and at 60% water holding capacity for 60 days. The released CO2 was trapped in NaOH solution and determined after 5, 10, 15, 20, 30 and 60 days. The trapped CO2 was precipitated with 1.0 mol L−<sup>1</sup> BaCl2 and then titrated with standard 0.1 mol L−<sup>1</sup> HCl to quantify the released CO2. Cumulative mineralized CO2-C was used for comparison between aggregate size classes and fertilization treatments.

The PLFAs were extracted and identified according to Bligh and Dyer [18]. Fatty acid methylesters were analyzed using an Agilent 6850 system (Agilent Technologies, Palo Alto, CA, USA), equipped with a DP-5MS capillary column (25 m × 200 μm, i.d. 0.33 μm). Assignment of microbial categories to the various PLFA biomarkers was based on Luo et al. [19].

#### *2.3. Data Analyses*

The aggregate stability was determined by mean weight diameter (MWD), and MWD was calculated with the following equation:

$$\text{MMWD} = \sum (\mathbf{\hat{x}\_i} \mathbf{W\_i}) \tag{1}$$

Of which Xi was the mean diameter of the class (mm), and Wi was the proportion of each aggregate class in relation to the weight of the soil samples.

Total respiration was calculated with the CO2-C produced from each aggregate class (CO2Agg) and the mass of each aggregate class (MAgg) in the total soil mass (SM), according to Xie et al. [15]:

$$\text{Total respiration (mg C kg}^{-1}\text{)}=\text{CO}\_2\text{ Agg} \times \text{MAgg/SM} \times 100\tag{2}$$

The contribution of organic C mineralization of each aggregate fraction to bulk soil was calculated with the ratio of cumulative CO2-C production of each aggregate multiplied by the mass of each aggregate fraction/total mineralization of organic C in bulk soil.

#### *2.4. Statistical Analysis*

For each of these variables, a mean value was obtained from the results for the three composite samples, significant differences between the means were identified by performing a one-way analysis of variance (ANOVA), and the least significant difference (LSD) was computed to compare the means of the above variables (*p* < 0.05) using the statistical package SPSS 19.0. Statistical correlations between mineralization of organic C and microbial community were assessed with correlation analysis, and Pearson correlation coefficients were presented. Prior to analysis, data were tested to ensure that they met homogeneity of variance and normality requirements. All the figures were accomplished with SigmaPlot 12.0 (Systat Software Inc., London, UK).

#### **3. Results**

#### *3.1. Distribution of Soil Aggregates*

Soil aggregate distribution was significantly affected by the bio-fertilizer application. As the dominant fraction, >2 mm aggregate accounted for 34.90% in CF treatment and for 48.94% in BF20 treatment. The proportion of 0.25–2 mm aggregate ranged from 26.08% (CF) to 36.60% (BF20) (Table 1). For the >2 mm aggregate, the highest value was recorded in BF20, followed by BF30, BF50, BF10 and CF. Compared with CF, the increases of >2 mm aggregate

were 18.70%, 40.22%, 35.97% and 32.75% in BF10, BF20, BF30 and BF50, respectively (*p* < 0.05). No significance was observed among BF20, BF30 and BF50 in >2 mm aggregate. Additionally, bio-fertilizer significantly increased the content of 0.25–2 mm aggregate, while BF20, BF30 and BF50 showed no difference between each other. Differently, the distribution of 0.053–0.25 mm aggregate showed a trend of CF > BF10 > BF20 = BF30 = BF50. Bio-fertilizer significantly decreased the content of <0.053 mm aggregate, and no difference was observed among the four bio-fertilizer treatments. Moreover, bio-fertilizer significantly increased MWD and the increase of MWD in BF20 was the highest (35.6%). *Trichoderma* bio-fertilizer promoted the formation of macro-aggregate (>0.25 mm) in soil but decreased proportions of micro-aggregate (0.053–0.25 mm) and clay (<0.053 mm), and then increased soil aggregate stability.

**Table 1.** Distribution of soil aggregate and soil aggregate stability assessed by dry sieving method.


Notes: Values are means ± standard deviation of three replicate plots. Different lower case letters in the same column indicate significant differences between treatments for the same aggregate fraction (*p* < 0.05). CF: 100% chemical fertilizer N; BF10: chemical fertilizer N supplemented with 10% bio-fertilizer N; BF20: chemical fertilizer N supplemented with 20% bio-fertilizer N; BF30: chemical fertilizer N supplemented with 30% bio-fertilizer N; BF50: chemical fertilizer N supplemented with 50% bio-fertilizer N.

#### *3.2. Organic C in Bulk Soil and Aggregate Fractions*

As shown in Figure 1, the aggregate-associated organic C was significantly affected by the bio-fertilizer application and aggregate size. Compared with CF, organic C in aggregates of >2, 0.25–2, 0.053–0.25 and <0.053 mm were significantly increased by 11.82%, 11.98%, 10.86% and 11.81% in BF10, respectively. The organic C in aggregates of >2, 0.25–2, 0.053–0.25 and <0.053 mm in BF20 was 18.15%, 18.60%, 17.75%, and 18.46% higher than that in CF, respectively. The increases of organic C in >2, 0.25–2, 0.053–0.25 and <0.053 mm aggregate fractions were 17.85%, 18.08%, 17.15% and 18.04% in BF30, and the increases were 17.99%, 18.31%, 17.15% and 17.74% in BF50, respectively, relative to CF. No difference in aggregate-associated organic C was observed in >2, 0.25–2 and 0.053–0.25 mm aggregate fractions among BF20, BF30 and BF50. Obviously, the organic C in <0.053 mm aggregate in BF50 was significantly lower than that in BF20.

#### *3.3. Mineralization of Organic C in Aggregate*

The mineralization dynamics of aggregate-associated organic C were noticeably affected by *Trichoderma* bio-fertilizer application, especially in the first 15 days (Figure 2), and then the mineralization rates were relatively low and remained stable. The BF50 treatment yielded the highest mineralization rate, and the lowest value was recorded in CF treatment across aggregate fractions and bulk soil. The >2 mm aggregate produced the highest amount of CO2-C, and the aggregate of <0.053 mm yielded the lowest value.

**Figure 1.** Organic C in different aggregate classes in the five treatments. Different lowercase letters indicate significant differences between treatments, and different uppercase letters indicate differences between aggregate size classes at *p* < 0.05. Error bars represent standard deviation of means (*n* = 3). CF: 100% chemical fertilizer N; BF10: chemical fertilizer N supplemented with 10% bio-fertilizer N; BF20: chemical fertilizer N supplemented with 20% bio-fertilizer N; BF30: chemical fertilizer N supplemented with 30% bio-fertilizer N; BF50: chemical fertilizer N supplemented with 50% bio-fertilizer N.

After 60 days of incubation, the cumulative mineralization of aggregate-associated organic C in the five treatments differed from each other (Figure 3a). The >2 mm aggregate shared the highest CO2-C among the five treatments, while 0.25–2 and 0.053–0.25 mm aggregate fractions produced similar amounts of CO2-C. The CO2-C yielded from the <0.053 mm aggregate was the lowest among the four aggregate fractions. Compared with CF, the highest values of CO2-C were recorded in the BF50 treatment across aggregate fractions. Increases of CO2-C in BF50 were 121.00% (>2 mm), 21.97% (0.25–2 mm), 35.18% (0.053–0.25 mm) and 77.66% (<0.053 mm) relative to CF, respectively. The bulk soil respired 130.28, 188.35, 244.37, 288.16 and 314.59 mg C kg−<sup>1</sup> soil in CF, BF10, BF20, BF30 and BF50, respectively (Figure 3b). Further analysis showed that the mineralization rate of organic C was generally similar for CF and BF10, while a significantly lower value was observed in BF20 relative to CF (Table 2). There was no difference between BF20, BF30 and BF50 in aggregate fractions of 0.25–2, 0.053–0.25 and <0.053 mm. However, the highest mineralization rate was observed in BF50 treatment in bulk soil.

**Figure 2.** The dynamics of organic C of mineralization in >2 mm (**a**), 0.25–2 mm (**b**), 0.053–0.25 mm (**c**), <0.053 mm aggregate fractions (**d**) and bulk soil (**e**) for the five treatments. Values are means of three replicates. CF, BF10, BF20, BF30 and BF50 are defined as in Figure 1.

**Figure 3.** Cumulative mineralization of organic C in different aggregate fractions (**a**) and bulk soil (**b**) under the five treatments. Different lowercase letters indicate significant differences among treatments at *p* < 0.05. CF, BF10, BF20, BF30 and BF50 are defined as in Figure 1.



Notes: Different lowercase letters within a column indicate significant differences between treatments, and different uppercase letters in the same row indicate differences between aggregate fractions at *p* < 0.05. CF, BF10, BF20, BF30 and BF50 are defined as in Table 1.

#### *3.4. Microbial Community Composition in Aggregate*

*Trichoderma* bio-fertilizer significantly increased the PLFAs of total microbes, bacteria, fungi, actinomycetes, G+ bacteria and G− bacteria across aggregate fractions and bulk soil (Figure 4). Total PLFA and bacteria PLFA in each aggregate fraction and bulk soil increased with increasing bio-fertilizer application rate, and BF50 yielded the highest PLFAs. Compared with CF, increases of bacteria PLFA in >2 mm fraction were 79.7%, 130.0%, 141.0% and 148.5% in BF10, BF20, BF30 and BF50; increases of fungi PLFA were 108.4%, 163.3%, 172.9% and 195.2%. Actinomycetes PLFAs were increased by 20.9%, 20.4%, 29.6% and 54.4% in BF10, BF20, BF30 and BF50, respectively, relative to CF. Increases of G+ PLFA reached 34.0%, 44.1%, 54.1% and 61.1% in BF10, BF20, BF30 and BF50, respectively, relative to CF. Additionally, G− PLFAs were enhanced by 22.5%, 34.3%, 38.8% and 45.5% in BF10, BF20, BF30 and BF50, respectively, in comparison with CF. The PLFAs in 0.25–2, 0.053–0.25 and <0.053 mm aggregate fractions showed a similar trend to that in >2 mm aggregate fraction.

**Figure 4.** The soil microbial community in >2 mm (**a**), 0.25–2 mm (**b**), 0.053–0.25 mm (**c**), <0.053 mm aggregate fractions (**d**) and bulk soil (**e**) for the five treatments. Error bars represent the standard deviation of the mean (*n* = 3). Different lowercase letters above bars for each PLFA indicate significant differences at *p* < 0.05. CF, BF10, BF20, BF30 and BF50 are defined as in Figure 1.

The soil microbial community diversity was significantly affected by *Trichoderma* biofertilizer application and aggregate fractions (Table 3). No significance of Shannon index *H* was observed among treatments in >2 and 0.25–2 mm aggregate fractions. Compared

with CF, BF30 and BF50 significantly decreased *H* by 11.72% and 10.94% in 0.053–0.25 mm aggregate and decreases of *H* were 8.00% and 4.80% in <0.053 mm in BF30 and BF50 treatments, respectively. Bio-fertilizer significantly decreased values of *H* in bulk soil, and the lowest *H* value was observed in BF30. Bio-fertilizer significantly increased the value of fungi PLFA/bacteria PLFA, and the highest value was recorded in BF30 in >2 mm aggregate fraction. Differently, BF20 yielded the highest value of fungi PLFA/bacteria PLFA in 0.25–2, 0.053–0.25 and <0.053 mm aggregates and bulk soil. BF10, BF20 and BF50 significantly increased values of G+ PLFA/G− PLFA, and the increases were 9.71%, 6.80% and 8.74%, respectively, relative to CF. No difference in G+ PLFA/G− PLFA was observed between CF and BF30 in >2 mm aggregate. Notably, decreases of G+ PLFA/G− PLFA in BF20 were the most in 0.25–2, 0.053–0.25 and <0.053 mm aggregate fractions and bulk soil relative to CF.

**Table 3.** The diversity of soil microbial community in aggregate fractions and bulk soil in different treatments.


Notes: Values are means ± standard deviation of three replicate plots. Different lowercase letters within a column for the same index indicated significant differences among treatments at *p* < 0.05. CF, BF10, BF20, BF30 and BF50 are defined as in Table 1.

As shown in Table 4, aggregate-associated organic C was significantly positively correlated with bacteria PLFA (*R =* 0.865, *p* < 0.01) and fungi PLFA (*R =* 0.881, *p* < 0.01), but was significantly negatively correlated with G+ PLFA/G− PLFA (*R =* −0.450, *p* < 0.05) and *H* (*R =* −0.665, *p* < 0.05). A significant positive correlation was observed between aggregate distribution and fungi PLFA (*R =* 0.141, *p* < 0.05). However, there was no correlation between mineralization of aggregate-associated organic C and microbial community.

**Table 4.** Relationships between mineralization of organic C and microbial community.


Notes: \* indicated significance at *p* < 0.05; \*\* indicated significance at *p* < 0.01.

#### **4. Discussions**

*4.1. Trichoderma Bio-Fertilizer Changed Aggregate Distribution and Its Associated Organic C Content*

Generally, both wet and dry sieving methods are used to assess soil aggregate conditions. However, dissolved organic matter in soil would lose, and the habitat of microbes would be damaged when aggregates are sieved by the wet sieving method [13]. Therefore, it would be more convictive to assess the change of aggregate-associated organic C in soil using the dry sieving method.

Our results showed that *Trichoderma* bio-fertilizer improved the aggregation capacity of macro-aggregates (>0.25 mm), consistent with the results of Zhu et al. [9]. Wang et al. [20] also found that the proportion of 0.20–2.00 mm aggregate fraction increased, but the <0.002 mm aggregate fraction decreased after bio-fertilizer application, bio-fertilizer significantly promoted the formation of large aggregates. This might be due to the fact that amounts of fungi introduced by bio-fertilizer promoted the decomposition of soil organic matter and the formation of organo-mineral materials [2,21]. Higher soil microbial biomass after bio-fertilizer application might also promote soil aggregation [9]. Negatively charged polysaccharides and polyuronic acid, along with the growth of bacteria, were beneficial for clay bounding in soil [22], and cementation of organic matter in soil contributed to the formation of macro-aggregates. Increases of macro-aggregates protected labile C from microbial mineralization, which in turn promoted stabilization of soil aggregate [23]. Therefore, a decrease in physical protection caused by the destruction of macro-aggregates had negative effects on labile C, causing organic C mineralization [24].

Application of bio-fertilizer markedly increased soil organic C stock in bulk soil and almost each aggregate fraction, and macro-aggregate contained higher organic C than micro-aggregate. However, an opposite trend was observed by Xie et al. [15], who conducted a field experiment on Anthrosols. Diverse results among studies might be attributed to variations in binding agents during soil aggregation. Generally, increases in aggregate-associated organic C are in line with classes of aggregate fractions since organic matter is the major binding agent [25]. The organic C content was significantly higher in bio-fertilizer treatments than in CF treatment across aggregate fractions, and the BF20 treatment shared higher values. When organic material was applied, plant growth and soil microbes were both promoted [9,26], which in turn increased aggregate-associated organic C. The photosynthetic activities of algae species after bio-fertilizer application would also contribute to increasing aggregate-associated organic C [27]. Consistent with Yilmaz and Snmez [28], our results also showed that bio-fertilizer application increased organic C across aggregate fractions relative to the chemical fertilizer alone treated soil. In addition, aggregate-associated organic C was at a similar level among bio-fertilizer treatments. However, Wang et al. [20] observed that increases of organic C in macroaggregate were higher than that in micro-aggregate after bio-fertilizer application. The diverse results suggested that soil type would be an inevitable factor in the sequestration of organic C in aggregates.

Most studies reported a linear relationship between organic C and organic materials addition rates [6,29,30]. However, in our study, organic C in aggregate fractions and bulk soil significantly increased when the bio-fertilizer application rate was below 20 t ha−1, while the increase was not significant when the bio-fertilizer application rate was above it. Thus, there was a threshold effect of the bio-fertilizer application on aggregate-associated organic C in the given soil.

#### *4.2. Mineralization of Organic C in Aggregates*

It is known that the rate of organic C respiration is an effective index to assess the capability of soil C sequestration. In the present study, significantly higher CO2-C in >2 mm aggregate was recorded. The mineralization rate of organic C in >2 mm aggregate was also significantly higher across aggregate fractions in bio-fertilizer-treated soil. These results were in line with that of Mustafa et al. [8], who recorded higher C mineralization per unit of organic C in macro-aggregate. Similarly, Kan et al. [6] found that organic C in macro-aggregate was the main source of mineralized C, of which >2 mm aggregate contributed 38.2–43.6% to the cumulative mineralization. This might be due to the fact that macro-aggregates dominated the bulk soil (>2 mm aggregate occupies 34.90–48.94% and 0.25–2 mm aggregate occupies 26.08–36.60% in the present study). Organic C in macroaggregate, mainly originated from plant materials, was more labile. Additionally, larger

pores in macro-aggregate increased the transposition of oxygen and microbial activities, which promoted the process of organic C mineralization in macro-aggregate. In the case of Anthrosols, Xie et al. [15] found greater CO2-C produced from micro-aggregate than macro-aggregate, suggesting that mineralizable C in micro-aggregate was much higher than that in macro-aggregate. However, Rabbi et al. [31] showed that there was no difference in organic C mineralization between macro-aggregate and micro-aggregate in Oxisols. Thus, mineralization of aggregate-associated organic C might be a better reflection of C sequestrated in aggregates. Generally, organic C in micro-aggregate was relatively stable for soil in which organic C dominated soil aggregation [32].

When compared with CF, the cumulative mineralization of organic C in aggregates and bulk soil significantly increased with increasing bio-fertilizer application rates. This might be due to the fact that bio-fertilizer application enhanced amounts of available nutrients and promoted microbial activities, which both contributed to the mineralization of organic C [33,34]. Similarly, Xie et al. [15] showed that manure application significantly increased the mineralization of aggregate-associated organic C, and that proportions of organic C mineralized were also similar in the surface layer. Notably, BF20 treatment shared a significantly lower rate of organic C respiration across aggregate fractions (except for <0.053 mm aggregate). The result might indicate that the bio-fertilizer application rate was an important factor in regulating the mineralization of aggregate-associated organic C. The decreased C mineralization in bio-fertilizer application treatments, especially in BF20, might be explained by the dilution effect. Namely, the significant increase of organic C mineralization was diluted by the large increase of aggregate-associated organic C. Generally, organic C in aggregate fractions and bulk soil in the BF20 treatment was more resistant to mineralization, which would be beneficial for soil C sequestration.

#### *4.3. Bio-Fertilizer Alter Microbial Community in Aggregates*

PLFA can well reflect changes in soil microbial communities influenced by fertilization. Our study showed that bio-fertilizer significantly increased total PLFA and each microbial group PLFA across aggregates and in bulk soil. Amounts of organic material introduced into the soil would provide extra nutrients and energy to soil microbes, which might directly increase the PLFAs of soil microorganisms within aggregates [35]. This was confirmed by the significant positive correlations between PLFAs and aggregate-associated organic C. In line with the results of Jiang et al. [36], our results also showed that total PLFA, bacteria PLFA and fungi PLFA in macro-aggregate were greater than that in micro-aggregate. It might be that levels of organic C and soil microbes were relatively low in micro-aggregates compared to those in macro-aggregates and that the turnover rate of C was slow. Differently, Wang et al. [37] showed that PLFAs of bacteria and fungi in macro-aggregate were significantly increased by manure application while no significant difference was observed in micro-aggregate. Due to the heterogeneity of soil properties, differences in microbial community composition across aggregate fractions would appear. Aggregates with larger particle sizes supported higher nutrient levels, which were conducive to microbial colonization [2,38].

Generally, bio-fertilizer high in C/N ratio would be more beneficial to fungal growth [39], which led to increased fungi PLFA/bacteria PLFA. In addition, significant changes in bacterial PLFA and fungal PLFA represented variations in soil microbial community, though bacteria were the dominant group among different treatments. It was believed that the higher ratio of fungi PLFA to bacteria PLFA indicated a more stable soil ecosystem and better soil quality [40]. The higher fungi PLFA/bacteria PLFA in BF20 treatment suggested that appropriate bio-fertilizer application helped maintain the stability of the soil ecosystem. Our results showed that bio-fertilizer significantly increased contents of >0.25 mm aggregates and its fungi PLFA, suggesting that fungi contributed to the formation of macroaggregate and soil aggregate stability via filamentous growth and excretion production [41]. Thus, bio-fertilizer could increase the content of macro-aggregate by promoting fungal growth, and fungi were important factors for bio-fertilizer to promote the formation of

macro-aggregates. Additionally, fungi could decompose foreign nutrients through hyphal movement and had high assimilation efficiency of C source, while bacteria did not have this advantage. However, extra organic material usually promoted the growth of G− bacteria, which preferred plant-derived carbon sources, while the G+ community usually participated in organic matter and litter decomposition [42]. A lower G+ PLFA/G− PLFA meant a better soil nutritional condition [43]. Thus, decreased G+ PLFA/G− PLFA after bio-fertilizer application indicated that the soil environment shifted to more eutrophic conditions. In addition, the fact that the ratio of G+ PLFA/G− PLFA was negatively correlated with aggregate-associated organic C also indicated the improved soil quality after bio-fertilizer application. Therefore, the bio-fertilizer application was beneficial to improving soil nutrient status and ecological buffer capacity of large aggregates in the farmland ecosystem. Though the mineralization of aggregate-associated organic C was positively correlated with fungi, *Trichoderma* was not distinguished from fungi in the present study. The abundance of specific microorganisms was also not assessed, which played vital roles in organic carbon turnover. Thus, further studies focusing on the effects of bio-fertilizer on soil microbial community composition are needed, which may help to develop an environmental bio-fertilizer application pattern in agricultural production.

#### **5. Conclusions**

Bio-fertilizer considerably increased proportions of macro-aggregate, organic C sequestration and microbial PLFAs across aggregate fractions and bulk soil. The increase of MWD in BF20 reached 35.6%, significantly higher than other treatments. Bio-fertilizer increased the fungi PLFA/bacteria PLFA but decreased G+ PLFA/G− PLFA, and the BF20 shared the greatest changes. Increases of aggregate-associated organic C of >2, 0.25–2, 0.053–0.25 and <0.053 mm in BF20 were 18.15%, 18.60%, 17.75% and 18.46%, respectively, relative to CF. However, the cumulative mineralization of organic C was relatively low in BF20. Thus, the promotion of organic C stock was generally higher in BF20, while the proportion of organic C mineralization was relatively low across aggregate fractions. Correlation analysis showed that microbial community in aggregate was correlated with increases of soil C, of which bacteria and fungi contributed more than actinomycetes. This study highlighted the vital role of *Trichoderma* bio-fertilizer in regulating C mineralization at the aggregate level and provided scientific bases for bio-fertilizer application in Shajiang Calci-Aquic Vertosol. However, we did not determine the distinct keystone taxa and their co-occurrence patterns. Further study concerning keystone taxa at the aggregate level after bio-fertilizer application is needed.

**Author Contributions:** Writing—original draft preparation, L.Z. and Y.C.; writing—review and editing, M.C., C.S. and T.L.; visualization, Y.Z.; funding acquisition, L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was financially supported by the Scientific and Technological Project of Henan Provincial Science and Technology Department of China (222102320276), the National Key Research and Development Program of China (2018YFD0300704).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available upon request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Carbon Storage Potential of Agroforestry System near Brick Kilns in Irrigated Agro-Ecosystem**

**Nayab Komal 1, Qamar uz Zaman 1,\*, Ghulam Yasin 2, Saba Nazir 1, Kamran Ashraf 3, Muhammad Waqas 1, Mubeen Ahmad 1, Ammara Batool 1, Imran Talib <sup>1</sup> and Yinglong Chen 4,5,\***


**Abstract:** The current study was conducted to estimate the carbon (C) storage status of agroforestry systems, via a non-destructive strategy. A total of 75 plots (0.405 ha each) were selected by adopting a lottery method of random sampling for C stock estimations for soil, trees and crops in the Mandi-Bahauddin district, Punjab, Pakistan. Results revealed that the existing number of trees in selected farm plots varied from 25 to 30 trees/ha. Total mean tree carbon stock ranged from 9.97 to 133 Mg C ha<sup>−</sup>1, between 5–10 km away from the brick kilns in the study area. The decreasing order in terms of carbon storage potential of trees was *Eucalyptus camaldulensis* > *Syzygium cumin* > *Popolus ciliata* > *Acacia nilotica* > *Ziziphus manritiana* > *Citrus sinensis* > *Azadirachtta Indica* > *Delbergia sisso* > *Bambusa vulgaris* > *Melia azadarach* > *Morus alba*. Average soil carbon pools ranged from 10.3–12.5 Mg C ha−<sup>1</sup> in the study area. Meanwhile, maximum C stock for wheat (2.08 <sup>×</sup> <sup>10</sup><sup>6</sup> Mg C) and rice (1.97 <sup>×</sup> <sup>10</sup><sup>6</sup> Mg C) was recorded in the cultivated area of Tehsil Mandi-Bahauddin. The entire ecosystem of the study area had an estimated woody vegetation carbon stock of 68.5 Mg C ha−<sup>1</sup> and a soil carbon stock of 10.7 Mg C ha<sup>−</sup>1. These results highlight that climate-smart agriculture has great potential to lock up more carbon and help in the reduction of CO2 emissions to the atmosphere, and can be further used in planning policies for executing tree planting agendas on cultivated lands and for planning future carbon sequestration ventures in Pakistan.

**Keywords:** agroforestry; brick kilns; carbon emissions; climate change; carbon sinks; carbon stock

#### **1. Introduction**

Pakistan is predicted to be among the ten countries most affected by climate change, according to the 2019 Global Climate Risk Index [1]. Global climate change, by increasing the amount of greenhouse gases (GHG) in the atmosphere, is causing severe environmental and climatic effects. Carbon dioxide (CO2) is one of the most commonly highlighted greenhouse gasses. Global climate change and the increase in the trend of CO2 emissions are a growing concern today [2]. Pakistan is confronting this powerful danger to social, environmental, and economic development [3]. The impacts of climate change can be categorized into extreme and non-extreme types [4,5]. The World Bank [6] recognizes five foremost factors through which climate change will affect agricultural production: change in precipitation pattern and temperature, climatic variability, CO2 fertilization and surface water runoff. Reilly et al. [7] found that higher rainfall results in reduction of yield.

**Citation:** Komal, N.; Zaman, Q.u.; Yasin, G.; Nazir, S.; Ashraf, K.; Waqas, M.; Ahmad, M.; Batool, A.; Talib, I.; Chen, Y. Carbon Storage Potential of Agroforestry System near Brick Kilns in Irrigated Agro-Ecosystem. *Agriculture* **2022**, *12*, 295. https:// doi.org/10.3390/agriculture12020295

Academic Editor: Rosa Francaviglia

Received: 18 January 2022 Accepted: 16 February 2022 Published: 18 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Variations in climate cause many people to move into poverty and food insecurity [8]. In Pakistan, the general origin of air pollution is customarily untidy industrial buildings of brick formation, present in the peri-metropolitan and rural regions [9]. Inferior quality fuels, including corncobs, rubber tires, rice straw, bagasse, rice husk, coal, oil and wood, used in these brick kilns produce fly ash particles that deposit on nearby plants affecting their photosynthetic potential [6]. Brick kilns release over 1072 million tons of CO2 emissions into the atmosphere per year, making 2.7% of total emissions. Most of the brick kilns in the rural areas in Pakistan use conventional technology that is very dangerous from an environmental aspect. According to the Punjab Disaster Management Authority (PDMA), 37.4% of brick kilns have moved to zig-zag technology in Punjab alone [1]. All the types of fuel utilized in this type of kiln cause a high concentration of pollutants in gaseous form in the air, with destructive effects on the atmosphere, plants and people [10,11]. Islam et al. [12] estimated that the soil close to the brick kiln was reduced in quantity when compared with the same soil further from brick kiln, showing a variation in agricultural production. The mean values of total nitrogen, available phosphorus and sulfur were significantly less in the soil samples close to the brick kiln (0.05%, 12.4, and 8.36 ppm, respectively) than those in the soil further from the brick kiln (0.06%, 24.6, and 11.7 ppm, respectively).

Agriculture is a key economic sector that contributes 21% to the gross domestic product (GDP), employs 45% of the total workforce and contributes about 60% to exports [13]. Changing climatic variables, particularly temperature and rainfall, will introduce several challenges to agriculture in the future. Changes in the frequency and intensity of droughts, flooding, and storm damage are anticipated [13]. Agroforestry could provide adaptation to this climate change [14] by protecting crops from temperature elevation. On the other hand, this causes a decrease in soil evaporation, wind seeds, and transpiration of crops. The carbon (C) sequestration above ground could easily be increased by planting trees and this also increases carbon in the soil on the land where crops are cultivated [15,16]. The general sequestration of carbon due to such actions has been assumed to be 9, 21, 50, and 63 Mg C ha−<sup>1</sup> in temperate, sub-humid, semiarid, and humid regions, respectively [17]. The planting of trees together with crops has many advantages, involving higher soil richness, limitation of soil erosion, lower water logging, decreased fermentation and eutrophication of streams and rivers, enhancement of local biodiversity, and reduction of pressure on common forests for fuel [18–20]. At a global scale, unproductive croplands of about 630 million ha could be used for agroforestry as part of an ecological engineering tool to potentially sequester 586,000 Mg C year−<sup>1</sup> by 2040. Moreover, in present national and global monitoring protocols for carbon, there is a need to include agroforestry in C stocks to estimate the share of this abandoned pool in a precise way. Agroforestry systems are considered as a carbon sink to sequester CO2, so there is a need to evaluate the carbon sequestration potential of agroforestry systems in order to decrease emissions near brick kilns on irrigated land. The main objectives of the current investigation were (1) to assess and quantify the potential of agroforestry systems in C stock at various distances from brick kilns in an irrigated agro-ecosystem, (2) to discover differences between agroforestry systems and the carbon capturing capacity of crops, trees and soil organic carbon (SOC).

#### **2. Materials and Methods**

#### *2.1. Study Locations and Sampling Methodology*

The present study was conducted in the Mandi-Bahauddin district, central Punjab, Pakistan (Figure 1). According to the population survey in 2017, Mandi-Bahauddin is the 41st city by population. This city is located between the river Chenab (39 km in south) and the river Jhelum (12 km in north). It lies between 30◦8 to 32◦40 north latitude, and 73◦36 to 73◦37 east longitude. It has a total area of 2673 km2. Major rivers are River Jhelum, with Rasool Barrage & Chenab and the Qadirabad Barrage. According to the Pakistan Bureau of Statistics (2019), the total population of Mandi-Bahuddin is 198,609 with a population density of 30 km2. The annual change in population is 1.68%. The elevation of this district is 220 m above the sea level. The weather conditions experience a maximum average

temperature of 45 ◦C and minimum average temperature of 3 ◦C, with average rainfall of 388 mm and wind direction mostly from north to east (Pakistan meteorological department, PMD). According to the US Department of Agriculture classification system, the soil type of the study area falls in the category of loam to clay loam soil with organic matter ranging from 0.50–1.01%. Most farmers in these areas use conventional practices (chemical fertilizer) for crop production instead of using organic methods for crop production to improve soil health. According to the agricultural department, the main crops in study area are wheat, rice and sugarcane and the total cultivated area is 214,348.83 ha. Maize and pulses are also grown in smaller areas.

**Figure 1.** Map of the study area and sampling points.

#### *2.2. Above and Below Ground Trees' Biomass Carbon Estimation*

Field visits were carried out from November 2020 to May 2021 at regular intervals for the collection of data in 3 selected Tehsils, i.e., Mandi-Bahauddin, Phalia, Malakwal. A total of 75 quadrate plots of 0.405 ha = 1 acre each with agroforestry practices were selected and a tree inventory such as diameter at breast height in cm (DBH) and height (m) was taken for each plot (Table 1). By adopting a lottery method of random sampling, the sample collection sites were selected on the basis of brick kiln location at 0, 5 and 10 km distance. Mean values of tree height and diameter at breast height of potential tree species at various distances from the brick kilns are depicted in Table 2.

**Table 1.** Number of plots measured from selected Tehsils in Mandi-Bahauddin sampling.


**Table 2.** Diameter at Breast Height (DBH) and Height of Potential Tree Species at Various Distance from Brick Kilns in the Agroforestry System of Mandi-Bahauddin district (Mean Values).


*Acacia nilotica* [Fabaceae], *Delbergia sisso* [Fabaceae], *Melia azedarach* [Meliaceae], *Citrus reticulate* [Rutaceae], *Popolus ciliate* [Salicaceae], *Eucalyptus camaldulensis* [Myrtaceae], *Melia azedarach* [Meliaceae], *Populus deltoides* [Salicaceae], *Syzygium cumini* [Myrtaceae], *Ziziphus mauritiana* [Rhamnaceae], *Azadirachtta indica* [Meliaceae], *Morus alba* [Moraceae] and *Bambusa vulgaris* [Poaceae] were the most commonly planted tree species in and along the farm fields in all three tehsils.

#### *2.3. Total Carbon Stock Estimate*

Allometric equations from the literature were used for the estimation of tree biomass and where appropriate corrected for log bias (Table 3). In case of non-availability of allometric equations, 26% of the above ground biomass was assumed as below ground biomass [21–23]. Next, biomass of individual trees was scaled to biomass/plot, biomass/hectare, and carbon stock/hectare. Contents of carbon were measured from biomass by presumption that the dry mass contains 48.1% of carbon [21,24]. The total tehsil tree carbon was estimated by multiplying carbon amount per hectare from sampled plots by the total area of the tehsil.


**Table 3.** Allometric equations for the calculation of above and below ground biomass.

AGB = aboveground biomass; B = belowground biomass; D = tree diameter at 1.3 m (cm); H = total tree height (m), BA = individual tree basal area (cm2).

#### *2.4. Soil Sampling Collection*

From 0, 5 and 10 km away from the brick kilns, soil was sampled in a random subset of plots to represent the major tree and crop combinations. Soil samples (*n* = 420) were collected at a depth of 0–30 cm near the base of a randomly selected tree, from the four cardinal directions. Samples were stored in polythene bags and analyzed at the Soil and Water Testing Laboratory for Research, Bahawalpur. A 100 cm3 stainless-steel cylinder was used to measure soil bulk density. After being air-dried and passed through a 2 mm sieve, organic carbon was measured using the Walkley–Black method [33]. To calculate the soil carbon per hectare, the values of bulk density, soil depth, and percentage of organic carbon were then multiplied [34].

#### *2.5. Crop Carbon Stock Determination*

Wheat and rice plants were manually harvested to the required depth, sun dried and then weighed with a spring balance to measure above and below ground biomass per plot, converted into Mg ha−1. From each plot at different locations, a different number of tillers of wheat and rice crops were selected randomly from an area of 1 m2. Above ground and below ground C stock in wheat and rice samples was determined by multiplying the respective above ground and below ground biomass with carbon conversion factor of 0.45, as explained by Prommer et al. [35]. Next, biomass of individual crop was scaled to biomass/hectare, carbon stock/hectare, and finally on a tehsil basis.

#### *2.6. Data Analysis*

Collected data was analyzed using "Statistix 8.1" and "Statistix 10" statistical packages. Data regarding analytical analysis of AGB, BGB and SOC were analyzed using descriptive statistics. Graphical work was performed using Microsoft Office software (Version, 2016; Microsoft Corporation, Albuquerque, NM, USA).

#### **3. Results**

#### *3.1. Tree Abundance and Distribution*

Variations in trees' abundance and distribution were observed in the Mandi-Bahauddin district. All the tehsils showed variations in trees' abundance; the maximum percentage of trees species was shown in Mandi-Bahauddin (48%), followed by Malakwal (34%) and Phalia (18%) (Figure 2A). All the values regarding the tree's distribution in Mandi-Bahauddin varied based on brick kiln distance. Mandi-Bahauddin showed variations in tree distribution, with maximum percentage (15%) noticed for *E. camaldulensis* (15%) and *V. nilotica* (15%) followed by *C. sinensis* (14%), *M. azadarach* (14%), *D. sisso* (12%) and *P. ciliate* (10%), while the lowest tree distribution was recorded in *S. cumin* (3%), *Z. manritiana* (3%) and *B. vulgaris* (2%) (Figure 2B). Phalia showed variations in tree distribution, the maximum percentage of tree species being *E. camaldulensis* (19%) and *V. nilotica* (17%) followed by *C. sinensis* (14%), *M. azadarach* (12%), *D. sisso* (10%) and *P. ciliate* (10%), while the lowest tree distribution was noticed in *A. indica* (6%), *S. cumin* (2%), *Z. manritiana* (2%) and *B. vulgaris* (2%) (Figure 2C). Malakwal showed variations in tree distribution, the maximum percentage of trees being *V. nilotica* (15%) and *E. camaldulensis* (14%), followed by *C. sinensis* (13%), *M. azadarach* (13%), *D. sisso* (10%) and *P. ciliate* (10%), while the lowest tree distribution was found in *A. indica* (6%), *S. cumin* (5%), *Z. manritiana* (3%) and *B. vulgaris* (3%) (Figure 2D).

#### *3.2. Carbon Stock of Trees*

The above ground and below ground carbon stock and total carbon stock potential of trees in Mandi-Bahauddin district were different, which describes the variations in carbon stock capacity of potential species of all three tehsils of the district (Table 4). Maximum C stock was recorded in *E. camaldulensis* (3951 Mg C ha−1), followed by *S. cumin* (282 Mg C ha−1), *P. ciliate* (75.9 t/ha), and minimum C stock was observed in *C. sinensis* (17.3 Mg C ha<sup>−</sup>1) in Mandi-Bahauddin. In Phalia, maximum C stock was shown in *A. nilotica* (18.5 Mg C ha<sup>−</sup>1), followed by *S. cumin* (16.8 Mg C ha<sup>−</sup>1) and *B*. *vulgaris* (2.5 Mg C ha<sup>−</sup>1) and minimum C stock was observed in *M. alba* (0.02 Mg C ha−1). Maximum C stock was shown in *P. ciliata* (44.5 Mg C ha<sup>−</sup>1), followed by *A. nilotica* (21.2 Mg C ha<sup>−</sup>1), *E. camaldulensis* (20.9 Mg C ha<sup>−</sup>1), and minimum C stock was observed in *M. azadarach* (0.05 Mg C ha<sup>−</sup>1) in Malakwal. Maximum C stock was shown in the trees of Mandi-Bahauddin followed by Malakwal, and minimum C stock was shown in Phalia. The descending order in terms of carbon sequestration potential of trees in Mandi-Bahauddin was *S. cumin* > *E. camaldulensis* > *P. ciliate* > *B. vulgaris* > *A. nilotica* > *Z. manritiana* > *M. azadarach* > *A. Indica* > *D. sisso* > *M. alba* > *C. sinensis*. The descending order in terms of carbon sequestration potential of trees in Tehsil Phalia was *E. camaldulensis* > *A. nilotica* > *D. sisso* > *Z. manritiana* > *A. Indica* > *P. ciliate* > *S. cumin* > *M. alba* > *M. azadarach* > *C. sinensis* > *B. vulgaris*. The decreasing order in terms of carbon stock potential of trees in Malakwal was *E. camaldulensis* > *P. ciliate* > *C. sinensis* > *A. nilotica* > *Z. manritiana* > *A. Indica* > *S. cumin* > *B. vulgaris* > *D. sisso* > *M. alba* > *M. azedarach* (Table 4).

**Figure 2.** Tree species abundance in Mandi-Bahauddin district (**A**) and tree species distribution expressed as a percentage of total basal area for agroforestry plots in tehsils Mandi-Bahauddin (**B**), Phalia (**C**) and Malakwal (**D**) in Punjab, Pakistan.



#### *3.3. Total Tree Carbon Stock in Mandi-Bahauddin District*

All the values regarding total tree carbon stock were different at various distances from the brick kilns in an agroforestry system in various tehsils of Mandi-Bahauddin district (Table 5). Maximum tree carbon stock (133 Mg C ha−1) was noticed in Mandi-Bahauddin followed by Malakwal (62.6 Mg C ha−1), while the minimum was recorded in Phalia (9.97 Mg C ha<sup>−</sup>1) in (Table 5).


Mean data from 3 locations of brick kilns.

#### *3.4. Total Soil Organic Carbon Stock*

Notable variations were observed for soil organic carbon at various distances from the brick kilns in the agroforestry system of various tehsils in Mandi-Bahauddin district (Table 6). Maximum total organic carbon of soil was noticed in Malakwal 10 km away from the brick kilns compared with control (near brick kilns). The highest stock of C in soil (15.30 Mg C ha<sup>−</sup>1) in Mandi-Bahauddin was noticed near the brick kiln, and minimum stock of soil C (9.06 Mg C ha−1) was measured 5 km away from the brick kiln, and an average organic C stock of (11.6 Mg C ha−1) was noticed. In Phalia maximum measured soil C stock (13.58 Mg C ha−1) was shown near the brick kiln and minimum soil C stock (11.6 Mg C ha−1) was shown 5 km away from the brick kiln, and average organic carbon of (12.5 Mg C ha−1) was noticed. In Malakwal, maximum soil C stock (15.3 Mg C ha−1) was measured 10 km from the brick kiln and minimum soil C stock (3.12 Mg C ha<sup>−</sup>1) was shown 5 km away, with an average organic C stock of (10.3 Mg C ha<sup>−</sup>1) (Table 6).

**Table 6.** Average Soil Organic Carbon Stock.


#### *3.5. Total Organic Carbon Stock in Potential Crops*

All the values regarding the carbon stock in Triticum aestivum and Oryza sativa at various distances were different in the agroforestry system of cultivated area for various tehsils of Mandi-Bahauddin district (Table 7). Maximum carbon stock was recorded in Triticum aestivum in Tehsil Mandi-Bahauddin (2.08 × 106 Mg C) followed by Malakwal (1.85 × <sup>10</sup><sup>6</sup> Mg C). Minimum carbon stock was noticed in wheat crop in Phalia (1.83 × 106 Mg C). For rice crop, maximum carbon stock was noticed in Mandi-Bahauddin (1.97 × 106 Mg C), followed by Phalia (1.83 × <sup>10</sup><sup>6</sup> Mg C), while minimum carbon stock in rice was recorded in Tehsil Malakwal (1.48 × <sup>10</sup><sup>6</sup> Mg C) (Table 8).



ABG = Above ground biomass; BGB = Below ground biomass; Mean data from 3 locations of brick kilns.


**Table 8.** Measured Carbon stock in *Oryza sativa*.

ABG = Above ground biomass; BGB = Below ground biomass; Mean data from 3 locations of brick kilns.

#### **4. Discussion**

This study provides the estimation of potential C pools for agroforestry systems in relation to brick kilns in Mandi-Bahauddin district, Pakistan. The agroforestry systems observed in the research site are crucial for the livelihoods of local farmers, as they have both commercial and subsistence production values [36,37]. In addition to timber from tree species, field crops such as wheat also have good market prospects at local and international level, as these are widely used for consumption [38,39]. Although not intrinsically C dense compared to systems such as forests or intensively managed agroforestry systems and pastures, C storage in agricultural farms can be increased by 20.4 to 21.4 t C ha−<sup>1</sup> globally [40,41] through the incorporation of long-living, deep-rooted trees [42,43]. While climatic attributes are consistent across the district where sampling was performed, the amount of C sequestered varied, largely depending on the tree species' distribution, density of tree species, basal area of tree, age of tree, area under crops and distance of agroforestry system from the brick kiln, emphasizing the importance of management decisions in determining carbon stocks. For example, district total mean tree carbon stock was the lowest (9.97 Mg ha−1) in Phalia, 5 km away from the brick kiln and the highest (133 Mg ha−1) in Mandi-Bahauddin, 10 km away from the brick kiln, among all tehsils in the study area. This appears to be related primarily to the level of tree stocking in the district, with Mandi-Bahauddin having the highest average basal area of all tehsils (8.9 m2 ha−1), and Malakwal having the lowest tree basal area (3.2 m2 ha<sup>−</sup>1). The current study revealed that maximum carbon sequestration was noticed in *E. camaldulensis* while the minimum was observed in *M. alba.* Among various tree species, the difference in biomass might be due to numerous aspects, i.e., number of trees/ha, tree age, quality and location of site, cultural practices, techniques and system of planting and ecological conditions in that area [21–24]. The average stored C stock of agroforestry systems in our study area was 11.46 Mg C ha−<sup>1</sup> in soil, 68.5 Mg C ha−<sup>1</sup> in trees, 60.4 Mg C ha−<sup>1</sup> in wheat and 54.2 Mg C ha−<sup>1</sup> in rice. The carbon stocks of simple systems (i.e., combination of single trees with cash crop or grass) and mixed-tree systems are similar to C stock on agroforestry systems reported by other studies in Indonesia [44,45]. Tree biomass accumulation representing the value of tree basal area has a correlation with C stock value [46]. The tree species having a higher basal area have the capacity for higher biomass accumulation which results in higher C stocks [47]. Our findings reflect that the plots near the brick kilns showed more growth of trees and potential crops (wheat and rice) in the study area, due to shifting of brick kilns to zig-zag technology that reduced dust emission. Similar patterns were reported in other studies [21,36,37,48,49].

In a terrestrial ecosystem, soil is a very notable system for CO2 mitigation. Many agroecologists have revealed that the soil carbon pool has prime importance in agroforestry systems [50,51]. The outcomes of the current investigation regarding SOC supports the hypothesis that SOC contents are maximum at 0–15 cm layers, with increased higher buildup of tree litter [52]. Higher soil carbon contents in agroforestry systems largely depend on the quality and amount of biomass input by non-tree and tree components of the system [53,54]. Moreover, vegetation detritus and litter from pruning under proper agroforestry management return a greater amount of organic carbon to the soil [55,56] (Figure 3). In this study, the average SC pools ranged from 10.3 to 15.3 Mg C ha−<sup>1</sup> 10 km away from the brick kilns. Kimaro et al. [57] reported that a significant variation was noticed in carbon sequestered

by an agroforestry system with legume trees, compared with a mono-crop system. Carbon storage in plants can be high in complex agroforestry systems and productivity of field crops depends on several factors such as age, structure and the way the systems are managed [58,59]. The results are comparable with the findings of other studies [42,60–62], which reported that agroforestry can store carbon in the range of 12–228 Mg ha<sup>−</sup>1. In our study it was noticed that maximum tree C stock potential (133 Mg C ha−1) was observed in Mandi-Bahauddin followed by Malakwal (62.6 Mg C ha<sup>−</sup>1), 10 km away from the brick kilns, while the minimum value was observed in Phalia (9.97 Mg C ha−1), 5 km away from brick kilns. In our research, the carbon sequestration potential of an agroforestry system (soil, trees and crops) showed differential response in relation to brick kilns (5 and 10 km away from brick kilns). Similar observations were reported by Gera et al. [63] who claimed that the variations in the carbon sequestration potential relate to the mean annual increment, which varies with site, age, density and plantation, as well as with the quality of planting stock.

**Figure 3.** Schemes of carbon sequestration in agroforestry.

Various factors such as species, land use type, cultural practices and CO2 supply play important roles in C stock and C sequestration rate [64]. Most farmers planted trees on their farmland for a short rotation. After regular intervals, harvesting of cultivated trees results in loss of C, but C is again stored when the harvested wood is converted into plywood, packaging materials, poles and the manufacturing of furniture [65]. Moreover, more accumulation of biomass was observed at the early stages and this decreased with the passage of time and age [66]. Tree stem sequesters the carbon for a longer time after felling as compared to the carbon stored in leaves and branch biomass [59]. These results highlight both the current and potential carbon sequestration potential of agroforestry in Pakistan, and can be further used in devising strategies for implementing tree planting programs on agricultural land and designing future carbon sequestration projects in Pakistan.

#### **5. Conclusions**

For the development of better management strategies, understanding of the influence of potential trees on farmlands in linear plantations is critical. Our intensive sampling in three tehsils showed that agroforestry systems in Punjab, Pakistan, currently store moderate amounts of carbon in plants and soil. In the agroforestry system, the increased soil organic carbon was due to litter fall and gave higher monetary returns in terms of more C stock 5 and 10 km away from the brick kilns. The decreasing order in terms of C stock potential of crops was wheat > rice. Variation was observed for C stock in crops (wheat and rice) at various distances from the brick kiln at the sampling sites. More C stock in the crops was noticed in Mandi-Bahauddin due to the maximum cultivated area, followed by Phalia and Malakwal. Crops had the maximum potential to store carbon dioxide at different sampling points, compared to trees and soil. The descending order in terms of C stock potential of agroforestry was crops > trees > soil. The findings of this study suggest that planting of tree species along with farm crops is a sustainable way to mitigate climate change by sequestering large amounts of carbon from the atmosphere. However, future studies should be conducted to highlight more indicators associated with the operation of the brick kilns. Given appropriate incentives, Punjab's farmers could help Pakistan meet its commitments to the Paris Climate Accord through reasonable changes in tree planting on existing agroforestry systems.

**Author Contributions:** Conceptualization, Q.u.Z., G.Y. and Y.C.; methodology, N.K., I.T. and S.N.; software, G.Y.; formal analysis, N.K.; investigation, Q.u.Z.; resources, N.K.; data curation, K.A.; writing—original draft preparation, N.K., M.A., M.W. and A.B.; writing—review and editing, Q.u.Z. and Y.C.; supervision, Q.u.Z.; project administration, Q.u.Z. and I.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not Applicable.

**Acknowledgments:** We are thankful to Umair Riaz (Scientific Officer), Soil and Water Testing Laboratory for Research Bahawalpur, Punjab Government, Pakistan for providing the research facilities for soil analysis.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

