Genome-wide association studies reveal novel QTLs, QTL-by-environment interactions and their candidate genes for tocopherol content in soybean seed

Abstract

Genome-wide association studies (GWAS) is an efficient method to detect quantitative trait locus (QTL), and has dissected many complex traits in soybean [Glycine max (L.) Merr.]. Although these results have undoubtedly played a far-reaching role in the study of soybean biology, environmental interactions for complex traits in traditional GWAS models are frequently overlooked. Recently, a new GWAS model, 3VmrMLM, was established to identify QTLs and QTL-by-environment interactions (QEIs) for complex traits. In this study, the GLM, MLM, CMLM, FarmCPU, BLINK, and 3VmrMLM models were used to identify QTLs and QEIs for tocopherol (Toc) content in soybean seed, including δ‐Tocotrienol (δ‐Toc) content, γ‐Tocotrienol (γ‐Toc) content, α‐Tocopherol (α‐Toc) content, and total Tocopherol (T-Toc) content. As a result, 101 QTLs were detected by the above methods in single-environment analysis, and 57 QTLs and 13 QEIs were detected by 3VmrMLM in multi-environment analysis. Among these QTLs, some QTLs (Group I) were repeatedly detected three times or by at least two models, and some QTLs (Group II) were repeatedly detected only by 3VmrMLM. In the two Groups, 3VmrMLM was able to correctly detect all known QTLs in group I, while good results were achieved in Group II, for example, 8 novel QTLs were detected in Group II. In addition, comparative genomic analysis revealed that the proportion of Glyma_max specific genes near QEIs was higher, in other words, these QEIs nearby genes are more susceptible to environmental influences. Finally, around the 8 novel QTLs, 11 important candidate genes were identified using haplotype, and validated by RNA-Seq data and qRT-PCR analysis. In summary, we used phenotypic data of Toc content in soybean, and tested the accuracy and reliability of 3VmrMLM, and then revealed novel QTLs, QEIs and candidate genes for these traits. Hence, the 3VmrMLM model has broad prospects and potential for analyzing the genetic structure of complex quantitative traits in soybean.

Introduction

Soybean [Glycine max (L.) Merr.] is an important crop, and provided a great source of protein, oil, vitamin, and other nutrients for humans around the world. As one of the functional nutrients of soybean, tocopherol (Toc) has strong antioxidative capabilities and benefits to human health. It can scavenge free radicals in the body and increase immune function (Meagher et al., 2001; Kumar et al., 2009). According to the chemical structure, Tocs are composed of four members: α-tocopherol (α-Toc), β-tocopherol (β-Toc), γ-tocopherol (γ-Toc), and δ-tocopherol (δ-Toc) (Wan et al., 2008; Rozanowska et al., 2019; Barouh et al., 2022). Among them, α-Toc has the highest activity (Shaw et al., 2016). Edible oil is one of the main sources of Toc (Packer and Fuchs, 1993). As the most widely produced vegetable oil in the world, soybean oil has the highest total-Toc content, however, γ-Toc in soybean oil accounts for more than 70%. Although γ-Toc has antioxidant and other physiological activities, α-Toc is more excellent (Bramley et al., 2000). Hence, elevating the α-Toc content and total-Toc content in soybean genetics is important for quality improvement.

The Toc content of soybean seed is a typical quantitative trait, and it is difficult to breed this target trait of soybean variety using traditional breeding. This requires a lengthy selection process (Britz et al., 2008; Seguin et al., 2010). As an ancient tetraploid plant (Blanc and Wolfe, 2004), the soybean owing to its large and complex genome background brings great challenges and difficulties in genetic improvement (Young and Bharti, 2012; Tian et al., 2020; Lemay et al., 2022).

Genome-wide association studies (GWAS) is a powerful genomics tool, and it can base on natural populations to detect quantitative trait locus (QTL) underlying complex quantitative traits (Burton et al., 2007; Hamblin et al., 2011). GWAS has the advantage of high-resolution and high-throughput, thus, this method for analysis provides great convenience for the study of genetic variation in soybean (Anderson et al., 2020). Since the first GWAS conducted in soybean until now, almost all the important agronomic traits have been covered and dissected (Zhou et al., 2015; Fang et al., 2017). And yet, different GWAS models yield different GWAS results when we owe high-quality genotype and phenotype data (Chatterjee et al., 2013). Therefore, selecting the most suitable model for GWAS analysis can increase the accuracy to identify QTLs.

The general linear model (GLM) (Price et al., 2006), the mixed linear model (MLM) (Yu et al., 2006), and the compressed mixed linear model (CMLM) (Zhang et al., 2010) are single-marker genome-wide scan models, and these models can comprise a one-dimensional genome scan by testing one marker at a time. Among them, CMLM is frequently used in the genomic dissection of soybean quantitative traits (Jing et al., 2018; Zhao et al., 2019; Sui et al., 2020). However, single-marker genome-wide scan models require Bonferroni correction and multiple tests (Wang et al., 2016). Bonferroni correction is a stringent criterion, although greatly reduced false positive rates, many important loci associated with the target traits were missed (Zhang et al., 2019). With the rapid development of statistical methods, several multi-locus GWAS approaches have been developed to improve the power of QTL detection (Segura et al., 2012; Wen et al., 2018). Such as the Bayesian-information and linkage disequilibrium iteratively nested keyway (BLINK) (Huang et al., 2018), and the fixed and random model circulating probability unification (FarmCPU) (Liu et al., 2016). The obvious advantage of these methods is not a Bonferroni correction, they can reduce the amount of calculation and improve the accuracy.

