Full-length genes of 30 fatty acid desaturase (GmFADs) were predicted from the genome of Glycine max var. Williams 82 (Glycine max Wm82.a4.v1) using HMM search, and were listed in Supplementary Table S1 . The identified proteins corresponding to the GmFAD gene family had an amino acid sequence length from 235 aa (GmFAD6.2) to 453 aa (GmFAD7-2). The relative molecular mass of these proteins varied between 26,800.2 Da (GmFAD6.2) and 51,550.41 Da (GmSLD1.2). The isoelectric points (pI) of the proteins were distributed within a range of 5.94 (GmFAB2.1) to 9.51 (GmFAD5.1), with the majority possessing pI values above 7. The lipolysis index values ranged from 79.26 (GmSACPD) to 94.62 (GmSLD1.4). The instability coefficients of proteins were found to vary between 30.72 and 48.67, with GmFAD2.3 protein being the most stable and GmFAD6.1 protein the least stable among these members identified. Out of the 30 family members, 21 members were classified as hydrophilic proteins, whereas the others were hydrophobic ( Supplementary Table S2 ).

The subcellular localization analysis predicted that proteins GmFAD2, GmDES, and GmSLD were localized on endoplasmic reticulum, whereas GmFAD6 and GmFAB2 proteins localized on chloroplasts. The members of GmFAD5 and GmFAD3 along with their isozymes GmFAD7/FAD8 proteins were found in both endoplasmic reticulum and chloroplasts as detailed in Supplementary Table S2 . These results suggest that the subcellular location diversity of these proteins could be related to multiple functions.

To determine the evolutionary relationships of FAD proteins among A. thaliana (At), rice (Os), and soybean (Gm), a phylogenetic tree was constructed using the neighbor-joining (NJ) method with p-distance model using amino acid sequences of GmFADs (30), AtFADs (24), and OsFADs (18) ( Supplementary Table S3 ). As depicted in Figure 1A , all FAD proteins were classified into two distinct clusters: soluble (FAB2) and membrane-bounding FAD proteins (ADS, SLD, DES, FAD6, FAD2, and FAD3/FAD7/8).

FAB2 subfamily represents the only soluble desaturase identified so far. In this study, both A. thaliana and rice were found to contain seven homologous proteins, respectively ( Figure 1A ). In contrast, it was observed that soybean possesses four FAB2, namely three GmFAB2 and one GmSACPD, and among them, two alleles GmFAB2.1 and GmFAB2.2 were clustered together, when GmFAB2.3 aligned with GmSACPD in a different branch ( Figure 1B ). Compared FAB members in three species, except for OsFAB2.1, AtFAB2.1 and OsFAB2.1, most FAB members of soybean and A. thaliana were closely clustered together but OsFAB in rice was located at the base position of this subfamily. This suggests that FAB2 in rice may have evolved earlier and could have been divided into monocotyledonous and dicotyledonous taxa.

As shown in Figure 1A , proteins encoding membrane-bound desaturases were categorized into six distinct subfamilies. Among them, the ADS subfamily was further classified into two groups, composed of nine ADS from A. thaliana along with their two homologous FAD5 from soybean. Interestingly, GmFAD5.1 and GmFAD5.2 were closely grouped with AtADS3 and AtADS5 as a small branch, indicating that GmFAD5 might evolved from AtADS. Furthermore, no members from rice were identified in the ADS subfamily, indicating that it might have undergone evolutionary loss in rice or were transmitted to dicotyledons (A. thaliana and soybean) following their divergence from the last common ancestor.

The SLD subfamily comprised six GmSLD, two AtSLD, and one OsSLD, all of which encode sphingolipid Δ7 desaturases. Except for OsSLD1 was grouped with both GmSLD1.4 and GmSLD1.5 as a separate branch, the SLD members from each plant species were distinctly clustered together, indicating a closer genetic relationship among them. The DES subfamily, which is responsible for sphingolipid Δ4 desaturases, had the members including OsDES1, AtDES1, GmDES1.1 and GmDES1.2, with OsDES1 in rice localized at the basal position of the subfamily. This suggests GmDES1 in soybean as well as AtDES1 in A. thaliana had more close relationship as the dicotyledonous plants, but they diverged from OsDES1 from the last common ancestor.

