Identification of Additive–Epistatic QTLs Conferring Seed Traits in Soybean Using Recombinant Inbred Lines

Abstract

Seed weight and shape are important agronomic traits that affect soybean quality and yield. In the present study, we used image analysis software to evaluate 100-seed weight and seed shape traits (length, width, perimeter, projection area, length/width, and weight/projection area) of 155 novel recombinant inbred soybean lines (NJRISX) generated by crossing “Su88-M21” and “XYXHD”. We examined quantitative trait loci (QTLs) associated with the six traits (except seed weight per projection area), and identified 42 additive QTLs (5–8 QTLs per trait) accounting for 24.9–37.5% of the phenotypic variation (PV). Meanwhile, 2–4 epistatic QTL pairs per trait out of a total of 18 accounted for 2.5–7.2% of the PV; and unmapped minor QTLs accounted for the remaining 35.0–56.7% of the PV. A total of 28 additive and 11 epistatic QTL pairs were concentrated in nine joint QTL segments (JQSs), indicating that QTLs associated with seed weight and shape are closely related and interacted. An interaction was also detected between additive and epistatic QTL pairs and environment, which made significant contributions of 1.4–9.5% and 0.4–0.8% to the PV, respectively. We annotated 18 candidate genes in the nine JQSs, which were important for interpreting the close relationships among the six traits. These findings indicate that examining the interactions between closely related traits rather than only analyzing individual trait provides more useful insight into the genetic system of the interrelated traits for which there has been limited QTL information.

Introduction

Soybean (Glycine max [L.] Merr.) is widely cultivated and consumed, accounting for 70% of protein meal and 28% of vegetable oil sources globally in 2018 (SoyStat, 2019). As yield is the most important target trait for improvement of this crop, its closely related trait – seed size including seed weight and shape (volume) – have been widely investigated (Fang et al., 2017; Yu et al., 2017). Seed size, which is measured as 100-seed weight (100-SW), is a fitness trait that is critical for adaptation to a particular environment (Tao et al., 2017). Seed shape (volume) traits have also been the focus of improvement by farmers and soybean breeders (Garg et al., 2017); round seed varieties are more suitable for sowing and mechanical handling, and large-seed soybeans are more attractive for direct consumption (e.g., as edamame) while small-seed varieties are consumed as sprouts. Thus, 100-SW and seed shape not only affect production and processing, but also influence growers’ preferences in the propagation of a cultivar (Maughan et al., 1996).

While 100-SW of soybean can be easily measured, seed shape has not been well studied because the traits are ill-defined and their measurements are tedious and inaccurate. Imaging technology has recently been applied to the study of crop phenotypes (Price et al., 2011). For instance, a computer image-based software has been developed that can accurately measure soybean seed morphology traits including seed length (SL), seed width (SW), seed perimeter (SP), and seed projection area (SA) (Ding et al., 2019), which can be used to calculate seed length-to-width ratio (SLW) and seed weight per projection area (SWA). This procedure is simple, accurate, and has high throughput compared to manual measurements using Vernier calipers.

Natural selection of larger seeds in soybean has resulted in an accumulation of minor QTLs (Liu et al., 2007), and QTL mapping has provided insight into these evolutionary changes (Salas et al., 2006; Xu et al., 2011; Hina et al., 2020). Early, composite interval mapping (CIM) using Windows QTL Cartographer (Wang et al., 2006) or inclusive composite interval mapping (ICIM) using QTL IciMapping (Li et al., 2007) were applied to QTL mapping. However, the CIM and ICIM methods can only detect additive QTLs and additive–epistatic QTLs, while providing no QTL × environment information. In contrast, the mixed-model-based composite interval mapping (MCIM) of QTL Network can detect additive, epistatic, and QTL × environment interactions, thus providing more detailed QTL information (Yang et al., 2008).

Linkage mapping has been widely used in detecting quantitative trait loci (QTL) of soybean 100-SW (Stombaugh et al., 2004; Kato et al., 2014; Yang et al., 2019). To date, 304 QTLs have been identified for 100-SW in soybean; most are minor QTLs and the candidate genes have yet to be validated (Karikari et al., 2019). There are 29 QTLs for SL mapped to 13 chromosomes, 25 for SW on 13 chromosomes, and 18 for SLW on 11 chromosomes in SoyBase¹.

Besides, genome-wide association studies (GWAS) has been also widely used to analyze soybean 100-SW and seed traits for natural populations (Zhang et al., 2016b; Li et al., 2019; Zhao et al., 2019; Qi et al., 2020). For example, a new seed size locus SW9-1 was found through the GWAS study for 100-SW, SL and SW (Li et al., 2019). GWAS is powerful in detecting additive QTLs both for natural population and bi-parental population, but if epistatic QTLs are involved, new GWAS procedure is not available yet. Therefore, in the present study, the linkage mapping procedure MCIM of QTL Network will be considered for a recombinant inbred line (RIL) population.

Given the lower marker density in early genetic linkage maps, chromosome regions potentially harboring 100-SW and seed shape QTLs were too broad and often overlooked, leading to imprecise mapping that was not useful for identifying candidate genes for marker-assisted breeding. Recent advances in sequencing technology have led to the discovery of new markers and the generation of high-density molecular genetic linkage maps to detect 100-SW QTLs for annotated candidate genes in soybean (Karikari et al., 2019). Thus, a high-density molecular genetic map is a basic requirement for QTL fine mapping and candidate gene discovery (Zhang et al., 2016a; Wei et al., 2019).

