Genome-wide association and RNA-seq analyses reveal a potential gene related to linolenic acid in soybean seeds

Determining fatty acid content of different soybean cultivars

A total of 1,060 soybean lines from China were collected to determine their fatty acid content. To obtain more accurate phenotypic data for the GWAS, soybean lines were planted in three regions of China representing typical climates (subtropical, tropical, and temperate zones). The seeds were planted at the following locations: 2019–2020 on 124.65°E, 43.14°N; in 2020–2021 on 114.03°E, 38.29°N; in 2021–2022 on 113.01°E, 24.17°N. After harvesting, the seeds of different soybean lines were dried and the fatty acid content was measured. Five fatty acids (linolenic acid, stearic acid, oleic acid, palmitic acid, and linoleic acid) in soybean cultivar seeds were determined using matrix-assisted laser desorption/ionization time-of-flight imaging mass spectrometry (MALDI-TOF IMS). The correlation coefficients of fatty acids were calculated using SPSS software (SPSS Inc., Chicago, IL, USA). Each calculation was repeated thrice to ensure accuracy.

Genotyping and GWAS analysis

We used the CTAB method to extract genomic DNA from fresh leaves of each soybean line (Porebski, Bailey & Baum, 1997). The Beijing Biomarker Biotechnology Co. Ltd. (Beijing, China) provided the sequencing services. A total of 2,935,226 SNP markers with a minor allele frequency (MAF) > 0.06, obtained from genotyping, were screened for GWAS using the fastlmmc model. GWAS was used to analyze the association between the genotype and phenotype datasets. The Haploview 4.2 software was used to construct a Manhattan map. The threshold value was set to –log(p) > 4.20. The plink2 software was used to calculate the LD decay distances. The NR database and Swiss-Prot were used to search for candidate genes in the 650 kb genomic region of the important SNPs.

Plant transformation

The GmWRI14 gene from the soybean cultivar 010a (approval number 2012010) was cloned into the BamHI-SacI site of plasmid pTF101 named pTF101- GmWRI14-Flag, which was induced by the CaMV35S promoter, and the target gene was terminated by the NOS terminator. The recombinant plasmid was transformed into calli using Agrobacterium tumefaciens strain LBA4404 (Holsters et al., 1978). Individual T₀ plant lines were established in a greenhouse, and three independent transgenic soybean lines with the same plasmid were harvested: GmWRI14-1, GmWRI14-2, and GmWRI14-3.

The generation and identification of GmWRI14 gene-edited lines

