Interplay between nitrogen fertilizer and biological nitrogen fixation in soybean: implications on seed yield and biomass allocation

Legumes rely on soil mineral nitrogen (N) and biological N fixation (BNF). The interplay between these two sources is biologically interesting and agronomically relevant as the crop can accommodate the cost of BNF by five non-mutually exclusive mechanisms, whereby BNF: reduces shoot growth and seed yield, or maintains shoot growth and seed yield by enhanced photosynthesis, or reduced root:shoot ratio, or maintains shoot growth but reduces seed yield by reducing the fraction of shoot biomass allocated to seed (harvest index), or reducing concentration of oil and protein in seed. We explore the impact of N application on the seasonal dynamics of BNF, and its consequences for seed yield with emphasis on growth and shoot allocation mechanisms. Trials were established in 23 locations across the US Midwest under four N conditions. Fertilizer reduced the peak of BNF up to 16% in applications at the full flowering stage. Seed yield declined 13 kg ha−1 per % increase in RAUR6. Harvest index accounted for the decline in seed yield with increasing BNF. This indicates the cost of BNF was met by a relative change in dry matter allocation against the energetically rich seed, and in favor of energetically cheaper vegetative tissue.

effect on soybean seed yield. This conclusion is, however, largely based on generic trials where coarse fertilization regimes were established to shift the contribution of mineral N and BNF. In contrast, a full-N treatment devised with a careful experimental protocol to ensure an ample N supply during the entire crop season increased soybean seed yield by 11% in relation to unfertilized controls, with a range from no effect for stressful environments (ca. 2500 kg ha −1 ) but increases of 900 kg ha −1 in high potential environments (ca. 6000 kg ha −1 ) 16 . The goal of this study was to investigate the effect of fertilizer N application on BNF and its implications for soybean seed yield and seed protein concentration. We tested the hypothesis that the cost of N fixation is mediated by reduced biomass, reduced allocation to seed captured in the HI 21 , and reduced concentration of protein and oil in seed. Quantification of these effects will provide insights of BNF impact on crop C and N economy, and will contribute to explain the apparent inconsistency in soybean seed yield responses to N fertilization.

