Genome-wide identification and expression analysis of diacylglycerol acyltransferase genes in soybean (Glycine max)

Sequence retrieval and data analysis

To identify all putative DGAT proteins in soybean, G. max genome data were downloaded from the latest NCBI database (http://www.ncbi.nlm.nih.gov/). Data were collected as previously described in Zhao (2018). According to Arabidopsis DGAT1 and DGAT2 proteins in the Pfam database (http://pfam.xfam.org/), DGAT proteins should have a PF03982 domain and WS/DGAT should have Pfam: PF03007 and PF06974 domains. These protein domains were used as queries in Hmmer3.0 with the E-value set at 1E-005. Based on the amino acid sequence in peanuts, the DGAT3 (cytoDGAT) protein was searched and screened via BLAST in the soybean database. The exon/intron organization for G. max DGAT genes was defined by loading its coding sequence (CDS) and the corresponding genomic sequence into Phytozome13 (Phytozome). Molecular weights and isoelectric points of candidate GmDGAT proteins were analyzed using ExPASy (http://web.expasy.org/protparam/). Predicting subcellular localization of DGAT genes in soybean was performed using Cell-PLoc 2.0 (Chou & Shen, 2010) (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/), DeepLoc (Armenteros et al., 2017), and WoLF pSORT (www.genscript.com/wolf-psort.html) (Horton et al., 2007).

Conserved domain, motif, CDS, and phylogenetic analyses

To identify the conserved domains, motifs, and CDS in the GmDGAT proteins, we used the Pfam (http://pfam.xfam.org/) and MEME (https://meme-suite.org/meme/) databases and Gene Structure Display Server 2.0 (GSDS2.0, http://gsds.cbi.pku.edu.cn/), respectively. Using MEGA X (https://www.megasoftware.net/), we constructed a phylogenetic tree with the Jones–Taylor–Thornton model; bootstrap tests were performed using 1,000 replicates for statistical reliability.

Cis-element, chromosomal localization, and collinearity analyses of GmDGAT genes

The 2-kb region upstream of the start codon of the DGAT genes served as the promoter sequences that were submitted to PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) for predicting cis-acting elements in the promoter region. The chromosomal location of DGAT genes was obtained from Phytozome (https://phytozome-next.jgi.doe.gov/). MCScanX was used to determine the syntenic relationships among DGAT gene family members. The TBTools software was used to describe the distribution of the identified cis-elements, chromosome location map, and collinearity analysis results.

Estimating Ka/Ks ratios for duplicated gene pairs

In order to understand the mechanisms of DNA sequence evolution, reconstruct phylogenic trees, and identify protein coding exons, nonsynonymous (amino acid replacing) and synonymous (silent) substitution rates between protein coding sequences were used to determine whether there was selection pressure acting on this protein coding gene (Yang & Nielsen, 2000; Nekrutenko, 2002). Nonsynonymous (amino acid substitution) and synonymous (silent) substitutions are called Ka and Ks, respectively. Ka reflects nonsynonymous substitutions per nonsynonymous site, and Ks reflects synonymous substitutions per synonymous site. The Ka/Ks ratio (denoted as ω) is widely used as an estimator of selective strength for DNA sequence evolution, with ω >1 indicating positive selection, ω <1 indicating purifying (negative) selection, and ω close to 1 indicating neutral mutation (Chattopadhyay et al., 1998). Data are processed according to the protocol described by Aylward (2018), calculated the Ka/Ks ratios of GmDGATs.

GmDGAT expression pattern analysis

Based on the soybase database’s (https://soybase.org/soyseq/) published sequencing results of the soybean whole transcriptome, we screened relevant data and analyzed GmDGAT expression in different tissues from the expression data available for 51,529 soybean genes. We used Heatmap in the R software (version 2.15.2) to construct a heat map.

RT-PCR analysis of GmDGAT genes

The soybean variety Williams 82 (w82), used in this experiment, was provided by the Institute of Plant Sciences of Jilin University and planted at its variety trial station (43°55′22.63″N, 125°16′3216′33.00″E). Roots, stems, leaves, flowers, and seeds of soybean plants were collected 10, 20, 30, 40, 50, and 60 days after flowering (DAF) for subsequent analysis. Tissue samples were frozen in liquid nitrogen and stored in refrigerators at −80 °C. The E.Z.N.A.® plant RNA Kit (Omega Biotechnology, Norcross, GA, USA) was used for extracting total RNA from all collected samples according to the manufacturer’s instructions. The ReverTra Ace PCR RT Master Mix with gDNA Remover (TOYOBO, Japan) was used to reverse-transcribe RNA into first-strand cDNA. Subsequently, RT-PCR was conducted using QuantStudio™ 5 Real-Time PCR software (Foster City, CA, USA) and KOD SYBR qPCR Mix (TOYOBO) following the protocol proposed by the manufacturer. Three biological and three technical replicates were used for each treatment in the RT-PCR analysis. Soybean GmActin was used as a reference gene. The 2^−ΔΔC(t) method was used to determine the relative gene expression levels of 10 randomly selected GmDGAT genes. The primers used for RT-PCR are listed in Table S1.

Identification of DGAT proteins in soybean (G. max)

A total of 26 DGAT putative proteins were identified: three proteins containing the MBOAT (PF03062) protein domain, 10 proteins containing the DGAT (PF03982) protein domain, 12 proteins containing the WS/DGAT (PF03007 and PF06974) protein domains (Fig. 1Bi), which were identified in soybean in the latest version of the G. max genome assembly (Wm82.a4.v1). Using the peanut cytoDGAT sequence as a reference query, the sequence alignment program BLASTP identified one cytoDGAT protein (Glyma.13G118300) in soybean.

Among the identified DGAT proteins, three were assigned to the DGAT1 subfamily because they contained the membrane-bound o-acyltransferase (MBOAT) domain. In turn, 10 genes were assigned to the DGAT2 subfamily as they contained the monoacylglycerol acyltransferase domain; these genes were further classified into DGAT2a and DGAT2b according to the similarity of their encoded amino acid sequences and motifs. Meanwhile, the WS/DGAT subfamily identified 12 genes containing the acyl-CoA wax-alcohol acyltransferase domain; these genes were further classified into WS/DGATa, WS/DGATb, and WS/DGATc. CytoDGAT was found to contain only one gene.

We evaluated the basic properties of these DGAT proteins, including gene name, gene ID, CDS length, isoelectric point, molecular weight, and number of exons (Table 1). The amino acid (aa) length of DGATs varied from 249 to 704 aa and the isoelectric point (pI) ranged from 6.06 to 10.03. There were similarities in basic gene and protein properties within as well as among groups.

The subcellular localization prediction results showed that GmDGATs were centrally localized in the chloroplast, followed by the cell membrane, cytoplasm, nucleus, mitochondria, peroxisome, and golgi apparatus (Table 1). GmDGAT2a5, GmDGAT2b3, GmDGAT2b4, and GmDGAT2b5 failed to predict their subcellular location.

Conserved domains and phylogenetic analysis

Multiple alignment analyses were conducted to determine the classification and evolutionary relationships of DGAT genes in soybean, and an unrooted maximum-likelihood phylogenetic tree was constructed using the MEGA 5.10 software with the amino acid sequences of the 26 soybean DGAT proteins. Based on the phylogenetic tree (Fig. 1Ai), soybean DGAT genes were divided into four subfamilies: DGAT1, DGAT2, WS/DGAT, and cytoDGAT. Three DGAT proteins were included in the DGAT1 protein subfamily. The members of the DGAT2 subfamily were divided into two subgroups, each possessing five members of DGAT2a and DGAT2b. In turn, the WS/DGAT subfamily members were divided into three subgroups that possessed 7, 3, and 2 members from WS/DGATa, WS/DGATb, and WS/DGATc, respectively.

To understand the structural diversity of DGAT, the conserved domains (Fig. 1Bi) were analyzed with Pfam. All DGAT1 subfamily members contained the MBOAT domain, and all DGAT2 subfamily members contained the DGAT domain. Furthermore, the DGAT2b family included the Abhydrolase_6 domain, DGAT2b1 contained the Hydrolase_4 domain, and DGAT2b2 and DGAT2b5 contained the Abhydrolase_1 domain. Additionally, the WS/DGAT subfamily contained the WES_Acyltransf and WS_DGAT_C domains (Fig. 1Bii). In contrast, cytoDGAT did not contain any of the conserved domains but contained a unique 2Fe-2S_Thioredoxin domain. The conserved motifs of the DGAT genes were analyzed using the MEME motif search tool, which revealed 15 significant (e-value <e-100) conserved motifs (Fig. 1Aii). Different subgroups contained different motifs, including numbers and composition; however, the same group in the phylogenetic tree had the same conserved motif structures. Motif8 existed in all DGAT1 subfamily members. Subgroup DGAT2a had motif7, motif10, motif12, and motif14, and subgroup DGAT2b had motif11, motif13, and motif15; the WS/DGAT subfamily contained motif8 and motif9. In contrast, the cytoDGAT subfamily had none of these motifs.

The exon/intron organization (Fig. 1Aiii) showed similar numbers of introns and exons in one subgroup. All members of DGAT1 had 16 exons. The DGAT2a members had 9 exons, and DGAT2b subgroup members had 13–14 exons. The number of introns in the WS/DGAT subfamily varied from three to seven. WS/DGATa and WS/DGATb members mostly contained five exons, WS/DGATa1 contained three exons, WS/DGATc had six and seven exons, and CtyoDGAT had two exons.

Chromosomal localization and collinearity analysis

The chromosomal location of the DGAT genes was plotted according to the physical positions of the DGAT genes on soybean chromosomes (Fig. 2) (Kim et al., 2017). The 26 DGAT genes were unevenly distributed on 12 of the 20 chromosomes in soybean. Chromosomes 9, 13, and 16 had four genes; chromosome 19 had three genes; chromosomes 6, 12, and 18 had two genes; and chromosomes 3, 7, 11, and 17 had only one gene. To understand the evolution and amplification of DGAT genes in the soybean genome, we performed collinearity analysis. All segmentally duplicated gene pairs belonged to the same phylogenetic subfamilies, likely because, as an ancient tetraploid, soybean shows many chromosomal rearrangements, high genomic duplication, and multiple gene copies (Schmutz et al., 2010). DGAT1a, DGAT1b, and DGAT1c, located on chromosomes 9, 13, and 17, respectively, were collinear.

Synteny and collinearity is used to clarify angiosperm evolution and reveal patterns of differential gene loss after genome duplication and different gene retention associated with evolution of some morphological complexity (Tang et al., 2008). Furthermore, DGAT2a1, DGAT2a2, DGAT2a3, and DGAT2a4 showed collinearity, whereas DGAT2a5 was only collinear with DGAT2a4. In turn, DGAT2b1 and DGAT2b2 were collinear, and DGAT2b 4 was collinear with DGAT2b 3 and DGAT2b 5. Furthermore, WS/DGAT a6 and WS/DGAT a7 were collinear, and WS/DGAT b1 was collinear with WS/DGAT a1 and WS/DGAT a3. The results of collinearity analysis show that DGAT subfamilies evolved independently. Meanwhile, gene duplication and orphan genes occur in subfamilies. A collinear relationship was also found between WS/DGAT b2 and WS/DGAT b3. Due to the chromosome doubling event in soybean evolution, chromosome mapping and collinearity can help us analyze which genes are generated by chromosome reduplication and which genes are generated by DNA fragment doubling. Genes with collinearity may have come from the same original gene during gene evolution.

Selective pressure on DGAT genes

The replacement times of non-synonymous replacement sites (Ka) to the replacement times of synonymous site (Ks). The Ka/Ks ratio is usually used to insight the pressure on the selective evolution of each gene in the whole gene family. In this study, 35 paralogous gene pairs were determined among the DGAT genes using the Clustal W software (Table 2). Based on 6.161029 synonymous mutations for soybean per synonymous site per year, the Ka, Ks, and Ka/Ks values varied from 0.165 to 0.539, 0.138 to 4.24, and 0.164 to 4.412 Mya, respectively. Most gene pairs showed an estimated Ka/Ks >1, suggesting they are under significant positive selection pressure. Moreover, members of duplicated gene pairs seemingly indicated positive selection. For pairs DGAT2a3–DGAT2b2, WS/DGATa6–WS/DGATc2, DGAT2a4–DGAT2b5, DGAT2a3–DGAT2b5, WS/DGATa3–WS/DGATc2, WS/DGATb2–WS/DGATc2, and WS/DGATa5–WS/DGATa6, the Ka/Ks ratio was estimated at <1, suggesting purifying selection. Our syntenic analysis revealed that the soybean DGAT gene family expanded through both segmental and tandem duplications. Furthermore, most orthologous genes are distributed on different chromosomes (Gm09 contains four genes: DGAT1C, DGAT2a3, WS/DGAT6, and WS/DGAT12, and Gm13 contains four genes: DGAT1A, cytoDGAT, WS/DGAT8, and WS/DGAT10); similarly, most paralogous genes are distributed on different chromosomes (WS/DGATa3–WS/DGATa2 on Gm06, WS/DGATa1–WS/DGATa4 on Gm12, DGAT2a4–DGAT2a5 on Gm16, DGAT2b3–DGAT2b4 on Gm16, and WS/DGATa5–WS/DGATa7 on Gm18), indicating that most duplications of the DGAT gene family occur with genome replication.

Promoter analysis

Positions of promoters of 26 DGAT genes are shown in Fig. 3A. Some diverse distribution patterns of cis-acting elements were observed in the promoter region of GmDGAT genes, indicating that the DGAT gene family of soybean participates in various biological processes. The analysis of the promoter regions of the 26 DGAT genes revealed 25 cis-acting elements that were classified into three divisions (Fig. 3B): hormone response, biotic and abiotic response, and development and metabolism. Light-responsive elements occur in all DGAT gene promoters. Elements responsive to abscisic acid and light were found in the promoters of the DGAT1 subfamily. Except for DGAT2a1 and DGAT2b1, the DGAT2 subfamily members showed enriched light-responsive and MYB-binding cis-elements in their promoters. Anaerobic induction elements were found in the WS/DGATa1 gene promoters. Interestingly, although cytoDGAT does not share protein domains, it does share promoter features with other DGAT genes.

Spatial DGAT gene expression patterns

The expression patterns of soybean DGAT genes in various tissues (root, stem, leaf at different times, flower, seed (stages 1–9), nodule, and shoot) were analyzed from public databases and were found to exhibit evident tissue specificity (Fig. 4). The DGAT1 subfamily members are primarily expressed during the intermediate stages of grain development; meanwhile, the DGAT2 subfamily members are primarily expressed in the late stage of grain development. WS/DGAT subfamily members are expressed in flowers, leaves, and stems as vegetative growth proceeds. In turn, the cytoDGAT subfamily members are primarily expressed late during grain development.

In the DGAT1 gene subfamily, DGAT1a showed higher expression than DGAT1b or DGAT1c, and was expressed in most tissues. The expression of DGAT1b, similar to that of DGAT1c, was very low in most tissues, suggesting that DGAT1a is the most active gene of the DGAT1 subfamily. Meanwhile, in the DGAT2 subfamily, DGAT2a1, DGAT2a3, and DGAT2a4 showed high expression levels in the roots, stems, and flowers, respectively. In turn, DGAT2a1 showed high expression in the early stage of seed formation, whereas DGAT2a4 showed high expression in the late stage of seed formation. DGAT2b2, DGAT2b4, and DGAT2b5 were highly expressed in various tissues. According to the correlation between DGAT2b1 and DGAT2b2, the latter was considered the most important gene among them. As for DGAT2b3, DGAT2b4, and DGAT2b5, the highest level of expression was observed in the mature stage of seed formation; hence, the important gene in this cluster could not be inferred. The expression of WS/DGATa1, WS/DGATa2, WS/DGATa3, WS/DGATa4, WS/DGATa5, and DGAT2a7 in various tissues was very low, whereas that of WS/DGATa6 was very high in stems and leaves.

Considering that a waxy layer covers the outermost surface of the primary aerial plant tissues to reduce water loss, WS/DGATa6 was inferred to be the important gene in this cluster. WS/DGATb2 was highly expressed in the stems, leaves, and flowers, and it was the most highly expressed of the three genes in the WS/DGAT b cluster. Expression of WS/DGATc1 in all tissues analyzed was higher than that of WS/DGATc2. The 12 genes of the WS/DGAT subfamily were highly expressed in stems, leaves, and flowers, suggesting the involvement of this family in the increased production of wax esters in these tissues. Furthermore, the cytoDGAT gene showed a very high level of expression in the mature stage of seed formation and the highest expression level among all gene families during seed formation at stages R7 and R8. Finally, the expression of WS/DGATA1, WS/DGATA2, and WS/DGATB1 in various tissues was very low; therefore, we speculate that these may be pseudo or redundant genes.

Expression of DGAT genes

To further investigate the possible functions of GmDGATs, qRT-PCR was used to measure the expression patterns of nine genes highly expressed in subfamilies, in four tissues (root, leaf, stem, and flower), and six stages of soybean embryo (Fig. 5). Nine GmDGATs were expressed in the ten tissues, with strong tissue-specific expression patterns. Expression analysis of the GmDGAT gene family (Fig. 5) showed that the GmDGAT1 subfamily members are gradually expressed in maturing pods; in turn, CtyoDGAT subfamily members are mainly expressed in the stem meristems, leaves, and mature pods, while GmDGAT2a1 is expressed in leaves, 10 days after flowering (DAF) pods, and 20 DAF pods. As for GmDGAT2b1, this gene is expressed in flowers, stems, and leaves, and the GmWS-DGAT subfamily members are expressed in roots, flowers, stems, and leaves.