100-seed weight and seed shape are heritable traits conferred by both major and minor genes (Zhao et al., 2019). For example, the GmGA20OX gene affecting SL, SW, and 100-SW encodes an enzyme involved in gibberellin synthesis (Lu et al., 2017), while the Arabidopsis homolog of GmCYP78A10, regulates SL, SW, seed thickness, and seed weight in soybean (Wang et al., 2015). The SoyWRKY15a gene associated with soybean seed volume and weight was identified through a combination of RNA sequencing and QTL mapping (Gu et al., 2017). However, the up-/downstream relationships of these genes and the mechanisms by which they regulate seed traits during plant development are unclear. The aim of the present study was to identify a QTL system for seed weight and morphology including additive, epistatic, and QTL × environment interactions of 100-SW, SL, SW, SP, and SA as well as SLW and SWA using an enhanced high-density genetic linkage map for a population of the RIL NJRISX. From the QTL system, candidate genes were annotated and validated using published RNA expression datasets. As seed weight and shape traits are interrelated, we speculated that their genetic constitutions may be somewhat overlapping, therefore, the genetic relationships among these traits were also examined.

Materials and Methods

Plant Materials and Field Experiments

A population of the RIL NJRISX was established from 155 F₂-derived homozygous lines obtained by crossing “Su 88-M21” and “XYXHD”. The large-sized round seed of “Su 88-M21”, small-sized oval seed of “XYXHD” and the RIL population (NJRISX) derived from the cross were provided by the National Soybean Improvement Center of Nanjing Agricultural University. RILs along with the parents were tested in three environments: Jiangpu Experimental Station of Nanjing Agricultural University, Nanjing, Jiangsu Province (latitude 33°03′ N; longitude 118°63′ E) in June 2017 (17JP); Wanjiang Station of Nanjing Agricultural University, Dangtu County, Anhui Province (latitude 32°87′ N; longitude 117°56′ E) in June 2017 (17DT); and the same location in June 2019 (19DT). Each line was planted in a single row plot (length × width, 1 × 0.5 m) in a randomized complete block design with three replicates. Uniform agronomic practices were used in the experiments.

Measurement of 100-Seed Weight and Seed Shape Traits

After harvesting, seeds were dried to a uniform moisture. Diseased, insect-infested, and physically damaged seeds were removed. Seed weight was measured using an electronic balance with 0.001-g accuracy. SL, SW, SP, and SA were measured from images using a high-speed camera (Model eloam-S1500A2; Shenzhen E-Loam Technology Co., Shenzhen, China). About 120 seeds per plot were collected and distributed evenly in the middle of the light-emitting diode backlight board of the equipment to obtain clear and complete images of the seeds. The board was calibrated and adjusted to an appropriate brightness according to the indoor light level. The equipment parameters were as follows: brightness, 64 cd/m²; contrast, 15; hue, 0; saturation, 43; clarity, 100; gamma, 100; white, balance 4600 auto, backlight contrast 0, power line frequency (anti-flicker) 50 Hz, focus 65, exposure −6; and scan size, 640 × 480 mm. The acquired images were processed using previously developed software (Ding et al., 2019). The image background was removed based on the “Otsu” threshold method to obtain the binary image of soybean seeds, and the adhesive soybean seeds in the binary map were segmented and counted based on the watershed transformation method. The sum of white pixels in each connected domain and the correction formula based on the freeman chain code algorithm is used to calculate seed area and perimeter, respectively. The second-order statistical moment is used to obtain the main axis direction of the soybean seeds and then the SL and SW was calculated by the extreme difference of the boundary point. Besides, this software can process hundreds of photos at a time and export the data to EXCEL. Data on the exact number of seeds on the backlight board, SL, SW, SP, and SA were directly obtained using the computer software and 100-SW was converted from seed number and seed weight values; SLW = SL/SW and SWA = 100-SW/SA were calculated. All seven seed weight and shape traits were measured in three replicates under three environmental conditions.

Before measuring the experimental seeds, in order to confirm that images obtained of a seed lot from a given plot were reproducible, two seed samples of the 138 lines were obtained from a single replicate and SL, SW, SP, and SA were measured under the same conditions using the same instrument. The correlation coefficients between the two measured values of SL, SW, SP, and SA were 0.95, 0.92, 0.93, and 0.94, respectively, indicating good consistency between the measurements from the same seed lot (Supplementary Table S1) and validating the utility of the imaging procedure used in this study.

Statistical Analysis of Phenotypic Data

Data were analyzed using Excel 2016 software (Microsoft, Redmond, WA, United States). Analysis of variance (ANOVA) and correlation analysis were performed using SAS v9.4 software (SAS Institute, Cary, NC, United States). Heritability (h²) was calculated as

where σ²_g, σ²_gy, and σ²_e are genotype, genotype × environment interaction, and error variance estimated from the expected mean squares in ANOVA, respectively; n is the number of environments; and r is the number of replicates. The genotypic coefficient of variation (GCV) was calculated as GCV = σ_g/μ, where μ is the mean value of the RIL population.

Specific-Locus Amplified Fragment Sequencing (SLAF-seq) and Genetic Linkage Map Construction

In 2015, two parents and 155 progeny soybean lines were planted at Jiangpu Experiment Station. Genomic DNA was extracted from young leaves and SLAF-seq (Sun et al., 2013; Zhang et al., 2016a) was performed by the Biomarker Technologies Corporation (Beijing, China). Briefly, the soybean genome sequence (G. max, Wm82.a1. v1) was used as a reference to predict digestion sites. RsaI and HaeIII were used to digest the genomic DNA and the obtained fragments (SLAF tags of 364–414 bp) were processed for target selection. After passing library quality inspection, sequencing was performed with the HiSeq 2500 system (Illumina, San Diego, CA, United States). Rice (Oryza sativa)² was processed in the same manner and served as the control for library construction and sequencing to determine whether the enzymes used in this experiment were activated or there were other quality problems in library construction. A total of 281.93 M reads were generated for the two parents and 155 recombinant inbred lines. After discarding low-quality reads, 247,821 SLAF tags were screened; of these, 207,180 and 212,958 were identified from female parent (Su 88-M21) and male parent (XYXHD), respectively, with sequencing depths of 41.56 and 44.04 fold, respectively. There were 150,515 SLAFs in the RILs with 7.75-fold coverage on average, corresponding to 1,473,406 reads. The reads of each identified sample were analyzed by clustering for SLAF tag screening. After removing low-quality reads, 52,988 tags that matched the parental separation pattern (Sun et al., 2013) were identified as polymorphic in the whole RIL population; these were screened according to SLAF tag filtering rules (Sun et al., 2013), yielding 9625 markers for linkage analysis.