Results
Effect of N fertilizer on N fixation. Data for this study were collected from 23 different locations across the US Midwest during the 2016 growing season (Table 1; Fig. 1). Table 2 shows the relative abundance of ureides in R6 (full seed stage; RAU R6 ) for each location and treatment, Fig. 2 illustrates the seasonal dynamics of the relative abundance of ureides (RAU) for crops grouped in high, medium and low BNF, and Table 3 summarizes the parameters of the curves. The fitted model (equation 2) returned R 2 between 0.62 and 0.87, with P < 0.001 in all cases.
The RAU R6 ranged between 42 to 93% ( Table 2) and responded to all three sources of variation: treatment, location and their interaction (P < 0.001). Fertilizer reduced RAU max in all BNF groups, but not where BNF was already low in the control treatment (Table 3). In the high BNF group, RAU max dropped from 90% in controls to 75% in their fertilized counterparts (averaging all N treatments; Table 3); the reduction in RAU was larger when the application of N was delayed from vegetative to reproductive stages. In the medium group, RAU max declined from 84% in controls to 68% in the V4 (fourth-leaf) application and 74% in both sowing and R2 (full flowering) applications, with no clear effect of N application timing.
The reduction in RAU max can be a consequence of a shorter time to peak RAU, a reduced rate or a combination of both effects. Combination of both traits contributed to reduced N fixation in the medium BNF group, as peak RAU and maximum rate were attained earlier in N-fertilized crops. In the low BNF group, N fertilizer reduced RAU rate but time of peak RAU was not affected. Reduction of the area under the curve (AUC) and the time when RAU reached 50% was observed from high to low BNF groups; however, treatments did not affect AUC (Table 3). Fertilization treatments affected the dynamics of RAU rate. Both the timing of peak rate (t m ), and the timing when rate became negative (t max ) were delayed from low to high BNF groups ( Fig. 2d-f). For the high and medium groups, reproductive N treatment was the most effective reducing both t m and t max hence contributing to an overall reduction of RAU max .
Effect of N fixation on seed yield, biomass, harvest index, protein and oil concentration. Seed yield ranged from 3151 to 7175 kg ha −1 , seed protein concentration from 31.9 to 41.8 g 100 g -1 , and seed oil concentration from 16.7 to 23.9 g 100 g −1 ( Table 2). For these traits, ANOVA showed location was a significant source of variation (P < 0.001), with no effect of treatment and its interaction with location (P > 0.05). Total biomass ranged from 6093 to 11376 kg ha −1 showing differences only among locations (P < 0.001). Harvest index ranged from 0.37 to 0.56, and was affected by both treatment (P < 0.05) and location (P < 0.001).
In this experiment, MGs were allocated to locations for agronomic relevance (Table 1; Fig. 1). The dominant effect of location on crop traits is therefore confounded with crop phenology. For instance, the range of thermal time to R6 was 909 to 1733 °Cd. We thus fitted bilinear models to account for the effect of phenology on crop traits (Fig. 3a,c,d), and regressed the residuals against RAU R6 (Fig. 3b,d,f; see section on materials and method). Analysis of residuals showed seed yield declined at 13 kg ha −1 per % of RAU R6 for the whole data set and 10 kg ha −1 per % of RAU R6 (p < 0.05, Fig. 3b) for top yielding crops (0.99 quartile). Under stressful conditions leading to low HI, the rate of decline in HI with RAU R6 computed as the slope of the 0.01 quartile, was 2.5 times larger than the rate for the pooled data (Fig. 3d). Oil seed concentration was negatively associated with RAU R6 (p < 0.05; Fig. 3f) and crops with higher oil concentration (0.99 quartile) were more responsive to RAU R6 . After removing the effects of phenology, biomass and seed protein concentration (not shown) did not relate to RAU R6 (both p > 0.05).
The association between BNF traits, soil attributes, seed yield, biomass, HI, seed protein and oil concentrations, adjusted by the effects of phenology, were explored using principal component analysis (PCA). Results from the ordination analysis were presented in a bi-plot showing the two first principal components ( Fig. 4) where angles between variable vectors denote the level of association among them (i.e., acute angles denote positive associations and obtuse angles negative associations between variables). Pearson correlation analysis complements the associations depicted on the PCA analysis (Supplementary Table S1). Data points from the same MG grouping together showed similarities for the variables of the PCA, reinforcing the dominant influence of phenology previously observed (Fig. 3a,c,e). The MGs, from the shortest to the longest, grouped in the biplot along the first principal component. For instance, data from groups 4.4 and 3.2 were associated with higher AUC when compared to 0.3 MG.
Soil attributes were positively correlated in the first principal component, discriminating the locations with shorter maturity groups and higher values in soil attributes. Interestingly, AUC was negatively correlated with OM, pH, clay and sand percentage, and the RAU R6 and rate were positively correlated with the OM and sand (Supplementary Table S1). Seed yield correlated positively with both HI and biomass and negatively with RAU max , RAU R6 , and protein concentration. The AUC correlated positively with RAU max and t max and negatively with the maximum RAU rate.

