Full-Season Cover Crops and Their Traits That Promote Agroecosystem Services

Non-marketable crops are increasingly being used as a tool to promote agroecosystem services and sustainable agriculture. Nevertheless, crops vary greatly in the traits by which they capture resources and influence the local ecosystem. Here we report on the traits and associated soil microbial communities that relate to aboveground biomass production, nutrient capture, weed suppression, erosion control and building particulate organic matter of 22 different full-season cover crops. All agroecosystem services were positively correlated with maximum canopy height and leaf area. Rooting density was positively associated with indices of bacterial diversity. While some legumes produced the greatest standing N and P in aboveground biomass, they were also poor at capturing soil nitrate and promoted high levels of potential plant fungal pathogens. Conversely, Brassicaceae crops had the lowest levels of potential plant fungal pathogens, but also suppressed saprophytic fungi and rhizobia. Thus, not all crops are equal in their ability to promote all agroecosystem services, and while some crops may be ideal for promoting a specific agroecosystem service, this could result in a trade-off with another. Nonetheless, our study demonstrates that plant functional traits are informative for the selection of crops for promoting agroecosystem services.


Introduction
The growth of agricultural production has increased dramatically over the past century with the aid of intensifying management strategies that include increased use of synthetic fertilizers, pesticides and reduction of traditional small scale integrative cropping systems [1,2]. However, over the last few decades, there has been a growing realization that intensive agricultural practices aimed at achieving higher yields also have undesirable long-term trade-offs that compromise local biodiversity and the ability of agroecosystems to sustain other desirable services [3][4][5][6][7]. These agroecosystem services include ecosystem properties such as efficient nutrient cycling in soils, building and maintaining soil organic matter and weed control, all of which support agricultural production and sustainability to benefit human societies. Consequently, there is now a need to develop practices for the ecological intensification of cropping systems to "leverage natures technologies" [8][9][10].
Full-season cover crops, more aptly referred to as 'agroecological service crops', are crops used in rotations to sustain, recover, or enhance desirable, but non-marketable, ecological attributes of agroecosystems [11]. Although typically non-marketable, cover crops are a versatile agricultural tool applied to provide desirable agroecological services that might be lacking in an agroecological system. These service crops are often chosen for the benefit they might provide in promoting the yield of a following cash crop in a rotational sequence [12,13] Simultaneously, agricultural service crops have the potential to mitigate the environmental cost of agricultural activities [11,14]. At present a large variety of plant species of varying characteristics have been identified as potential cover crops. A multitude of ecosystem properties can be achieved using plant specific functional traits. For instance, legumes can provide organic sources of N through their N-fixing abilities, potentially displacing part of a crop's fertilizer N requirements [15][16][17][18]. Members of the Brassicaceae family (e.g., mustards) are known for their potential to suppress soil pathogens [19][20][21], and various crops in this family are used for their abilities to suppress invasive agricultural weeds [22][23][24][25]. The high variety and flexibility of cover crop applications provide the potential to allow them to be tailored to suit various agricultural production systems and to target specific agroecosystem services. The use of functional traits and identities has been proposed for selecting cover crops tailored to target and enhance specific ecosystem services [26][27][28][29].
Plants have evolved a wide variety of traits to capture resources, defend against predators and compete with neighboring species [30,31]. These plant 'functional traits' can also provide insights into how plants influence and respond to the local ecosystem beyond their taxonomic identity [32][33][34][35][36]. Specifically, leaf area, leaf mass per area, maximum canopy height and rooting traits are often traits associated with nutrient capture and competition against neighboring species [37][38][39][40][41], defense against herbivores [42][43][44][45], promotion of desirable soil microbial communities [46,47] and enhancement of soil erosion control [48][49][50]. Thus, it has been proposed that the use of plant functional traits could be insightful for predicting how crops may promote agroecosystem services [34,51]. However, the use of functional traits of crops to predict their influence on various agroecosystem services is still in its infancy and lacks the full breadth of potential crops; therefore, this needs further empirical study.
Here we assessed the effects of 22 full-season cover crops (agricultural service crops) that represent a wide variety of ecological attributes and functional traits on six ecosystem services: crop biomass production, nutrient capture (standing N and P and soil ammonium and nitrate at the end of the growing season), weed suppression (weed cover and biomass), soil building (soil particulate organic matter N and C, active carbon and soil aggregate stability), erosion control potential (ground cover and rooting density) and soil health (using various soil fungal and bacterial indices). Because crops vary greatly in their life strategies and functional traits, we expect that not all crops are equal in their abilities to promote various ecosystem services. This means that it is unlikely that a single crop species will be able to simultaneously maximize multiple ecosystem services, and, therefore, specific crops are likely needed to target specific ecosystem services. Under this expectation we aim to (1) identify cover crops and their associated traits that maximize specific agroecosystem services and (2) assess the effects of crop functional traits on the soil microbial communities they promote that are associated with enhancing these ecosystem services.