Markers with high collinearity were removed for map construction, leaving 5351 SLAF markers. To obtain high-quality molecular tags, modified logarithm of odds scores between the tags were calculated and used for linkage grouping. HighMap software (Liu et al., 2014) was used to construct a genetic map for each linkage group. The software uses an efficient maximum likelihood estimation method to correct label classification based on the layout results, and after multiple cycles of layout correction–layout, a high-quality map is obtained. Map quality was evaluated in terms of co-linearity, genetic relationships, and monomer sources. The linkage map (Supplementary Figure S1) was drawn using the R-based software LinkageMapView (Ouellette et al., 2018).

Mapping QTLs Conferring 100-Seed Weight and Seed Shape Traits

The MCIM function of QTL Network v2.0 was used to detect additive QTLs, additive × additive epistatic QTL pairs, and additive QTL × environment and epistatic QTL pair × environment interactions. The critical F value of MCIM was calculated with 1,000 permutation tests. The QTL effects were estimated by using the Monte Carlo Markov Chain method with 20,000 Gibbs sampler iterations and candidate interval selection; putative QTL detection and QTL effects were calculated with an experiment-wise type I error under a = 0.05 (Yang et al., 2007; Xing et al., 2012). The CIM method of Windows QTL Cartographer v2.5 software was used to perform CIM scanning of chromosomes to identify additive QTLs of soybean seed traits in different environments and verify the QTL mapping results of QTL Network v2.0 software. For CIM, the LOD significance threshold determined empirically using 1,000 permutation tests. Neighboring QTLs of different seed traits within the same support interval were grouped into a joint QTL segment (JQS). The 1-logarithm of the odds support (confidence) intervals were calculated using the QTL network procedure (Yang et al., 2007), which is defined by points on the genetic map at which the likelihood ratio has fallen from the maximum by a factor of 10 (Lander and Botstein, 1989).

Positive and negative additive effects indicated that the alleles were from “Su 88-M21” and “XYXHD” for all seed traits, respectively. Genetic contribution of the collective unmapped minor QTLs is equal to total genetic contribution minus variation explained by all detected additive and epistatic QTLs (Korir et al., 2011). Random error variation is equal to total phenotypic variation minus total genetic contribution, and variation explained by all detected additive QTL × environment and epistatic QTL× environment interactions (Xing et al., 2012).

Annotation of Candidate Genes Conferring 100-Seed Weight and Seed Shape Traits

Candidate gene annotation was carried out for JQSs based on physical locations in SoyBase³. Gene ontology (GO) annotation v1.1 was downloaded from SoyBase⁴ and GO classification was based on clusterProfiler package (Yu et al., 2012) in R software (p value < 0.01, q value < 0.05) to identify terms related to seed weight and shape traits in soybean. To determine whether these genes are expressed in seeds, gene expression data of 14 soybean tissues at seven seed development stages (seed_10DAF, seed_14DAF, seed_21DAF, seed_25DAF, seed_28DAF, seed_35DAF, and seed_42DAF, where DAF stands for days after flowering) and 7 other tissues (young_leaf, flower, one.cm.pod, pod.shell.10DAF, pod.shell.14DAF, root, and nodule) (Severin et al., 2010) were downloaded from SoyBase³. Expression data from cultivated soybean were used as an approximation for analyzing the parents in the present study. The candidate genes were classified according to Protein Class using Protein ANalysis THrough Evolutionary Relationships (PANTHER)⁵ and annotated based on the National Center for Biotechnology Information (NCBI)⁶ and UniProt Protein⁷ databases to determine gene function.

Results

Phenotypic Variations (PVs) in 100-Seed Weight and Seed Shape Traits

The seed of the maternal parent “Su 88-M21” is round and large, while that of the paternal parent “XYXHD” is flat and small (Figure 1H); as these seeds differ significantly in terms of 100-SW, SL, SW, SP, SA, SLW, and SWA, it was possible to establish a RIL population with potential for genetic variation in these seven seed weight and shape traits (Table 1). Accordingly, the frequency distributions all showed a large variation (Figure 1 and Table 1), while transgressive segregation was observed for SL, SP, and SA in both directions (Figures 1B,D,E). In the joint ANOVA for multiple environments, there existed significant differences among Lines and Line × Env. interactions for all the traits except no significant differences among lines in SWA (Table 2). The phenotypic data from multiple environments were used to estimate heritability. The heritability of the seven traits ranged from 60.3 to 88.2% (Table 1). These results showed that further QTL constitution analysis for the traits except SWA would be meaningful.

FIGURE 1

Of the seven seed traits, 100-SW, SL, SW, SP, and SA are of the first order and can be directly measured, while SLW and SWA are second-order traits that are calculated from first-order traits. To clarify their relationships, a correlation analysis was performed for these traits. The correlation coefficients of the five first-order traits ranged from 0.78 to 1.00, showing that they are closely related (Table 3). Notably, the correlation coefficient between SP and SA was approximately 1.00; additionally, the high correlation between 100-SW and the other four first-order traits (0.88–0.94) implied that they have a common genetic basis. For the two second-order traits, SLW – which is related to seed shape – was not correlated with 100-SW or SW, while SWA – related to seed volume weight – was not correlated with SLW, their other correlations with the other first rank traits were not high, therefore, are not attractive traits. In this case, the second-order trait SWA was neglected and excluded in the QTL mapping analysis.