Following a previous protocol, the novel CRISPR/Cas9 Vector pGES201 was used for soybean genome editing (Do et al., 2019). CRISPRdirect (http://crispr.dbcls.jp) was used to design the target site sequence based on the exon sequence of the GmWRI14 gene in the soybean genome. The GmWRI14 gene from soybean cultivar 010a (approval number 2012010) was cloned into the BamHI-SacI site of plasmid pTF101 named pTF101-GmWRI14-3×Flag, which was induced by the pM4 promoter, and the target gene was terminated by the NOS terminator. The CRISPR/Cas9 expression vector was constructed using the method described by Bai et al. (2020). The CRISPR/Cas9 expression vector was transformed into Agrobacterium tumefaciens LBA4404 using Agrobacterium-mediated soybean cotyledon node genetic transformation method (Holsters et al., 1978). The EasyPure® Plant Genome DNA Kit (TransGen, Beijing, China) was used according to the manufacturer’s instructions to extract total genomic DNA from the seed samples of the putative mutants. The T₀–T₂ generation seeds and control seeds of the transgenic plants were planted in a greenhouse under a 16/8-h light/darkness cycle. When the T₁–T₂ generation plants grew the first triple-compound leaf, they were screened for glyphosate resistance. The genomic DNA of genetically modified soybeans and controls was extracted using Kangwei Century Company’s new plant genome DNA extraction kit, and PCR amplification was performed using DNA as a template to detect positive plants. Positive plants were selected and the PCR products were sent to a sequencing company for sequencing. DSDecodeM was used to analyze the sequences of the T3 generation plants to characterize CRISPR/Cas9-induced mutations. Successfully edited strains were identified using sequence peaks and were aligned to the reference sequence. Heterozygous mutants exhibited overlapping peaks near the target site, whereas homozygous mutants exhibited a single peak at the target site. Three independent gmwri14 soybean mutant lines containing the same plasmid were harvested. They are gmwri14-1, gmwri14-2, and gmwri14-3.

The RNA-seq data analysis

To analyze the transduction pathways induced by GmWRI14, the total RNA of each sample (GmWRI14 transgenic soybeans, CK, and gmwri14 mutant) was isolated from the seeds for RNA sequencing, with three biological and three technical replicates, using the Super total RNA® extraction kit (TaKaRa, San Jose, CA, USA). A total of 2 μg RNA was obtained from each plant sample for the RNA sequencing. RNA quality was assessed using Nanodrop 2000c (Thermo Fisher Scientific, Waltham, MA, USA). All transcriptome data were analyzed and calculated using methods outlined in previous studies (Zhong et al., 2011; Sultan et al., 2012; Kumar et al., 2012). The RNA-sequencing depth was 10×, and the 2^−ΔΔCt method was used for the fold change of gene expression (Livak & Schmittgen, 2001). The original data were deposited in the China National GenBank database (accession number PRJNA435633). Each sample had an average of 5.68 GB of data for a total of 102.11 GB of sequencing data. We used bioinformatics software (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) to improve the data quality, resulting in 41.66 GB of clean data. The clean data of all samples reached 6 GB, and the base percentage of Q30 was more than 93%. The statistical power of this experimental design was calculated in cluster Profiler Power with a q value $\le$ 0.01. A KEGG statistical analysis (Kyoto Encyclopedia of Genes and Genome) and Gene Ontology analysis (http://www.geneontology.org/) were performed to identify the functions of the differentially expressed genes (DEGs).

Gene expression analysis by qRT-PCR

After 25 days of seedling growth, RNA samples were extracted from the seeds of soybean lines with low LA content and high LA content using an RNA extraction kit (TaKaRa, San Jose, CA, USA). Using overexpressed GmWRI14 soybeans, gmwri14 mutants, and control soybean seeds as materials, total RNA was extracted and reverse transcribed into cDNA at 10, 11, 12, 13, and 14 days after podding, using 24-h sampling (T0, T4, T8, T12, T16, T20, and T24 soybean seeds, respectively). qRT-PCR primers for bZIP54, GmFAD3C and GmFAD3B, were used to perform a 24-h expression regulation analysis of the candidate target genes of GmWRI14. qRT-PCR was performed with three biological replicates, according to the MIQE guidelines (Bustin et al., 2009). Real-time fluorescence analysis was performed using a LightCycler (Roche, Basel, Switzerland). The lectin gene (Gene ID: 100775957) was used as an internal control (Qin et al., 2022). All primers used are listed in the additional files (Table S1).

Yeast one-hybrid (Y1H) assay

To test whether GmWRI14 could affect the expression of FAD3s by binding to the promoter region of FAD3s, both the FAD3s (GmFAD3B and GmFAD3C) and the GmbZIP54 promoter regions were amplified from soybean cultivar 010a (Approval number 2012010). The Y1H assay was performed according to the manufacturer’s instructions (Zhao et al., 2019).

Identification of candidate gene GmWRI14 by GWAS

The LA content in 1,060 soybean lines collected from different regions of China between 2019 and 2021 ranged from 1.01% to 15.92%. These 3-year measurements suggest that there is a wide variation in LA content in the soybean population, the distribution is continuous, and the LA content conforms to a normal distribution (Fig. 1A). A total of 2,923,211 SNP markers were obtained from the dataset (Fig. 1B), of which 899,218 SNPs had an MAF > 6%. To identify potential candidate genes associated with LA, a GWAS was used to screen candidate genes using the fastlmmc model. These SNP markers were evaluated for their association with LA content. We identified 17 candidate genes associated with LA content from to 2019–2020, ten candidate genes from to 2020–2021, and ten candidate genes from to 2021–2022. Ten SNPs were associated with LA, all of which were located on chromosome 4. Glyma.04G127300.1, a candidate gene for LA, was detected by GWAS in each of the 3 years (Fig. 1C), indicating that Glyma.04G127300.1 is a candidate gene for LA. The results of the bioinformatics analysis showed that Glyma.04G127300.1 belongs to the plant WRI1 family containing two AP2/EREB domains; therefore, we named it GmWRI14. The expression of the GmWRI14 gene in the flowers, internodes, stems, and roots of the soybean plants was investigated using qRT-PCR. GmWRI14 was expressed in different organs of soybean plants but was mainly expressed in soybean seeds, with low levels of GmWRI14 expression detected in the other organs (Fig. 1D). The qRT-PCR results showed that the GmWRI14 gene expression was highly tissue-specific in soybeans. The raw qRT-PCR data are available in File S1.

The GmWRI14 gene from the soybean cultivar 010a (Approval number 2012010) was cloned into the BamHI-SacI site of the plasmid pTF101, named pTF101- GmWRI14-Flag (Fig. S1). GmWRI14 transgenic soybean and gmwri14 mutant were then produced (Figs. 2A–2D), and Southern blotting was used to detect GmWRI14 expression (Fig. S2), and MALDI-TOF MS was used to analyze the LA distribution in soybean seeds (Fig. S3). MALDI-TOF IMS results showed that the total fatty acid content of GmWRI14 transgenic soybean increased from 41.32% to 56.02% (Fig. 2E); however, the LA content of GmWRI14 transgenic soybean decreased from 12.24% to 6.12%. MALDI-TOF IMS results showed that, compared to soybean receptors (CK), the linolenic acid content in the gmwri14 mutant significantly increased from 21.01% to 23.92% (Fig. 2E). These results indicate that the expression of GmWRI14 and its ability to regulate specific genes are important for determining LA content in soybean seeds.

Overexpression of GmWRI14 suppresses GmFAD3s and bZIP54 expression in soybean seeds

RNA-seq was used to verify the genes related to fatty acid synthesis that were induced by GmWRI14 in transgenic soybean seeds and control plants (CK, GmWRI14 transgenic lines, and gmwri14 mutant). RNA extracted from the seeds was used for RNA-seq to analyze the differentially expressed genes (DEGs) between the three varieties (CK, GmWRI14 transgenic lines, and gmwri14 mutant). The RNA-seq results showed that the expression of more than 50 genes related to fatty acid synthesis changed significantly. The KEGG and GO analyses revealed that the 891 identified DEGs were mainly involved in fatty acid and linolenic acid metabolism (Fig. 3A). Weighted correlation network analysis (WGCNA) was performed to identify GmWRI14 co-expressed genes using the WGCNA R package. The correlation coefficient between any two genes was calculated using WGCNA to analyze the expression patterns of different genes. The square of the threshold of the correlation coefficient was 0.8, and 15,134 potential genes were transferred to the co-expression network to explore the relationship between phenotypes and candidate genes. A gene co-expression correlation matrix was defined by calculating the dissimilarity coefficients of different nodes, and a hierarchical clustering tree was constructed. Abnormal samples were removed based on clustering tree results. Different colors were used to distinguish the 13,211 genes, with gray representing the genes that could not be classified into any module (Fig. 3B). The 13,211 genes were divided into 45 distinct modules (Fig. 3B), and the turquoise, yellow, and light-yellow modules had the highest correlations with LA content in the module–trait correlation analysis (Fig. 3C). RNA-seq showed that DEGs related to linolenic acid metabolism were differentially expressed (more than 2.5 times), namely GmFAD3C and GmFAD3B (Fig. S4). qRT-PCR showed that the expression of GmWRI14 was significantly increased in the seeds of GmWRI14 transgenic soybean lines (Fig. 3D), and the expression of GmbZIP54 in GmWRI14 transgenic soybean seeds on different days was also measured by qRT-PCR (Fig. 3E), with GmbZIP54 expression ranging from 2.45 ± 0.52 to 8.44 ± 0.22. The expression of GmFAD3C and GmFAD3B was significantly decreased in the seeds of GmWRI14 transgenic soybean lines (Table S2), with GmFAD3C expression ranging from 22.75 ± 0.26 to 26.12 ± 0.26 and GmFAD3B expression ranging from 4.12 ± 0.15 to 9.23 ± 0.13 in GmWRI14 transgenic seeds (Figs. 3F and 3G). The expression levels of GmFAD3B and GmFAD3C were 24% and 31% higher, respectively, in control seeds than in the GmWRI14 transgenic soybean lines (Table S2). Notably, GmFAD3B and GmFAD3C expression remained virtually unchanged in the stems and leaves of GmWRI14 transgenic soybeans, with GmFAD3B and GmFAD3C expression in seeds ranging from 11.25 ± 0.41 29.22 ± 0.25. qRT-PCR results suggested that the GmWRI14 gene inhibited the expression of GmFAD3s (GmFAD3C and GmFAD3B) in soybean seeds (Table S2).

Disruption of GmWRI14 altered total fatty acid content, and expression of GmFAD3s increased

MALDI-TOF MS was used to profile fatty acid species in the mature seeds of gmwri114 mutant and control soybean lines. There was a significant difference between the gmwri114 and wild-type soybeans (P < 0.01). The LA content increased from 11.31% to 20.59% in the gmwri14 mutant, and owing to the reduction in total fatty acids, the size and quality of the gmwri14 mutant seeds were also reduced. The RNA-seq data showed that in gmwri14 mutant soybean seeds, the DEGs related to fatty acid synthesis from the two varieties (CK and GmWRI14 mutant) changed significantly (Fig. 4A), with a threshold of >2.5-fold change indicating significance. In the gmwri14 mutant, qRT-PCR results indicated that the lack of GmWRI14 induced the transcriptional activity of GmFAD3s or GmbZIP54 genes (Figs. 4B–4D). GmWRI14 is a transcriptional inhibitor that negatively regulates the expression of FAD3s. The abundance of FAD3s genes is closely related to linolenic acid content in different plant seeds (Flores et al., 2008; Hu et al., 2006). Therefore, the GmWRI14 gene likely regulates linolenic acid synthesis by regulating FAD3s genes. The qRT-PCR results suggested that disruption of GmWRI14 enhanced the expression of GmFAD3s and GmbZIP54 in soybean seeds.

GmWRI14 can inhibit GmFAD3s in the presence of GmbZIP54

To investigate whether the GmWRI14 gene binds to the promoter regions of GmFAD3B, GmFAD3C, and GmbZIP54 genes in soybean, a yeast one-hybrid assay (Y1H) was performed (Fig. 5A). The results showed that the GmWRI14 gene activated the GmbZIP54 promoter region, but not the promoter regions of the GmFAD3s genes. GmbZIP54 directly bound to the promoter region of GmFAD3s genes (Fig. 5A). When GmFAD3s were co-infiltrated with GmbZIP54, the promoter regions of the GmFAD3s (GmFAD3B and GmFAD3C) were significantly activated (Figs. 5B and 5C), with the promoter regions of the GmbZIP54 were significantly suppressed. The raw data for the dual-luciferase assay are shown in Table S3. Together, these data suggest that GmWRI14 inhibits the expression of GmFAD3s genes by interacting with GmbZIP54, leading to decreased linolenic acid production, and that strong inhibition of FAD3s can potentially be achieved by the combined action of GmbZIP54 and GmWRI14.