Experimental Design
We chose 22 crops to trial as full-season service crops (see Table 1) in order to cover a wide variation in plant functional traits. We categorized the crops by plant functional group as grasses, forbs, legumes and those belonging to the Brassicaseae family. We separated the Brassicaseae family out from the forbs as they are non-mycorrhizal plants with antimicrobial properties that are often used as a biofumigant against soil pathogens [19][20][21]. The experiment was located at the Fredericton Research and Development Centre, Fredericton, New Brunswick, Canada. The crops were sown as monocultures in 4 m × 4 m plots and replicated four times. Replicate plots were randomly arranged into four randomized complete blocks split between two separate fields for a total of 88 plots. Block 1 and 2 occurred in field 1 (45 • 55 12 N, 66 • 36 15 W) and blocks 3 and 4 were located in field 2 (45 • 54 51 N, 66 • 36 25 W). These two fields differed in their soil properties and these two fields were used to provided robustness in our results against background soil properties and histories. Field 1 was characterized by higher levels of phosphate (P 2 O 4 = 417 ppm) and potash (K 2 O = 92.5 ppm) and lower nitrate (NO 3 < 5 ppm), calcium (Ca = 558 ppm) and soil organic matter (2.5%) compared to field 2 (P 2 O 4 = 169 ppm; K 2 O = 45 ppm; NO 3 = 6.2 ppm; Ca = 1212 ppm, organic matter = 3.2%). The soil pH was 6.1 and 6.7 in field 1 and 2 respectively. Prior to our experiment, both fields were a perennial red clover (Trifolium pratense) and timothy (Phleum pratense) grassland. In 2019, pots were prepared by disc tilling the soil on 3 and 29 May for blocks 1 & 2 in Field 1 and 30 May and 4 June for blocks 3 & 4 in Field 2. All plots were then harrowed on June 17 following the addition of calcium ammonium nitrate at a rate of 34 kg N ha −1 and K Mag at a rate of 22 kg K ha −1 . Plots were sown on June 18 using a Wintersteiger seeder with a depth of 1-3" and 6" row spacing. Seeding rates were based on recommended rates for each of the crops (Table 1). No further fertilizer was applied, and no pesticides or herbicides were applied throughout the growing season in order to characterize the ability of the cover crop to promote ecosystem services without additional inputs.

Data Collection
All measurements taken, units and associated ecosystem services are listed in Table 2. Throughout the growing season, we measured the maximum canopy height, the ground cover of the crops, and the ground cover of weeds at 21, 30, 44 and 51 days after planting (DAP). Canopy height was determined using the average of three measures along a transect through the center of the plot. The spring relative growth rate (RGR) of the crops was then calculated as the log of the canopy height at 21 DAP per day. The maximum canopy height is the maximum height the crops achieved throughout the growing season. Ground cover of crops and weeds was visually estimated using a modified Londo decimal scale [52], where 0.5 <10% cover, 1 = 10-25%, 2 = 25-40%, 3 = 40-55%, 4 = 55-70%, 5 = 70-85%, and 6 = 85-100% cover. The average cover throughout the season was used as crop cover and ground cover. The ground cover of crops post-harvest/mowing re-growth was also recorded in the fall at 115 DAP as a measure of winter cover. All further mentioned data was collected by block after 55 DAP. Plant biomass was collected from a 0.25 m × 0.25 m quadrat that was randomly placed at least 0.25m from the plot edge. Plant biomass was harvested at 3 cm above the soil surface and sorted by sown crop and weeds (pooled). From each sown crop a sample of five leaves were taken to calculate leaf mass, leaf area and visually score the % leaf damage (in 10% ranges) of leaf mining, chewing and sucking insects as defined in [53]. However, these herbivory data were often 0, and on average, the total herbivory damage was <5%. For these reasons we omitted these herbivory measures from further analyses. Leaf area was measured using the WinFOLIA software and scanner (Regent Instrument, Québec, QC, Canada) on these leaves as well. Crop, weed and leaf samples were all dried at 65 • C for at least 48 h and weighed to determine their aboveground biomass. Leaf mass was used to calculate the leaf mass per area (LMA). The dried biomass of the crops was finely ground to <2 mm and used to assess the %N content of the crops using the prescribed method of the Vario MACRO elemental analyzer (Elementar, Langenselbold, Germany) and the % P content measured by the Olsen P method. The standing N and P of the crop in g m −2 was then calculated as the product of the % N or P and the biomass of the standing crop. Soil cores were also taken at the time of biomass sampling using a slide hammer soil sampler (7 cm diameter by 15 cm deep) from three randomly selected points in the plot, which were subsequently pooled and homogenized. A 500 mg sample of fresh soil was used to extract DNA for next-generation sequencing to characterize the fungal and bacterial communities (described below). The remaining sample was air dried and sieved through a 2 mm sieve to quantify soil ammonium and nitrate (KCl extracted) [54,55] and the particulate organic matter (POM) of carbon and nitrogen, which was measured in soil collected in a 53 µm sieve and run on an Vario MACRO elemental analyzer (Elementar, Langenselbold, Germany). Aggregate stability was determined by wet sieving using 4.0 g of air-dried soil aggregates of 1-2 mm size by shaking using an Eijkelkamp wet sieving apparatus. Samples were placed into a 250 mm sieve, gently moistened, and repeatedly immersed for 3 min in water. Aggregate fragments that passed through the sieve were filtered, dried, and weighed. Particles remaining on the sieve were repeatedly immersed in a 2 g/L NaOH dispersing solution for intervals of 5 min until there were only sand particles remaining. Soil aggregate stability was calculated as the percent of mass aggregates remaining minus the sand [56]. Active carbon (POX-C) was determined with the permanganate oxidizable carbon method [57] in duplicate samples using 2.5 g of air-dried soil mixed with 0.02 mol L −1 KMnO 4 , then shaken for 2 min at 240 rpm and allowed to settle. A 0.5-mL aliquot of supernatant was diluted in 49.5 mL of deionized water, and absorbance was measured at 550 nm on a Biochrom Libra S60 Spectrophotometer (Biochrom Ltd., Cambridge, UK). The absorbance of four standard solutions were used (0.00005, 0.0001, 0.00015 and 0.0002 mol/L KMnO 4 ). At the time of plant and soil sampling, we also measured rooting density (root mass per volume of soil) by taking an additional 7 cm diameter and 15 cm deep soil core around the focal crop at the center. Roots were washed clean of soil and dried at 65 • C for at least 48 hrs. Fields 1 and 2 were flail mowed at 66 and 80 DAP, respectively, and left as a green manure over winter.