Genetic Linkage Map and Genetics of 100-Seed Weight and Seed Shape Traits

A molecular genetic linkage map with 5351 SLAF-seq markers was constructed that spanned 3046.52 cM with an average intermarker interval of 0.57 cM. Chr9 harbored the most markers at 500 and Chr11 had the fewest at 80; the latter spanned the shortest distance at 106.38 cM, whereas markers on Chr9 covered the largest distance at 199.24 cM (Supplementary Figure S1 and Supplementary Table S2). By comparing the physical positions of markers on each chromosome, a consistent relationship between physical and genetic distances was observed on all chromosomes except Chr11 and Chr17, indicating that the genetic map was of good quality (Supplementary Figure S2).

A total of 42 additive and 18 epistatic QTL pairs were detected for the six seed traits, accounting for 24.9–37.5% and 2.5–7.2% of the PV, respectively; thus, a large part of the genetic variance (heritability minus total QTL contribution, 35.0–56.7%) was not explained by these QTLs, and were instead attributed to a collection of undetected minor QTLs (Table 4). In addition to these QTLs, the PV was explained by additive QTL × environment interaction (1.4–9.5%), epistatic QTL × environment interaction (0.4–0.8%), and random error (2.4–14.0%). Thus, the largest portion of the genetic variation was explained by additive QTLs, with epistatic QTL pairs accounting for only a small part of the genetic variation in the six seed traits. Identification of the unmapped minor QTLs, which collectively accounted for a relatively large portion of the genetic variation, depends on improvements in the precision and sensitivity of mapping procedures. In addition, there were both additive and epistatic QTL × environment interactions, but these explained only a small part of the PV. Overall, the genetic and gene × environment components of the six traits were similar.

Additive and Epistatic QTLs Conferring 100-Seed Weight and Seed Shape Traits

Table 5 shows information on each of the additive QTLs of the six seed traits. A total of 42 QTLs distributed on 13 chromosomes were identified by MCIM of the QTL network, of which 25 were also identified by CIM using Windows QTL Cartographer. For 100-SW, eight QTLs were identified on six chromosomes, each accounting for 1.8–8.2% of the PV (Figure 2 and Table 5). Two QTLs, q100SW-6-1 and q100SW-19-1 interacted with the environment and contributed 3.9 and 2.5% PV, respectively. Three QTLs have been previously reported in the literature (Supplementary Table S3). Seven QTLs on seven chromosomes were identified for SL, accounting for 3.3–6.8% of the PV; two have been previously reported. For SW, five QTLs on five chromosomes were identified, accounting for 2.3–8.4% of the PV; one QTL has been previously reported. For SP, eight QTLs on eight chromosomes were identified, accounting for 2.1–7.8% of the PV. Eight QTLs on eight chromosomes were identified for SA, accounting for 2.2–8.4% of the PV. For SLW, six QTLs on six chromosomes were identified, each accounting for 3.6–10.1% of the PV; two have been reported in the literature. Some of the identified additive QTLs interacted with each other to form significant epistatic QTL pairs involving all traits. There were three epistatic QTL pairs for 100-SW, two for SL, two for SW, four for SP, four for SA and three for SLW, with contributions to PV ranging from 0.6 to 3.0% for a single epistatic QTL pair (Figure 2 and Table 6) and 2.5–7.2% for a single trait. These low rates indicated that there was epistasis among genes governing seed weight and shape traits, although epistatic QTL pairs accounted for just a small part of the PV. Epistatic QTL pairs also interacted with the environment, with phenotypic contributions ranging from 0 to 0.7%; they also explained 0.4–0.8% of the PV for a single trait, indicating that they were relatively stable across environments.

FIGURE 2

The major QTL qSLW-2-1 was located at 52.2–54.9 cM on Chr2, accounting for 10.1% of the PV, and was detected by Cartographer under all three environmental conditions. These seed trait additive QTLs interacted with the environment and their phenotypic contribution ranged from 0 to 3.9%; this is not a large proportion, it indicated that the identified QTLs were relatively stable.

JQSs Related to 100-Seed Weight and Seed Shape Traits

A QTL of a trait may be located in the same chromosomal region as QTL(s) of other trait(s). These QTLs may be either different loci or the same locus (referred to as JQSs) due to random shifting. Of the 42 identified QTLs, 28 were located in nine JQSs on nine chromosomes (Table 7).

The nine JQSs were designated as JQS-1, -4, -6, -7, -11, -12, -15, -17, and -19 (Table 7) and were composed of 2, 3, 5, 2, 4, 4, 2, 2, and 4 QTLs, respectively, on Chr 1, 4, 6, 7, 11, 12, 15, 17, and 19, respectively. JQS-6 harbored five QTLs, the most among JQSs (Figure 2). The 28 QTLs conferring different seed traits constituted nine QTL segments.

Joint QTL segment-1 was located at 73.0–74.1 cM on Chr1, covering 1.1 cM and containing two QTLs: q100SW-1-1 and qSP-1-1. JQS-4 was located at 62.5–64.9 cM on Chr4, covering 2.4 cM and containing three QTLs: q100SW-4-1, qSL-4-1 and qSA-4-1. JQS-6 was located at 40.3–42.3 cM on Chr6, covering 2.0 cM and containing five QTLs. JQS-7 was located at 102.8–105.4 cM on Chr7, covering 2.6 cM and containing two QTLs. JQS-11 was located at 22.8–26.4 cM on Chr11, covering 3.6 cM and including q100SW-11-1, qSW-11-1, qSP-11-1, and qSA-11-1. JQS-12 was located at 66.2–66.5 cM on Chr12, covering 0.3 cM and containing four QTLs. JQS-15 was located at 2.3–13.6 cM on Chr15, covering 11.3 cM and containing two QTLs. JQS-17 was located at 161.2–164.9 cM on Chr17, covering 3.7 cM and containing two QTLs. JQS-19 was located at 15.6–17.2 cM on Chr19, covering 1.6 cM and containing four QTLs (Table 7). Among the nine JQSs, JQS-1, -4, -6, -11, -12, and -19 harbored 100-SW and seed shape QTLs whereas only the latter were present in JQS-7, -15, and -17 (Table 7 and Supplementary Table S4).

