Site Description
To examine the potential for a cover crop to improve soil health and plant production on dredged sediments, we manipulated the presence of a cover crop in plots at the Great Lakes Dredged Material Center for Innovation (GLDMCI) located in Toledo, Ohio (41.6700354oN, -83.5029711oW; Rúa et al. 2023). The GLDMCI was created between 2016 and 2017 by pumping ~ 8,800 m3 of dredged sediments from the Toledo Harbor and the Maumee River into four ~ 1 ha plots (Hull and Associates 2018). Dredged sediments in each plot were allowed to dry and acclimate to the on-land environment for two years prior to the experiment (Rúa et al. 2023).
Plot Preparation and Planting
In November 2018, the two available 1 ha plots were mowed and disked to remove existing vegetation. Then one plot was broadcast seeded with the common Midwestern cover crop winter rye (Secale cereale; The CISCO Companies, Indianapolis, IN, USA; ‘cover crop’) while the second plot was left fallow (‘control’). Approximately 43 days following the last freeze (CFAES 2022), both plots were tilled such that the existing vegetation was sowed into the dredged material prior to sowing rows of corn (variety: W2903DP; Wellman Seeds Inc., Delphos, OH, USA) via tractor on 15 May 2019. To eliminate weedy vegetation, halfway through the growing season (8 July 2019), both plots were sprayed with the herbicide Glyphosate with a 53.8% concentration at a rate of 0.0004 L/m2 (1.5 qt per acre; Buccaneer 5; CommoditAg, Effingham, IL, USA). To mirror local farming practices, plots were also fertilized with 48% ammonium sulfate (AMS; CommoditAg, Effingham, IL, USA) on 8 July 2019 at a rate of 0.0002 kg/m2 into 75.7 L of water per 0.4057 ha (1.5 lbs into 20 gallons of water per acre).
Soil Sample Collection
To assess differences in soil properties due to the use of a cover crop, we collected soil samples for analyses of physical, biological, and chemical properties from 10 sampling locations per plot at three collection times throughout the experiment: following site prep but prior to cover crop planting (12 November 2018, ‘pre-cover crop’), after completion of the cover crop life cycle / prior to corn planting (12 April 2019; ‘post-cover crop’) and following final corn harvest (12 October 2019; ‘post-harvest’). Soil samples from 10 haphazardly chosen locations were collected per plot (20 samples total), and each location was marked for future sampling. Samples were transported within six hours of collection to the Rúa lab at Wright State University in Dayton, OH, for further analyses.
Upon arriving at the Rúa lab, a bag containing ~ 750 g of soil from each sampling location (n = 20) was stored at -20 oC for up to three months prior to analysis for physical and chemical properties. In addition, one 2 mL tube per sampling site (n = 20) was stored on dry ice during transport and then stored at -20 oC for up to six months until analyzed for enzymatic activities. Finally, for post-harvest sampling only (12 October 2019), we also collected one soil core (11.5 cm depth x 10.5 cm diameter) from each sampling location and stored the core at room temperature for 18 hours until processing for bulk density measurement.
Soil Physical and Chemical Properties
The soil physical properties, including soil texture, BD, and gravimetric moisture content, were quantified to assess the effects of cover crop application. To assess soil texture, we measured the % sand, % silt, and % clay within each soil sample using the LaMotte soil texture kit (LaMotte, Chestertown, MD, USA). Gravimetric water content (‘moisture’) was quantified using the mass ratio such that the samples were weighed wet, dried for 48 hours at 105 oC, reweighed dry, and then the ratio of dry to weight was taken to determine the percent of water in each sample (Woods et al. 2019).
To measure the BD of soils, we multiplied the ratio of the volume of the soil core (‘field moist volume’) to the weight of the soil in the core (‘field moist mass’) by the ratio of the subsample wet mass to the subsample dry mass from the moisture calculation to calculate the bulk density of the sample (Eq. 1; Onufrak et al. 2019). As soil cores were not used for the first and second collection times, we acquired field moist volume using a 250 mL beaker as the ‘core’, filled it with the soil sample, took the final weight (‘field moist mass’) and used these values to calculate the bulk density of the sample (Eq. 1). To determine any differences between the core and modified methods, we ran the modified method using the same soils from the core method for three samples per treatment of the third collection time and compared the values obtained per method. The side-by-side comparison resulted in a difference of 0.2652 g/mL between the methods that was then used as a correction factor for the first and second collection times.
$$g/{ml}^{-1}=\left(\frac{field moist mass \left(g\right)}{field moist volume \left(mL\right)}\right)*\left(\frac{subsample dry mass \left(g\right)}{subsample wet mass \left(g\right)}\right)$$
Equation 1. Soil bulk density calculation.
We measured several soil chemical properties, including total N, total P, K, magnesium (Mg+), calcium (Ca+), CEC, pH, and % OM content. Soil pH was assessed using a 1:2 soil to water ratio and was measured using a Fisherbrand™ accumet™ AB15 Basic and BioBasic™ pH / mV / oC meter (ThermoFisher, Waltham, MA, USA; Woods et al. 2019). Three composite samples from each sampling location per plot (Sample A = 1, 2, 3; Sample B = 4, 5, 6; Sample C = 7, 8, 9, 10) for each collection time (n = 9) were sent to A & L Great Lakes Laboratories (Fort Wayne, IN, USA) for CEC, K, Mg+, Ca+, and OM analysis following EPA recommended laboratory protocols. Total N and total P were measured by the Midden Lab at Bowling Green State University (Bowling Green, OH, USA) using alkaline persulfate digestion with automated colorimetric detection of N and P using an AQ2 + discrete chemical analyzer (Seal Analytical, Mequon, WI, USA), as recommended by the EPA and US Geological Survey (Francy et al. 2020).
Bulk Soil Enzyme Activity
We quantified activities of six extracellular enzymes, including the fluorometric enzymes of alkaline phosphatase (AKP), β-glucosidase (BG), and leucine aminopeptidase (LAP), as well as the colorimetric enzymes of peroxidase (PER), polyphenol oxidase (PPO), and urease (UR). All enzyme assays followed protocols modified from Woods et al. (2019) and Sinsabaugh et al. (2000) and used a homogeneous soil slurry with 0.25 g dry weight of soil sample and 31.5 mL of 50 mM sodium acetate buffer (pH 5.6), except for AKP which required 31.5 mL of NaOH Tris buffer (pH 11.0). Soil slurries (‘soil homogenates’) were stored in the dark at 4 oC for up to one week prior to analysis. All enzymes were assayed in triplicate alongside a buffer blank, substrate blank, soil homogenate blank, and soil homogenate and substrate mix using a 96-well plate. Fluorometric enzyme analyses also included a standard curve for each sample.
For AKP and BG assays, the fluorescence indicator Methylumbelliferone was used in the form of the substrates 4-Methylumbelliferyl phosphate and 4-Methylumbelliferyl beta-D-glycopyranoside. For LAP, assays were conducted with the substrate L-Leucine-7-amid-4-methylcoumarin hydrochloride and 7-Amino-4-methylcoumarin. Once all reagents and samples were added to the 96-well plate, the prepared plates were left in the dark to incubate at room temperature (22 oC) for the following times: AKP, nine hours; BG, one to eight hours; and LAP, 96 hours. Following incubation, a NaOH stop solution was added to all wells, and the samples were incubated at room temperature for 10 minutes before reading at 360 nm excitation and 450 nm emission on a BioTek Synergy HT microplate reader (Agilent, Santa Clara, CA, USA). Using the dilution factor (DF = mL volume of buffer / g of soil sample), extinction coefficient \((\text{Ɛ}\); slope of the sample’s standard curves), and the calculated net fluorescence units (NFU; Eq. 2), all fluorometric enzyme activities were expressed as µmol h− 1 g− 1 using Eq. 3.
$$NFU =SH -SHB-B*\left(\frac{slope of standard with sample}{slope of standard with buffer}\right)$$
Equation 2. Net Fluorescence Unit (NFU) calculation for fluorometric enzymatic activities where ‘SH’ is the soil homogenate + substrate fluorescence, ‘SHB’ is the soil homogenate blank fluorescence, and ‘B’ is the substrate blank fluorescence.
$${\mu }\text{m}\text{o}\text{l} {h}^{-1}{g}^{-1} =\frac{\left(\frac{\text{N}\text{F}\text{U}}{\frac{\text{Ɛ}}{0.25}\text{m}\text{L}}\right)}{\begin{array}{c}\left(incubation time*\left(\frac{1}{DF}\right)*0.200 mL\right)\end{array}}$$
Equation 3. Fluorometric enzymatic activity calculation using the Net Fluorescence Unit (NFU) and Dilution Factor (DF) for alkaline phosphatase, β-glucosidase, and leucine aminopeptidase.
The colorimetric enzymes PER and PPO were tested using 25 µM 3,4-Dihydroxyphenylalanine with the addition of 3% hydrogen peroxide for PER. Plates were incubated in the dark for 24 hours at room temperature (22 oC). The plates were read at 450 nm emission on a Molecular Devices Corporation SpectraMax 190 microplate reader (Molecular Devices Corporation, Sunnyvale, CA, USA), and activity in µmol h− 1 g− 1 was calculated using net absorbance units (NAU = soil homogenate absorbance – soil homogenate blank absorbance – substrate blank absorbance) and the extinction coefficient \(\left(\text{Ɛ}\right)\) 1.8446 for PER and 2.4942 for PPO (Eq. 4).
$${\mu }\text{m}\text{o}\text{l} {h}^{-1}{g}^{-1} =\frac{\text{N}\text{A}\text{U} *Buffer volume \left(ml\right)}{\left(\text{Ɛ}\text{*}0.200 \text{m}\text{L}*Incubation Time \right(hr\left)*sample \right(g)}$$
Equation 4. Colorimetric enzyme activity calculation using the Net Absorbance Unit calculation (NAU) for peroxidase, polyphenol oxidase, and urease.
UR activity was tested using 10 µl of a 400 mM Urea substrate with an incubation time of between two and seven hours at room temperature (22 oC) in the dark before adding 40 µl of both salicylate and cyanurate acid and a one-hour incubation time after the addition. We read the plates at 610 nm on a Molecular Devices Corporation SpectraMax 190 microplate reader (Molecular Devices Corporation, Sunnyvale, CA, USA). We used net absorbance units (NAU = soil homogenate absorbance – soil homogenate blank absorbance – substrate blank absorbance) and the extinction coefficient \(\left(\text{Ɛ}\right)\) of 0.2403 to calculate UR activity in µmol NH4 g− 1 h− 1 (Eq. 4).
Plant Responses
We assessed the response of corn to the cover crop treatment by measuring plant height, leaf count, reproduction status, photosynthetic efficiency, total (above and below ground) biomass, and final yield for 10 focal plants located at the previously sampled soil locations in each plot (total of 20 plants). We measured leaf count and plant height weekly throughout the growing season, starting from plant emergence on 24 May 2019 and ending with the final harvest on 12 October 2019. Plant height (cm) was assessed using a tape measure from the ground to the top of the arch of the highest fully formed leaf (APHIS USDA 2020). Leaf count was measured following the guidelines from the USDA (2020) by counting only the fully formed leaves with an arch. Leaf count was also used to indicate vegetative growth stages (i.e., for V6, the corn plant had six fully developed leaves). Reproductive growth stages of tassel (VT) and kernel development (R1-R6) were assigned when applicable (Nafziger 2017).
Photosynthetic efficiency, which evaluates the efficiency of a plant’s photosystem II to capture sunlight and transfer energy, was measured twice during the growing season, once during the vegetation stage (21 June 2019, day 63 of growth) and again during kernel production (31 August 2019, day 110 of growth). We measured photosynthetic efficiency (FVM) using an OS-30 chlorophyll fluorometer (Opti-Sciences, Hudson, NH, USA) following the procedure outlined in Friedman et al. (2020). Briefly, we placed a closed clip on two different leaves, facing the top of the leaf, and allowed the clipped area to equilibrate for 30 minutes prior to light exposure. We then exposed the leaf enclosed by the clip to a low light emission followed by a high light emission and recorded the FVM for each clip per plant. Photosynthetic efficiency was then calculated by averaging the FVM of both clipped leaves for each plant.
To assess yield, we quantified number of ears, kernel count per ear and per plant, the weight of kernels per ear (g/ear) and per plant (g/plant), and total plant biomass (g). Above and below ground biomass were measured twice during the experiment, once to mimic silage harvest (31 August 2019) and once at final harvest (12 October 2019). For silage harvest, plants of similar height and leaf count located near the focal plants were used, while the final harvest plants were the focal plants themselves. For both harvest collections, plants were hand harvested, and the roots were disconnected from the stalk and cleaned to remove excess soil. All plant material was then dried for 48 hours at 60 oC, and the above and below ground plant material was weighed separately (g). Ears were counted on each corn plant at final harvest. Next, kernels were hand removed from the cob, counted for each ear, and totaled for each plant. Finally, kernels were weighed (g) per ear and per plant.
Statistical Analysis
To determine if the application of a cover crop improved dredged sediment soil property and corn responses compared to a plot left fallow, we ran univariate and multivariate statistical analyses. All statistical analyses were performed in the statistical programming environment R, version 4.2.1 (R Core Team 2022). Data were normalized through log transformation as needed to match model assumptions. All data visualizations were created using the ggplot2 (Wickham et al. 2022) and ggbiplot (Vu 2011) packages.
We created separate linear models for each soil property metric to quantify changes in soil physical properties, chemical properties, and enzymatic activities due to adding a cover crop compared to a plot left fallow. We created linear models using the lm() function in base R for % OM, CEC, Ca+, K, and Mg+ since these properties were assessed using only three composite samples per plot per collection period. We created linear mixed effects models using the lme() function in the nlme package (Pinheiro et al. 2021) for % sand, % silt, % clay, and pH and the log transformed data for BD, soil moisture, and activities for AKP, BG, LAP, PER, PPO, and UR with a random effect in each model for unique plant ID. All soil property models were created such that each soil property metric was a function of the interaction of treatment (cover crop vs. control) and collection time [pre-cover crop (12 November 2018), post-cover crop (12 April 2019), post-harvest (12 October 2019)]. We then tested each model with an analysis of variance (ANOVA) using the anova() function in base R. In cases with significant relationships between response variables and interaction terms, we performed post-hoc analyses using the emmeans() function in the emmeans package (Lenth et al. 2022), where all pair-wise comparisons a Tukey’s p-value adjustment method. To determine if soil property responses differed between groupings of treatments (cover crop vs. control), we used separate principal component analyses (PCAs) for physical, chemical, and enzyme properties were also created with the prcomp() function in base R to condense changes in properties into two linear principal components. We performed an analysis of similarity (ANOSIM) for all PCAs using the adonis2() function in vegan package (Oksanen et al. 2022).
To quantify changes in corn vegetative development between plots with or without a cover crop application, we created separate linear models for plant stage, height, leaf count, relative growth rate, and root to shoot ratio. We evaluated differences in development time by week due to cover crop use by creating a linear model using the lm() function in base R, using week of plant data collection as a function of the interaction between treatment (cover crop vs. control) and plant stage (VE - R6).
To assess differences in plant height over time, we created a linear model using the lme() function in the nlme package (Pinheiro et al. 2021), where height was a function of the interaction of treatment (cover crop vs. control) and week of data collection with unique plant IDs was a random effect. We then tested our model with repeated measures ANOVA using the anova() function in base R, with post-hoc analyses using the emmeans() function in the emmeans package (Lenth et al. 2022) where all pair-wise comparisons used a Tukey’s p-value adjustment method. We created a generalized linear mixed effects model with a Poisson distribution using the glmer() function in the lme4 package (Bates et al. 2022) to evaluate leaf count as a function of the interaction between treatment (cover crop vs. control) and week, with individual plant ID as the random effect. We tested our model with an ANOVA using the Anova() function in the car package (Fox et al. 2022) and post-hoc analyses for significant interactions using the emmeans() function in the emmeans package (Lenth et al. 2022) where all pair-wise comparisons used a Tukey’s p-value adjustment method.
To quantify changes in plant matter allocation, we calculated the ratio of roots (below ground biomass) to shoots (above ground biomass). We then created a linear model using the lm() function in base R, where the root to shoot ratio was a function of treatment (cover crop vs. control) and was evaluated for both silage and final harvests. We then tested our models with an ANOVA using the anova() function in base R. Significant interactions were subjected to post-hoc analyses using the emmeans() function in the emmeans package (Lenth et al. 2022) where all pair-wise comparisons used a Tukey’s p-value adjustment method.
Finally, we calculated the relative growth rate (RGR) by week for each plant by taking the height of the current week, subtracting it from the height in the previous week, and dividing that value by the height in the previous week (Eq. 5). We then evaluated changes in weekly RGR with a linear mixed effects model using the lme() function in the nlme package (Pinheiro et al. 2021), where RGR was a function of the interaction between treatment (cover crop vs. control) and current week with individual plant ID as the random effect. We then tested our model with an ANOVA using the anova() function in base R and significant interactions were subjected to post-hoc analyses using the emmeans() function in the emmeans package (Lenth et al. 2022), where all pair-wise comparisons used a Tukey’s p-value adjustment method.
$$RGR for Week n=\frac{Height (Week n+1)-Height \left(Week n\right)}{Height \left(Week n\right)}$$
Equation 5. Calculation for weekly relative growth rate (RGR).
We assessed photosynthetic efficiency (FVM) in both vegetative (V5-V7) and reproductive stages (R3) by creating a linear mixed effects model using the lme() function in the nlme package (Pinheiro et al. 2021), where photosynthetic efficiency (FVM) was a function of treatment (cover crop vs. control) and plant stage (vegetative vs. reproductive) with the individual plant ID as a random effect. We again tested our model with an ANOVA using the anova() function in base R and significant interactions were subjected to post-hoc analyses using the emmeans() function in the emmeans package (Lenth et al. 2022), where all pair-wise comparisons used a Tukey’s p-value adjustment method.
We quantified differences between corn reproductive development when planted on a plot with a cover crop application compared to a plot left fallow by assessing final ear count, kernel count per plant, and kernel weight per plant. For final ear count, we created a generalized linear model with a Gaussian distribution that better captured ears where only half the kernels developed, using the glm() function in base R, where ear count was a function of treatment (cover crop vs. control). Final kernel count was assessed with a generalized linear model with a Poisson distribution using the glm() function in base R, with kernel count as a function of treatment (cover crop vs. control). Differences in final kernel weight were assessed with a linear model using the lm() function in base R, where kernel weight was a function of treatment (cover crop vs. control). All corn reproductive models were tested with an ANOVA using the anova() function in base R. We then performed post-hoc analyses using the emmeans() function in the emmeans package (Lenth et al. 2022), where all pair-wise comparisons used a Tukey’s p-value adjustment method for all significant interactions.
Furthermore, since corn yield can encompass both kernel production and the whole plant (silage harvest), yield was also captured with above ground, below ground, and total plant biomass measurements. To evaluate differences in biomass, we created three separate linear models using the lm() function in base R, with biomass (above ground, below ground, or total) as a function of treatment (cover crop vs. control) and these models were used for both silage and final harvest. All biomass models were tested with an ANOVA using the anova() function in base R, and any significant interactions were tested using post-hoc analyses using the emmeans() function in the emmeans package (Lenth et al. 2022) where all pair-wise comparisons used a Tukey’s p-value adjustment method.
Finally, we assessed the relationship of crop responses associated with production (biomass, kernel count) with individual soil properties. We created linear models using the lm() function in base R, where total biomass was a function of a single soil property (Physical: % sand, % silt, % clay, % OM, BD, moisture; Chemical: pH, CEC, N, P, K, Mg+, Ca+; Enzyme Activity: AKP, BG, LAP, PER, PPO, and UR) and treatment (cover crop vs. control). We used generalized linear models using the glm() function in base R for kernel count with a Poisson distribution, where kernel count was a function of a single soil property. We then tested each model with an ANOVA using the anova() function in base R, and any significant interactions were assessed using the emmeans() function in the emmeans package (Lenth et al. 2022), where all pair-wise comparisons used a Tukey’s p-value adjustment method.