Both FAD6 and FAD2 subfamilies encode Δ12 desaturases but were differently localized in the subcellular organelles. In the plastidial FAD6 subfamily, GmFAD6.2 notably occupied at a basal position within this subfamily. In the microsomal FAD2 subfamily, three rice OsFAD2 were clustered together, while six GmFAD2 and one AtFAD2 grouped in another branch, indicating this family has diverged into monocotyledonous and dicotyledonous groups from the same ancestor. Additionally, in the FAD3/FAD7/FAD8 subfamily, eight FADs from soybean (four GmFAD3, two GmFAD7 and two GmFAD8), three from Arabidopsis, and four from rice, were clustered together and encoded their corresponding microsomal or plastidial ω3 desaturases.

As shown in Figure 2 , different GmFAD genes were located in different soybean chromosomes. Except for the chromosomes 4, 5, 6, 12, and 16, the identified 30 GmFAD genes in soybean were found to distribute on the other 15 of the 20 soybean chromosomes. The maximum number of GmFAD genes were located on chromosome 2 as four genes. There were three genes located on the chromosomes 7, 14, and 18, respectively, while only one or two genes were located on other chromosomes. In addition, different gene members of each subfamily were mostly distributed in different chromosomes except that GmFAD2.4 and GmFAD2.5 was clustered in a small region of chromosome 19. This suggests that GmFAD genes are broadly dispersed throughout the soybean genome, and they might originate from diverse ancestors.

The diversity and expansion of gene families often arise from crosstalk and segmental duplication events, and gene duplication serves as a key mechanism for enhancing plant genetic diversity and the generation of novel genes. In this study, covariance analysis based on MCScanX was performed to investigate gene duplication events in the GmFAD gene family, and results showed that a large number of GmFAD genes had covariance between and within soybean chromosomes, and the most GmFAD genes with covariance were localized on the chromosome 19. As shown in Table 1 , both tandem and segmental duplication can be observed in the GmFAD family. Except for GmFAD2.4 and GmFAD2.5 on the chromosome 19 were caused by tandem duplication, a total of 29 segmental duplication events were identified ( Table 1 ). Among them, the ω-3 desaturase subfamily (GmFAD3/GmFAD7/GmFAD8) had the maximum duplicated gene pairs as 13 pairs followed by the GmSLD subfamily with seven pairs, while the ω-6 desaturase subfamilies (GmFAD2 and GmFAD6) had five duplicated gene pairs. The subfamilies of GmFAB2, GmDES and GmFAD5(ADS) had less duplicated gene pairs than other subfamilies in our study. Thus, the role of segmental duplication events plays a prominent role in increasing the genetic diversity of soybean GmFAD gene families rather than tandem duplication.

To assess the selection of duplicated GmFAD gene pairs, we calculated the substitution ratios between non-synonymous and synonymous (Ka/Ks) based on a whole genome analysis of gene duplication ( Table 1 ). In soybean, most Ka/Ks ratios for the GmFAD gene replication pairs were found to be less than 1, indicating that these genes are relatively conserved during the evolutionary process and have been subjected to purifying selection pressure to maintain gene functions and species stability. However, three gene pairs, such as GmFAD7-2/GmFAD3.4, GmFAD2.1/GmFAD2.2, and GmFAD2.2/GmFAD2.5, exhibited a Ks value greater than 1, and it indicates that a positive selection pressure occurs to produce new protein functions aiming to promote gene evolution and the adaptive changes of species.