Among the nine JQSs, JQS-11 and -17 were independent of the others, whereas the remaining seven JQS interacted with other JQSs (epistasis) or else harbored epistatic QTL. Of these JQSs, JQS-19 interacted with JQS-6 and JQS-12 for multiple traits. Although there were no interactions among the other JQSs, interactions were observed at the individual QTL level. Thus, a single or multiple QTLs in a JQS may interact with a single or multiple QTLs in another JQS, while some QTLs are independent and do not interact with other QTLs. Among the 28 QTLs on nine JQSs, eight were independent; 11 epistatic QTL pairs were located on JQSs; three epistatic QTL pairs partly located on JQS, e.g., qSL-6-1 located on JQS, while qSL-19-1 isolated from JQS; and the remaining four epistatic QTL pairs were isolated from JQSs (Figure 2 and Tables 6, 7).

Functional Annotation of Candidate Genes in the Nine JQSs

Annotation of candidate genes in the nine JQSs was carried out using SoyBase⁸. JQS-1, -4, -6, -7, -11, -12, -15, -17, and -19 contained 6, 8, 57, 26, 48, 37, 149, 95, and 72 genes, respectively, for a total of 498 genes. In the GO enrichment analysis of 498 genes, in the molecular function category, most genes were enriched in terms related to activity and binding factors such as transcriptional regulator activity, ribonucleoside binding, and purine nucleoside binding (Figure 3A). In the cell component category, genes were mostly enriched in plasma membrane, organelle envelope, and envelope (Figure 3B). In the biological process category, the genes were enriched in multicellular organismal process, multicellular organism development, response to chemical, and other terms related to growth, development, and reproduction (Figure 3C).

FIGURE 3

Previously published gene expression data for soybean (Severin et al., 2010) were used for gene annotation. Of the 498 genes, we selected 294 located in the nine JQSs with expression in seed tissues; 99 were included in PANTHER protein classes, including 16 that were related to seed weight and morphology. Additionally, two genes were identified by analyzing protein functions in the gene database combined with literature searches (Supplementary Table S4). Ultimately, 18 genes that were mostly related to seed weight and shape were selected as candidate genes for further analysis, including 1, 1, 1, 1, 3, 1, 3, 5, and 2 candidate genes located in JQS-1, -4, -6, -7, -11, -12, -15, -17, and -19, respectively. Six of the genes were related to ubiquitin protein ligase; four to protein class of ribosomal protein; four to protein phosphatase; one (Glyma19g37910) to basic leucine zipper transcription factor; and one (Glyma01g26950) to tubulin; and one each encoded RING-H2 finger protein ATL52 and auxin response factor 18 (Glyma04g15820 and Glyma07g06060, respectively) (Supplementary Table S4). The 18 genes have been shown to be related to seed size and shape traits in multiple plant species including soybean; in particular, Glyma06g22900, and Glyma19g37910 have high expression in different seed tissues and low expression in other tissues (Severin et al., 2010), suggesting that they are important for seed development (Figure 3D and Supplementary Table S4). Thus, genes with similar functions in seed development are distributed in nine JQSs, which could account for the close relationship between seed weight and shape traits. In addition, some JQSs had a single candidate gene that conferred multiple traits whereas others had multiple candidate genes conferring several traits, indicating that some genes are pleiotropic and that multiple candidate genes may exist within a JQS.

Discussion

Genetic Basis of 100-Seed Weight and Seed Shape Traits

The results of this study provide an outline of the genetic structure of the RIL population NJRISX. We identified 42 additive QTLs that contributed 24.9–37.5% to the PV in seed weight and shape traits (100-SW, SL, SW, SP, SA and SLW), as well as 2–4 of 18 epistatic QTL pairs per trait for all six seed traits that contributed 2.5–7.2% of the PV. The remaining PV (35.0–56.7%) was explained by unmapped minor QTLs. The 28 additive and 11 epistatic QTL pairs were located in nine JQSs, suggesting that they were closely related and interacted. In addition, their interaction with the environment contributed 1.4–9.5% and 0.4–0.8% to the PV, respectively, although not large but significant.

Additive QTL is an important genetic component of seed weight and shape traits. Many QTLs associated with soybean seed weight and size have been identified (304 for 100-SW, 29 for SL, 25 for SW, and 18 for SLW in SoyBase⁹), but few have been examined in detail. Of the 42 QTLs detected in the present study, 8 were harbored in known loci (Supplementary Table S3) and the other 34 were identified for the first time; in 15 of these QTLs, each accounted for >5% of the PV.

Unmapped minor QTLs accounted for a considerable portion of the total PV, implying that more QTLs might be discovered by examining more lines in the population or more markers by improving the efficiency of the mapping procedure. Previous studies have demonstrated that restricted two-stage multi-locus genome-wide association study using single nucleotide polymorphism linkage disequilibrium block markers showed a superior performance to CIM, MCIM, joint inclusive CIM, and MLM-GWAS in mapping QTLs associated with days to flowering in soybean (Li et al., 2017; Pan et al., 2018). However, this procedure can identity additive QTLs and gene × environment QTLs but not epistatic QTL pairs. For a comprehensive analysis of the genetic structure of seed traits, both procedures can be used on the same dataset or a new mapping procedure can be applied.

Correlation Among 100-Seed Traits and JQS Properties

