High Seeding Rates and Low Soil Nitrogen Environments Optimize Weed Suppression and Profitability in Organic No-Till Planted Soybean

No-till planting crops into rolled-crimped cover crops can improve soil health while reducing labor and fuel requirements compared with traditional tillage-based production. However, little information is available to help farmers optimize the management of organic no-till planted crops. Weed suppression, crop yield, and profitability were assessed across soybean [Glycine max (L.) Merr.] seeding rates and soil nitrogen environments in an experiment conducted at two sites in central New York. Soybeans were no-till planted into rolled-crimped cereal rye (Secale cereale L.) at 0, 185,000, 371,000, 556,000, and 741,000 seeds ha−1. Three rates (0, 63, or 125 kg ha−1) of sodium nitrate (15-0-2) were applied across seeding rates to create different soil nitrogen environments. When pooled over sites, the lowest weed biomass occurred at the highest soybean density in the lowest soil nitrogen environment. An interaction was observed between soybean seeding rate and nitrogen treatments on weed communities. Soybean yield increased asymptotically with crop density and was not affected by nitrogen or site treatments. When pooled over nitrogen treatments and sites, partial returns to the soybean seeding rates were maximized at $2,238 ha−1 with 527,800 seeds ha−1. Results suggest that crop density is an important lever for optimizing weed suppression and crop yield in organic no-till soybean, and that managing for low soil nitrogen conditions may further enhance weed suppression while maintaining high yields.