Recently, a novel model was presented, named 3V multi-locus random-SNP-effect mixed linear model (3VmrMLM) (Li et al., 2022a). It is a multi-marker genome-wide scan model, this model not only provides high QTL detection power and sensitivity, at the same time, but it can also detect the QTL-by-environment interaction (QEI) and the QTL-by-QTL interaction (QQI). In this study, based on 23,149 SNPs and 175 soybean germplasms, we used six models (including 3VmrMLM, BLINK, FarmCPU, GLM, MLM, and CMLM) and conducted GWAS of individual and total-Toc content across three environments. The aim of this study is to reveal novel QTLs and QEIs of soybean Toc content and screen candidate genes.

Materials and methods

Plant materials, field trials, and phenotypic evaluation

The material used in this study included 175 diverse soybean accessions (Table S1), which encompassed most of the northeast regions of China and other countries. These materials were collected from the Chinese National Soybean GeneBank (CNSGB) and can represent the genetic diversity inside and outside of China. In this study, all experimental materials were planted at Harbin (117°17′E, 33°18′N), Liaoning (41°48′N, 123°25′E), and, Jilin (124°82′E, 43°50′N) in 2021. The field trials used a single-row plot (3 m-long rows and spaced 0.65 m) and were arranged in a randomized complete block design with three replicates per test environment. After full maturity, mature kernels of 10 randomly selected plants in each line were collected and used for evaluation of individual and total Toc content. The soybean seed Toc extraction and measurement were performed according to previous reports (Ujiie et al., 2005).

DNA isolation and sequencing

The genomic DNA of each sample from 175 tested accessions was isolated from young leaf was isolated by the method of CTAB (Han et al., 2015), and simplified-sequenced via specific locus amplified fragment sequencing (SLAF-seq) (Sun et al., 2013). The digest enzyme group of MseI (EC: 3.1.21.4) and HaeIII (EC: 3.1.21.4) (Thermo Fisher Scientific Inc, Waltham, MA, USA.) were used to obtain more than 50,000 sequencing tags, each 300-500 bp in length. The obtained markers were evenly distributed in unique genomic regions of the 20 soybean chromosomes. The short oligonucleotide alignment program 2 software (SOAP2) was used to align the raw paired-end reads to the soybean reference genome. Based on over 58,000 high-quality SLAF labels from each test sample, raw reads from the same genomic location were used to define SLAF groups. Genotypes were considered heterozygous if the minor allele depth or total allele depth of the sample was greater than 1/3 (Han et al., 2016).

Population structure evaluation and linkage disequilibrium analysis

The principle component analysis (PCA) was performed using the genome association and prediction integrated tool (GAPIT) R package to analyze the population structure of the natural panel (Lipka et al., 2012). The linkage disequilibrium (LD) parameter (r²) for estimating the degree of LD between pair-wise SNPs (MAF ≥ 0.05 and missing data ≤ 10%) was calculated by TASSEL 5.0 (Bradbury et al., 2007). Unlike GWAS, missing SNP genotypes were not classified as major alleles prior to LD analysis. Parameters in the program included MAF (≥ 0.05) and completeness (> 80%) for each SNP.

Genome-wide association studies