The five first-order traits (100-SW, SL, SW, SP, and SA) were closely related with correlation coefficients ranging from 0.78 to 1.00, which was supported by the mapping results. Firstly, QTLs conferring different traits were located in the same segment, and could be the same QTL/gene with pleiotropic functions. Of the 42 QTLs, 28 were present in nine JQSs and 6, 3, 5, 7, and 7 QTLs were associated with 100-SW, SL, SW, SP, and SA, respectively. The number of JQSs shared between any two traits was highly consistent with their correlation coefficients (Table 3); that is, the number of shared JQSs forms the genetic basis of closely related seed traits. Secondly, along with additive QTLs, epistatic QTL pairs of traits were located on the same set of two JQSs, including JQS-6/-19, and JQS-12/-19; the parallel QTL interactions increased the correlation between the two traits. Additionally, additive effects of QTLs in the same JQS were in the same direction (positive or negative), suggesting a pleiotropic effect of the same QTL/gene or an aggregation of multiple QTLs/genes occurring in the same direction in the related traits. In other words, the allele effects contribution of the different QTLs in the same JQS are from a same parent. The number of QTLs and traits in the nine JQSs in the present study varied, indicating that each JQS has unique characteristics and functions. Genes with similar functions in seed development were separated in the JQSs while each JQS also harbored QTLs of genes with similar function, which could underlie the close association between seed traits.

The average segment length of overlapped confidence intervals in the nine JQSs was 3.2 cM (rang: 0.3–11.3 cM) – e.g., 11.3 cM for JQS-15 and 1.6 cM for JQS-19. The length depends on the number of QTLs linked at the same location (i.e., the density of markers in the segment). It is expected that by using higher-resolution genome-wide markers, more JQSs can be identified. From the present results, the JQSs can be grouped into three types based on whether individual QTLs in the JQS interact with others. (i) QTLs in the JQS all have epistatic interaction effects and interact with QTLs within the same JQS, resulting in a parallel relationship between traits. For example, every QTL in JQS-6 and -19 except those for SL (there was maybe a random shifting for qSL-19-1) had significant interactions with each other that appeared as interactions between the two JQSs, suggesting that the two segments had two interacting genes when they were actually the same QTL/gene with pleiotropic functions. JQS-12 and -19 also belong to this class. (ii) Some QTLs in JQSs (e.g., JQS-4, -7, and -15) show epistatic interactions but these do not occur in parallel between traits; or else different QTLs in same JQS interact with different JQSs, such as JQS-1. This class of JQSs contains QTLs with distinct properties (i.e., different QTLs for seed weight and seed shape traits) and hence, different genes. (iii) There is no epistatic interaction of QTLs in the JQS (e.g., JQS-11 and -17).

The above JQS results were observed in a set of mutually related traits, and would not be revealed if only a single trait was involved. When QTL mapping was performed for closely correlated traits, QTLs of different traits aggregated on neighboring or overlapping segments on the same chromosome (Liu et al., 2017), which are referred to as QTL hotspots (Zhang et al., 2019) or QTL clusters (Cai and Morishima, 2002). Because of limited marker density, QTL clusters were previously defined as regions with different QTLs located in the same or adjacent segments (Xu et al., 2011). While QTL clusters or hotspots have no defined thresholds, JQS is defined as a group of QTLs linked within confidence intervals; thus, using multiple related traits in JQSs can yield more precise and accurate results in the identification of QTLs/candidate genes. However, additional studies are needed to evaluate whether QTL grouping according to the support interval criterion of Lander and Botstein (1989) is appropriate.

Conclusion

In this study, we examined the QTL system of six seed weight and shape traits (100-SW, SL, SW, SP, SA, and SLW) in the soybean RIL population NJRISX, and identified 42 additive QTLs and 18 epistatic QTL pairs accounting for 24.9–37.5% and 2.5–7.2% of the PV, respectively, with the remaining part of the PV (35.0–56.7%) attributable to unmapped minor QTLs. Notably, 28 additive QTLs and the 11 epistatic QTL pairs were concentrated in nine JQSs, indicating that seed weight and shape QTLs are closely related and interact; moreover, additive QTLs and epistasis QTL pairs interaction with the environment made a small but significant contribution to the PV (1.4–9.5% and 0.4–0.8%, respectively). Thus, the JQS is important for interpreting the close relationships among the six seed weight and shape traits, especially the five first-order traits. Our findings provide a whole picture of the genetic structure of soybean seed traits and demonstrate that examining a group of closely related traits can be more informative than analyzing individual traits.

References

• 1

  CaiH. W.MorishimaH. (2002). QTL clusters reflect character associations in wild and cultivated rice.Theor. Appl. Genet.1041217–1228. 10.1007/s00122-001-0819-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  ChengZ. J.ZhaoX. Y.ShaoX. X.WangF.ZhouC.LiuY. G.et al (2014). Abscisic acid regulates early seed development in Arabidopsis by ABI5-mediated transcription of SHORT HYPOCOTYL UNDER BLUE1.Plant Cell261053–1068. 10.1105/tpc.113.121566

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  DingQ.XuW.LiM.WangX.LuW.GaiJ.et al (2019). Measurement method of soybean seed morphological parameters based on watershed and statistical moment.Soybean Sci.38960–967. 10.11861/j.issn.1000-9841

  • CrossRef
  • Google Scholar