INTRODUCTION
Constraints to organic soybean production include fuel use, repair, and labor costs from soil tillage and cultivation for weed management (McBride and Greene, 2015). Organic rotational no-till production is an alternative management system for soybeans that typically has lower labor and fuel inputs (Ryan, 2010;Mirsky et al., 2012) and enhanced soil health benefits (Crowley et al., 2018) compared to traditional tillage-based organic management. No-till production of organic crops can be accomplished by planting into cover crop mulches (Vincent-Caboud et al., 2019).
To supplement the weed suppressive effects of cover crop mulches, various physical and cultural weed management tactics can be used to further reduce weed abundance (Mirsky et al., 2012). High-residue inter-row cultivation is an option when soybean are planted in wide rows (Zinati et al., 2017); However, efficacy of high-residue inter-row cultivation can vary with equipment, soil conditions, and mulch levels, and might serve best as a rescue treatment rather than a planned operation (Keene and Curran, 2016). Conversely, increasing the seeding rate of no-till planted soybean is a proactive weed management strategy that suppresses weeds through faster canopy closure (Weiner, 1990;Schwinning and Weiner, 1998), increased light interception (Steckel and Sprague, 2004), and greater crop-weed size asymmetries (Bastiaans et al., 2008). Recommended soybean seeding rates in tilled organic production are 370,500 seeds ha −1 (Cox et al., 2019) but higher seeding densities are desirable for organic no-till planted soybean production. In organic notill planted soybean, weed suppression increased with soybean seeding rates between 195,000 and 914,000 seeds ha −1 (Liebert and Ryan, 2017). In that research, profitability was maximized at 646,000 seeds ha −1 at one site and 728,000 seeds ha −1 ac at another site, suggesting that high soybean seeding rates are a viable weed management strategy in an organic no-till system.
Soil nitrogen availability is known to affect several ecological processes that are relevant to organic no-till soybean production systems including weed-crop competition and biological nitrogen fixation. In conventional conservation tillage systems, starter nitrogen fertilization can provide a yield advantage to soybean (Touchton and Rickerl, 1986;Osborne and Riedell, 2006) but results are inconsistent (Mendes et al., 2003;Hungria et al., 2006). However, in organic cropping systems, where herbicides are prohibited, increased nitrogen availability can encourage weed growth and weed-crop competition. As an agroecological weed management strategy, Ryan and Peigné (2017) hypothesized that reduced plant-available nitrogen can increase the competitive ability of legume crops against weeds and promote biological nitrogen fixation, enhancing legume crop yields. In a theoretical context, low nitrogen environments likely promote legume dominance because of their reduced dependence on soil nitrogen relative to many non-legume weeds (Tilman, 1982). Supporting research in soybean production systems has found greater biological nitrogen fixation in low compared to high soil nitrogen environments (Schipanski et al., 2010). Mulching with cereal rye reduces soil nitrogen availability because the high carbon to nitrogen ratio of cereal rye residue encourages nitrogen immobilization as it decomposes (Wells et al., 2013;Ketterings et al., 2015;Williams et al., 2018). Cover crop mulches also lower soil temperature (Teasdale and Mohler, 1993a), which can decrease nitrogen mineralization (Miller and Geisseler, 2018). Despite the potential for a low nitrogen environment to reduce weed competitiveness, nitrogen fertilization has been proposed in organic no-till systems with low soil nitrogen to encourage cereal rye growth (Mirsky et al., 2013). However, in previous research, Ryan et al. (2011a) demonstrated that nitrogen fertilization increased cereal rye biomass, but the resulting increase in rolled cereal rye mulch did not improve weed suppression. To understand the effects of reduced soil nitrogen availability from cereal rye mulch on weed management and crop yields, it is important to test the effect of different soil nitrogen environments under a uniform cereal rye mulch.
An experiment was conducted to test the effects of crop density and soil nitrogen environment on weed community structure and crop yields. Partial returns to the soybean seeding rates were also assessed to determine the economic optimum no-till planted soybean seeding rate in the three nitrogen environments. We hypothesized that adding nitrogen would increase weed biomass, which would reduce yields and increase the economic optimum seeding rate.
Seeding rate treatments were applied in 36.5 by 6 m strips at both sites. Each seeding rate main plot was split into 12.2 by 6 m thirds, and 0, 63, or 125 kg of nitrogen ha −1 was applied within 24 h of soybean planting by broadcasting sodium nitrate (15-0-2, North Country Organics, Bradford, VT, USA). The added nitrogen treatments were used to create different soil nitrogen environments; however, soil nitrogen was not measured directly. This design was replicated in four blocks at each site. The soybeans did not receive additional fertilization or weeding throughout the season. At the Geneva site, Pyganic R (5% pyrethrins, McLaughlin Gormley King Company; Minneapolis, MN, USA) was applied at 4.2 L ha −1 on 19 July 2019 to mitigate a Japanese beetle (Popillia japonica, Newm.) infestation. This OMRI approved insecticide was applied from the field margin with a Cannon airblast sprayer (Jacto; Tualatin, OR, USA) to avoid damaging the soybean canopy.
To confirm the establishment of our treatments in the field, we assessed mulch biomass and soybean populations. Immediately after soybean planting, rye biomass was sampled with a 0.25 m 2 quadrat in each seeding rate main plot (n = 20 site −1 ). Soybean emergence was assessed at growth stages V1 (24 June 2019 in Aurora; 25 June 2019 in Geneva) and R8 (8 to 9 October 2019 in Aurora; 10 to 11 October 2019 in Geneva) by counting soybean plants within a one-meter row length in each sub-plot (n = 60 site −1 ). A comparison of soybean populations in growth stage V1 and R8 allowed for quantification of soybean self-thinning.

Weed and Soybean Sampling
Weed biomass was sampled when soybean reached the R6 growth stage (15 to 17 September 2019 in Aurora; 19 to 20 September 2019 in Geneva). One 0.5 m 2 quadrat was sampled in each subplot; all weeds taller or wider than 5 cm were sampled at the soil surface and identified at the species level. Sampled biomass was dried at 55 • C for 2 weeks before weighing. A second sampling to assess soybean yield occurred when soybeans were at growth stage R8 (8 to 9 October 2019 in Aurora; 10 to 11 October 2019 in Geneva). During this harvest, soybean was sampled along a 1 m row length and seed pods were removed from the soybean by hand. Seed pods were dried at 55 • C for two weeks and beans were weighed after shelling. Soybean yields were adjusted to 13% moisture for all analyses.

Statistical Analysis
All statistical analyses were done in R version 3.6.1 (R Core Team, 2019). Linear mixed-effects models were fit with the "lmerTest" package (Kuznetsova et al., 2017) and analyzed with least-squares means comparisons using the "emmeans" package (Lenth et al., 2019). Non-linear models were fit with "nlme" and compared with F-tests. Multivariate weed community analyses were done with PERMANOVA tests using the vegan package (Oksanen et al., 2020).
Assumptions of linearity were assessed by plotting fitted against residual values. Non-linear functions were confirmed by first fitting local polynomial regressions and verifying that their shape was similar to the non-linear equation to be fit. After fitting the non-linear model, the accuracy of the calculated parameters was confirmed by plotting the function through the data. For PERMANOVA tests of the weed community, the homogeneity of multivariate groups was confirmed by assessing their beta dispersion (Anderson, 2006).