In total, 23,149 polymorphic SNP markers and 175 tested accessions were used to perform GWAS, it was performed using six models, including three single-locus model: MLM, GLM, CMLM, and three multi-locus models: FarmCPU, BLINK, 3VmrMLM. Among these, the GLM, MLM, CMLM, FarmCPU, and BLINK models were implemented with the R package “GAPIT” and visualization used scripts from the R package “qqman” (https://cran.r-project.org/package=qqman) and “CMplot “ (https://github.com/YinLiLin/R-CMplot).

The significant threshold value for the association between SNP and traits were determined by -log10 (P) ≥ 4, which is equivalent to P ≤ 0.0001, for MLM, GLM, CMLM, FarmCPU, and BLINK. The R software IIIVmrMLM (Li et al., 2022b) of the 3VmrMLM method (Li et al., 2022a) was downloaded from GitHub website (https://github.com/YuanmingZhang65/IIIVmrMLM). In this study, we used the single environment and multiple-environment methods to identify QTLs and QEIs. The significant threshold value was determined by LOD score ≥ 4.

Prediction of candidate genes

Candidate genes located in the 200-kb genomic region (100 kb upstream and 100 kb downstream) of each significant or suggested QTL then identified and annotated the candidate genes with the soybean reference genome (Wm82.a2.v1, http://www.soybase.org) (Cheng et al., 2017). The gene ontology (GO) enrichment analysis of candidate genes using the online tool (https://www.soybase.org/goslimgraphic_v2/dashboard.php). In addition, the whole genome and QEIs candidate genes among soybean relatives were compared using OrthoVenn2 (https://orthovenn2.bioinfotoolkits.net/task/create) (Xu et al., 2019).

Association analysis of candidate genes

Genome resequencing data were used to select the SNP variations within candidate genes. These SNP were located in exonic, intronic regions, upstream and downstream regions. Then, we combined the phenotype values of 56 soybean germplasms in three environments, these soybean germplasms were selected from the 175 diverse soybean accessions (Table S1) (including 9 high and low individual and total Toc germplasms), using the general linear model (GLM) in TASSEL 5.0 to identify SNPs of candidate genes that related to individual or total Toc content (Bradbury et al., 2007). Significant SNPs associated with the target trait were claimed when the test statistic was P < 0.01.

Haplotype analysis

The haplotypes were classified based on all of the SNPs with an MAF >0.05 in each candidate gene. Best linear unbiased predictors (BLUP) value were calculated using the “Phenotype” (https://cran.r-project.org/package=Phenotype) in R package. For each Toc component, haplotypes containing 18 soybean germplasms accessions were used for comparative analysis. One-way ANOVA and Two-tailed unpaired t -test were used to compare the differences in TC-BLUP value among the haplotypes. Finally, we compared the individual or total Toc content among these different haplotypes.

RNA-Seq data analysis of candidate genes

For candidate genes expression pattern analysis, first, we performed a differential expression pattern analysis at different tissues by downloading the RNA expression data from the plant public RNA seq database (PPRD) (http://ipf.sustech.edu.cn/pub/soybean/), which integrated all publicly available RNA-Seq soybean libraries (4,085) (Yu et al., 2022). Then, we also analyzed the expression of candidate genes in the development stage (R6) at different germplasms using the transcriptome data (unpublished data) from our laboratory. Additionally, we constructed a heat-map plot, and it was performed using the R package pheatmap (Kolde, 2012).

Quantitative real−time PCR (qRT−PCR)

Total RNA was isolated using the RNAprep pure Plant Kit (DP432, Tiangen). First-strand cDNA was synthesized from total RNA using TIANScript RT kits (KR104, Tiangen). And qRT-PCRs were performed using SYBR Green (FP205, Tiangen) reagents on an ABI 7500 fast real-time PCR platform. All qRT-PCRs were performed in three independent repeats, and the relative levels of transcript abundance were calculated using the 2^−ΔΔCT method (Livak and Schmittgen 2001). The GmActin4 (Glyma.12G063400) was used as an internal control for data normalization. Primer sequences for candidate genes were obtained from the qPrimerDB database (Table S2) (Lu et al., 2018).

Statistical analysis

Descriptive statistical analysis of phenotypic data including mean, minimum, maximum, coefficient of variation (CV), heritability, skewness, and kurtosis was performed using IBM SPSS statistics 25.0 (SPSS, Chicago, USA). One-way ANOVA with Dunnett’s multiple comparisons test and unpaired two-tailed t-test were performed using GraphPad Prism 9.4.1.

Results

Statistical and variation analysis of Toc content

Statistical analysis showed a wide range of phenotypic variations in the levels of the individual and total Toc content of the 175 soybean accessions from Harbin, Liaoning, and Jilin in 2021 (Table 1). The coefficient of variation (CV%), skewness, and kurtosis of Toc content of the association panel are also presented in Table 1. The CV varied a lot among different Toc content, especially the α-Toc content under three locations were observed from 35.21% to 44.9%, but all Toc content was no significant skewness or kurtosis (Figure 1). These results showed that Toc content was mainly influenced by genetic factors with less effect by environmental factors. Therefore, the tocopherol content of soybean in this study was appropriate for GWAS.

Figure 1

SNP genotyping, linkage disequilibrium estimating, and population structure for the GWAS panel

The genotyped samples included 175 soybean germplasms (including landraces and elite cultivars). The genomic DNA of these 175 accessions was sequenced using SLAF-seq. A total of 23,149 high-quality markers (MAF ≥ 0.05, missing data ≤ 10%) were identified from 153 million paired-end reads with 45 bp-read lengths and the sequencing depth was about 6.5 fold. The number of SNPs varied across the 20 soybean chromosomes. The highest number of SNPs was observed in Chr.18 (1732) and the lowest was detected in Chr.11 (685) (Figure 2A).

Figure 2

We assessed the mapping power of GWAS by the average distance of LD decay. The mean LD decay of the population was estimated at 97466 bp, when r² dropped to 0.2 (Figure 2B). Then, all 23,149 SNPs were used for scanning the population stratification of association panels through the principal component (PC), and evaluation of the variation of the first 10 PCs analysis revealed an inflection point at PC3, which demonstrated that the first 3 PCs dominated the population structure on the association mapping (Figures 2C, D). Additionally, a lower level of genetic relatedness among the 175 tested accessions based on pairwise relative kinship coefficients was observed (Figure 2E).

Quantitative trait locuss associated with Toc content by GWAS

GWAS was conducted using GLM, MLM, CMLM, FarmCPU, BLINK, and 3VmrMLM models. All of which accounted for kinship and population structure. First of all, we used different thresholds of significance (by -log10 (P) or LOD score= 3, 4, 5, 6, 7, 8, and 9) for testing six GWAS models and counted the number of QTLs detected (Figure 3A).Then, when -log10(P) ≥ 4 as significant thresholds, a total of 86 QTLs significantly associated with individual and total Toc concent in soybean seeds were detected via GLM, 18 QTLs were detected by MLM, 41 QTLs by CMLM, 41 QTLs by BLINK, and 34 QTLs by FarmCPU (Figure 4A, Figures S1–S5 and Tables S3–S7). Among them, only 4 QTLs were co-detected by all six models (Figure 3B). Furthermore, the largest number of QTLs were detected with the 3VmrMLM model. Among them, the single-environment method detected 101 QTLs (Figure S6, Table S8), the multiple-environments method detected 57 QTLs (Figure S7, Table S9), and 13 QEIs (Figure S8, Table S10). Among them, 11 QTLs were co-detected by single-environment and multiple-environment method (Figure 3C). The results showed that the number of QTLs detected by 3VmrMLM are more abundant and stable under different significance thresholds.

Figure 3

Figure 4

Finally, the QTLs, which were repeatedly detected in multiple GWAS models, were selected as reliable QTLs—group I. As shown in Figure 3B, Table 2, 19 QTLs were co-detected by at least three times or at least two models, which were distributed among 24 genomic regions in 14 chromosomes. Among these, 9 QTLs (rs9337368, rs1834346, rs17125409, rs330000, rs9782629, rs19530677, rs5680781, rs17266245, and rs53062844) were located in genomic regions or QTLs reported by previous studies, confirming the accuracy of QTL detection. We regard the remaining 15 QTLs as the novel QTLs (rs39895210, rs2960931, rs19310064, rs31044180, rs7543892, rs4992837, rs14593163, rs24979561, rs588498, rs19962490, rs6204830, rs8720462, rs37558520, rs34774232, and rs35815938). Moreover, a total of 161 QTLs were identified by 3VmrMLM (Figure 3A), in order to test the reliability of the 3VmrMLM model, we selected the QTLs only detected in 3VmrMLM. 9 QTLs (detected by at least two times) were repeatedly detected as specific QTLs—group II (Table 3), which were distributed among 9 genomic regions in 8 chromosomes. rs41784197 was located in genomic regions or QTLs reported by previous studies. Again, we regard the remaining 8 QTLs as the novel QTLs (rs7167202, rs9140707, rs18105573, rs2669053, rs40595691, rs43000771, rs5779917, and rs46814888).

Prediction of candidate genes for Toc content in soybean seeds

Based on annotations for the soybean reference genome in SoyBase, we further predicted candidate genes within the 200-kb flanking regions of the novel QTLs. In two group novel QTLs, a total of 248 genes were obtained (Table S11). And a total of 134 genes were obtained in QEIs (Table S12). Then, we used GO annotation to perform enrichment analysis for group I and group II genes. The results categorized as molecular function, cellular component, and biological process, were shown in Figure 4. Both group I and group II candidate genes are involved in a variety of functions, such as carbohydrate metabolic process, translation, protein binding, cytoplasm component, DNA binding, and so on.

Comparative genome analysis

In order to predict the authenticity of the QEIs, firstly, we selected four closely related species, Glyma_max, Vigna_radiate, Vigna_augularis, and Phaseolus_vulgaris,for comparative genomic analysis. A total of 12847 core gene clusters were found in the four species, and 1197 gene clusters were unique to Glyma_max (Figure 5A), specific genes clusters account for 5.4% (1197/22159). Then, we used candidate gene of QEIs for comparative genomic analysis, 12 gene clusters were unique to candidate gene of QEIs (Figure 5B), specific genes clusters account 9.23% (12/130), this result shown that these QEIs have more abundant specific genes. As shown in Figure 5C, these specific genes are involved in various biological processes, metabolic processes, response to stimulus, etc. More detailed statistics on the number of shared gene clusters are shown in Figure 5D. Figure 6E is count of proteins by type of cluster.

Figure 5

Figure 6

Gene-based association analysis of candidate genes

Two groups of candidate gene association analysis were performed using the GLM model with the TASSEL, using the genome resequencing of 56 germplasms (including 9 high and low individual and total Toc germplasms). A total of 4537 SNPs with MAF ≥ 0.05 were identified among 248 candidate genes. Among them, a total of 50 SNPs from 11 candidate genes were found to reach the threshold with -log10(P) ≥ 2.0 (Table S13), of these, 4 SNPs are located in upstream regions, 10 SNPs are located in intronic regions, 26 SNPs are located in exonic ;regions, and 10 SNPs are located ;in downstream regions. Those SNPs are considered to be significantly associated with individual and total Toc concentrations in soybean seeds. Among these genes, 4 candidate genes from group I and 7 candidate genes from group II. These genes can be considered potential candidate genes for individual and total Toc-related. For example, as shown in Figure 6A, the significant SNPs correlated to α‐Toc and δ‐Toc on basis of association analysis for two candidate genes were respectively identified (Glyma.17G188700 and Glyma.20G235100 were shown in Figure 6A, others were shown in Figure S9).

Haplotype analysis of candidate genes

For the haplotype analysis, first, all the SNP markers within each gene are used to construct haplotypes. Then, we performed one-way ANOVA with TC-BLUP values of each soybean accession. The results are shown in Table 4, each gene contains haplotypes that are significant differences from TC-BLUP values. In addition, 14 haplotypes of 11 candidate genes respectively conferred an increased individual and total Toc content in soybean seeds (Glyma.17G188700 and Glyma.20G235100 were shown in Figure 6B, others were shown in Figure S10). Therefore, these haplotypes are beneficial and can be adjusted for individual and total Toc content in soybean seeds.

RNA-Seq data analysis of candidate genes for Toc content in soybean

In order to confirm the possible effect of candidate genes in the regulation of Toc content, we firstly used PPRD to analyze the expression patterns of 11 candidate genes in different tissues. The result showed that all candidate genes were expressed in soybean seed (Figure S11), and Glyma.10G171600 is most abundantly expressed in seed compared with other tissues. Then, for the 11 candidate genes of 56 soybean germplasms at the development stage (R6), RNA-Seq data analysis was done. The result showed that the expression levels of the 11 candidate genes in low and high Toc content germplasms were different. Among them, Glyma.17G188700 can regulate α-Toc content in soybean seeds. The range of the expression levels of Glyma.17G188700 in higher α-Toc germplasms was much higher than those of lower. Other genes regulate Toc content as shown in Figure 7. Interestingly, Glyma.01G054800, Glyma.09G032100, and Glyma.10G171600 can regulate both the γ-Toc and Total-Toc content. Glyma.09G032100 in higher γ-Toc and total-Toc germplasms were much higher than those expression levels of lower. However, Glyma.01G054800, and Glyma.10G171600 in higher γ-Toc germplasms have higher expression levels, but in higher total-Toc germplasms have lower expression levels. Moreover, these candidate genes results of qRT-PCR are consistent with the RNA-seq data (Figure S12).

Figure 7

Discussion

As one of the vitamin E family members, Toc plays a crucial role for humans, plants, and animals (Bramley et al., 2000). For humans, daily Toc supplementation can decrease the risk for cancer and cardiovascular disease (Shaw et al., 2016). For plants, Toc can protection of chloroplasts from photooxidative damage (Munne-Bosch and Alegre, 2002). For animals, Toc must be added to animal feed to improve and maintain growth and health (Pinelli-Saavedra et al., 2008). Soybean is a major crop used worldwide as a source of food, oil, and animal feed. Soybean oil compared to other oil crops contains a higher total Toc content, but γ-Toc comprises 70% (Park et al., 2019). The physiological activity of γ-Toc was lower than that of α-Toc (Wan et al., 2008). Therefore, increasing the α-Toc and total Toc content in soybean seeds is important to improve the nutritional variety and feed quality of soybean. However, the genetic background of Toc content is complex quantitative inheritance. The reason why quantitative traits are complex is that they are controlled by unequal polygenes and are susceptible to environmental influences. In this study, individual and total Toc content of 175 soybean accessions were evaluated. The results showed that the Toc content of tested germplasms was relatively stable to the environment, and Toc content had a wide range of variation among the different germplasms.

GWAS has been widely used in the mining of QTLs in most crops including soybean. It is a method to identify the genetic variation among the natural populations to establish genetic markers based on linkage disequilibrium (LD) (Yano et al., 2019; Xiao et al., 2022). How improve the power of GWAS has been a major challenge for the last decade. In recent years, a variety of new methods have been proposed, with the rapid development of computing technology and sequencing technology (Wang et al., 2016; Huang et al., 2018; Xiao et al., 2021; Li et al., 2022a). Although this propelled much of the practicability of GWAS, it is particularly important to select the appropriate sequencing method and suitable model for improving the positioning efficiency according to the research needs (Liu et al., 2017; Kim et al., 2021). For this study, we adopted six models (GLM, MLM, CMLM, BLINK, FarmCPU, and 3VmrMLM), to conduct GWAS of Toc content in soybean seeds. And the results were divided into two groups, revealed a total of 23 novel QTLs, other QTLs were located in the regions of QTLs in previous studies or overlapped our previous GWAS studies, and these known QTLs are all covered by 3VmrMLM.

3VmrMLM is a new algorithm, different from other algorithm, the 3VmrMLM use single-marker genome-wide scanning to select potentially associated markers and uses empirical Bayes and the likelihood ratio test in a multi-locus model to identify significant QTLs, this undoubtedly improves its detection capability (Li et al., 2022a). Additionally, it can be simultaneously estimated in a vector manner that QEI and QQI effects. Although the QQI detection in this study did not achieve good results, the 3VmrMLM still showed better detection ability than the GLM, MLM, CMLM, BLINK, and FarmCPU, indicating a more reliable tool for complex trait dissection.

In soybean and other plants, only a few definite genes have been characterized, associated with an individual or total Toc. Among them also includes most of the key enzyme genes (Dwiyanti et al., 2011; Zhang et al., 2013). To accurately screen candidate genes, we selected a total of 248 genes within the 200-kb flanking regions of the 23 novel QTLs and using a gene-based association by the GLM method, a total of 11 genes were finally determined to be significantly related to individual or total Toc in soybean seeds. Moreover, almost all these genes have beneficial haplotypes. Glyma.06G038000 encoded alpha/beta-Hydrolases superfamily protein. Glyma.01G054800 encoded plant proteins of unknown function, Glyma.03G186500 encoded a WD-40 repeat family protein, Glyma.20G235800 encoded a WD40 repeat-like superfamily protein, Glyma.03G186200 is a RAB GTPase homolog C2A, Glyma.10G171600 encoded a RAB GTPase homolog A5A, Glyma.17G188700 encoded transposas, Glyma.09G032100 encoded a myb domain protein, Glyma.20G235100 encoded an indeterminate domain protein, Glyma.20G235400 encoded a P-loop containing nucleoside triphosphate hydrolases superfamily protein. Of these genes, Glyma.01G054800 and Glyma.10G171600 are the most special, and these two genes are higher expressed in higher γ-Toc content germplasms, but lower expressed in higher total-Toc content germplasms. The soybean oil contains a higher proportion of γ-Toc, this is very different from the other oil crops (Cahoon et al., 2003). Therefore, we conclude that the Glyma.01G054800 and Glyma.10G171600 inhibited the transformation of α-Toc and δ-Toc, resulting in the excessive accumulation of γ-Toc, while the total-Toc content decreased. This requires further experiments to prove. The precise functions and mechanisms of 11 candidate genes will be planned in future studies.

In general, the 3VmrMLM algorithm achieved good results in the GWAS. In this study, Toc content in soybean seed in group I QTLs, 10 known QTLs are all covered by 3VmrMLM. The results of GO enrichment analysis showed that group I; and group II candidate genes had similar GO biological process terms. for the 11 candidate genes finally identified in this study, 7 genes were alone identified by the 3VmrMLM. All candidate genes were able to detected by the 3VmrMLM. In addition, a higher percentage of the Glyma_max specific genes have also been found in candidate genes near QEIs by comparative genomic analysis. These results have preliminarily determined the detection efficiency of the 3VmrMLM algorithm. Thus, we hope that using 3VmrMLM could be used to dissect more important complex quantitative traits in the future, and this algorithm is advantageous to promoting the development of soybean breeding.

Funding

This study was financially supported by National Key Research and Development Project of China (2021YFF1001204), the Chinese National Natural Science Foundation (31971967, 31871650), National Key Research and Development Program of China (2021YFD1201604, 2019YFD1002601), the Youth and Middle-aged Scientific and Technological Innovation Leading Talents Program of the Crops (2015RA228), the National Ten Thousand Talent Program (W03020275), Postdoctoral Scientific Research Development Fund of Heilongjiang Province (LBH-Z15017, LBH-Q20004), Program on Industrial Technology System of National Soybean (CARS-04-PS06).

Acknowledgments

This study was conducted in the Key Laboratory of Soybean Biology of the Chinese Education Ministry, Soybean Research & Development Center (CARS) and the Key Laboratory of Northeastern Soybean Biology and Breeding/Genetics of the Chinese Agriculture Ministry.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

• 1

  AlbertE.KimS.Magallanes-LundbackM.BaoY.DeasonN.DaniloB.et al. (2022). Genome-wide association identifies a missing hydrolase for tocopherol synthesis in plants. PNAS119 (23), e2113488119. doi: 10.1073/pnas.2113488119

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AndersonR.FernandezC.YuanY.GoliczA.EdwardsD.BayerP. (2020). Method for genome-wide association study: A soybean example. Method Microbiol.2107, 147–158. doi: 10.1007/978-1-0716-0235-5_7

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BarouhN.Bourlieu-LacanalC.Figueroa-EspinozaM. C.DurandE.VilleneuveP. (2022). Tocopherols as antioxidants in lipid-based systems: The combination of chemical and physicochemical interactions determines their efficiency. Compr. Rev. Food Sci. F.21 (1), 642–688. doi: 10.1111/1541-4337.12867

  • CrossRef
  • Google Scholar

• 4

  BlancG.WolfeK. H. (2004). Widespread paleopolyploidy in model plant species inferred from age distributions of duplicate genes. Plant Cell.16 (7), 1667–1678. doi: 10.1105/tpc.021345

  • CrossRef
  • Google Scholar

• 5

  BradburyP. J.ZhangZ.KroonD. E.CasstevensT. M.RamdossY.BucklerE.SBucklerE.S (2007). eTASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics23 (19), 2633–2635. doi: 10.1093/bioinformatics/btm308

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BramleyP.ElmadfaI.KafatosA.KellyF.ManiosY.RoxboroughH.et al. (2000). Vitamin e. J. Sci. Food Agric.7), 80, 913–938. doi: 10.1002/(SICI)1097-0010

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BritzS. J.KremerD. F.KenworthyW. J. (2008). Tocopherols in soybean seeds: genetic variation and environmental effects in field-grown crops. J. Am. Oil Chem. Soc85 (10), 931–936. doi: 10.1007/s11746-008-1286-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BurtonP. R.ClaytonD. G.CardonL. R.CraddockN.DeloukasP.DuncansonA. (2007). Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature447, 661–678. doi: 10.1038/nature05911

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  CahoonE. B.HallS. E.RippK. G.GanzkeT. S.HitzW. D.CoughlanS. J. (2003). Metabolic redesign of vitamin e biosynthesis in plants for tocotrienol production and increased antioxidant content. J. Nat. Biotechnol.21, 1082–1087. doi: 10.1038/nbt853

  • CrossRef
  • Google Scholar

• 10

  ChatterjeeN.WheelerB.SampsonJ.HartgeP.ChanockS. J.ParkJ. H. (2013). Projecting the performance of risk prediction based on polygenic analyses of genome-wide association studies. Nat. Genet.45 (4), 400–405e4053. doi: 10.1038/ng.2579

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  ChengY.MaQ.RenH.XiaQ.SongE.TanZ.et al. (2017). Fine mapping of a phytophthora-resistance gene RpsWY in soybean (Glycine max l.) by high-throughput genome-wide sequencing. Theor. Appl. Genet.130 (5), 1041–1051. doi: 10.1007/s00122-017-2869-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  DwiyantiM.S.YamadaT.SatoM.AbeJ.KitamuraK. (2011). Genetic variation of γ-tocopherol methyltransferase gene contributes to elevated α-tocopherol content in soybean seeds. BMC Plant Biol.11, 152. doi: 10.1186/1471-2229-11-152

  • CrossRef
  • Google Scholar

• 13

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

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  HamblinM.BucklerE. S.JanninkJ.-L. (2011). Population genetics of genomics-based crop improvement methods. Trends Genet.27, 98–106. doi: 10.1016/j.tig.2010.12.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  HanY.ZhaoX.CaoG.WangY.LiY.LiuD.et al. (2015). Genetic characteristics of soybean resistance to HG type 0 and HG type 1.2.3.5.7 of the cyst nematode analyzed by genome-wide association mapping. BMC Genomics16, 598. doi: 10.1186/s12864-015-1800-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  HanY.ZhaoX.LiuD.LiY.LightfootD. A.YangZ.et al. (2016). Domestication footprints anchor genomic regions of agronomic importance in soybeans. New Phytol.209, 871–884. doi: 10.1111/nph.13626

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  HuangM.LiuX.ZhouY.SummersR. M.ZhiwuZ. (2018). BLINK: a package for the next level of genome-wide association studies with both individuals and markers in the millions. Gigascience8, 1–12. doi: 10.1093/gigascience/giy154

  • CrossRef
  • Google Scholar

• 18

  JingY.ZhaoX.WangJ.TengW.QiuL.HanY.et al. (2018). Identification of the genomic region underlying seed weight per plant in soybean (Glycine max l. merr.) via high-throughput single-nucleotide polymorphisms and a genome-wide association study. Front. Plant Sci.91392. doi: 10.3389/fpls.2018.01392

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  KimM.-S.LozanoR.KimJ.BaeD.KimS.ParkJ.-H.et al. (2021). The patterns of deleterious mutations during the domestication of soybean. Nat. Commun.12 (1), 97. doi: 10.1038/s41467-020-20337-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  KoldeR. (2012) Pheatmap: Pretty Heatmaps. Available at: https://CRAN.Rproject.org/package=pheatmap (Accessed Aug 15, 2022).

  • Google Scholar
  • View reference in article

• 21

  KumarV.RaniA.DixitA. K.BhatnagarD.ChauhanG.S (2009). Relative changes in tocopherols, isoflavones, total phenolic content, and antioxidative activity in soybean seeds at different reproductive stages.J. Agr Food Chem.57 (7), 2705–2710. doi: 10.1021/jf803122a

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  LemayM. A.SibbesenJ. A.TorkamanehD.HamelJ.LevesqueR. C.BelzileF. (2022). Combined use of Oxford nanopore and illumina sequencing yields insights into soybean structural variation biology. BMC Biol.20 (1), 53. doi: 10.1186/s12915-022-01255-w

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  LiH.LiuH.HanY.WuX.TengW.LiuG.et al. (2010). Identification of QTL underlying vitamin e contents in soybean seed among multiple environments. Theor. Appl. Genet.120 (7), 1405–1413.

  • Google Scholar
  • View reference in article

• 24

  LiH.WangY.HanY.TengW.ZhaoX.LiY.et al. (2016). Mapping quantitative trait loci (QTLs) underlying seed vitamin E content in soybean with main, epistatic and QTL x environment effects. Plant Breed135 (2), 208–214. doi: 10.1111/pbr.12346

  • CrossRef
  • Google Scholar

• 25

  LipkaA.TianF.WangQ.PeifferJ.LiM.BradburyP.et al. (2012). GAPIT: genome association and prediction integrated tool. Bioinformatics28 (18), 2397–2399. doi: 10.1093/bioinformatics/bts444

  • CrossRef
  • Google Scholar

• 26

  LiuH.CaoG.WuD.JiangZ.HanY.LiW. (2017). Quantitative trait loci underlying soybean seed tocopherol content with main additive, epistatic and QTL × environment effects. Plant Breeding.136(6), 924–938. doi: 10.1111/pbr.12534

  • CrossRef
  • Google Scholar

• 27

  LiuX.HuangM.FanB.BucklerE. S.ZhangZ. (2016). Iterative usage of fixed and random effect models for powerful and efficient genome-wide association studies. PloS Genet.12 (2), e1005767. doi: 10.1371/journal.pgen.1005767

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  LivakK. J.SchmittgenT. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2^-ΔΔCT method. Methods.25, 402–408. doi: 10.1006/meth.2001.1262

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  LiM.ZhangY. W.XiangY.LiuM. H.ZhouY. H. (2022b). IIIVmrMLM: the r and c++ tools associated with 3VmrMLM, a comprehensive GWAS method for dissecting quantitative traits. Mol. Plant15 (8), 1251–1253. doi: 10.1016/j.molp.2022.06.002

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  LiM.ZhangY. W.ZhangZ. C.XiangY.LiuM. H.ZhouY. H.et al. (2022a). A compressed variance component mixed model for detecting QTNs and QTN-by-environment and QTN-by-QTN interactions in genome-wide association studies. Mol. Plant15 (4), 630–650. doi: 10.1016/j.molp.2022.02.012

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  LuK.LiT.HeJ.ChangW.ZhangR.LiuM.et al. (2018). qPrimerDB: a thermodynamics-based gene-specific qPCR primer database for 147 organisms. Nucleic Acids Res46 (D1), 1229–1236. doi: 10.1093/nar/gkx725

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  MeagherE. A.BarryO. P.LawsonJ. A.RokachJ.FitzGeraldG. A. (2001). Effects of vitamin e on lipid peroxidation in healthy persons. JAMA285 (9), 1178–1182. doi: 10.1001/jama.285.9.1178

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  Munne-BoschS.AlegreL. (2002). The function of tocopherols and tocotrienols in plants. J. Crit. Rev. Plant Sci.21, 31–57. doi: 10.1080/0735-260291044179

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  PackerL.FuchsJ. (1993). Vitamin e in health and disease. J. Crc Press.

  • Google Scholar
  • View reference in article

• 35

  ParkC.DwiyantiM. S.NaganoA. J.LiuB.YamadaT.AbeJ. (2019). Identification of quantitative trait loci for increased alpha-tocopherol biosynthesis in wild soybean using a high-density genetic map. BMC Plant Biol.19 (1), 510. doi: 10.1021/jf100455f

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  Pinelli-SaavedraA.Calderón de la BarcaA. M.HernándezJ.ValenzuelaR.ScaifeJ. R. (2008). Effect of supplementing sows' feed with alpha-tocopherol acetate and vitamin c on transfer of alpha-tocopherol to piglet tissues, colostrum, and milk: aspects of immune status of piglets. Res. Vet. Sci.85 (1), 92–100. doi: 10.1016/j.rvsc.2007.08.007

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  PriceA. L.PattersonN. J.PlengeR. M.WeinblattM. E.ShadickN. A.ReichD.. (2006). Principal components analysis corrects for stratification in genomewide association studies. Nat. Genet.38 (8), 904–909. doi: 10.1038/ng1847

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  RozanowskaM.EdgeR.LandE. J.NavaratnamS.SarnaT.TruscottT. G. (2019). Scavenging of retinoid cation radicals by urate, trolox, and α-, β-, γ-, and δ-tocopherols. IJMS.20 (11), 2799. doi: 10.3390/ijms20112799

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  SeguinP.TremblayG.PageauD. (2010). Soybean tocopherol concentrations are affected by crop management. J. Agric. Food Chem.58 (9)5495–5501. doi: 10.1021/jf100455f

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  SeguraV.VilhjálmssonB.PlattA.KorteA.SerenÜ.LongQ. (2012). An efficient multi-locus mixed-model approach for genome-wide association studies in structured populations. Nat. Genet.44, 825–830. doi: 10.1038/ng.2314

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  ShawE. J.KakudaY.RajcanI. (2016). Effect of genotype, environment, and genotype×environment interaction on tocopherol accumulation in soybean seed. J. Crop Sci.56, 40–50. doi: 10.2135/cropsci2015.02.0069

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  ShawE.RajcanI. (2017). Molecular mapping of soybean seed tocopherols in the cross AC Bayfield X OAC Shire. Plant Breeding136, 83–93. doi: 10.1111/pbr.12437

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  SuiM.JingY.LiH.ZhanY.LuoJ.TengW.et al. (2020). Identification of loci and candidate genes analyses for tocopherol concentration of soybean seed. Front. Plant Sci.11. doi: 10.3389/fpls.2020.539460

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  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 (3), e58700. doi: 10.1371/journal.pone.0058700

  • CrossRef
  • Google Scholar

• 45

  TianD.WangP.TangB.TengX.LiC.LiuX.et al. (2020). GWAS atlas: a curated resource of genome-wide variant-trait associations in plants and animals. Nucleic Acids Res.48 (D1), D927–D932. doi: 10.1093/nar/gkz828

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  UjiieA.YamadaT.FujimooK.EndoY.KitamuraK. (2005). Identification of soybean varities with high levels of α-tocopherol content. Breed Sci.55 (2), 123–125. doi: 10.1270/jsbbs.55.123

  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  WangS.-B.FengJ.-Y.RenW.-L.HuangB.ZhouL.WenY.-J.et al. (2016). Improving power and accuracy of genome-wide association studies via a multi-locus mixed linear model methodology. Sci. Rep.6, 19444. doi: 10.1038/srep19444

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  WanJ.ZhangW.JiangB.GuoY.HuC. (2008). Separation of individual tocopherols from soybean distillate by low pressure column chromatography. J. Am. Oil Chem. Soc85 (4), 331–338. doi: 10.1007/s11746-008-1198-x

  • CrossRef
  • Google Scholar

• 49

  WenY.ZhangH.NiY.HuangB.ZhangJ.FengJ. (2018). Methodological implementation of mixed linear models in multi-locus genome-wide association studies. Brief. Bioinf.19, 700–712. doi: 10.1093/bib/bbw145

  • CrossRef
  • Google Scholar

• 50

  XiaoQ.BaiX.ZhangC.HeY. (2022). Advanced high-throughput plant phenotyping techniques for genome-wide association studies: A review. J. Adv. Res.35, 215–230. doi: 10.1016/j.jare.2021.05.002

  • CrossRef
  • Google Scholar

• 51

  XiaoJ.ZhouY.HeS.RenW.-L. (2021). An efficient score test integrated with empirical bayes for genome-wide association studies. Front. Genet.12. doi: 10.3389/fgene.2021.742752

  • CrossRef
  • Google Scholar

• 52

  XuL.DongZ.FangL.LuoY.WeiZ.GuoH.et al. (2019). OrthoVenn2: a web server for whole-genome comparison and annotation of orthologous clusters across multiple species. Nucleic Acids Res.47, W52–W58. doi: 10.1093/nar/gkz333

  • CrossRef
  • Google Scholar

• 53

  YanoK.MorinakaY.WangF.HuangP.TakeharaS.HiraiT.et al. (2019). GWAS with principal component analysis identifies a gene comprehensively controlling rice architecture. Proc. Natl. Acad. Sci.116 (42), 21262–21267. doi: 10.1073/pnas.1904964116

  • CrossRef
  • Google Scholar

• 54

  YoungN.BhartiA. (2012). Genome-enabled insights into legume biology. Annu. Rev. Plant Biol.63, 283–305. doi: 10.1146/annurev-arplant-042110-103754

  • CrossRef
  • Google Scholar

• 55

  YuJ.PressoirG.BriggsW. H.VrohB. I.YamasakiM.DoebleyJ. F.et al. (2006). A unifed mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat. Genet.38, 203–208. doi: 10.1038/ng1702

  • CrossRef
  • Google Scholar

• 56

  YuY.ZhangH.LongY.ShuY.ZhaiJ. (2022). Plant public RNA-seq database: a comprehensive online database for expression analysis of ~45 000 plant public RNA-seq libraries. Plant Biotechnol. J.20 (5), 806–808. doi: 10.1111/pbi.13798

  • CrossRef
  • Google Scholar

• 57

  ZhangZ.ErsozE.LaiC.TodhunterR.TiwariH.GoreM. (2010). Mixed linear model approach adapted for genome-wide association studies. Nat. Genet.42, 355–360. doi: 10.1038/ng.546

  • CrossRef
  • Google Scholar

• 58

  ZhangY.JiaZ.DunwellJ. (2019). Editorial: the applications of new multi-locus GWAS methodologies in the genetic dissection of complex traits. Front. Plant Sci10. doi: 10.3389/fpls.2019.00100

  • CrossRef
  • Google Scholar

• 59

  ZhangL.LuoY.ZhuY.ZhangL.ChenR.XuM.et al. (2013). GmTMT2a from soybean elevates the α-tocopherol content in corn and arabidopsis. Transgenic Res.22 (5), 1021–1028. doi: 10.1007/s11248-013-9713-8

  • CrossRef
  • Google Scholar

• 60

  ZhanY.LiH.SuiM.ZhaoX.JingY.LuoJ.et al. (2020). Genome wide association mapping for tocopherol concentration in soybean seeds across multiple environments. Ind. Crops Products.154, 1–15. doi: 10.1016/j.indcrop.2020.112674

  • CrossRef
  • Google Scholar

• 61

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

  • CrossRef
  • Google Scholar

• 62

  ZhouZ.JiangY.WangZ.GouZ.LyuJ.LiW.et al. (2015). Resequencing 302 wild and cultivated accessions identifies genes related to domestication and improvement in soybean. Nat. Biotechnol.33 (4), 408–414. doi: 10.1038/nbt.3096

  • CrossRef
  • Google Scholar