• 4

  FangC.MaY. M.WuS. W.LiuZ.WangZ.YangR.et al (2017). Genome-wide association studies dissect the genetic networks underlying agronomical traits in soybean.Genome. Biol.18:161. 10.1186/s13059-017-1289-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  FunatsukiH.KawaguchiK.MatsubaS.SatoY.IshimotoM. (2005). Mapping of QTL associated with chilling tolerance during reproductive growth in soybean.Theor. Appl. Genet.111851–861. 10.1007/s00122-005-0007-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  GargR.SinghV. K.RajkumarM. S.KumarV.JainM. (2017). Global transcriptome and coexpression network analyses reveal cultivar-specific molecular signatures associated with seed development and seed size/weight determination in chickpea.Plant J.911088–1107. 10.1111/tpj.13621

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  GuY.LiW.JiangH.WangY.GaoH.LiuM.et al (2017). Differential expression of a WRKY gene between wild and cultivated soybeans correlates to seed size.J. Exp. Bot.682717–2729. 10.1093/jxb/erx147

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  HinaA.CaoY.SongS.LiS.SharminR. A.ElattarM. A.et al (2020). High-resolution mapping in two RIL populations refines major “QTL Hotspot” regions for seed size and shape in soybean (Glycine max L.).Int J Mol Sci.21:1040. 10.3390/ijms21031040

  • CrossRef
  • Google Scholar

• 9

  KarikariB.ChenS.XiaoY.ChangF.ZhouY.KongJ.et al (2019). Utilization of interspecific high-density genetic map of RIL population for the QTL detection and candidate gene mining for 100-seed weight in soybean.Front. Plant Sci.10:1001. 10.3389/fpls.2019.01001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  KatoS.SayamaT.FujiK.YumotoS.KonoY.HwangT. Y.et al (2014). A major and stable QTL associated with seed weight in soybean across multiple environments and genetic backgrounds.Theor. Appl. Genet.1271365–1374. 10.1007/s00122-014-2304-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  KorirP. C.QiB.WangY.ZhaoT.YuD.ChenS.et al (2011). A study on relative importance of additive, epistasis and unmapped QTL for aluminium tolerance at seedling stage in soybean.Plant Breed.130551–562. 10.1111/j.1439-0523.2011.01862.x

  • CrossRef
  • Google Scholar

• 12

  KurodaY.KagaA.TomookaN.YanoH.TakadaY.KatoS.et al (2013). QTL affecting fitness of hybrids between wild and cultivated soybeans in experimental fields.Ecol. Evol.32150–2168. 10.1002/ece3.606

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  LanderE. S.BotsteinD. (1989). Mapping mendelian factors underlying quantitative traits using RFLP linkage maps.Genetics121185–199.

  • Google Scholar

• 14

  LiH. H.YeG. Y.WangJ. K. (2007). A modified algorithm for the improvement of composite interval mapping.Genetics175361–374. 10.1534/genetics.106.066811

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  LiJ.ZhaoJ.LiY.GaoY.HuaS.NadeemM.et al (2019). Identification of a novel seed size associated locus SW9-1 in soybean.Crop J.7548–559. 10.1016/j.cj.2018.12.010

  • CrossRef
  • Google Scholar

• 16

  LiS.CaoY.HeJ.ZhaoT.GaiJ. (2017). Detecting the QTL-allele system conferring flowering date in a nested association mapping population of soybean using a novel procedure.Theor. Appl. Genet.1302297–2314. 10.1007/s00122-017-2960-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  LiuB.FujitaT.YanZ. H.SakamotoS.XuD.AbeJ.et al (2007). QTL mapping of domestication-related traits in soybean (Glycine max).Ann. Bot.1001027–1038. 10.1093/aob/mcm149

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  LiuD.MaC.HongW.HuangL.LiuM.LiuH.et al (2014). Construction and analysis of high-density linkage map using high-throughput sequencing data.PLoS One9:e98855. 10.1371/journal.pone.0098855

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  LiuN.LiM.HuX.MaQ.MuY.TanZ.et al (2017). Construction of high-density genetic map and QTL mapping of yield-related and two quality traits in soybean RILs population by RAD-sequencing.BMC Genomics18:466. 10.1186/s12864-017-3854-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  LuX.XiongQ.ChengT.LiQ. T.LiuX. L.BiY. D.et al (2017). A PP2C-1 allele underlying a quantitative trait locus enhances soybean 100-seed weight.Mol. Plant10670–684. 10.1016/j.molp.2017.03.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  MaughanP. J.MaroofM. A. S.BussG. R. (1996). Molecular-marker analysis of seed-weight: genomic locations, gene action, and evidence for orthologous evolution among three legume species.Theor. Appl. Genet.93574–579. 10.1007/BF00417950

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  OuelletteL. A.ReidR. W.BlanchardS. G.BrouwerC. R. (2018). LinkageMapView—rendering high-resolution linkage and QTL maps.Bioinformatics34306–307. 10.1093/bioinformatics/btx576

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  PanL.HeJ.ZhaoT.XingG.WangY.YuD.et al (2018). Efficient QTL detection of flowering date in a soybean RIL population using the novel restricted two-stage multi-locus GWAS procedure.Theor. Appl. Genet.1312581–2599. 10.1007/s00122-018-3174-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  PriceC. A.SymonovaO.MileykoY.HilleyT.WeitzJ. S. (2011). Leaf extraction and analysis framework graphical user interface: segmenting and analyzing the structure of leaf veins and areoles.Plant Physiol.155236–245. 10.1104/pp.110.162834

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  QiZ.SongJ.ZhangK.LiuS.TianX.WangY.et al (2020). Identification of QTNs controlling 100-seed weight in soybean using multilocus genome-wide association studies.Front. Genet.11:689. 10.3389/fgene.2020.00689

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  SalasP.Oyarzo-LlaipenJ.WangD.ChaseK.MansurL. (2006). Genetic mapping of seed shape in three populations of recombinant inbred lines of soybean (Glycine max L. Merr.).Theor. Appl. Genet.1131459–1466. 10.1007/s00122-006-0392-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  SchruffM. C.SpielmanM.TiwariS.AdamsS.FenbyN.ScottR. J. (2005). The AUXIN RESPONSE FACTOR 2 gene of Arabidopsis links auxin signaling, cell division, and the size of seeds and other organs.Development133251–261. 10.1242/dev.02194

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  SegamiS.KonoI.AndoT.YanoM.KitanoH.MiuraK.et al (2012). Small and round seed 5 gene encodes alphatubulin regulating seed cell elongation in rice.Rice5:4. 10.1186/1939-8433-5-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  SeverinA. J.WoodyJ. L.BolonY. T.JosephB.DiersB. W.FarmerA. D.et al (2010). RNA-Seq atlas of Glycine max: a guide to the soybean transcriptome.BMC Plant Biol.10:160. 10.1186/1471-2229-10-160

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  SongX. J.HuangW.ShiM.ZhuM. Z.LinH. X. (2007). A QTL for rice grain width and weight encodes a previously unknown RING-type E3 ubiquitin ligase.Nat. Genet.39623–630. 10.1038/ng2014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  SoyStat (2019). Soystats. Available online at: http://www.soystats.com/(accessed January 25, 2020).

  • Google Scholar