Soil Fungal and Bacterial Community Characterization
Soil DNA was extracted using the FastDNA SPIN Kit for Soil (MP Biomedicals) and quantified on the Qubit 4 fluorometer with the high sensitivity sdDNA assay kit. DNA samples were sent to Genome Quebec for PCR amplification using the primers ITS1F [58] and ITS2 [59] to target the fungal community (ITS) in 25 µL PCR reactions using Qiagen HotStart Taq with a 52 • C annealing temperature for 33 cycles. The primers 341F and 805R [60] were used to target the bacterial (16S rRNA gene) community in 25 µL PCR reactions with the New England Biolabs Q5 HiFi polymerase with an annealing temperature of 60 • C for 25 cycles. Amplicons were sequenced using the Illumina MiSeq PE250 platform.The QIIME2 platform [61] v2020.8 was used to filter 16S and ITS sequencing reads to amplicon sequence variants (ASVs) and classify them taxonomically; DADA2 [62] was used for the filter step to trim, denoise, merge (minimum overlap of 12nt), and remove chimeric sequences, and the pseudo-pooling parameter was applied due to working with soil (diverse) samples. 16S reads were trimmed at the front ends according to primer length and at the tail ends once mean quality score declined below 35 (resulting in 1nt position trimmed). The ITS reads were first trimmed using Cutadapt [63] to remove primer sequences and read-through, before DADA2 was applied. The resulting 16S ASVs were assigned taxonomically using SILVA [64] v132 99% identity majority-taxonomy-alllevels taxonomy file and the corresponding rep-set-16S-only sequence file. The ITS ASVs were assigned to fungal taxonomy using UNITE v8 [65] developer dynamic-delineation taxonomy and sequence files. Mitochondria and chloroplast sequences were also removed. Sequence reads were not rarified [66]. Additionally, the sequencing depth was unrelated to the identity of the crop or their functional grouping (ITS: F 21, 56 = 1.35, p = 0.183 and F 3, 75 = 0.57, p = 0.635; and for 16S: F 21, 65 = 0.67, p = 0.845 and F 3, 83 = 1.91, p = 0.135 for crop identity and functional group respectively). Sequencing depth was found to be spatially dependent (total ITS and 16S rRNA gene ASV reads varied among blocks: F 3, 71 = 5.28, p = 0.002 and F 3, 80 = 22.67, p > 0.001 respectively and 16S rRNA gene read depth changed significantly with the spatial coordinates of the plot depending on the field: F 2, 80 = 7.21, p = 0.001). These spatial structures of the experiment are incorporated into our ANOVA models (see Section 2.4 Data Analyses below).
We assessed the abundance of key functional guilds of microbes in each sample using the log(x + 1) cumulative sum of sequence reads of taxa within a guild. Fungal taxa were assigned to trophic guilds using the FunGuild data base [67] and the "FUNGuildR" R package (https://github.com/brendanf/FUNGuildR/ (accessed on 4 February 2021)). The key functional fungal guilds were identified as arbuscular mycorrhizal fungi (AMF, phylum Glomeromycota), fungi identified as saprotrophs and plant pathogens. For bacteria, we identified plant symbiotic N-fixing Rhizobia as genera belonging to the family Rhizobiaceae and the genus Bradyrhizobium [68]. Nitrifying bacteria were characterized as those belonging to the Family Nitrosomonadaceae specifically the genera Nitrolancea, Nitrosospira, Nitrospira, Nitrosomonas, Nitrosococcus, Nitrosolobus, Nitrosovibrio, Nitrobacter, Nitrococcus, Nitrotoga, Nitrospina and Nitrolancetus [69][70][71]. Finally, we also assessed the abundance of Cyanobacteria (phylum) due to their important role in agroecosystems as free-living soil N-fixing bacteria [72][73][74]. We did not characterize bacterial denitrifiers or plant pathogens, as these functions are not specific to particular genera but are rather unique to specific bacterial taxa among various genera; our data does not have the necessary taxonomic resolution for this characterization. We also calculated the richness and evenness (Pielou's J') using the ASV log(x + 1) count data for both bacterial and fungal ASVs.

Data Analyses
All analyses were performed using R version 3.6.1 (R Core Team). All data were assessed using mixed-effects models that included the field (1 or 2) as a random intercept and the plant functional group identity or cover crop identity (22 levels) as a fixed term. The spatial correlation (x, y position of plots within a field) of residuals was also assessed by comparing models that included a spatial autocorrelation structure using either an 'AR1 or 'Spat' function vs. the model that excluded it using the maximum likelihood estimation (ML) with the R package 'nlme'. In nearly all cases the models that included a spatial autocorrelation provided a better model fit to the data and were included in the model. Homoscedasticity in the residuals was visually assessed and crop and weed biomass, standing N and P and maximum canopy height were improved by a square-root transformation. POM-C, crop and weed % cover, and rooting density were improved by log transformation.
We grouped plant and soil measurements into six ecosystem services: productivity, weed suppression, nutrient capture, erosion control, soil building and soil health (see Table 2). Soil building, we define as the levels of particulate organic matter C and N, active carbon and soil aggregate stability [75]. Although soil erosion was not directly measured in our experiment, we infer greater erosion control from greater root mass density, growing season and winter ground cover [76,77]. While soil health is a multifaceted term that may include various biotic and abiotic components, it is acknowledged that soil microbial diversity is a key component of soil health as greater soil microbial diversity and relative abundance of key desirable fungal and bacterial functional guilds, and fewer plant pathogens, underpin greater nutrient and carbon cycling and plant productivity [78][79][80][81][82][83].
Further, we wanted to group the above and belowground responses into separate nonoverlapping ecosystems services to avoid redundancy in data use and to avoid inherent correlations among the ecosystem service groupings. We first inspected the correlations among variables that were to represent each agroecosystem service. Since our index of soil health comprises 10 different fungal and bacterial indices, we also assessed the variation among crops using principal component analysis (PCA). The individual response measures representing each agroecosystem service were then z-transformed (unit variance and mean of 0) so that they are on a common standardized scale and where then averaged [81,84,85]. Weed biomass, weed cover, soil NH 4 , NO 3 and plant fungal pathogens were considered measures that indicate an undesirable condition and were multiplied by -1 prior to averaging to invert values such that greater positive values indicate a more desirable condition. While greater plant available soil NH 4 and NO 3 may be desirable at the beginning of the growing season, we associated the reduced nutrient crop uptake of these forms of N as an agroecological disservice as it may lead leaching of soil NH 4 and NO 3 in the fall post-harvest and/or spring during rain fall events and run-off [86][87][88]. To make our composite ecosystem service indices more intuitive for interpretation we scaled each ecosystem service between 0 and 1 where 0 is the observation with the lowest score and 1 is the maximum score for each ecosystem service. Ecosystem services were also assessed using mixed effects models as mentioned above for individual measures. We also correlated (Spearman) the above soil microbial indices, used as indicators of soil health, with plant functional traits and individual ecosystem services to provide insights as to whether certain crops traits are able promote key guilds of soil microbes and their link to ecosystem services. To assess the correlations and trade-offs among ecosystem services and how they group out with different crops we performed a PCA on the six ecosystem services.

Aboveground Productivity and Functional Traits
Aboveground biomass of the crops varied significantly among the functional grouping of the crops and their identity ( Table 2). Crops that produced a substantially greater aboveground biomass than other crops were buckwheat, oats and field peas ( Figure 1a). In contrast, the legumes were on average the least productive, while the forbs were the most productive followed by the Brassicaceae and grasses (Figure 1a). Functional traits of the crops were all highly dependent upon the crop identity and the functional group to which it belonged ( Table 2). Leaf area was greatest in the sunflower, followed by the two Brassicaceae oilseed and tillage radish (Figure 1b). Leaf mass per area (LMA) was greatest in the grasses, particularly Teff, sorghum-sudan grass as well as sunflowers (Figure 1c). Maximum canopy height was greatest in the forbs, particularly buckwheat, hemp var. anka, and sunflowers. Maximum canopy height was the lowest in the legumes; however, field peas, faba beans and hairy vetch achieved the greatest canopy height of the legumes (Figure 1d). The spring relative growth rate (RGR) of the plants during the first three weeks post seeding was greatest in oats, followed by buckwheat, field peas and annual ryegrass. The spring RGR was lowest for the legumes, apart from field peas, hairy vetch and faba beans (Figure 1e). Crop aboveground biomass production was significantly positively correlated with greater leaf area, maximum canopy height, rooting density and spring RGR (Figure 1f).  Table 2. Significant correlations (Spearman) between crop biomass production and plant traits are shown (f). The association with rooting density is not shown as it was unrelated to crop biomass (r = 0.148, p = 0.169).

Nutrient Capture
The standing N and P in the crop varied significantly among the crops, but not their functional groupings ( Table 2). The standing N was greatest in the field pea, followed by buckwheat and oats, and on average was greatest in the legumes, although not significantly (Figure 2a). The standing P was greatest in the field pea, followed by buckwheat and oats (Figure 2b). The ammonium (NH 4 ) concentration in the soil at the end of the growing season did not vary significantly among crops or their functional groups ( Table 2) but was observed to be highest in plots of oilseed radish, faba beans and sunflowers (Figure 2c).
However, plant available soil nitrate (NO 3 ) concentration at the end of the growing season varied significantly among crops and their functional groupings (Table 2), where plots with field peas, hairy vetch, and the legumes in general had higher levels of soil NO 3 at the end of the growing season (Figure 2d). The ecosystem service of nutrient capture, represented by the average of the standing N and P and the reduction of the mobile forms of plant available N (NH 4 and NO 3 ), did not significantly vary among crop functional groups, but did significantly vary among crops. The field pea scored the highest for nutrient capture followed by buckwheat and oats (Figure 2e). Crops with larger leaf area, maximum canopy height and a faster spring RGR were all significantly correlated to greater nutrient capture (Figure 2f).  Significance is indicated by † p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001.

Weed Suppression
The biomass and cover of weeds varied significantly among plant functional groups and the identity of the crop ( Table 2). The legumes resulted in the greatest, and the Brassicaceae and forbs the least, weed biomass ( Figure 3a) and cover (Figure 3b). However, there was considerable variation among the crops with the Brassicaceae oilseed radish, the forb buckwheat, the grass oats, and the legume field pea having the lowest weed biomass and cover (Figure 3a,b). Consequently, weed suppression (the inverse of weed biomass and cover) also varied significantly among plant functional groups and the identity of the crops (Table 2) and was greatest in the Brassicaceae and lowest in the legumes with the greatest suppression with oilseed radish, buckwheat, oats, and field pea (Figure 3c). Weed suppression was positively correlated with greater leaf area, LMA, maximum canopy height and the spring RGR of the crops (Figure 3d).  Table 2. Significant correlation (Spearman) between weed suppression and plant functional traits are shown in (d). Note weed suppression was not correlated with rooting density (r = 0.141, p = 0.192). Note, LMA = leaf mass per area and RGR = relative growth rate.

Erosion Control
The average crop cover throughout the growing season, winter cover, rooting density and the specific root length all varied significantly among crops and the four plant functional groups ( Table 2). The ground cover throughout the growing season was greatest in the oilseed radish, followed by buckwheat, oats, annual ryegrass and field peas (Figure 4a). Overall, the legumes provided the least growing season cover while the Brassicaceae and forbs provided some of the greatest (Figure 4a). However, winter cover was the greatest in the legumes and particularly in hairy vetch and red clover (Figure 4b). Annual ryegrass was also able to provide a high level of ground cover by post-harvest re-growth prior to winter (Figure 4b). Rooting density was greatest in the crops with large taproots. Brassicaceae had the highest rooting density due to the oilseed and tillage radishes, as well as sunflowers, which also had a high rooting density (Figure 4c). The scaled average of all these indicators representing the ecosystem service erosion control varied significantly among both crops and functional groups ( Table 2). Oilseed and tillage radishes as well as annual ryegrass scored the highest for potential erosion control. Erosion control was positively correlated with greater leaf area, LMA, maximum canopy height and the spring RGR of the crops (Figure 4e).  Table 2. Note, LMA = leaf mass per area and RGR = relative growth rate.

Soil Building
The soil particulate organic matter (POM) carbon and nitrogen, as well as the active carbon (POX-C), did not vary significantly among crops or their functional groupings (Figure 5a-c, Table 2). Aggregate stability did vary significantly among crops with red clover, oats and phacelia having the highest percentages of aggregate stability, and hemp and annual ryegrass having the lowest (Figure 5d). The average of these measures, indicating the soil building properties of the crops, also did not vary significantly among crops or their functional groups (Figure 5e, Table 2). However, the composite of these measures representing the agroecosystem service of soil building was significantly positively correlated with greater leaf area and maximum canopy height (Figure 5f).  Table 2. Significant correlation (Pearson's) between soil building and plant functional traits are shown in (f). Note soil building was not correlated with leaf mass per area, spring relative growth rate or rooting density (r = −0.095, p = 0.392; r = 0.090, p = 0.403; r = −0.006, and p = 0.958 respectively).

Soil Health
With our experimental design of four replicates spread over two fields, we were not able to detect statistically significant differences in fungal richness and evenness among crops or their functional grouping (Table 2). However, fungal richness and evenness were lowest in soils conditioned by the Brassicaceae family (Figure 6a,b). The greatest fungal richness occurred with oats, followed by alfalfa and field peas (Figure 6a). The abundance of AMF sequence reads in the soils did not vary significantly among crops or functional groups, but the highest levels were observed with the legume galega and the forbs sunflower and hemp var. ferimon (Figure 6h). The abundance of saprotrophic fungi and plant pathogenic fungal guilds varied significantly among plant functional groups ( Table 2) and were lowest in the Brassicaceae family (Figure 6i,j). The abundance of fungal saprotrophs also varied significantly among crops. Oilseed and tillage radishes resulted in a significantly lower abundance of saprotrophic fungi than most other crops (Figure 6d). Fungal richness was positively correlated with a greater spring RGR, and fungal saprotrophs were negatively associated with crops with a larger leaf area, LMA and taller canopy (Table 3).  Figure S1 for the variation among crops that are associated with the different microbial community characteristics.
Bacterial richness and evenness both varied significantly among the crop functional groups (Table 2). Bacterial richness and evenness were greatest in soils conditioned by grasses and lowest in soils conditioned by legumes (Figure 6f,g). However, bacterial evenness also varied significantly among crops (Table 2). Of the legumes, galega, white clover and faba beans had the greatest bacterial evenness, which was similar to that of the grasses (Figure 6g). The abundance of rhizobia detected in the bulk soil varied significantly among both plant functional groups and individual crops ( Table 2). Rhizobia were least abundant in the soils conditioned by the Brassica family and greatest in the forbs and legumes (Figure 6h). Specifically, field peas, crimson clover, red clover and alfalfa had the greatest abundance of Rhizobia, as well as the forb phacelia (Figure 6h). Nitrifying bacteria only varied significantly among plant functional groups, where grasses had the highest, and legumes the least, abundance of nitrifying bacteria (Table 2, Figure 6i). Cyanobacteria did not vary significantly among crops or their functional groups but were most abundant in soils conditioned by tillage radishes, oats, faba beans and hemp var. ferimon and least abundant in soils conditioned by buckwheat (Table 2, Figure 6e). Bacterial evenness was significantly correlated with greater LMA, while Rhizobia were marginally negatively associated with greater LMA in crops (Table 3). Cyanobacteria were significantly positively correlated with a greater rooting density ( Table 3). All fungal indices were significantly positively correlated with each other, with the exception of AMF and plant pathogens, and all bacterial indices were positively correlated with each other, with the exception of rhizobia and Cyanobacteria (see Table S1 and Figure S1a). The averaging of the standardized values of these soil community indices as a measure of soil health varied significantly among crop functional groups (Table 2), where the brassicas scored the lowest for promoting soil health and grasses the highest ( Figure S1b). Soil health was marginally associated with a greater rooting density (Table 3).

Correlations among Ecosystem Services
The agroecosystem services soil building, nutrient capture, weed suppression and erosion control were all positively correlated with biomass production (Figure 7). Intriguingly, soil health was negatively associated with all agroecosystem services, with the exception of crop nutrient capture (Figure 7 and see Figure S2). By assessing the independent correlations between agroecosystem services with specific soil microbial community characteristics that made up our index of soil health, we found that soil building was significantly negatively associated with greater abundance of sequence reads of AMF, saprotrophic fungi, rhizobia, nitrifying bacteria and Cyanobacteria as well as bacterial richness (Table 3). Fungal richness was significantly positively correlated with a greater nutrient capture (Table 3). Fungal evenness and AMF were significantly negatively correlated with greater crop productivity as well as greater weed suppression (Table 3). Weed suppression was also significantly negatively correlated with greater bacterial richness and nitrifying bacteria (Table 3).

Discussion
In this study, we assessed the effects of a wide variety of crops on six ecosystem services in order to (1) identify full season cover crops that maximize specific ecosystem services and (2) assess effects of crop functional traits on the soil microbial communities that are associated with enhancing these ecosystem services. As expected, we found that cover crops vary greatly in their traits' influence on ecosystem services. These results show that not all crops are equal in their ability to promote various ecosystem services, as some crops may promote particular agroecosystem services, but may have a potential trade-off with another agroecosystem service. Beyond the identity of the cover crops and their functional grouping, we found that the traits of crops can also be good predictors of their effect on ecosystem services. Specifically, leaf area and maximum canopy height were significantly positively correlated with all ecosystem services, with the exception of soil health, even if we were unable to detect statically significant effects of individual crops or their functional groupings on ecosystem services. For instance, even though measures of soil building did not show any statistically significant variation among crops or their functional grouping, we did find that soil building was significantly positively associated with crops with a greater leaf area and maximum canopy height. Further, plant functional traits that were linked with various ecosystem services were also linked to various soil microbial community characteristics that we used as indicators of soil health. These results show that the functional traits of crops may be considered as an informative tool as to the potential of a crop to promote a desired agroecosystem service that may be equally as informative as their identity or functional group assignment.

Productivity, Nutrient Capture and Weed Suppression
Crops that exhibited a more rapid growth rate 20 days after planting (spring RGR), achieved a greater maximum canopy height and produced a larger leaf area were all associated with crops that produced greater biomass, higher nutrient capture and were more effective in suppressing weeds. These results are expected as multiple studies have demonstrated the competitiveness of a plant for capturing resources over neighboring species is determined by growth rate, leaf area and canopy height [40,[89][90][91]. Weed suppression was greatest with oat paralleling a previous study where oat was also the most weed suppressive compared to wheat and barley [91]. Buckwheat was also particularly effective in suppressing weeds in our study, which was also among the crops with the greatest canopy height, leaf area and spring RGR. Buckwheat has often been shown to be a successful crop for suppressing weeds due to aggressive growth [92,93], as observed here, and potential allelopathic effects of buckwheat to inhibit root growth of competing weeds [94][95][96]. While the oilseed radish was not one of the most productive crops, it was in the top three for weed suppression in our study. Oilseed and forage radish have been known to suppress weed populations in the spring under low fertilizer inputs due to its aboveground competitiveness [97,98]. While in our study oilseed radish was not necessarily among the crops with the greatest spring RGR it was among the top three for greatest leaf area and the greatest % ground cover throughout the growing season demonstrating that the weed suppressive abilities of oilseed radish are due to aboveground light competition. Hemp has also been known to have weed suppressive properties due to its rapid and tall growth [99]. Here we also found hemp var. ferimon also ranked within the top five crops for weed suppression even when given minimal fertilizer input, but the variety anka ranked lower indicating that the weed suppressive abilities of a full season cover crop can be variety specific.
Legumes were particularly effective at capturing more nutrients from the environment due to their symbiotic reliance on rhizobia and arbuscular mycorrhizal fungi [100,101]. As a green manure, the high biomass producing legumes, such as field pea in our study, can be particularly effective for naturally building up soil fertility and organically bound nutrients [102,103]. However, we found that the residual NO 3 remaining in the soil at the end of the growing season was highest with legumes. This means that while legumes are highly efficient at capturing N from the environment, these crops may not be ideal as a high fertilizer N use efficiency crop for capturing excess mobile forms of N in the soil because legumes obtain most of their N through N 2 fixation [104,105].

Soil Building and Erosion Control
Soil particulate organic matter is an important component of soil organic carbon and nutrient cycling as it provides a pool of organically bound nutrients and is a primary resource for nutrient mineralization by soil microbes [106][107][108]. While there was little significant variation among crops and their functional grouping on soil building (POM-C and N) after a single growing season, we did observe that a greater leaf area and maximum canopy height, and thus greater crop biomass production, were significantly associated with greater soil POM. This is likely due to greater photosynthetic capacity of taller plants with larger leaves to accumulate greater soil POM-C and N that may have become deposited in the soil through root exudates and allocation of photosynthetically derived carbon to soil microbes [109][110][111]. However, building soil POM-C and N pools may take multiple years to accumulate [112].
Soil erosion control potential was greatest with the two tillage and oilseed radishes due to their large taproot that resulted in a high rooting density. Such a rooting structure can break up compacted soil providing greater vertical transmission of precipitation into the soil, thus reducing potential surface run off related erosion [48][49][50]. Annual ryegrass also scored high for potential erosion control as it provided a high level of vegetative ground cover throughout the growing season and regenerated vegetative growth postharvest prior to winter. Previous studies have also shown that the use of ryegrass as a cover crop can substantially reduce soil erosion and sediment loss over various slope and rain intensities [113]. While the radish crops and annual ryegrass were not the greatest biomass producing crops, we did find that crop biomass production was also positively related to erosion control and weed suppression, likely due to biomass production, weed suppression and erosion control being all associated with greater leaf area, maximum canopy height and a more rapid spring growth rate. While the same crop management was used for all crops in this study to avoid bias, it should be noted that the potential of soil erosion control for a given crop can be improved by optimizing seeding rates, date and field operations such as tillage and fertilization. Also, the erosion control potential was inferred based on measurement of indirect parameters (surface cover and root density) in this study. Thus, the erosion control abilities of the potential crops identified here needs to be further confirmed with measured water erosion in plot or field scale studies to establish ways to enhance erosion control with these crops.

Soil Health: Linking Microbial Communities, Agroecosystem Services and Crop Traits
We found the Brassicaceae significantly reduced the abundance of potential plant fungal pathogens. Brassicaceae, and more specifically mustards, have been deemed ideal for biologically suppressing plant pathogens such as plant parasitic nematodes and fungal pathogens due to their production of glucosinolates and isothiocyanates [114][115][116][117][118][119][120][121][122]. However, it has also been shown that isothiocyanates produced by Brassica kaber and Brassica nigra inhibit mycorrhizal fungi [123]. The invasive Alliaria petiolata (garlic mustard) has also been shown to use allelopathy to suppress AMF that native competitor species depend upon for enhanced growth [124]. While in our study the Brassicaceae did not significantly reduce AMF, they were also associated with a strong reduction in the abundance of saprotrophic fungi and plant symbiotic Rhizobia. These results show that while Brassicaceae can suppress potential plant fungal pathogens, they may also suppress other functionally important soil microbes that are needed for promoting soil carbon and nutrient cycling.
The grater bacterial richness and evenness in grasses that we observed is likely due to the greater fine rooting structure of grasses and high root turnover rates that provide resources and habitat for bacterial communities [125][126][127]. In line with this root trait driven microbial community composition, we also found that greater rooting density was positively associated with greater bacterial richness, evenness and abundance of nitrifying bacteria. While members of the Brassicaceae are known for their anti-fungal properties mentioned above, these crops here did not seem to suppress characteristics of the bacterial community, and in particular tillage radish. This may seem counter intuitive, but it has been observed in another short-term cover crop study that oilseed radish can promote microbial biomass [128], yet others have observed no effect of this crop species on fungal abundance [129]. Taken together along with our findings, it would seem that tillage and oilseed radish may have strong selective effects on promoting associated soil bacterial and or suppressing fungal community characteristics that seems not well understood. Arbuscular mycorrhizal fungi are well known to be important for soil resource acquisition and growth for the majority of vascular land plants and have been considered an important feature of agricultural soils [130][131][132]. While these are functionally important plant symbionts in terrestrial ecosystems these fungi do not distinguish between agricultural crops and weeds and these fungi can interconnect both intra-and inter-specific plant species to mediate the allocation of resources among competing plants [133][134][135][136]. In our study the abundance of AMF was negatively associated with the biomass of the sown crops and weed suppression. This would suggest that AMF abundance in the soil supported the productivity of weeds over the productivity of the sown crop in our system. While it has been proposed that AMF may be able to suppress agricultural weeds [137,138], it should be noted many crops, such as the Brassicaceae used in this study, are non-mycorrhizal plants [139] and many weed species are also mycorrhizal host plants from which they likely benefit from the association. Therefore, our results show that AMF may benefit the productivity and abundance of agricultural weeds over the promotion of sown crops. This does not mean that AMF have undesired impacts in agricultural systems, but rather that these agriculturally important fungi are indiscriminate between promoting desired crops and undesired weeds. Yet, whether and how AMF facilitate the coexistence or exclusion of agricultural weeds has not been well researched.
While we could not detect a statistically significant association between Rhizobia and AMF with nutrient capture across our wide variety of cover crops, we did observe that nutrient capture was significantly positively associated with greater fungal richness. Greater soil fungal diversity has been associated with reduced nutrient losses from the soil through leaching and N 2 O emissions and greater litter decomposition and plant nutrient uptake [80][81][82][83][140][141][142][143]. This infers that crops that promote greater fungal richness in soils may also enhance more efficient nutrient cycling processes that crops are able to capture. In our study these crops were oat and the legumes alfalfa and field pea. Legumes are well known to produce strigolactones that are plant hormones that stimulate fungal hyphal development, particularly in AM fungi, and can alter fungal community composition [144,145]. But such fungal promoting properties of legumes can also promote the development of plant fungal pathogens in cropping rotations [146,147]. While nutrient capture in our study was only significantly positively related to greater fungal richness in the bulk soil, a greater fungal richness was also strongly related to greater abundance of potential plant fungal pathogens further suggesting potential trade-offs between the benefits and undesirable consequences of legumes as cover crops.
Greater biomass production had a marginally non-significant negative association with Cyanobacteria, but unlike AMF the Cyanobacteria were unassociated with the abundance of weeds. This seems counter intuitive as Cyanobacteria are frequently considered as plant growth promoting bacteria [72][73][74]. In our study, the negative association between Cyanobacteria and crop biomass production can be attributed to buckwheat as it was one of the most productive crops having the lowest abundance of Cyanobacteria. Omitting buckwheat from the biomass-Cyanobacteria correlation could explain away the significant association (r = −0.17, p = 0.132). This negative effect of buckwheat on Cyanobacteria may be due to buckwheat's production of antimicrobial peptides that can inhibit bacterial and plant pathogenic fungal growth [148]. The lack of variation in the relative abundance of Cyanobacteria reiterates current knowledge that Cyanobacteria, while perhaps playing an important role in our soils, are generally ubiquitous in soils and there is not much known at to their functional role in agroecosystems [74,149].
It was unexpected that we found many desirable characteristics of the soil microbial community indicating soil health, such as microbial richness, saprotrophs, and AM fungi, to be negatively associated with soil building (see Table 3). However, this finding may also make sense in light that the POM-C and N in soils are the resources that many soil microbes utilize and deplete through mineralization processes [80,108,150]. Thus, in soils with a greater fungal evenness, AM fungi, bacterial richness, nitrifying bacteria, and Cyanobacteria our results inherently suggest that as there are more microbes present that are involved in the breakdown of organically bound resources that could have lead to the lower levels of POM-C and N by the end of the growing season. However, recently it was observed that greater soil C was related to greater fungal and bacterial biomass which in turn is negatively related to greater richness [151]. This indicates that soils with greater C are dominated by a few highly effective microbes in capturing resources resulting in lower richness, while lower levels of soil C drives greater facilitation and niche partitioning through specialization resulting in greater microbial richness. This may explain the negative POM-microbial diversity relationships in our study where more POM-C rich soils had reduced bacterial richness and fungal evenness.

Conclusions
Our results show that nearly all agroecosystem services measured here could be related to cover crops with a more rapid, taller growth and larger leaf area. This, in turn, produced more biomass, generating a greater standing N and P and weed suppression, which is in line with previous studies showing that agroecosystem services are often a function of cover crop biomass production [11][12][13][14]. However, we also show that not all crops are equal in their ability to promote all ecosystem services, and while some crops may be ideal for promoting a specific agroecosystem service, it could result in a trade-off with another. For instance, legumes may be ideal as a green manure due to their high standing N and P content but are relatively poor at capturing soil NO 3 and may promote potential plant fungal pathogens. Conversely, while members of the Brassicaceae family are particularly efficient at suppressing potential plant fungal pathogens, this may have undesirable effects on suppressing fungal saprophytic fungi and rhizobia for enhancing soil carbon and nutrient cycling. Thus, careful selection of the identity and functional traits of cover crops needs to be considered to target particular agroecosystem services. While different crops promote different ecosystem services, this may indicate that mixtures of species that differ in functional traits may complement each other to enhance particular agroecosystem services or enhance multiple ecosystem services simultaneously. However, whether mixtures of crops that are functionally different can enhance single or multiple ecosystem services through functional complementarity requires further detailed research, as competition among crops and lower densities of individual crops within mixtures may respectively interfere or dilute their abilities to enhance an ecosystem service.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agriculture11090830/s1: Table S1. Pearson's correlations among microbial community characteristics; Figure S1. Principal coordinate ordination of the ten soil microbial community characteristics and their average scores indicating the ecosystem service soil health; Figure S2. Principal coordinate ordination of the six ecosystem services.