Mulch Establishment, Soybean Density, and Self-Thinning
To quantify mulch biomass, cover crop samples were collected from 0.5 m 2 quadrats at the block level (n = 4 site −1 ) just prior to rolling. Mulch biomass was compared between the two sites with t-tests. Variation of mulch biomass within each site was described by fitting a linear model to mulch biomass as a function of experimental block.
Soybean plant density was assessed through linear mixedeffects models that described soybean density as a function of site and the interaction of seeding rate, nitrogen treatment, and soybean growth stage. Experimental block and seeding rate were nested random effects to account for variability within each field and the split-plot design of the experiment. In these models, data from the 0 seeds ha −1 treatment were removed because there was no variation in soybean density at this seeding rate. To address violations of the assumption of normality, soybean density data were log-transformed.

Weed Biomass
The effect of soybean density on weed suppression was modeled with an inverse hyperbolic function where weed biomass (W b ; kg ha −1 ) was driven by soybean density (S d ; plants ha −1 ) (Liebert and Ryan, 2017). (1) The W 0 parameter represents the amount of weed biomass (kg ha −1 ) in the absence of soybean plants and i (ha plant −1 ) is the hyperbolic shape term, describing the reduction in weed biomass across soybean density where the greater the i parameter, the steeper the initial slope of the hyperbola. The effect of site and nitrogen treatment was elucidated by comparing a reduced model, which pooled data across all sites and nitrogen treatments to semi-reduced and full models. The semi-reduced model allowed the initial amount of weed biomass (W 0 ; kg ha −1 ) to vary as a function of site or nitrogen treatment level. In the full model, both the initial weed biomass (W 0 ; kg ha −1 ) and shape (i; ha plant −1 ) parameters differed as a function of site or nitrogen treatment levels. Following the method described by Weisberg (1985) and used by Ryan et al. (2011b), Teasdale et al. (2006), and Seefeldt et al. (1995), the reduced, semi-reduced, and full models were compared with F-tests. Model selection was based on F-test results, which allowed us to choose the simplest model that best described weed biomass. All weed biomass models considered field block as a random effect that was nested within each experimental site. To get model convergence, two outliers had to be removed from the data set.

Weed Communities
Weed communities were assessed with a PERMANOVA on the Bray-Curtis distance matrix of weed species biomass.
PERMANOVAs were sequential and modeled weed communities as a function of the interaction between seeding rate and nitrogen treatment. The weed community data were modeled separately within each site to maintain homogenous group variance across nitrogen and seeding rate treatments. Two missing experimental units at the Geneva site were imputed as the mean biomass of each weed species in the treatment replicates at the site to balance the data. To account for the field blocks and splitplot design of the experiment, permutations were constrained to field blocks and the nested effect of blocks and soybean seeding rate treatments.
The 0 soybean seeds ha −1 treatment was not included in the multivariate weed community analysis because weed community data from the 0 seeds ha −1 treatment had differing variation compared to the other seeding rates, violating the PERMANOVA assumption of equal multivariate group variance. Furthermore, omitting the 0 seeds ha −1 treatment allowed the multivariate analysis to report the effect of changing soybean seeding rates rather than the differences in the weed community arising from the presence and absence of soybean.
The multivariate weed community analysis was supplemented with linear mixed-effects models that described the amount of total annual and perennial weed biomass as a function of seeding rate and nitrogen treatment in the two sites. In these mixed effects models, block and soybean seeding rate were nested random effects. Annual and perennial biomass data for the two sites were square root transformed to conform to linear model assumptions. Data from the 0 seeds ha −1 treatment, was not included in these models to align with multivariate analysis and emphasize the effect of changing soybean seeding rate on annual and perennial weed biomass. For all linear mixed-effects models, ANOVA and least-square means tests were used to assess main and treatmentlevel effects, respectively. Dominant weeds were described by calculating the average cumulative biomass of each weed species in the three nitrogen treatments at the two sites. Weed species that compromised up to 90% of the cumulative weed biomass were reported. A 90% threshold was chosen because species often compromised <5% of the cumulative weed community biomass at higher thresholds.