Discussion
Biological N fixation in soybean has been quantified at different scales, from field to country [22][23][24] . The range of RAU in unfertilized controls in our study ranged from 48 to 93%. This compares with an average of 60% of N derived from BNF for the US 23 . In Argentina, BNF in 86 location-years averaged 60% and ranged from 12 to 90% 24 . In Brazil, measurements in 6 environments returned an average of 81% and a range from 69 to 94% 22 . All these studies reflected a similar BNF ceiling around 90%, comparable to the maximum recently reported by Ciampitti and Salvagiotti 25 . Nitrogen fertilization reduced BNF and increased seed yield by enhancing C allocation to seed. Phenology was the main source of variation in seed yield and traits associated with BNF clustered with maturity group (Figs 2 and 4). Soybean maturity group influences not only phenology but also growth, and allocation of biomass and nitrogen 26 . Variation in BNF with MG relates to both the duration of the reproductive period when BNF and biomass growth rate peak, and the delay in the exponential phase of BNF [27][28][29] . In our study, application of N reduced peak and altered dynamics of RAU during the season (Fig. 2), but phenology masked the association between N fixation and seed yield. After removing the dominant effect of phenology, crop yield declined with increasing BNF (Figs 3 and 4). This effect of maturity group has not been considered in previous comprehensive studies 20,30 . Reduced HI was the primary driver of the reduction in seed yield with increasing RAU R6 , with an additional weak but significant reduction in oil seed concentration. Low yielding environments showed a steeper decline of HI with increasing RAU R6 (Fig. 3d), highlighting the interaction with overall environmental conditions affecting dry matter allocation; stress during reproduction often reduces HI 31,32 .
Changes in plant C allocation in association with BNF have been reported at different levels of organization and time scales. Reduction in BNF associated with phosphorus deficiency altered short-term allocation of C in lupin (Lupinus luteus), reducing photosynthesis:respiration ratio, and increasing the ratio between growth respiration and maintenance respiration 33 . Likewise, low magnesium supply altered carbohydrate allocation in soybean, increasing sucrose and starch allocation to leaves that later limited nodule growth 34 . Decreases in biomass allocation in seeds for chickpea (Cicer arietinum L.) and pea (Pisum sativum L.) were reported with increasing BNF 7,35 . Re-analyzing the data of Sadras et al. 7 , where 20 chickpea varieties were grown in 8 environments, HI declined witn BNF at a rate of 0.0022 units per %BNF (Supplementary Fig. S1), in comparison with 0.0011 units per %RAU R6 for soybeans in our study; similar to our trial, the decline in HI with BNF was larger for stressed chickpea crops.
Generically, plants require 1 g of glucose to produce either 0.33 g of lipid, 0.40 g of protein, or 0.83 g of carbohydrates 36 . Reducing oil concentration in seed is therefore an energetically effective way to meet the cost of BNF, as found in this study (Fig. 3e). This is in contrast to previous studies where seed protein concentration was reduced and oil concentration did not change in response to BNF 37,38 .
Our findings are in contrast to other studies where N fertilizer reduced BNF but did not increase soybean seed yield 39,40 . For example, Santachiara et al. 39 found no seed yield response in heavily fertilized crops (600 kg N ha −1 spread over the season) that reduced BNF to 16% in comparison to 69% in unfertilised controls. However,  Santachiara et al. 39 neither report seed yield in equivalent glucose nor changes in protein and oil concentrations. Results from these experiments can be influenced from soil variables influencing BNF activity. In our study, long-term stable soil attributes were included in the analysis. Organic matter was positively correlated with RAU R6 and the maximum RAU rate but negatively with the AUC (Supplementary Table S1). The negative association between AUC and soil organic matter might be attributed to the soil N mineralized from N organic fraction during the season 41 . Collino et al. 24 compared average BNF in soybean production systems of Argentina and Brazil, and attributed the lower BNF in the former to better soil fertility. Of the five putative mechanisms to account for the metabolic cost of BNF, enhanced sink-driven photosynthesis 13 and reduced root:shoot ratio are the remaining hypotheses to explain the lack of seed yield response with reduced BNF in Santachiara et al. 39 . Implicit in the photosynthesis hypothesis is that soybean seed yield is sink-limited; although unlikely, it requires further research. Enhancing photosynthesis by increasing atmospheric   Table 3. Parameters ± standard error of the evolution curves of RAU fitted for crops with high BNF (>75th quartile), medium BNF (25th-75th quartiles), and low BNF (<25th quartile). RAU max is the peak RAU during the growing season (V2 to R7) reached at t max , thermal time t 0.5 when 50% of RAU max is reached, the maximum rate reached at t m , and the area under the curve AUC. Syx is the standard deviation of the residuals of the fitted curve, and R 2 is the coefficient of determination of the fitted curve (all p < 0.001). a Different letters represent differences between parameters of the same group according to the AICc comparison.  42 reported increases in nodule activity for shadowed white clover (Trifolium repens L.) plants but not in their non-stressed counterparts. It is likely that effects of photosynthetic rates on nodule activity depend on reserve carbohydrates 43 , suggesting a link with the differential trends in HI between favourable and poor environments observed in this study. An alternative, less explored mechanism for the maintenance of seed yield in crops relying on N fixation relative to fertilized crops is the reduction in root:shoot ratio; reduced root:shoot ratio is a generic response of plants to high availability of soil N 15 . 44 and is a major source of oil and protein worldwide. Improving BNF can be achieved by breeding and selection targeting the plant, the N-fixing bacteria, and better matching plant and bacteria [45][46][47] . Selection for maintenance of BNF in dry soil has been proposed to improve seed yield of soybean under drought 48,49 . Sinclair et al. 49,50 combined ureide concentration in petioles and acetylene reduction activity to test this proposition. Selected lines were compared with high-yielding commercial cultivars under broad environmental conditions. Two lines were identified that outperformed commercial checks under water deficit, but trade-offs were apparent under high yielding conditions. In this context, the trade-off between BNF and seed yield mediated by HI needs attention. Solving this trade-off needs quantification of the costs (seed yield reduction), agronomic and environmental benefit of BNF. Selection for high biomass partitioning to seed in genotypes growing under low concentration of soil nitrate represents a possible breeding strategy as higher rates of BNF are expressed. In both soybean 51 and common bean (Phaseolus vulgaris L.) 52 , sensitivity of N fixation to soil mineral N varies with genotype. In alfalfa (Medicago sativa L.), selection for BNF improved plant growth 53 .