Analysis of the exon/intron structures of 30 GmFAD genes identified in soybean ( Figure 3A , Supplementary Table S4 ), showed that the gene structure of different GmFAD subfamily was generally variable. Both GmFAD3/FAD7/FAD8 and GmFAD6 gene subfamilies had more than 6 exons, especially GmFAD6.1 possessed ten exons. GmFAD5 subfamily had five exons as compared to GmDES and GmFAB2 subfamilies which contained two to three exons. In addition, members in the GmFAD and GmFAB subfamilies had different number of exons, for example, GmFAD6.2 and GmFAD6.1 had four and ten exons, respectively. Interestingly, GmFAD2 and GmSLD had the simplest structures with only one exon. Furthermore, most GmFAD genes within the same subfamily exhibited high similarity in their exon/intron patterns as compared to those in different subfamilies, thus the exon/intron distribution could provide strong supports for the phylogenetic classification of GmFAD proteins above.

Furthermore, the gene motifs were analyzed in this study, and total 20 motifs were identified in GmFAD genes ( Figure 3B ). Evaluation of motifs 1-20 using the Pfam database (http://pfam.xfam.org) showed that motifs 2, 3, 6, 10, 11, 12, 16, 19, and 20 corresponded to membrane-FADS-like superfamily structural domains, while other three protein conserved structural domains PLN02505, Δ6-FADS-like, and PLN02598 were also identified as the membrane-FADS-like superfamily. Motifs 2, 3, and 10 were observed in all members of the ω-3 subfamilies (GmFAD3/GmFAD7/GmFAD8), while two unique motifs 16 and 20 were found in the GmFAD7 and GmFAD8 subfamilies. All the GmFAD2 subfamily members contained motifs 4, 6, and 11 except that GmFAD2.5 lacked motif 4. Motif 17 was present only in the GmFAD6 subfamily, but this motif is not functional. In the GmFAB2 subfamily, motifs 5, 14, 15 and 20 corresponded to the acyl-[acyl-carrier protein] desaturase structural domain (PLN00179), whereas Motif 11 was a domain of unknown function (with PLN00179 domain) and only existed in the branch of GmFAB2.3 and GmSACPD as compared to GmFAB2.1 and GmFAB2.2. The members of GmFAD5(ADS), GmDES, and GmSLD subfamilies comprised two, two, and four motifs, respectively. Notably, motif 9, as a component of the Cyt-b5 domain, was uniquely identified in the GmSLD subfamily. In short, members with similar conserved motifs and gene structures clustered together in the GmFAD gene family. The motif distribution patterns of each subfamily are similar, whereas differences between subfamilies may be related to subfamily functional convergence.

To elucidate the three-dimensional conformation of the GmFAD proteins, homology modeling method was adopted and homology modeling was performed using Swiss-Model online software. As shown in Figure 4 , GmFAD primarily consisted of α-helixs, irregular coils, and β-turns, and proteins clustered into the same clade exhibited analogous 3-dimensional (3D) structures. GmFAD3.3 within the ω-3 desaturase subfamilies and GmFAD2.6 within the ω-6 desaturase subfamilies exhibited distinct characteristics compared to GmFAD proteins from other subfamilies, indicating that there exists a balance between conservation and divergence within the GmFAD proteins.

Promoters located in the upstream regions of genes are crucial in gene expression regulation involving in plant growth and development as well as environmental adaptation. To understand the roles of GmFAD genes and the precise regulation of gene expression, cis-acting elements within the promoter regions were analyzed ( Figure 5 , Supplementary Table S5 ). The results showed that the promoter region of the GmFAD genes contained three distinct types of cis-elements. The first type of cis-elements was related to plant growth and development, including light responsiveness, circadian control, the differentiation of fenestrated chloroplasts, and cell cycle regulation. The second type was responsible for stress responses, including those triggered by methyl jasmonate (MeJA), abscisic acid (ABA), salicylic acid (SA), GA, growth hormone, low-temperature, drought, defense and stress responsiveness. The third category is associated with specific biological processes, such as those involved in endosperm expression, zein metabolism regulation, flavonoid biosynthesis regulation, and cell cycle regulation.

As illustrated in Table 2 , the promoters of most GmFAD genes contained cis-elements related to phytohormone responses, including ABRE, TGA, P-box, TATC-box, GARE-motif, CGTCA-motif, TGACG, AuxRR-core, and TGA-box, which correspond to ABA, SA, GA, MeJA, and Auxin signaling, respectively. Additionally, TC-rich repeats, which often function on the regulation of defense and stress responses, have been identified within the promoter sequences of 10 GmFAD genes (GmFAD3.3, GmFAD7-2, GmFAD2.1, GmFAD2.2, GmFAD2.5, GmDES1.2, GmSLD1.1, GmSLD1.4, GmSLD1.5, and GmSLD1.6). Moreover, the distribution patterns of TC-rich repeats within the GmFAD2 subfamily were similar to those in the GmSLD subfamily, but it is notably absent in the GmFAD6, GmFAD5, and GmFAB2 subfamilies.

All members of the GmFAD gene family possess abundant light-responsive cis-elements at the start codon upstream, with the detection of 17 such regulatory elements, such as Box4, G-Box, and GT1-Motif and others, suggesting that light signals play crucial roles in the accumulation of soybean fatty acid desaturases. Overall, these findings illuminate the potential functions of GmFAD genes in facilitating plant growth and development, enhancing stress responses, and modulating hormone signaling.

To gain a deeper understanding of the functional pathways associated with GmFAD genes, we conducted Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis to explore the potential functions of GmFAD. The results showed that GmFAD was mainly enriched in three GO terms: “molecular function”, “cellular component” and “biological process” ( Figure 6A ). Notably, the GmFAD gene family were significantly implicated in oxidoreductase activity, facilitating the transfer of electrons between paired donors and binding or reduction of molecular oxygen, and GmFAD genes also exhibited acyl-acyl-carrier-protein desaturase activity. Furthermore, GmFAD genes showed significant enrichment in membrane components, particularly endoplasmic reticulum membrane network, suggesting their potential involvement in the stability of membrane assembly. Additionally, both GO and KEGG annotations revealed that GmFAD genes were significantly involved in the biosynthetic and metabolic pathways of lipids and fatty acids ( Figure 6B ). In summary, GmFAD gene family participate in the biosynthesis and metabolism of fatty acid and lipid, and they also are responsible for redox and desaturation activities.

To investigate the regulation of GmFAD genes in soybean seed resistance to F. fujikuroi, the dominant seed decay pathogen, we analyzed expression patterns of 16 GmFAD genes in soybean seeds of susceptible cultivar (CX12) and resistant cultivar (ND25) after F. fujikuroi inoculation ( Figure 7 ). As shown in Figure 7 , all representative GmFAD genes belonging to ω-3 and -6 desaturase subfamilies were significantly induced by F. fujikuroi infection but their expression patterns differed in the resistant ND25 and susceptible CX12. Most genes in the susceptible CX12 were significantly up-regulated at much earlier infection (12 hpi), whereas they were strongly induced at 48 hpi in the resistant ND25. In particular, relative expression level of five genes including GmFAD8.2, GmFAD3.3, GmFAD7-1, GmFAB2.2, GmFAD2.3 and GmFAD2.5 were as higher as 6.0 fold at 48 hpi as compared to 0 hpi. In contrast, expression of GmFAB2.3 and GmSACPD was inhibited at the early infection of F. fujikuroi (6 hpi) in both cultivars, but their expression were more rapidly recovered and reached the peak at 12 hpi in the resistant ND25 rather than those in the susceptible CX12. Additionally, GmFAD2.4 were dramatically up-regulated by F. fujikuroi at 6 hpi and then decreased after 24 hpi in the susceptible CX12, but it had different expression pattern in the resistant ND25. Expression of GmFAD2.1 and GmFAD7.2 were weakly induced in the resistant ND25 upon F. fujikuroi infection but they were up-regulated differently in the susceptible CX12. This suggests that GmFAD genes play important roles in soybean seed resistance to F. fujikuroi infection.

The genome data of Glycine max var. Williams 82 were retrieved from the Genome Warehouse (GWH) with the Phytozome v13 database (Annotation version: Glycine max Wm82.a4.v1, https://phytozome-next.jgi.doe.gov/) (Xu et al., 2019). The published protein sequences and gene sequences of 24 AtFADs were obtained from TAIR release 10 (http://www.arabidopsis.org) (Shaheen et al., 2023), while the published sequences of 18 OsFADs were downloaded from RGAP release 7 (http://rice.plantbiology.msu.edu) (Zhiguo et al., 2019).

The hidden Markov model (HMM) profiles for the FA_desaturase (PF00487), FA_desaturase 2 (PF03405), and TMEM189 (PF10520) domains (Cheng et al., 2022) downloaded from the Pfam protein family database (https://www.ebi.ac.uk/interpro/entry/pfam/) were searched against the soybean (Glycine max var. Williams 82) protein data using an e-value threshold of ≤1e ^-5. The GmFAD gene family members were then filtered to eliminate duplicates and identify potential gene family members. Subsequently, the protein physicochemical properties and subcellular localization of the GmFAD members were further analyzed using Expasy (https://www.expasy.org/) and Cell-PLoc (http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/).

We performed multiple alignments of FADs were conducted by Clustal W using full-length protein sequences from soybean (30), A. thaliana (24), and rice (18), respectively. A Neighbor-Joining (NJ) phylogenetic tree was constructed using the p-distance model by MEGA7 (http://www.megasoftware.net) (Kumar et al., 2016). The bootstrap support values were calculated from 1000 repeats. Subsequently, the evolutionary tree was refined and visualized through the iTOL online platform (https://itol.embl.de/).

The annotation file for the general feature format version 3 (GFF3) in the soybean genome database was utilized to pinpoint the chromosome locations of the GmFAD genes. Visualization of the chromosomal localization of the GmFAD genes were achieved using TBtools software, which was based on the starting position on the soybean chromosome (Chen et al., 2020). The matched sequences covered more than 80% of the length of the longer gene, exhibited over 80% similarity within their respective regions, and were products of a single duplication event (Jiang et al., 2013; Singh and Jain, 2015). To further assess the evolutionary pressure on the GmFAD gene family, the synonymous (Ks) and non-synonymous (Ka) substitution rates of the GmFAD gene pairs were computed using TBtools, along with the Ka/Ks calculator (Bailey et al., 2006; Kong et al., 2013; Kumar et al., 2016). In addition, to explore the evolutionary relationships within the soybean species, MCScanX and BLASTP were employed to detect gene pairs that were co-variantly associated with FAD members (Xue et al., 2020).

The soybean coding sequence and genome file were applied to explore the splicing phase of the GmFADs family (Zhang et al., 2021). The Gene Structure Display Server (GSDS, http://gsds.cbi.pku.edu.cn/).and TBtools software (Chen et al., 2020) were employed to map the distribution of introns, exons and non-coding regions within genes. NCBI Conserved Domains database (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) was used to identify the conserved structural domains of GmFAD gene family proteins, and MEME (https://meme-suite.org/meme/) to analyze the conserved motifs of FAD gene family proteins with the following parameter settings: the motif discovery was set to classical mode with a predicted motif count of 20, and each motif was allowed to occur 0 or 1 times (Bailey et al., 2006). The Pfam database was then utilized to evaluate the functions of the aforementioned motifs (Mistry et al., 2021). The conserved motifs and structural domains of the FAD gene family proteins were simultaneously visualized, and their interactions were analyzed using TBtools. The integrated structure of GmFAD was constructed using the Swiss-Model platform (https://swissmodel.expasy.org/) based on fragments of iterative templates (Wei et al., 2022), and subsequently refined and visualized by Chimera software to yield a three-dimensional structural model.

Sequences encompassing the 1500 bp region upstream of the start codon (ATG) for each GmFAD gene were retrieved from the soybean genome database using TBtools. Subsequently, the promoter sequences of these GmFAD genes were uploaded to Plant CARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to identify cis-regulatory elements (Lescot et al., 2002).

The functional annotation of genes was conducted using the eggNOG-mapper database (http://eggnog-mapper.embl.de/). To further understand the potential pathways that might be associated with the GmFAD genes, the TBtools eggNOG-mapper Helper tool was subsequently employed to systematically organize and process the results derived from the eggNOG-mapper. Following this, the GO-basic file and the KEGG-backend file were exported for GO enrichment analysis and KEGG Pathway enrichment, respectively (Wei et al., 2022).

The fungal isolate F. fujikuroi (No. S100) was isolated from the rot seeds of soybean in the fields and identified using sequence analysis of translation elongation factor 1 alpha (EF-1α) and DNA-directed RNA ploymerase II second largest subunit (RPB2) (Chang et al., 2020a). The cultivar CX12 exhibited susceptibility to F. fujikuroi, whereas the cultivar ND25 showed a high resistance to the fungus, and both these two cultivars were chosen to examine the expression of GmFAD genes following inoculation of soybean seeds with F. fujikuroi.

Spore suspensions of F. fujikuroi were prepared following the protocol by Chang et al. (2022). For sporulation, a mung bean liquid medium was prepared by boiling 30 g of mung bean in 1 L of sterilized water for 20 min, filtering the mixture with cheesecloth, and then autoclaving at 121℃ for 30 min (Lv and Sun, 2007). Disease-resistant and susceptible soybean seeds were inoculated with F. fujikuroi suspensions at a concentration of 1 × 10⁶ spores per milliliter, supplemented with 0.1% Tween 20. As a control (CK), soybean seeds were treated with an equal volume of mung bean liquid medium. The treated and control seeds were then arranged on water agar medium (WA), and were incubated in the dark at a constant temperature of 25 ℃ for varying durations of 0, 6, 12, 24, and 48 h, respectively. The experiments were conducted in triplicate. At specified intervals post inoculation, soybean seeds from both control and treatment groups were collected, immediately frozen in liquid nitrogen, and subsequently stored at -80°C for RNA extraction.

The expression profiles of twenty-six GmFAD genes were analyzed using qRT-PCR with specific primer pairs listed in Supplementary Table S6 . Total RNA was extracted from soybean seed samples using Fast Pure^® Universal Plant Total RNA Isolation Kit (Vazyme-Bio, Chengdu, China) and subsequently assessed with a NanoDrop fluorometer (Thermo Fisher Scientific, Stuttgart, Germany) to ascertain the concentration and quality of RNA. First-strand cDNA was synthesized from 2.5 μg RNA in a 20-μL reaction volume according to the instructions of BeyoRT™II First Strand cDNA Synthesis Kit (Beyotime-Bio, Shanghai, China). The qRT-PCR assay was conducted on the Chromo4 Real-Time PCR System (Bio-Rad, CA, USA), within a 10-μL reaction mixture containing 5 μL of the SYBR qPCR Mix (2×) (Vazyme-Bio, Shanghai, China), 1 μL of each primer (10 mM), 1 μL of template DNA (10 ng), and 2 μL of ddH₂O. The PCR thermal cycling conditions were set as follows: an initial heat activation at 95°C for 3 minutes, followed by 40 cycles of amplification at 10 s at 95°C for 10 s, annealing temperature for 30 s, 15 s at 95°C, 1 min at 60°C, and a final 15 s at 95°C for melting curve analysis. Relative gene expression levels were determined using the 2^-ΔΔCt method (Livak and Schmittgen, 2001) with the GmActin gene serving as an internal reference for normalization. Each sample was subjected to four technical replicates and three biological replicates. Data analysis was preformed using SPSS 24 software (SPSS Software Inc., Chicago, IL, USA), and statistical analysis was conducted with Analysis of variance (ANOVA) and Duncan’s multiple range test (DMRT). The results were visualized with GraphPad Prism 10 (GraphPad Software Inc., San Diego, CA, USA).