Soybean Yield
Soybean yield (Y c ; kg ha −1 ) was modeled as a function of soybean density (S d ; plants ha −1 ) with an asymptotic function (Liebert and Ryan, 2017). In the asymptotic function, a is the horizontal asymptote and represents a yield plateau, and b is the natural logarithm of the rate constant.
To elucidate the effect of site and nitrogen treatment on crop yield (Y c ; kg ha −1 ), a reduced model, where yield was only a function of soybean density (S d ; plants ha −1 ), was compared to a full model. In the full model, the maximum yield parameter (a) varied as a function of site or nitrogen treatment level. Comparisons between the reduced and full models were accomplished through F-tests (Weisberg, 1985;Seefeldt et al., 1995;Teasdale et al., 2006;Ryan et al., 2011b) using the ANOVA function in R (R Core Team, 2019). Field blocks were accounted for as a random effect that was nested within each experimental site in all yield models.

Partial Returns to Soybean Seeding Rates
Partial returns were calculated to find optimal seeding rates for the organic no-till planted soybeans. For each experimental unit, the partial return was calculated as where P r denotes partial return (dollars ha −1 ); Y c is soybean yield (kg ha −1 ); M p is the average market price of organic feed-grade soybean ($0.69 kg −1 ) and was calculated from biweekly farmgate market prices across the United States from 2018 to 2019 (USDA-AMS, 2020); S c the price of seed used in the experiment ($0.000321429 seed −1 ); S r is soybean seeding rate (seeds ha −1 ). Partial returns across soybean seeding rates were described by a quadratic function (Liebert and Ryan, 2017;Purucker and Steinke, 2020): where P r denotes partial return (dollars ha −1 ) of the organic notill system, S r is the soybean seeding rate (seeds ha −1 ), a and b are constants of the quadratic equation. To assess whether the experimental site or nitrogen environment affected partial returns to the seeding rate gradient, F-tests were used to compare a reduced model to a full model, which allowed all formula parameters to vary as a function of site or nitrogen treatment level. In all partial return models, field blocks were accounted for as a random effect that was nested within each site.

Field Conditions and Soybean Density Gradient
During the growing season, the average monthly temperature from May through October in Aurora (16.7 • C) and Geneva (16.5 • C) was slightly lower than the 15-year average, except in July which was warmer at both sites ( Table 1). Cumulative monthly precipitation from May through October was slightly higher in Geneva (66.2 cm) than in Aurora (63 cm). At the two sites, May and June were wetter and September was drier than the 15-year average. Cereal rye biomass was greater at Geneva (6,560 kg ha −1 ) than Aurora (5,310 kg ha −1 ) (P < 0.01). Within each site, cereal rye mulch biomass did not differ among blocks (P > 0.05), indicating that cereal rye growth and consequently mulch was uniform across sites. Although there was reduced soybean establishment at high crop densities, the seeding rate treatments resulted in a crop density gradient ranging from 0 to 701,801 plants ha −1 (Figure 1; P < 0.001). Site and nitrogen treatments did not affect soybean density (P > 0.05) and we did not observe differences in crop lodging across treatments. Likewise, during the growing season, there was no self-thinning because soybean density was consistent at growth stages V1 and R8. FIGURE 1 | Linear model of the observed soybean density gradient during the beginning (soybean growth stage V1; R 2 = 0.97; P < 0.001) and end of the season (soybean growth stage R8; R 2 = 0.94; P < 0.001). The dashed line represents expected soybean density at 90% germination, the germination rate of the soybean (cv. "Viking 0.1518N") used in this experiment according to the information on the seed tag.