• 32

  StombaughS. K.OrfJ. H.JungH. G.ChaseK.LarkK. G.SomersD. A. (2004). Quantitative trait loci associated with cell wall polysaccharides in soybean seed.Crop Sci.442101–2106. 10.2135/cropsci2004.2101

  • CrossRef
  • Google Scholar

• 33

  SunX.LiuD.ZhangX.LiW.LiuH.HongW.et al (2013). SLAF-seq: an efficient method of large-scale de novo SNP discovery and genotyping using high-throughput sequencing.PLoS One8:e58700. 10.1371/journal.pone.0058700

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  TaoY.MaceE. S.TaiS.CruickshankA.CampbellB. C.ZhaoX.et al (2017). Whole-genome analysis of candidate genes associated with seed size and weight in sorghum bicolor reveals signatures of artificial selection and insights into parallel domestication in cereal crops.Front. Plant Sci.8:1237. 10.3389/fpls.2017.01237

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  TianS.WuJ.LiuY.HuangX.LiF.WangZ.et al (2017). Ribosomal protein NtRPL17 interacts with kinesin-12 family protein NtKRP and functions in the regulation of embryo/seed size and radicle growth.J. Exp. Bot.685553–5564. 10.1093/jxb/erx361

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  WangS. C.BastenC. J.ZengZ. B. (2006). Windows QTL Cartographer 2.5.Raleigh, NC: Department of Statistics, North Carolina State University.

  • Google Scholar

• 37

  WangX.LiY.ZhangH.SunG.ZhangW.QiuL. (2015). Evolution and association analysis of GmCYP78A10 gene with seed size/weight and pod number in soybean.Mol. Biol. Rep.42489–496. 10.1007/s11033-014-3792-3

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  WeiJ.ChenY.WangW. (2019). A high-density genetic linkage map and QTL mapping for sex and growth-related traits of large-scale loach (Paramisgurnus dabryanus).Front. Plant Sci.10:1023. 10.3389/fgene.2019.01023

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  XingG.ZhouB.WangY.ZhaoT.YuD.ChenS.et al (2012). Genetic components and major QTL confer resistance to bean pyralid (Lamprosema indicata Fabricius) under multiple environments in four RIL populations of soybean.Theor. Appl. Genet.125859–875. 10.1007/s00122-012-1878-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  XuY.LiH. N.LiG. J.WangX.ChengL. G.ZhangY. M. (2011). Mapping quantitative trait loci for seed size traits in soybean (Glycine max L. Merr.).Theor. Appl. Genet.122581–594. 10.1007/s00122-010-1471-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  YangH.WangW.HeQ.XiangS.TianD.ZhaoT.et al (2019). Identifying a wild allele conferring small seed size, high protein content and low oil content using chromosome segment substitution lines in soybean.Theor. Appl. Genet.1322793–2807. 10.1007/s00122-019-03388-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  YangJ.HuC.HuH.YuR.XiaZ.YeX.et al (2008). QTL Network: mapping and visualizing genetic architecture of complex traits in experimental populations.Bioinformatics24721–723. 10.1093/bioinformatics/btm494

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  YangJ.ZhuJ.WliiliamsR. W. (2007). Mapping the genetic architecture of complex traits in experimental populations.Bioinformatics231527–1536. 10.1093/bioinformatics/btm143

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  YuC.QuZ.ZhangY.ZhangX.LanT.AdelsonD. L.et al (2017). Seed weight differences between wild and domesticated soybeans are associated with specific changes in gene expression.Plant Cell Rep.361417–1426. 10.1007/s00299-017-2165-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  YuG.WangL. G.HanY.HeQ. Y. (2012). Clusterprofiler: an R package for comparing biological themes among gene clusters.Omics16284–287. 10.1089/omi.2011.0118

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  ZhangD.LiH.WangJ.ZhangH.HuZ.ChuS.et al (2016a). High-density genetic mapping identifies new major loci for tolerance to low-phosphorus stress in soybean.Front. Plant Sci.7:372. 10.3389/fpls.2016.00372

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  ZhangJ.SongQ.CreganP. B.JiangG. L. (2016b). Genome-wide association study, genomic prediction and marker-assisted selection for seed weight in soybean (Glycine max).Theor. Appl. Genet.129117–130. 10.1186/s12864-019-6009-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  ZhangX.HinaA.SongS.KongJ.BhatJ. A.ZhaoT. (2019). Whole-genome mapping identified novel “QTL hotspots regions” for seed storability in soybean (Glycine max L.).BMC Genomics20:499. 10.1186/s12864-019-5897-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  ZhaoX.DongH.ChangH.ZhaoJ.TengW.QiuL. (2019). Genome wide association mapping and candidate gene analysis for hundred seed weight in soybean [Glycine max (L.) Merrill].BMC Genomics20:648. 10.1186/s12864-019-6009-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar