Unraveling the genetic architecture for carbon and nitrogen related traits and leaf hydraulic conductance in soybean using genome-wide association analyses

Background Drought stress is a major limiting factor of soybean [Glycine max (L.) Merr.] production around the world. Soybean plants can ameliorate this stress with improved water-saving, sustained N2 fixation during water deficits, and/or limited leaf hydraulic conductance. In this study, carbon isotope composition (δ13C), which can relate to variation in water-saving capability, was measured. Additionally, nitrogen isotope composition (δ15N) and nitrogen concentration that relate to nitrogen fixation were evaluated. Decrease in transpiration rate (DTR) of de-rooted soybean shoots in a silver nitrate (AgNO3) solution compared to deionized water under high vapor pressure deficit (VPD) conditions was used as a surrogate measurement for limited leaf hydraulic conductance. A panel of over 200 genetically diverse soybean accessions genotyped with the SoySNP50K iSelect BeadChips was evaluated for the carbon and nitrogen related traits in two field environments (Athens, GA in 2015 and 2016) and for transpiration response to AgNO3 in a growth chamber. A multiple loci linear mixed model was implemented in FarmCPU to perform genome-wide association analyses for these traits. Results Thirty two, 23, 26, and nine loci for δ13C, δ15N, nitrogen concentration, and transpiration response to AgNO3, respectively, were significantly associated with these traits. Candidate genes that relate to drought stress tolerance enhancement or response were identified near certain loci that could be targets for improving and understanding these traits. Soybean accessions with favorable breeding values were also identified. Low correlations were observed between many of the traits and the genetic loci associated with each trait were largely unique, indicating that these drought tolerance related traits are governed by different genetic loci. Conclusions The genomic regions and germplasm identified in this study can be used by breeders to understand the genetic architecture for these traits and to improve soybean drought tolerance. Phenotyping resources needed, trait heritability, and relationship to the target environment should be considered before deciding which of these traits to ultimately employ in a specific breeding program. Potential marker-assisted selection efforts could focus on loci which explain the greatest amount of phenotypic variation for each trait, but may be challenging due to the quantitative nature of these traits.


Background
Soybean [Glycine max (L.) Merr.] seeds are an important source of protein and oil for a range of applications. Drought stress is the most important abiotic factor affecting soybean production, and can cause large decreases in yield [1]. Use of irrigation during drought stress could ameliorate this issue; however, less than 10% of U.S. soybean hectares are irrigated [2]. Therefore, the development of soybean cultivars that can withstand periods of drought stress is necessary to protect yield when water resources are limited.
Certain morphological and physiological traits could reflect the ability of soybean plants to better tolerate drought stress. Carbon isotope composition has been previously identified as a useful screening method to understand photosynthetic tradeoffs and water-saving capabilities of C3 plant species in certain environments [3][4][5][6][7]. C3 plants readily assimilate the 12 C isotope of carbon in photosynthesis, and therefore discriminate against the heavier 13 C isotope, which constitutes only around 1% of the atmosphere [4]. Carbon isotope composition can be expressed as either carbon isotope discrimination (Δ 13 C, CID) or carbon isotope ratio (δ 13 C). Carbon isotope composition has been used as an indirect method for selection of genotypes with improved productivity in droughtstressed environments. However, it should be noted that in some cases CID has not been a good indicator for drought tolerance or did not produce consistent genotypic rankings across environments [8][9][10].
Additionally, previous genome-wide association studies (GWAS) and quantitative trait locus (QTL) mapping studies have identified genomic regions controlling carbon isotope composition in soybean. In one of these studies, 373 diverse maturity group (MG) IV soybean genotypes were grown in four environments and 39 single nucleotide polymorphisms (SNPs) were identified with GWAS that had significant association with δ 13 C in at least two environments [11]. Another study using the same set of accessions and phenotypic data, but with2 0,000 additional SNP markers and a different GWAS model, found 54 environment-specific SNPs tagging 46 putative loci for δ 13 C [12]. Previous QTL mapping in soybean identified five loci controlling CID [13].
Soybean is a legume which uses a symbiotic association with bradyrhizobia to fix N 2 from the atmosphere. This nitrogen fixation provides a supply of nitrogen (N) to the plant that is used for growth and development, as well as providing nitrogen in the crop residue for subsequent crops when soybean is used in a crop rotation. However, symbiotic N 2 fixation can be affected by limited water availability, and certain soybean genotypes are more sensitive than others in regards to N 2 fixation during drought stress [14][15][16][17][18]. A previous simulation study that investigated the benefits of altered soybean drought traits found that sustained N 2 fixation during water deficits had the most consistent and greatest yield advantage compared to four other traits using 50 years of weather data across U.S. soybean growing regions [19].
Using a three-stage screening process, [20] identified eight soybean genotypes with superior N 2 fixation during water deficits. In addition, PI 471938 has been reported to have tolerant N 2 fixation as soil dries [21]. Differences in the amount of N present in leaf tissue have previously been used as a way to determine a soybean genotype's sensitivity to N 2 fixation during drought conditions, with lower foliar N concentrations having superior fixation during water deficits [14,17,18]. This could be due to genotypes with higher plant N concentrations under well-watered conditions being closer to a threshold N level in the plant that can trigger a negative feedback of nitrogen compounds decreasing N 2 fixation rate. In contrast, genotypes with lower plant N concentrations may continue to fix nitrogen during water deficits due to a lack of this feedback. Four QTLs for foliar N concentration were previously identified on Chr 13, 16, and 17 using a 'KS4895' × 'Jackson' RIL population [22].
Nitrogen isotope composition (δ 15 N) could be a useful evaluation tool given that 15 N is present at much greater levels in soil compared to the atmosphere [23][24][25]. The fraction of 15 N found in a soybean plant would be decreased if it is actively fixing N 2 from the atmosphere, and could be an indicator of how much nitrogen fixation is affected by drought stress [26]. A previous association mapping study using 373 soybean genotypes in MG IV found 19 and 17 SNP markers significantly associated with N concentration and the fraction of N derived from the atmosphere (Ndfa), respectively, that were found in at least two of the four environments tested [26].
Leaf hydraulic conductance is defined as the water flux through the leaf per unit water potential driving force, and is a measure of how readily water flows through the leaf [27]. Limited leaf hydraulic conductance is a trait related to soybean drought tolerance that results in conserved soil moisture for use during subsequent water deficits. According to previous research, decreased hydraulic conductance allows certain soybean plants, namely PI 416937, to conserve soil water and express a slow canopy-wilting phenotype in the field after extended periods with little to no precipitation [28]. Additionally, it was hypothesized that differences in hydraulic conductance were a result of different populations of aquaporins, water-conducting membrane proteins that are involved in water movement through cell membranes. It was suggested that these aquaporin populations could be differentiated due to differences in sensitivity to exposure to certain chemical inhibitors [29]. Subjecting de-rooted soybean shoots to a silver nitrate (AgNO 3 ) solution under high vapor pressure deficit (VPD) conditions resulted in some genotypes expressing a decreased transpiration rate, and it was hypothesized that this decrease in transpiration was a result of silver ions blocking silver-sensitive aquaporins. PI 416937, a slow-wilting genotype with low hydraulic conductance, exhibited an insensitivity to silver nitrate by not decreasing its transpiration rate when subjected to the inhibitor solution [30]. Given the possible relationship of the transpiration response to silver nitrate and hydraulic conductance, soybean genotypes could be characterized using this procedure to potentially differentiate aquaporin populations and identify drought tolerant germplasm. A previous QTL mapping study identified four QTLs explaining 17.7 to 24.7% of the phenotypic variation for the limited leaf hydraulic conductance trait using transpiration response to silver nitrate as the measurement for the trait [31].
In this study, a genetically diverse panel of over 200 soybean genotypes was evaluated for δ 13 C, δ 15 N, and foliar nitrogen concentration from leaf samples collected in two field environments. Additionally, this panel was evaluated for transpiration response to silver nitrate under high VPD conditions in a growth chamber. The objectives of this study were to identify genomic regions controlling these traits using genome-wide association analyses, validate genomic loci for these traits across environments or studies, and identify genotypes in the panel which have favorable breeding values for these traits.
Results δ 13 C, δ 15 N, and N concentration Carbon isotope composition (δ 13 C), nitrogen isotope composition (δ 15 N), and foliar nitrogen (N) concentration were evaluated in two field environments (GA-15 and GA-16). Based on the analyses of variance (ANOVA), genotypes, environments, and their interaction were statistically significant (p < 0.05) for all carbon and nitrogen related traits (Table 1). Genotype mean values within environments of δ 13 C ranged from − 29.97 to − 25.14‰ (Fig. 1), and had a correlation of r = 0.74 between environments. Broad-sense heritability of δ 13 C on an entry-mean basis for each environment was 61% (GA-15), 72% (GA- 16), and 62% across both environments (Table 2). δ 15 N had a correlation of r = 0.28 between environments, and ranged from − 1.23 to 4.50‰ based on mean genotype values within environments (Fig. 1). Heritability for δ 15 N was lower than for all other carbon and nitrogen related traits at 24% (GA-15), 40% (GA-16), and 17% across both environments (Both) ( Table 2). The range of leaf nitrogen concentrations observed for genotype means within environments was from 16.67 to 55.45 g kg − 1 , and the correlation between the two environments was r = 0.73. Broad-sense heritability for N concentration was between 63 and 73% ( Table 2).
In general, these carbon and nitrogen related traits had fairly strong relationships with one another. Using best linear unbiased predictors (BLUP) values calculated from across both environments, correlations between the carbon and nitrogen related traits were from r = − 0.52 to 0.71 (Table 3). The most negative correlation (r = − 0.52) was between δ 13 C and δ 15 N, and the most positive correlation (r = 0.71) was observed between δ 13 C and N concentration (Table 3).
PI 398823, a MG IV accession had the highest breeding value for δ 13 C using the sum across the two individual environments (Additional file 1). In addition, PI 416937, a slow-wilting check genotype, had a relatively high breeding value for this trait and ranked within the top 10% of genotypes tested (Additional file 1).
A MG VI accession from China, PI 567377B, had the most negative (favorable) breeding value for N concentration using the sum across both individual environments (Additional file 1). PI 471938, which was previously identified as a genotype possessing nitrogen fixation drought tolerance [21,33], had the 40th lowest breeding value for N concentration (Additional file 1).

Transpiration response to silver nitrate aquaporin inhibitor
Normalized decrease in transpiration rate (NDTR) values ranged from − 2.33 to 1.00 within individual replications (Fig. 2), and from − 0.99 to 0.48 based on genotype means. Genotype effects were statistically significant (p < 0.05) ( Table 1), and broad-sense heritability on an entry-mean basis was 17% ( Table 2). Using BLUP values across replications and environments, the relationships between NDTR in response to AgNO 3 and the carbon and nitrogen related traits were also evaluated ( Table 3). Silver nitrate NDTR was not correlated (r = − 0.02 to 0.05) with the previously described carbon and nitrogen related traits.
Twelve out of the 15 accessions with the most negative breeding values for transpiration response to AgNO 3 originated from China (Additional file 1). PI 416937 was previously identified as a genotype with a transpiration response that is relatively insensitive to silver nitrate [30], and ranked 123rd based on NDTR breeding values.

GWAS of carbon and nitrogen related traits
A total of 35 unique SNPs tagging 32 loci were identified either in individual environments or when using the BLUP calculated across both environments for δ 13 C (Additional file 2 and Table 4). Two SNPs for δ 13 C (ss715587736 and ss715587739) on Chr 4 were in the same genomic region, and were found in GA-15 and across both environments, respectively (Table 4). Of all other SNPs identified for δ 13 C, each SNP tagged a single genomic region, with the exception of two SNPs identified on Chr 4 and 16. The allelic effects across all significant (p < 0.0001; −log 10 (P) > 4) SNPs ranged from − 0.19 to 0.13 (Table 4), with all significant SNPs explaining a total of 29-44% of the variation, depending on the environment (Table 4).
For δ 15 N, 23 loci were identified in the GWAS (Additional file 2 and Table 4). Depending on the environment, 36 to 51% of the phenotypic variation for δ 15 N was explained by the significant (p < 0.0001; −log 10 (P) > 4) SNPs. The allelic effects ranged from − 0.14 to 0.11 for the SNPs significantly associated with δ 15 N (Table 4). One SNP (ss715635458) was found for δ 15 N both in GA-16 and using the across both environments BLUPs   (Table 4). All other SNPs identified tagged a single genomic region. Twenty seven SNPs tagging 26 loci were identified in the GWAS for nitrogen concentration (Additional file 2 and Table 4). One SNP (ss715610522) was identified in both an individual environment (GA-15) and with the BLUP value from across both environments (Table 4). All other SNPs tagged a single genomic region, except for two SNPs (locus 17) on Chr 13. Allelic effects for nitrogen concentration ranged from − 1.33 to 1.46 (Table 4). Phenotypic variation explained (R 2 ) across all significant SNPs for N concentration was 50, 35, and 21% for GA-15, GA-16, and across both environments (Both), respectively.

GWAS for transpiration response to silver nitrate aquaporin inhibitor
Nine SNPs tagging nine loci were significantly (p < 0.0001; −log 10 (P) > 4) associated with NDTR following silver nitrate treatment ( Fig. 3 and Table 5). Thirty one percent of the phenotypic variation for the trait was explained by these nine SNPs. The allelic effects for these significant SNPs ranged from − 0.04 to 0.03 (Table 5).

Candidate genes for carbon and nitrogen related traits
For every trait evaluated, candidate genes were identified within plus or minus 10 kb (approximately spans the mean distance between all markers) of the SNPs with the lowest p-value (highest -log 10 (P)) in each environment Table 3 Correlations among canopy wilting, carbon isotope composition (δ 13 C), nitrogen concentration, nitrogen isotope composition (δ 15 N), and normalized decrease in transpiration (NDTR) rate in response to silver nitrate (AgNO 3 )    Allelic effects were calculated by taking the difference in mean phenotypic value between the two alleles at a particular SNP, and the direction, negative or positive, of the allelic effect estimates are relative to the alphabetical order of the nucleotides at each particular marker f Environment written as location-year and using across environments data. Eight, six, and seven candidate genes were identified for δ 13 C, δ 15 N, and N concentration, respectively, near these most significant SNPs (Additional file 3).

Rationale for trait evaluation
In this study, a genetically diverse panel of over 200 soybean genotypes was evaluated for δ 13 C, δ 15 N, and nitrogen concentration from leaf samples collected in two different field environments. In addition, this panel was also evaluated for transpiration response to silver nitrate under high vapor pressure deficit conditions in a growth chamber. Using genome-wide association mapping, genomic regions were identified controlling each of these different drought tolerance related traits and the results were compared to previous mapping studies for these traits. In addition, genotypes in the panel were identified which possessed favorable breeding values for these drought tolerance related traits. Carbon isotope composition can relate to photosynthetic tradeoffs that result from variation in water-saving capabilities. Nitrogen fixation can be highly sensitive to drought stress [14][15][16], and above-ground measurements such as nitrogen concentration and nitrogen isotope composition might relate to nitrogen fixation rate and soybean drought tolerance [17,24,25]. The amount of 15 N found in a soybean plant would be decreased if it is actively fixing N 2 from the atmosphere, and lower N concentrations have been shown to correlate with superior fixation during water deficits. However, given the high protein content of soybean, and the amount of nitrogen required to produce protein in seed, lower N concentrations could well be a poor trait for a soybean genotype to possess. Water-transporting proteins called aquaporins are involved in water movement through cell membranes [34], and populations of aquaporins in soybean lines can vary as detected by transpiration response to chemical inhibitors such as silver nitrate [30,31,35]. It is hypothesized that insensitivity to silver nitrate is correlated with the limited leaf hydraulic conductance trait, a beneficial trait associated with water conversion and improved drought tolerance in certain environments [28,29]. All of these traits were evaluated in the current study in order to develop insight about the genetic architecture of these drought tolerance related traits and identify germplasm with favorable breeding values for these traits.  Values for δ 13 C were in a similar range to those observed in two previous carbon isotope association mapping studies [11,12] (Fig. 1). The range of values observed for nitrogen concentration was wider and concentrations were higher compared to those observed in a previous study [26]. Direct comparisons to [26] were not able to be made for δ 15 N due to differences in the units used for these measurements. Analyses of variance (ANOVA) showed that genotype, environment, and their interaction were statistically significant (p < 0.05) for all carbon and nitrogen related traits evaluated with the association panel (Table 1). Although these genotype-byenvironment interactions were significant (p < 0.05), correlations were generally high between the two environments. Correlations for δ 13 C and nitrogen concentration were all above r = 0.70 between the two environments tested, indicating the genotypes performed similarly across environments. The lowest correlation was for δ 15 N at r = 0.28, which suggests this trait could be subject to environmental influence, such as nitrogen levels in the soil. Heritability for δ 15 N was substantially lower and ranged from 17 to 40% (Table 2). This lower heritability for δ 15 N could potentially be explained by the fact that we did not adjust our values to a non-nodulating reference crop, and that these values are also affected by field variation in soil nitrogen concentration [36]. However, heritability estimates for all of these carbon and nitrogen related traits are comparable to the values observed in other studies [11,12,26].

Transpiration response to AgNO 3
Low or negative DTR to silver nitrate values (transpiration less affected by AgNO 3 ) have been previously correlated with limited leaf hydraulic conductance, which is a beneficial trait in certain drought stress environments [29]. Given the hypothesis that silver nitrate blocks only specific aquaporins and reduces transpiration, and that most previously reported DTR values were positive, we observed an unexpected distribution of NDTR values given that many of the genotypes we tested had negative non-normalized DTR (negative NDTR). This could indicate that silver nitrate blocked some aquaporins as expected, but in some genotypes this blockage resulted in a stimulus in the number or activity of other silverinsensitive aquaporins. However, this hypothesis needs further experimental investigation.
Analyses of variance found that genotype effects were statistically significant (p < 0.05) ( Table 1), and heritability for this trait was 17% (Table 2). This low heritability estimate could have been a result of a technical issue or that this phenotyping method may not be a reliable proxy for limited leaf hydraulic conductance, and would make it difficult for soybean breeders to make effective selection for this trait. One potential technical issue which could have explained the low heritability observed was variation in VPD throughout each experimental replication and between each of the eight replications, as well as VPD values lower than the desired 3.00 kPa for our protocol. As shown in Table 6, average VPD by replication ranged from 1.56 to 2.33 kPa. In addition, VPD variation within each replication was relatively stable, but varied by as much as 0.5 kPa during a single replication due to the size of the walk-in growth chamber and its ability to maintain the environmental settings we aimed to achieve in the protocol. While temperature remained relatively constant throughout the experiments, relative humidity (RH) was more variable, and was the primary driver in the varying VPD observed ( Table 6). Given the genotypes tested were a diverse panel from different maturity groups and geographic origins, there was some variation in the size of the plants as they were growing in the greenhouse in preparation for the experiments. This variation in size was accounted for in our DTR calculations, because each plant's difference in transpiration rate between water and silver nitrate solution was relative to itself. However, it is still worth noting that plant size differences could cause some degree of soil moisture deficit in the relatively small pots we used to grow the plants to V3-V4 stage, and may be another factor to explain the low heritability we observed. In addition, during the process of cutting the shoot of the soybean plants from the roots it is possible that some plants were embolized. However, as part of our protocol, we made a second cut underwater away from our initial cut to help potentially avoid this issue.

Comparison to previous mapping results for carbon and nitrogen traits
Given that FarmCPU uses the most significant markers as covariates in the GWAS model, SNPs are seldom identified within the same LD block for an environment-specific dataset. However, two genomic regions were found both in individual environments and when using the across both environments BLUP data for these carbon and nitrogen related traits. Significant (p < 0.0001; −log 10 (P) > 4) SNPs for carbon and nitrogen related traits were found on all 20 soybean chromosomes (Table 4). Previously identified QTLs for CID are numbered with their approximate physical positions on the SoyBase website (www.soybase.org). Locus 32 identified with GWAS for δ 13 C in the current study is found within the CID 1-5 QTL on Chr 19 identified in [13] (Table 4). A comparison of SNPs significantly associated with δ 13 C from two previous association mapping studies [11,12] and the current study was conducted (Fig. 4a). Two SNPs on Chr 6 and 11 from the current study are near significant markers identified in [12], and one SNP on Chr 13 and another SNP on 18 were found near the significant SNPs for δ 13 C in [11].
No QTLs for δ 15 N identified with linkage mapping are reported on the SoyBase website. One previous linkage mapping study for foliar nitrogen concentration identified four QTLs, of which one QTL on Chr 16 was 256 kb away from locus 21 identified in the current study [22]. A comparison of SNPs identified for nitrogen related traits in a previous association mapping study [11] and the current study was also performed (Fig. 4b). SNPs on Chr 9 and 15 were found in common for δ 15 N in the current study and nitrogen derived from the atmosphere (Ndfa) in [26]. No SNPs were within 1 Mb of previously identified genomic regions for nitrogen concentration. Additionally, when making comparisons only across studies and different nitrogen related traits, only two regions on Chr 15 and 16 had common SNPs within 1 Mb of each other. Within the current study only, two regions contained nitrogen related significant (p < 0.0001; −log 10 (P) > 4) SNPs within 1 Mb of each other on Chr 13 and 20 ( Table 4). The relatively small number of consistent associations across these studies could be due to differences in the maturity groups tested or the tissue collection method. However, the consistent QTLs and genomic regions across environments, studies, and traits, along with SNPs explaining a high amount of phenotypic variation in the current study could be useful as breeding targets for these carbon and nitrogen drought tolerance related traits.
Genetic mapping for transpiration response to AgNO 3 and proximity of identified regions to aquaporin gene models This is the first report of association mapping for this trait to the authors' knowledge in any crop species. A previous QTL mapping study for limited leaf-hydraulic-conductance identified QTLs on Chr 3, 5, 10, and 12 [31]. The locus identified on Chr 12 in the current study is located approximately 2 Mb away from the Chr 12 QTL from that previous study. A lack of overlap in the genomic regions observed in these two studies could be due to differences in the populations utilized for the mapping, and could also be affected by the low heritability for this trait ( Table 2). A search on Phytozome for gene models with a functional annotation which contained the word "aquaporin" was also conducted given the hypothesized relationship between this limited leaf hydraulic conductance trait and aquaporins, and found 88 gene models. The physical locations of these gene models and the loci identified in the current study with association mapping were compared (Fig. 4c). Three SNPs identified in the GWAS were within 1 Mb of four gene models with an aquaporin functional annotation. These regions could be further investigated to see how this trait relates to aquaporins.

Candidate genes at identified genomic regions for carbon and nitrogen related traits
A total of 21 gene models were identified near the most significant SNP across each trait and environment tested. A gene model located at locus 11 for carbon isotope composition, Glyma.10 g047500, is a protein phosphatase 2C family protein (Additional file 3). This gene family has been shown to function at the intersection of drought, oxidative, and heat shock stresses in tobacco [37]. The gene model Glyma.09 g043900 is a transducing/WD40 repeat-like superfamily protein located near locus 9 (ss715603834) associated with nitrogen isotope composition (Additional file 3). A report in Arabidopsis thaliana showed that a member of the WD40 gene family functions in drought stress tolerance by modulating nitric oxide accumulation and stomatal closure [38]. A C2H2-type zinc finger family protein gene (Glyma.12 g065800) located at locus 15 is associated with nitrogen content (Additional file 3). In rice, a zinc finger transcription factor, drought and salt tolerance (DST), was shown to play a role in stomata-regulated abiotic stress tolerance [39]. These gene models could be potential targets for understanding and improving these drought tolerance related traits given their relationship with drought stress tolerance response or enhancement.

Relationship between drought tolerance related traits
Another measurement related to soybean drought tolerance, canopy wilting, was added to the Table 3 correlation matrix using data from [32]. This additional data from the same field experiments provides another trait to compare to carbon and nitrogen related traits and NDTR in response to silver nitrate treatment. Canopy wilting and NDTR to silver nitrate had relatively low correlations with each of the other traits evaluated and with one another. A previous study also found that there was not a consistent relationship among genotypes within slow-or fast-canopy wilting groups and CID [40]. Drought tolerance is a complex, quantitative trait, so it is expected that multiple different traits and loci are responsible for soybeans' ability to withstand water deficit stress. Wm82.a2 version of the soybean genome. SNPs identified in GWAS from current study that met -log 10 (P) > 4 significance threshold are shown as larger circles for a) carbon isotope composition (δ 13 C), b) Nitrogen concentration and nitrogen isotope composition (δ 15 N), and c) normalized decrease in transpiration rate (NDTR) in response to silver nitrate treatment. Smaller circles represent SNPs identified in a) [11,12], b) [26] that were converted from version 1 to 2 physical positions of the soybean genome assembly, and c) location of gene models with the term "aquaporin" in their functional annotation from Phytozome v12.1. BARC_1.01_Gm20_46575262_G_A identified for nitrogen concentration in [26] does not have a perfect match in the version 2 assembly, and therefore was excluded from this comparison

Breeding implications
Many different genotypes were identified in the current study with favorable breeding values for drought tolerance related traits and could be utilized by breeders to improve soybean drought tolerance directly with forward breeding or be used as parents to create mapping populations to further understand the genetic architecture for these traits. Genotypes with positive breeding values for δ 13 C, negative breeding values for N concentration, negative breeding values for δ 15 N, and accessions with lower NDTR values and low negative breeding values could be candidate parents to use for drought-tolerance improvement in a soybean breeding program. However, the challenge as a breeder would be to determine which trait(s) to target given the quantitative nature of the genetic architecture for many traits that could lead to soybean drought tolerance improvement, and some of these traits could be associated with poor agronomic performance.
In addition, accessions in the current study often had favorable breeding values for certain traits, but then also had less favorable breeding values for other traits (Additional file 1). As a reference point, PI 416937, a genotype previously identified as possessing the slow canopy-wilting trait [41], was ranked as the 133rd best accession tested based on an overall median rank across breeding value ranks for canopy wilting, carbon isotope composition, nitrogen concentration, nitrogen isotope composition, and NDTR in response to silver nitrate (Additional file 1). It ranked 69th best for canopy wilting and 15th best for carbon isotope composition, but ranked 189th for nitrogen concentration, 140th for nitrogen isotope composition, and 123rd best for transpiration response to silver nitrate (Additional file 1). One hundred thirty-two accessions with overall median ranks lower than PI 416937 were identified in this research (Additional file 1). To make selections based on multiple traits an index accounting for trait heritability, economic importance, and genetic and phenotypic correlations among the traits would likely need to be employed with consideration for phenotyping costs and genotype by environment interactions for these traits. Ultimately, a breeder may need to weight traits according to which would provide the best drought tolerance in their given target environment, and then utilize the germplasm and genomic regions identified for that specific trait.

Conclusions
Genome-wide association analyses were conducted for δ 13 C, δ 15 N, and nitrogen concentration from two environments using over 200 genetically diverse soybean genotypes. Thirty two, 23, and 26 loci were identified for δ 13 C, δ 15 N, and nitrogen concentration, respectively. One locus detected with the GWAS for δ 13 C was co-located with a previously identified QTL for CID, and four SNPs were near SNPs found in previous association mapping studies. Two SNPs for δ 15 N were found in the GWAS near genomic regions identified in an association mapping study for nitrogen related traits. Nine SNPs tagging nine loci were identified with a GWAS approach for normalized DTR to silver nitrate, and three of the SNPs identified were found near four aquaporin related gene models. Breeding values calculated with the significant SNPs from the GWAS enabled the identification of accessions which possess favorable combinations of alleles for these drought tolerance related traits. The genomic regions and germplasm identified in this study, especially those found in common across environments, studies, and traits, can be used to understand the genetic architecture for these traits and by soybean breeders to improve drought tolerance.

Soybean populations
An association panel of 211 genetically diverse soybean genotypes was evaluated for transpiration response to a silver nitrate solution. The panel was previously described in [32], but with the addition of two lines and replacement of 10 other lines that did not produce enough seed for the field evaluations of drought tolerance related traits described in [32] and also in the current study. This panel was selected based on SoySNP50K genotype data to be genetically diverse, consisted mostly of maturity group (MG) VI-VIII plant introductions, and included drought tolerant and susceptible genotypes. One hundred ninety-five and 205 of the soybean genotypes described in [32] were evaluated in 2015 and 2016 in Athens, GA, respectively, for carbon and nitrogen related traits in the field. The majority of these lines had not previously been evaluated for drought-tolerance related traits, and are later maturing lines than those previously tested (MG IV) and used for association mapping of these traits [11,12,26].

Isotope analysis and sample collection
Leaf samples were collected from field plots of the association panel grown in Athens, GA in 2015 (GA-15) and 2016 (GA-16) and used for stable isotope analysis. More information about sowing dates, row spacing, and management of these plots can be found in [32]. Based on soil sample testing, no fertilizer was added to the field in 2015, and a 4-15-30 fertilizer was applied at a rate of 392 kg ha − 1 in 2016 prior to sowing. These plots were grown under rain-fed conditions and experienced intermittent drought stress periods in both years. In 2015, the leaf samples were collected on 23 September and on 12 September in 2016. All of the soybean genotypes in the panel were in reproductive growth stages (R3-R6) at the time of sample collection. Five leaves were randomly selected from each of the two-row plots at the third trifoliolate leaf below the top of the plants. These leaves were placed in seed envelopes, and stored in a − 20 0 C freezer until they could be processed at a later date. For isotope analysis, 100-150 samples were processed at one time by transferring leaf samples to 50 ml Falcon tubes and placing them in a lyophilizer for two days to freeze dry. The samples were then ground to a fine powder by placing 4.5 mm zinc plated BBs in the tubes and grinding them using a Geno/Grinder (SPEX SamplePrep, Metuchen, New Jersey, USA). Immediately before using this ground leaf tissue for isotope analysis, the tubes were placed in a drying oven to ensure all residual moisture was removed. In an effort to further keep out moisture, the Falcon tube caps were wrapped with Parafilm immediately after this second drying step.
Stable isotope analysis was then performed using a Carlo Erba NA1500 CHN combustion analyzer coupled to a Delta V isotope ratio mass spectrometer via the Conflo III open split interface. Three experimental replications of the dry leaf tissue of each genotype were analyzed at the Center for Applied Isotope Studies, University of Georgia, Athens, GA. A detailed protocol for the procedure can be found at http://sisbl.uga.edu/ ratio.html. The quantity of 13 C in the leaf samples was compared to a reference standard Pee Dee Belemnite, and these δ 13 C values were used for further analyses. δ 13 C was expressed in units per mil (‰) using the following equations [4]: The quantity of 15 N in the leaf samples was compared to air and expressed in units per mil (‰) according to the following equations: Nitrogen concentration was expressed as g kg − 1 .

Evaluation of response to silver nitrate inhibitor
Soybean plants for evaluation of transpiration response to silver nitrate were grown in a greenhouse at the University of Georgia in Athens, GA, USA under a 16 h day and eight hours night lighting regime. Three seeds of each genotype were sown in 32 oz. styrofoam cups using a Fafard 2B soil media (Sun Gro Horticulture, Agawam, MA, USA). Approximately 1.5 weeks after seedling emergence, the plants were thinned to one plant per cup and maintained under well-watered conditions by watering each pot twice daily until the soil reached water holding capacity. Once the soybean plants reached the V3-V4 growth stage (approximately four weeks after sowing), the tests for response to the silver nitrate inhibitor began [29].
The tests were conducted over two days. In the afternoon of the first day, the soybean plants were removed from their growing media in the greenhouse and de-rooted using clippers. A second cut on the stem was then made underwater adjacent (1-3 cm away) to the first cut using a razor blade. The remaining shoot was then placed in a 250 mL Erlenmeyer flask filled with deionized water and the mouth of the flask was sealed with Parafilm to avoid water evaporation. Plants in flasks were then placed in a walk-in Conviron growth chamber at approximately 20 0 C and 60% relative humidity (RH) overnight in dark conditions.
In the morning of day 2, the growth chamber settings were adjusted to turn the lights on, raise temperature to 30 0 C, and decrease RH to 30% to obtain a higher vapor pressure deficit (VPD) in the growth chamber. The observed VPD for the chamber was between 1.56-2.33 kPa across replications of the experiment ( Table 6). The plants were allowed to acclimate to the higher VPD condition for 60 min. Then, each flask/soybean was weighed inside the growth chamber using a balance with a resolution of 0.001 g in order by flask number. Sixty min after the first weighing, they were weighed again in the same order to determine the transpiration rate in water (TR W ). Each soybean shoot was then transferred to a 60 mL amber glass bottle containing a 200 μM solution of silver nitrate (AgNO 3 ) under semi-dark conditions. This AgNO 3 solution concentration was previously shown to best differentiate the transpiration response of drought tolerant versus susceptible soybean plants in [29]. Parafilm was again used to seal the mouth of the amber bottles to avoid evaporation and spilling of any chemical. Then, the plants were returned to the growth chamber and allowed to acclimate to the inhibitor treatment for 60 min. The amber bottles with shoots were then weighed for their initial weight in order by bottle number. After approximately 120-160 min, the bottles were reweighed in bottle order to determine the transpiration response to the silver nitrate inhibitor (TR I ). Differences in the amount of time that elapsed between weight measurements were accounted for in the TR W and TR I calculations by changing the denominator in increments of minutes. Decrease in transpiration rate (DTR, %) was then calculated as follows: Due to limitations in the size of the walk-in growth chamber and ability to weigh the flasks/bottles in an orderly and timely fashion, eight separate replications of this experiment were conducted (Table 6). Each replication consisted of the entire panel of 211 soybean genotypes, and the flask/bottle order was randomized for each replication. To account for small differences in the range of DTR among the eight replicate experiments due to plant size and environmental differences with each replication, the results were normalized against the genotype with the highest DTR value within each replication using the following equation: Genotype data and quality control The association panel was genotyped with the SoySNP50K iSelect BeadChip [42]. DNA extraction and genotyping procedures for this panel were conducted as described in [32]. A total of 42,079 genome-wide SNP markers resulted from the genotyping effort, with most marker data being downloaded from SoyBase [43]. Markers with minor allele frequencies (MAF) lower than 0.05 were eliminated leaving 35,262 SNP markers for the association analysis of transpiration response to silver nitrate. For the carbon and nitrogen related traits, 35,234 (Both), 35,101 (GA-15), and 35,219 (GA-16) markers were used after eliminating markers with MAF lower than 0.05. The number of markers varied, because certain SNPs with a MAF close to 0.05 were either included or excluded depending on the number of entries tested in the given environment. Physical positions are based on the Glyma.Wm82.a2 version of the soybean genome.

Statistical analyses
Analyses of variance (ANOVA) was conducted using PROC GLM in SAS version 9.4 (SAS Institute Inc., Cary, NC, USA). For the response variables relating to carbon and nitrogen traits, genotype was treated as a fixed effect, and environment, genotype-by-environment interaction, and replication within environment were random effects. For transpiration response to silver nitrate, a model was created with genotype as a fixed effect and replication as a random effect, with NDTR as the response variable. Broad-sense heritability was calculated on an entry-mean basis according to [44] with the variance components being calculated with PROC MIXED of SAS 9.4 using a model where all variables were treated as random.
Best linear unbiased predictors (BLUPs) were calculated from both across and within environments and used as the phenotype values for subsequent GWAS analyses. The BLUP calculations for carbon and nitrogen related traits across both environments were performed using JMP Pro (JMP®, Version 13, SAS Institute Inc., Cary, NC, USA). The model was built by treating genotype, environment, genotype-by-environment, and replication within environment as random variables using the Standard Least Squares personality and REML method. For individual environments for carbon and nitrogen related traits and transpiration response to silver nitrate, genotype and replication were used as variables and treated as random to calculate BLUPs.

Genome-wide association analyses
Fixed and random model Circulating Probability Unification (FarmCPU) was used to perform the genomewide association analyses for all traits evaluated [45]. FarmCPU is an R package that implements a multiple loci linear mixed model incorporating a modified mixed linear model that includes the most significant markers as covariates. It uses fixed and random effect models iteratively to help reduce potential confounding between the markers and kinship. This model has previously been successfully utilized in soybean genome-wide association analyses to identify genomic regions controlling canopy wilting [32,46], carbon and oxygen isotope ratios [12], and resistance to Sclerotinia sclerotiorum [47].
Manhattan plots were visualized with the 'qqman' [48] and 'CMplot' R packages using the p-values generated from the FarmCPU output. The significance threshold (p < 0.0001; −log 10 (P) > 4) was used to determine if SNPs were significantly associated with the traits of interest. This threshold is less stringent than a Bonferronicorrected threshold, but is more stringent than many other soybean GWAS studies using 50 K SNP genotyping data [12,46,49,50]. It is also near the point at which the p-values deviated from the linear expected pvalues in the quantile-quantile (QQ) plots (Additional file 2). Days to flowering (DTF) was recorded in both field environments as the number of days from sowing until 50% of the plants in a plot reached the first bloom (R1) growth stage. The carbon and nitrogen related traits evaluated had relatively strong correlations (data not shown) with DTF in both environments, so DTF was used as a fixed effect covariate, along with the first four genetic principal coordinates, in the GWAS to account for this correlation and population structure, respectively.
Haploview version 4.2 software [51] was used to calculate pairwise estimates of D′ and r 2 and estimate linkage disequilibrium (LD) blocks. Using D′ > 0.8 to extend the spine, LD blocks were identified by chromosome with the Solid Spine of LD option. These LD blocks were used to determine if significant (p < 0.0001; −log 10 (P) > 4) SNPs that are physically close (less than 1 Mb) were at the same locus (genomic region) controlling the trait of interest. Significant SNPs not part of the same LD block were deemed different loci controlling the trait.