Weed Suppression
Comparisons of weed biomass and soybean density in Aurora and Geneva show that weed suppression was not affected by site (P > 0.05). Consequently, weed biomass data from the two sites were combined when testing for the effect of nitrogen treatment on weed biomass (Figure 2). When pooled across sites, weed biomass ranged from 0.02 to 332.5 kg ha −1 (Figure 2). Weed suppression differed as a function of nitrogen treatment (P < 0.01). However, the similarity between the semi-reduced and full weed suppression models (P = 0.56) indicates that the effect of soil nitrogen environment on weed suppression was due to increasing weed biomass in the absence of soybeans, not an interaction between nitrogen treatment and soybean density. In the absence of soybeans, plots that were fertilized with 63 and 125 kg N ha −1 had 70% and 107% more weed biomass relative to plots that did not receive nitrogen, respectively (Figure 2). FIGURE 2 | Weed suppression across soybean density at growth stage V1 and 0, 63, and 125 kg N ha −1 . The inverse hyperbolas depict the semi-reduced model of weed suppression (Equation 1), where initial weed biomass (W 0 ) was allowed to vary across nitrogen treatments, but the shape parameter of the curves (i) was the same in all nitrogen treatments. All model parameters are described in Supplementary Table 1.
The weed suppressive effect of increasing soybean density diminished at the upper range of density (Figure 2). Previous research has also reported that increased soybean density reduces weed biomass following an inverse hyperbolic trend (Ryan et al., 2011b;Liebert and Ryan, 2017). Diminishing rates of weed suppression at high soybean densities align with findings that crop density promotes weed suppression through light interception (Steckel and Sprague, 2004). The plateau in weed suppression reflects maximal canopy closure, where increasing soybean density has an increasingly smaller effect on light transmittance.
The increase in weed biomass in heightened nitrogen environments indicates that altered nutrient availability is a critical component of weed suppression from rolled cereal rye mulches. Weeds often incur similar or greater benefits from nitrogen fertilization relative to crops (Liebman and Mohler, 2001;Liebman and Gallandt, 2002) because weeds tend to uptake nutrients more efficiently than crops (Di Tomaso, 1995;Blackshaw et al., 2003).

Weed Communities
Soybean seeding rate and nitrogen treatments had an interactive effect on weed community composition at the two sites ( Table 2; P < 0.05). The interaction between treatments suggests that changes in soybean density have differing effects on weed communities depending on soil nitrogen environment. Changes in soil nitrogen had larger effects on annual than perennial weed biomass (Figure 3). Past experiments also reported that The analyses were based off a Bray-Curtis distance matrix of weed species abundance.
To emphasize the effect of changing soybean density on weed community composition rather than the presence or absence of soybeans, data from the 0 seeds ha −1 was omitted. Terms were added sequentially to determine significance and reflect the main plot (seeding rate) and split-plot (nitrogen) treatment structure.
fertilization had a filtering effect on weed communities Storkey et al., 2010). Common ragweed (Ambrosia artemisiifolia L.) was the most or second most dominant weed in all soil nitrogen environments at the two sites, compromising 65-73% and 14-17% of the average weed biomass in Aurora and Geneva, respectively ( Table 3). Common ragweed was reported as the dominant weed in other research on organic no-till planted soybean (Mirsky et al., 2011), primarily because it emerges before roller-crimping cereal rye and survives the process (Mirsky et al., 2012). At both sites, nitrogen additions also stimulated the growth of some relatively nitrophilic weeds, including Panicum capillare L. in the 63 and 125 kg N ha −1 treatments at Aurora and Echinochloa crus-galli (L.) P. Beauv. in the 125 kg N ha −1 treatment at Geneva ( Table 3).

Crop Yield and Partial Returns to Soybean Seeding Rates
Soybean yield increased asymptotically with soybean density and did not differ between sites (P > 0.05; Figure 4). Furthermore, when the yield data were pooled across sites, soil nitrogen environment did not affect soybean yield (P > 0.05). Ryan and Peigné (2017) hypothesized that reduced plant-available nitrogen would increase the relative competitiveness of soybean against weeds, leading to greater yields. However, the neutral effect of soil nitrogen environment on soybean yields did not support expectations of increased yield in low soil nitrogen environments. Previous findings have also shown that spring nitrogen fertilization can have neutral effects on soybean yield (Mendes et al., 2003;Hungria et al., 2006). In our study, increased weed biomass (P < 0.01; Figure 2) may have offset a soybean yield response to added nitrogen.
Partial returns followed a quadratic trend across soybean seeding rates because of diminishing yield improvement at high soybean densities (Figure 5). Similar to soybean yield, partial returns to the soybean seeding rates did not differ between sites or FIGURE 3 | Annual and perennial weed biomass in Aurora and Geneva as a function of nitrogen treatment. Different letters denote significant pairwise comparisons (P < 0.05) within each site and were calculated from linear mixed effects models that assessed annual or perennial weed biomass as a function of seeding and nitrogen treatment with field block and soybean seeding rate as nested random effects. Sites were modeled separately and data from the 0 seeds ha −1 was omitted.  (Ellenberg, 1991) and are from Julve (1998). nitrogen treatments (P > 0.05) and profit optimization occurred at $2,238 ha −1 with 527,800 seeds ha −1 . Soybean seeding rates at the economic optimum are higher than recommended seeding rates in tilled organic (370,500 seeds ha −1 ; Cox et al., 2019) and conventional (321,000 seeds ha −1 ; Orlowski et al., 2012) systems. Higher recommended seeding rates compared to the tilled systems may be due to a greater reliance on cropdensity based weed management in the absence of tillage. The profit-optimizing seeding rate for our feed-grade soybeans is lower than recommended for food-grade no-till planted organic soybeans (646,000 and 728,000 seeds ha −1 ; Liebert and Ryan, 2017). The difference between the recommended no-till seeding  2). Site or nitrogen treatment levels are not shown because they did not affect yield (P = 0.21 and 0.80 the reduced compared to full model of site and nitrogen treatments, respectively). All model parameters are described in Supplementary Table 2. FIGURE 5 | Partial returns across soybean seeding rates. The quadratic curve depicts the reduced model of partial returns (Equation 4) where the shape terms of the quadratic model (a and b) did not vary across sites or nitrogen treatments. The arrow in the figure shows maximum yield. All model parameters are described in Supplementary Table 3. rates for organic feed-and food-grade soybean stems from the higher market price of food-grade grain, which justifies greater seeding rates.

CONCLUSIONS
High seeding rates increased weed suppression, yields, and profits in organic no-till planted soybeans. Partial returns to the soybean seeding rates were maximized at $2,238 ha −1 with 527,800 seeds ha −1 . The profit-maximizing seeding rates in this study were greater than soybean seeding rate recommendations for tilled organic (370,500 seeds ha −1 ; Cox et al., 2019) and conventional (321,000 seeds ha −1 ; Orlowski et al., 2012) soybean systems, respectively. Despite greater weed biomass, especially of annual species, in plots with added nitrogen, soybean yield was not affected by nitrogen treatment. Although a reduction in soybean yield from increased weed biomass was expected, the experiment was not designed to compare weed-crop competition across nitrogen treatments because weed abundance was confounded with soybean density, and yield was not assessed under weedfree conditions. For further development of organic no-till systems, future research should assess fertilization method and seeding rate effects on weed abundance and crop productivity. Fertilizer effects on weed species vary by product (Gaskell and Smith, 2007;Tang et al., 2014;Little et al., 2015Little et al., , 2021 and could influence weed community assembly, weed-crop competition, and crop yield. Common fertilizer application methods, like banding fertilizer along a crop row or injecting it into the soil near crops, can reduce weed biomass relative to broadcasting, which was used in this experiment (Di Tomaso, 1995;Blackshaw et al., 2004). Furthermore, while high soybean density suppressed weeds and enhanced yields, there is a greater risk of disease (Hwang et al., 2006;Chang et al., 2011;Swoboda et al., 2011) and crop lodging (Costa et al., 1980;Neugschwandtner et al., 2019) in dense crop stands. Understanding how different management decisions affect weed communities and competition, as well as crop productivity and crop susceptibility to other pests, will enable the development of robust weed management strategies for organic no-till systems.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
UM designed and conducted the experiment, did the statistical analysis, and wrote the manuscript. SP designed the experiment, helped coordinate field work, and helped with the analysis and manuscript. CP conducted the experiment and helped with the analysis and manuscript. RS and AD helped with the analysis and manuscript. MR designed and conducted the experiment and helped with the analysis and manuscript. All authors contributed to the publication and approved the submitted version.

FUNDING
This work was funded through the USDA NIFA Organic Transitions Program grant# 2018-51106-28775.