Conclusion
Seasonal characterizations over a wide range of agronomic and environmental conditions revealed that N application reduced maximum RAU at R6, particularly for late applications. At the crop level, soybean met the cost of BNF by a reduction in seed yield mediated by lower HI, particularly in stressful environments, and a secondary contribution of reduced seed oil concentration. The lower-level mechanisms underlying shifts in HI associated with BNF warrant further attention.

Methods
Experimental sites and treatments. During the 2016 growing season, soybean N fertilization studies were replicated in 23 sites across the US Midwest in a latitude range from 34°16′ N to 48°14′ N and from 90°29′ W 98°53′ W (Fig. 1, Table 1). Due to the range of latitudes between locations, the length of the growing season differs among sites. Thus, sowing dates and MG were considered following local management practices and recommendations which ranged from 0 to IV due to the large range of latitude in the locations 54 . The seeding rate of 300,000 seeds per hectare, targeted maximum seed yield. Crops were rainfed and received supplementary irrigation to avoid severe water stress (Table 1). On-site weather stations recorded daily temperature, precipitation, and relative humidity; the vapor pressure deficit (VPD, kPa) was estimated using the maximum daily temperature and relative humidity 55 . Soil parameters from every location are presented in Table 1. Data for percentages of clay and sand, organic matter (OM), soil pH and cation exchange capacity (CEC) was extracted from the California Soil Resource Lab (http://casoilresource.lawr.ucdavis.edu, accessed 11 June 2018) using latitude and longitude coordinates from each experiment. Four treatments were established: an unfertilized control, and 112 kg N ha −1 as urea (46-0-0 N-P-K,) at one of three stages: at sowing, at V4 (fourth-leaf), and between R2 (full flowering) to R3 (beginning of pod formation) 56 . The experimental design at each location was a randomized complete block with three replicates. Plot size was 8.4 m long by eight rows at 0.76 m row spacing. The supply of other nutrients was done with N-free fertilizer.
Phenology, biomass, seed yield, harvest index, seed protein and oil concentration. Plant development stages in relation with calendar time is usually referred as phenology. Phenological stages were recorded during the season following Fehr and Caviness 56 . Shoot biomass samples were collected at the R8 stage (full maturity) from 1.6 linear m and fractioned into stem, leaves, and seeds. The relative proportion of seeds to the total shoot biomass was quantified as the harvest index (HI) 57 . Variations on this ratio can be associated with the influence of the environmental effects on seed yield and biomass production. Seed yield was collected from two-central rows at maturity and adjusted to 13.5 g 100 g −1 seed moisture content. Seed samples were collected from harvest for oil and protein determination by near infrared (NIR) spectroscopy using a completely automated Fourier Transform-IR imaging Microscope (Hyperion 3000, Bruker Optics, Ettlingen, Germany) and a sample of >50 seeds. Seed protein and oil concentrations are reported on dry weight basis ( Time units used to measure the progress of RAU during the season was thermal units (degree-days; °Cd) to account for thermal differences in growing conditions and be independent from the temperature in which different developmental stages occurs. A degree-day is the result from every degree on the daily mean temperature above the base temperature 60 . Thus, cumulative thermal time was calculated as: where T max and T min is the maximum and minimum daily air temperature (°C), and Tb is the base temperature (8 °C) 61 . Changes in RAU during the growing season has been described as a sigmoidal pattern with a slow increase early in the season and a maximum attainable between R5 (beginning of seed filling) and R6 (full seed) 62 . However, owing to the large variation in genotypes and growing conditions, RAU at R6 varied widely. To account for this variation, we used RAU at R6 in unfertilized controls to split data into three groups: below the 25 th quartile (low BNF), between the 25 th -75 th quartiles (medium BNF), and above the 75 th quartile (high BNF). The low BNF comprised five sites with RAU below 72%; the medium BNF included twelve sites with RAU from 72 to 88%, and the high BNF group comprised of six sites with RAU above 88%. For the data combined for each group, the seasonal RAU evolution was described with the beta growth function 63 with three parameters: where t is the thermal time from V2 (vegetative leaf), RAU max is the maximum RAU at thermal time t max , and t m is the thermal time for maximum RAU growth rate. Biological meaning on parameters allowed us to make inferences on the magnitude of the N treatments on the RAU dynamics by statistical comparisons. Differences between parameters of equation (2) were tested using the Akaike's Information Criteria (AIC) by performing pairwaise comparisons of individual curves against a global fit. Maximum rate, t 0.5 , and the AUC where compared using the 95% confidence interval. Both RAU and thermal time were estimated through least squares mean analysis by fitting a mixed model with PROC MIXED procedure 64 (lsmeans statement) to adjust corrected means to the factors of the model. For this analysis, treatment and locations were considered as fixed factors, and block was nested within location as a random factor. The goodness of fit of the model was assessed with the coefficient of determination (R 2 ) and the standard deviation of residuals (Syx) 65 . Using equation (3), we derived three related traits: the AUC to integrate seasonal N fixation 66 normalized to the maximum of the data set; t 0.5 the thermal time when RAU is 50% of RAU max; and the maximum rate of RAU expressed in changes of % units of RAU per unit of thermal time (°Cd) 63 : The first derivative of equation (3) with respect to thermal time can be solved to calculate the RAU rate changes across the growing season 63 : Same approach has been utilized to describe other biological process such as N uptake rate in corn 67 or grain growth rate 68 . Analysis of treatment effects and associations between traits. Analysis of variance (ANOVA) was used to investigate effects of treatments on crop traits (seed yield, biomass, HI, RAU R6, seed protein and oil concentration). Sources of variation in ANOVA included N treatment, location, and their interaction as fixed factors, and block as a random effect nested within location; this analysis was implemented by using the R software (version 3.4.0, lme4 package, lmer function) 69 The effects of BNF on seed yield, biomass, HI, seed protein and oil concentration were analyzed in two steps. First, due to the geographical distribution of the experiments, responses on crop traits are confounded with the different duration of the developmental stages. Thus, effects of phenology were captured with non-linear models: where Y is the trait, X is the thermal time to R6 (°Cd), and a, b, TT o , and d, are parameters. Next, linear regressions and quantile regressions were fitted between residuals of these models and RAU R6 . This simple approach on the use of residuals allows to netting out 71 the effect of the phenology on the traits observed when are regressed against RAU R6 . There are relationships on other parts of the distribution of the response variable that can provide more complete view of the processes studied besides of the mean effect observed. Slopes from quantile regression analysis estimate the changes at the maximum and minimum response that can be missed when other regression methods are used 72 . Thus, regressions for 0.99 and 0.01 quantiles capture the boundaries of the relationships, and were fit in R (quantreg package 73 ). The rest of linear and non-linear regression analyses, computation of AUC, and estimation and comparison of parameters from equations (3) and (5) were performed using GraphPad Prism 74 . Principal component analysis (PCA) was used to analyze general associations between traits allowing the identification of any grouping association within the data set when environmental and crop attributes are analyzed together 75 . Data were classified according to MG, which in turn had a geographical correlation (Table 1). Traits included RAU max , t max , maximum rate of RAU evolution, seed protein and oil concentration, AUC, residuals from seed yield, biomass, and HI vs thermal time to R6 relationship, and soil attributes (clay, sand, organic matter, pH, and CEC). Principal component analysis was fit using the "FactoMineR" package in R 76 . Pearson correlation coefficients were calculated to complement associations found in the PCA.