Mapping additive and epistatic QTLs for forage quality and yield in soybean [Glycine max (L.) Merri.] in two environments

Abstract Soybean plants have high protein content and can be used as a supplementary source of high-protein feed. To map quantitative trait loci (QTL) for the content of crude protein (CP), neutral detergent fibre (NDF), acid detergent fibre (ADF) and dry weight of plant (DWP) in R2 stage of soybean, two recombinant inbred lines, RIL3613 and RIL6013, containing 134 and 156 RILs, derived from the cross of Dongnong L13 × Heihe 36 and Dongnong L13 × Henong 60, were planted for two consecutive years. Based on a simple sequence repeat (SSR) linkage map, QTLs of CP, NDF, ADF and DWP were mapped by interval mapping (IM) and inclusive composite interval mapping method (ICIM) using additive effect, epistatic effect and environmental interaction model. The variance components of genotype, environment and genotype × environment (G × E) interaction for quality and yield traits in the two RIL populations were significant under multiple environmental conditions. Eighteen additive effect QTLs on 10 of 20 soybean chromosomes explained 7.02%–15.67% and 2.13%–11.42% of the phenotypic variation in RIL3613 and RIL6013, respectively. Three epistatic QTL pairs related to CP and six ones for DWP were identified. Eight additive effect QTLs for CP, eight ones for NDF, three ones for ADF, and another three ones for DWP were identified by genotype × environment interaction analysis. One epistatic QTL for CP, 11 epistatic QTL pairs for ADF and 38 ones for NDF were identified. These results can provide better understanding of the genetic basis of soybean feed quality and yield.


Introduction
Soybean [Glycine max (L.)] is not only an important source of edible protein and oil, but also has rich nutritional value. Like most legumes, soybean feeds tend to be high in protein and low in fibre content, and are considered to have similar nutritional value compared to alfalfa [1]. Soybean as a feed crop has a long history of cultivation. Soybean was introduced to North America in 1765 and was grown primarily as a forage crop for the next 155 years. Starting from the 1920s, soybean has been grown as a seed crop mainly for seed protein and oil. With the growing demand for soybean oil and meal, in 1941, the acreage sown to seed soybeans in the United States exceeded that sown to forage soybeans for the first time [2]. The traditional soybean breeding target is mainly the seed yield and quality, which is reflected in the seed protein and oil content and other quantitative traits [3]. However, the progress made in the last two decades in improving the feed quality and biomass of soybeans indicates a growing interest in breeding forage soybeans [1,[4][5][6].
In the study of biomass per plant, shoot dry weight is one of the acceptable indicators. This method is often used to estimate plant yield, but it is also an accurate measure of plant biomass [7]. With the increasing demand for high quality feed, it is imperative to combine the dry weight of plant (DWP) per plant with crude protein content (CP), neutral degraded fibre content (NDF) and acid degraded fibre content Soybean; crude protein content; neutral detergent fibre content; acid detergent fibre content; dry weight of plant; QTL (ADF) as a comprehensive evaluation index to evaluate the genetic variation of soybean feed quality and yield for the improvement of soybean varieties. Significant genetic variation in feed quality traits had been detected in soybean [8][9][10]. With the objective to improve the feed quality and yield related traits in soybean efficiently, several studies on mapping quantitative trait locus (QTL) for feed quality and yield had been conducted based on a variety of genetic backgrounds, environments and statistical methods. Asekova et al. [1] used a set of 188 recombinant inbred lines (RILs) derived from a cross between wild soybean (PI 483463) and cultivated soybean (Hutcheson) to locate QTL controlling feed quality traits. Six CP QTLs, six NDF QTLs and two ADF QTLs were identified. These QTLs were distributed on 8 out of 20 soybean chromosomes, and some QTLs for CF and ADF on L chromosome were identified stably in multiple environments. The results indicated that the genomic region labelled BARC-029975-06765 may carry a gene that has a significant effect on the expression of feed quality traits. Santos et al. [11] used soybean varieties Bossier (high biological nitrogen fixation capacity) and Embrapa 20 (medium biological nitrogen fixation capacity) to obtain a RIL, and identified 2 DWP QTLs from it, explaining 15.4% of the total variation of DWP [12]. Zhang et al. [13] found eight QTLs controlling DWP using 152 RILs derived from a cross between Nannong94-156 and Bogao. A coincident QTL for SDW flanked by S506-Satt534 on linkage group B2-1 was detected across years at low phosphorus condition, which explained 41.4% of the phenotypic variation in combination. Up to date, 7 CP QTLs, 8 NDF QTLs, 2 ADF QTLs and 47 DWP QTLs have been integrated on the public genetic linkage map (http://www.soybase.org/).
The genetic effects include additive effects, epistasis and interactions of QTLs with the environment. Epistasis effects refer to additive × additive interaction components for RIL population, which may play an important role in the improvement of quantitative traits [14]. Under the assumption that there is no epistasis, QTL mapping genetic models will lead to biased estimation of QTL parameters, therefore, models containing the epistasis are proposed [12]. In quantitative genetics, additive × environment (AE) and epistasis × environment (AAE) are two major genetic components that contribute significantly to phenotypic variation of complex traits [15,16]. With advances in molecular techniques and statistical methods in recent decades, many studies on soybean have identified many pairs of epistatic QTLs [17][18][19][20][21][22][23]. However, the studies on mapping QTL for forage quality and yield traits have mainly focused on the additive effect and neglected AA and AAE [1,11,13].
With the objective to provides theoretical basis and to search technical support for forage soybean breeding, mapping QTLs for CP, NDF, ADF and DWP was conducted based on two soybean RIL populations derived from the crosses of Dongnong L13 × Henong 60 and Dongnong L13 × Heihe 36 in this study. . Both parents and recombinant inbred lines were planted in the field. The field experiment was carried out in a randomized block design and repeated 3 times. Each plot is 5 m long, the row spacing is 0.7 m, and the plant spacing is 0.06 m. The field management was similar as the general field cultivation. Three soybean plants from each plot of each RIL were collected at the R2 growth stage (full-bloom stage) [24] for forage quality and dry weight measurements from individual plants.

Measurement of forage traits
Each harvested plant is placed in a separate paper bag and placed in a 60 °C forced-air oven (YS-120203N, Vision Scientific Co.,106 Ltd) until constant weight. Dried samples were ground with a cyclone mill (Hood Mixer, AC 220 Super grinder) and passed through a 1-mm sieve for traits of forage quality and yield determination. CP was determined by Kjeldahl nitrogen determination method (Kjeltec 8400, FOSS), NDF and ADF were determined by Van Soest method (Fibertec 8000, FOSS).

Statistical analysis
For forage quality and yield traits between parents and inbred lines Microsoft Excel 2016 were used to perform statistical analysis on the mean value, SD, CV and other parameters of parental samples and RILs samples, and Graph Prism 8 was used for plotting. The SAS program (Version 9.3;SAS Institute, 2013) was applied to conduct analysis on of variance and estimate generalized heritability (h 2 ). The model of analysis of variance (ANOVA) included environment, genotype, genotype × environment interaction effect, and the following formula is used to estimate heritability.
For single environments: For the multi-environment average values: Where h 2 is broad-sense heritabilit,

QTL mapping
The SSR linkage map has been constructed in previous studies [25]. The total lengths of the SSR linkage maps for RIL3613 and RIL6013 were 2849.54 cM and 1886.8 cM, respectively, and the mean interval lengths were 21.92 cM and 16.13 cM, respectively. QTL mapping was performed for the average of quality and yield traits of each line in each environment by interval mapping (IM) and inclusive composite interval mapping (ICIM) [26,27]. Using the software QTL IciMapping v4.2, firstly, the ICIM-ADD and ICIM-EPI algorithms in BIP were used to map the additive-effect and epistatic-effect QTLs in single environment and the average of two environments, and then MET model was used to analyze the additive and epistatic effect QTLs in two environments jointly. Recombination frequencies between linked loci were transformed into distances (cM) using Kosambi's mapping function [28]. The existence of QTL was determined by setting the minimum LOD threshold to 2.5. The naming of QTL follows the QTL nomenclature of McCouch et al. [29].

Phenotypic analysis
CP, NDF, ADF and DWP of soybean were measured in two RIL populations and three parents for two environments. Among RIL3613 and RIL6013 populations, four traits showed significant variation among individuals (Table 1), and the values of each trait were distributed normally (Figures 1 and 2). Compared with the phenotype of the parents (Table 2), transgressive inheritance existed in this four traits. Furthermore, the skewness and kurtosis values of most of the traits of RIL3613 and RIL6013 were less than 1, indicating that the four traits of RIL3613 and RIL6013 conform to the quantitative genetic model [30]. Therefore, RIL3613 and RIL6013 are suitable for QTLs localization. The variance of genotype, environment and gene-environment interaction of the four traits was significant under two environments (p < 0.0001) ( Table 2). The generalized heritabilities of soybean forage quality and yield traits under single environment was significantly higher than those under two combined environments. These results indicated that the genetic basis of feed quality and yield traits of soybean were different in various environments and specific QTL would be detected in multiple environments.
In RIL6013 population, two QTLs with multiple effects simultaneously controlled NDF and ADF, and one QTL with multiple effects simultaneously controlled CP and NDF. Meanwhile, two QTL SSR intervals (sat_096-sat_289, Sat_413-Sat_160) simultaneously controlled NDF and ADF and one QTL SSR interval (satt546-staga002) simultaneously controlled CP and NDF. The additive effect of QTLs controlling CP was positive, while the additive effect of QTLs controlling NDF and ADF were negative, which means that the use of these QTLs may simultaneously increase CP and reduce ADF and NDF.

Additive × environment interaction QTL
The phenotypic data obtained in the two environments were used for QTL mapping of CP, ADF, NDF and DWP based on the multi-environment QTL mapping model, and a total of 22 additive × environment QTLs were detected (Table 5). Eight CP QTLs were located on six linkage groups, A2, D1A, F, K, M and N. Three ADF QTLs were mapped and distributed on three linkage groups, C2, I and M. Eight NDF QTLs were located on six linkage groups, A1, A2, B1, C1, C2 and M. Three DWP QTLs were located and distributed on three linkage groups, A1, A2 and F.
One epistatic effect QTL pair for CP was identified, qCP-Gm02-2~qCP-Gm19-1. The epistatic (additive × additive) contribution rate for phenotype variat i o n ( A AC P V ) w a s g r e a t e r t h a n t h e epistatic × environment interaction contribution rate for phenotype variation (AAECPV), indicating that CP was strongly controlled by the epistatic effects.
Thirty-eight epistatic effect QTL pairs for NDF were identified, and the PVE of AA was 6.514%~60.095%. The PVE range of AAE is between 0%~0.287%. For all NDF epistatic-effect QTL pairs, the PVE for the epistatic QTLs was greater than the PVE for the AAE interaction; that is, it was deeply impacted by the epistatic effects.
Among them, 15 pairwise interaction sites showed positive epistatic effects, while the remaining 23 pairwise interaction sites showed negative epistatic effects.

Associated populations can improve the efficiency of QTL mapping
In previous studies, a recombinant inbred line derived from two parents was often used for QTL mapping. The number of polymorphism markers between parents may be limited, leading to low marker density, fewer detected intervals and low efficiency of QTL detection. To overcome this limitation, plant breeders had adopted a multi-population improvement strategy, which has been used in rice, arabidopsis, maize and soybean [31][32][33][34][35][36]. In this study, we used two RIL populations with the same female parent (Dongnong L13) to locate QTLs for forage quality and yield traits (Tables  3 and 5). Ten additive effect QTLs were identified in the RIL3613 population, whereas 18 additive effect QTLs were identified in the RIL6013 population, and the intervals of these QTLs were not duplicated, indicating that the genetic mechanisms of CP, ADF, NDF and DWP in the two populations were different, and the detection efficiency could be improved by using two populations.

Multi-environment experiment could improve the detection efficiency of QTL for forage quality and yield
The quality of soybean forage is easily affected by environmental conditions [1]. This study found that in addition to forage quality, DWP was also easily affected by environmental conditions (Table 1). In order to overcome the influence of the environment on the results of QTL mapping, identify the stable expression of QTLs in different environments, and accurately estimate the location and function of QTLs, it is necessary to conduct analysis in multiple environments over many years. There are two methods to analyze the results of multi-environment experiments. One is to analyze the results of a single environment alone. The other is to combine the data collected from multiple environments [25,[37][38][39][40]. Compared with the single environmental analysis method, the statistical method combining multiple environments can accurately estimate the error, improve the detection ability of QTLs, and estimate the location and effect of QTLs more accurately [41]. In this study, both single environment (BIP function) and multiple environment (MET function) methods were used to map QTLs for three forage quality traits and one yield trait. Fourteen additive effect QTLs and nine epistatic effect QTL pairs physical regions 2 loD loD(aa) loD(aabye) pVe pVe(aa) pVe(aabye) add1 add2 addbyadd a1bye_01 a1bye_02 a2bye_01 a2bye_02 aabye_01 aabye_02 Quantitative traits are controlled by both genetic factors and environmental influences. The sensitivity of a certain QTL effect to the environment can be evaluated by comparing the phenotypic contribution rate of QTL main effect and the phenotypic contribution rate of genotype × environmental interaction effect. In this study, the additive phenotypic contribution rates of 18 additive effect QTLs were higher than that of additive × environmental interaction effects, and the phenotypic contribution rates of 50 epistatic effect QTL pairs were higher than that of epistatic × environmental interaction effects (Tables 5  and 6). These QTLs showed stability in multiple environments, indicating that their effects could be expressed stably in multiple environments, which is of great significance for breeding with multiple environment stability. On the contrary, the phenotypic contribution rate of additive × environment interaction was higher than that of additive effect in four QTLs, and there were no epistatic QTL pairs in which the phenotypic contribution rate of epistatic × environment interaction was higher than that of epistatic effect in this study. In addition, 11 QTLs were even only identified in a single environment ( Table 3). The QTLs that are only suitable for a single environment, are environment-specific QTLs, and have practical significance for breeding varieties in specific environments.

The epistasis effect QTL is an important genetic basis for forage quality and yield traits
The epistasis of QTL × QTL and its interaction with the environment play an important role in the genetic research of quantitative traits, and it is an important component of quantitative trait inheritance [42]. It is of great significance to study the epistasis effect of complex quantitative traits for genetic breeding. Asekova et al. [1] only studied the additive effect of QTLs, but did not study the epistatic effect of QTLs. Laurie et al. [43] divided epistatic effect QTL pairs into three types: significantly additive QTL × significantly additive QTL, significantly additive QTL × non-significantly additive QTL, and non-significantly additive QTL × non-significantly additive QTL according to the two QTLs interacting with each other. In this study, we mapped a total of 59 epistatic effect QTL pairs by the epistatic effect model, including 4 CP QTL pairs, 11 ADF QTL pairs, 38 NDF QTL pairs and 6 DWP QTL pairs. A total of 64 QTLs were identified. Among them, 6 QTLs were identified as significant additive QTLs by the additive effect model, whereas the other 58 QTLs were non-significant additive QTLs. Therefore, through epistatic effect analysis, QTL locus that could not be identified by additive effect model could be found, providing genetic basis for further exploration of forage quality and yield traits. Therefore, QTL analysis of additive effect and epistatic effect should be carried out simultaneously when analyzing the genetic basis of forage quality and yield traits, which can improve the efficiency of molecular assisted selection and genetic analysis.

Breeding value of co-localization with QTLs
Some QTLs are localized to the same genomic region, and this pleiotropy is common in most QTLs [44]. Most QTLs are aggregated, which is probably the result of close association of multiple genes. It is an important problem in genetic breeding to analyze whether a single gene affects multiple traits simultaneously. If pleiotropy is caused by the tight linkage of multiple genes, then the beneficial genes can be exploited by breaking the linkage. However, if a single effect controls multiple traits, it may not be possible to break the link between them because it carries some negative effects at the same time, thus limiting its effective use. In this study, some QTLs controlling multiple traits were identified. For example, in the Gm01 linkage group, qADF-Gm01-1, qNDF-Gm01-1 and qDWP-Gm01-1 share a genomic inter val (satt482-Sat_160). In the Gm02 linkage group, the gene interval staga002-sat_289 affects CP, ADF, NDF and DWP simultaneously. In the Gm11 linkage group, the gene interval satt197-sat272 affects CP, ADF and DWP simultaneously. In the Gm18 linkage group, qADF-Gm18-1 and qNDF-Gm18-2 overlap in position and coexist in the intergenic region (Satt610-Satt570). The QTLs controlling ADF and NDF in satt482-Sat160, staga002-sat_289 and satt197-sat272 are all negative additive effects, and CP are all positive additive effects. The alleles of the female parent Dongnong L13 promoted simultaneous reduction of ADF and NDF and increased CP. Fibre is an important forage component that affects animal digestibility and energy utilization. An increase in fibre content will reduce the digestibility of the forage, ultimately leading to a decrease in energy content. Fibre content is usually determined by measuring NDF and ADF [45,46]. With the increase of ADF, the digestibility of forage decreased [47]. Cell wall digestibility of grasses and legumes is often limited by the content and composition of lignin components. Since digestibility is negatively correlated with cell wall lignin and fibre content, reducing fibre content can help develop feed crops with higher digestibility [48]. Since the protein content corresponds to the digestible portion of dry matter in hay [49], increasing the crude protein content should improve the digestibility of hay. Therefore, using the gene regions of QTLs above can better achieve the purpose of improving the nutritional quality of forage soybean.

Comparison with previous studies
Comparing the QTLs identified in this study with the results of previous studies, it was found that the genetic regions of qDWP-Gm01-1 and qDWP-Gm02-4 overlap with the genetic interval of QTLs previously discovered by Zhang et al. [50]. The genetic interval of qDWP-Gm11-1 overlaps with the genetic interval of QTLs mapped by Zhang et al. [50] and Liang et al. [20]. The genetic interval of qCP-Gm19-3 overlaps with the QTL previously discovered by Asekova et al. [1]. The comparison with previous studies shows the authenticity of these QTLs. The other QTLs in this study are all reported for the first time and need further study in the future.

Conclusions
Eight CPs, eight ADFs, nine NDF and three DWP additive effect QTLs were identified in two soybean populations. Four epistatic effect QTL pairs for CP, 12 for ADF, 38 for NDF and 8 for DWP were also identified. Four QTLs were consistent with previous studies, and the rest were newly discovered QTLs. These results proved theoretical and practical significance for future maker assisted selection.

Disclosure statement
No potential conflict of interest was reported by the authors.

Data availability
The data that support the findings reported in this study are available from the corresponding author upon reasonable request.