How the Transition to Organic Grain Effects Biological Indicators of Soil Health

Treatments

The experiment includes three or four treatments (depending on the location) that follow a continuum of soil disturbance and cover crops. The treatments are: 1) Perennial alfalfa-grass hay, untilled; (UT) 2) Minimum-till corn-soybean-wheat rotation with precision-zoned high biomass diverse cover crops (MT); 3) Reduced-till corn-soybean-wheat with high-biomass cover crops (RT); and 4) Traditional full-tillage organic soybean-corn, with simple cover crops and input substitution (FT). 

Study sites

In January 2020, field-scale plots were developed in four different locations in Maryland, including both university experiment station land and commercial farms. Three locations are on two different soil types that reflect the coarser and finer textured soils found in Delmarva's grain cropping coastal plain area; the fourth has loamy soils in the hilly Piedmont region. In the commercial farms, three transition systems (MT, RT, and FT) are being compared, and four systems are being studied in experiment stations (NT, MT, RT and FT). The systems are being tested in randomized complete block experiments, with the treatments repeated four times on each farm. The individual treatment plots range in size from 30 ft x 300 ft to 60 x 150 ft, large enough to allow commercial scale farm equipment to be used. In each treatment, two representative and relatively uniform sampling areas for collecting soil and plant samples have been demarcated.

Sample collection

In early April 2020 the transition to organic farming began at all four sites by planting a uniformly-managed spring oats crop that was harvested in July, 2020. In July-August 2020, composite soil samples were collected in each half of each treatment plot using a 3.1 cm diameter x 30 cm long cores which were subdivided into 0-10, 10-20 and 20-30 cm depth increments. This set of soil samples represent the baseline soil condition before establishment of the four transition treatments. These baseline soil samples were tested for the physical and chemical properties. Subsamples of the baseline soil samples were saved unground and frozen at -20 C for the possibility of soil microbial DNA extraction. After beginning the NE-SARE funded project, soil samples were collected during the fall of 2022. For each treatment, ten cores were taken from a single experimental replicate and composited. Cores were taken at least 2 m from the edge of the plot and distributed throughout the whole plot. Cores were taken using a 3.1 cm soil probe. Debris was removed from the soil surface and the probe was inserted to a depth of 30 cm. The column of soil was divided into the depths of 0-10 cm, 10-20 cm, and 20-30 cm. The soil samples from each depth were composited for each treatment.  The core or loose soil was wrapped tightly in the bag, sealed and immediately placed on ice. Gloves, probe, and knife were wiped clean and sterilized with 70% isopropyl alcohol spray between treatments. Samples were transported to the laboratory in a cooler and frozen immediately upon arrival.  The baseline soil samples collected in fall 2020 and the samples from  2022 are being used in microbial analysis. 

Microbial analysis

The lead graduate student has been trained for DNA extraction and quantification at the UMD/ENST microbial ecology lab. Total microbial DNA from each frozen soil sample are being extracted using 250 mg of fresh soil mass using the DNeasy PowerLyzer PowerSoil Kit (QIAGEN Sciences Inc, Hilden, Germany) following the manufacturer’s recommendations. Subsequently, DNA extracted from each sample are being quantified using a Qubit 2.0 fluorometer (Invitrogen) and stored at -20 °C. When needed, DNA concentration for each extraction are being diluted to 5ng/μl for PCR amplification. Until February 2024, 20 soil samples have been extracted and frozen at -20 °C. 

After  DNA extraction: For prokaryotic communities, a fragment of the 16S ribosomal ribonucleic acid (rRNA) gene will be targeted using the primers 515F/806R and fungal internal transcribed spacer (ITS) gene regions will be targeted using the primer ITS1f/ITS2r (Coporaso et al., 2011; White et al., 1990). Illumina sequencing will be performed according to the 16S metagenomic sequencing library preparation manual (Part number 1504423 rev. B, Illumina, Inc.). PCR cleanup and indexing will be carried out using AMPure CP beads (Beckman Coulter, Pasadena, CA) and Nextera XT 96 index kit (Illumina, Inc.), respectively. Samples will be pooled; the amplicon library size will be verified using a Bioanalyzer 2100 DNA chip (Agilent Technologies). The library will be quantified using qPCR and the final library will be run on an Illumina MiSeq using a 600-cycle v3 cartridge.

The R package dada2 (version 1.14.1; Callahan et al., 2016) will be used for processing fastq files returned from Illumina Sequencing (filter and trim, dereplication, sample inference, merging pair-end reads, and removing chimeras). Prior to dada2 workflow, sequences containing ‘N’ values (unreadable bases) will be pre-filtered and the forward and reverse primers will be removed using cutadapt (version 2.9; Martin, 2011). Taxonomy will be assigned to the chimeric-free sequence table with the dada2 native implementation of the naïve Bayesian classifier method. The SILVA (version 132; Quast et al., 2013) and UNITE (release 10.10.2017; Nilsson et al., 2019) databases will be used as reference sequences for the identification of prokaryotic 16S rRNA genes and fungal ITS gene amplicon sequence variants (ASV), respectively. The ASV tables will be converted into phyloseq objects using the R phyloseq package (McMurdie and Holmes, 2013). 

Before performing microbial community analysis, features with ambiguous phylum annotation (phylum classified as ‘NA’) and phyla with low-prevalence will be removed. In addition, ASVs available in less than 10% of the total samples will be discarded; low-performing samples with total sequence reads less than 1000 will also be discarded. To account for differences in sequencing depth across samples, samples will be rarefied by randomly sub-sampling the minimum number of sequence reads present in the sample.  

Statistical analysis

For Objective (1), the Shannon diversity index will be calculated as a measure of α-diversity using the R phyloseq package. The Bray-Curtis dissimilarity (β-diversity) will be calculated after transforming the count data into relative abundances. Non-metric multidimensional scaling (NMDS) ordination plots based on Bray-Curtis dissimilarity will be constructed to visualize differences in prokaryotic and fungal community structure between the treatments. The effect of treatments on microbial community structure will be statistically tested using the permutational multivariate ANOVA (PERMANOVA) with 1000 permutations via the adonis function of the R vegan package (Oksanen et al., 2019). To check the assumption of PERMANOVA, permutation-based tests of multivariate homogeneity of group dispersions will be performed using the betadisper and permutest functions of the vegan package. If the overall treatment effect is found significant based on PERMANOVA, pairwise treatment comparisons will be performed using the pairwiseAdonis package in R (Martinez Arbizu, 2019). 

The associations between the most abundant taxa and the microbial community structure will be investigated using the envfit function of the R vegan package, and the taxa that shows significant correlations will be fitted to the NMDS plots. The VennDiagram package in R will be used to construct venn diagrams to depict the number of unique and shared ASVs between soil samples for each treatment (Chen, 2018). 

For objective Objective (2), the potential metabolic functionality of the prokaryotic community will be predicted based on 16S rRNA gene data using the t4f (Tax4Fun) function of the R themetagenomics package and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Woloszynek et al., 2017). The relative gene abundances will be calculated by summing all KEGG orthologs associated to a specific metabolic pathway (C cycling, N cycling, etc) (Lüneberg et al., 2018). Similarly, FUNGuild will be used for functional guild and trophic mode assignments to the ITS gene data for the fungal community (Nguyen et al., 2016). For simplicity, we will classify the assigned fungal trophic modes into four categories: Pathotrophs (plant pathogens, fungal parasites), Saprotrophs, Pathotroph-Saprotrophs, and others (animal pathogens, symbiotrophs). As recommended by FUNGuild developers, only trophic modes assigned with ‘highly probable’ or ‘probable’ confidence ranking will be reported. Kruskal-Wallis and Dunn’s test will be used to statistically test the effect of treatments on microbial diversity and potential microbial functionality of the soil microbiota using the R rstatix package (Kassambara. 2020). Significant differences will be determined at α value of 0.05, unless otherwise stated. Objective (3) The physico-chemical parameters of the soil (such as pH, organic carbon, active carbon, bulk density, structural stability) measured from the larger study will be correlated with the bacterial diversity (OTUs number – number of identified operational taxonomic units). Also, a Pearson’s correlation analysis will be done to determine the interaction between the soil physico-chemical properties and identified microbial taxonomic groups (e.g. Order, Phyla), using R software. 

The data on microbial diversity and metabolic functionality for each of the treatments will be compared across sites and over the time points to assess the consistency in response of microbes. Correlation of biological parameters with the physico-chemical properties will be used to understand the possible contrast in microbial response. A multivariate CLUSTER will be used to distinguish differences in the overall community structure involving tens to hundreds of microbial taxa.  Because of the complex variables involved, a quality index will be developed to simplify the tracking and comparison of changes in soil health among the four treatments. Several biological quality functions will be chosen (such as the microbial diversity, nitrifiers, denitrifiers, plant-growth promoting rhizobacteria, etc) and corresponding indicator properties will be assigned functional scores. The t-scores for each variable will be weighted and then averaged to give a relative quality index. 

The ongoing DNA extraction step is expected to complete in March 2024. This will be followed by amplification of the extracted samples using qPCR and sequencing. qPCR and sequencing steps are expected to be completed in May 2024. Data analysis and completion of final report is planned to be finished by August 2024.