Genome-wide identification and comparative analysis of diacylglycerol kinase (DGK) gene family and their expression profiling inBrassica napusunder abiotic stress

Background

Diacylglycerol kinases (DGKs) are signaling enzymes that play pivotal roles in response to abiotic and biotic stresses by phosphorylating diacylglycerol (DAG) to form phosphatidic acid (PA). However, no comprehensive analysis of theDGKgene family had previously been reported inB. napusand its diploid progenitors (B. rapaandB. oleracea).

Results

In present study, we identified 21, 10, and 11DGKgenes fromB. napus,B. rapa, andB. oleracea, respectively, which all contained conserved catalytic domain and were further divided into three clusters. Molecular evolutionary analysis showed that speciation and whole-genome triplication (WGT) was critical for the divergence of duplicatedDGKgenes. RNA-seq transcriptome data revealed that, with the exception ofBnaDGK4andBnaDGK6,BnaDGKgenes have divergent expression patterns in most tissues. Furthermore, someDGKgenes were upregulated or downregulated in response to hormone treatment and metal ion (arsenic and cadmium) stress. Quantitative real-time PCR analysis revealed that differentBnaDGKgenes contribute to seed oil content.

Conclusions

Together, our results indicate thatDGKgenes have diverse roles in plant growth and development, hormone response, and metal ion stress, and in determining seed oil content, and lay a foundation for further elucidating the roles ofDGKsinBrassicaspecies.

在植物、p等重要的信号分子hosphatidylinositol lipids phosphatidic acid (PA), diacylglycerol (DAG), and some lysophospholipids can be activated by distinct environmental stresses, triggering signal transduction cascades and endowing plants with stress resistance [1,2]。Among these molecules, PA is a key second messenger in the responses to various stresses. PA levels are typically low, accounting for 0.67% of total phospholipids, and PA formation depends on the perception of extracellular stimuli [3], such as drought [4], salinity [5], chilling damage [6], osmotic pressure [7], wounding [8], and pathogen attack [9], or of phytohormones such as ethylene [8], abscisic acid (ABA) [10], brassinolide (BR) [11], and gibberellic acid (GA) [12,13]。Plants’ ability to overcome such stress events depends on the role of PA in signal perception and transduction [2,14,15]。PA is a precursor to all phosphoglycerolipids and its bioaccumulation is pivotal for lipid metabolic flux and membrane composition [16]。PA is generated in the plasma membrane by two different biosynthetic pathways in eukaryotic cells. In the first pathway, phospholipase D (PLD) hydrolyzes phosphatidylcholine (PC) and phosphatidylethanolamine (PE) to generate PA directly. In the second pathway, phospholipase C (PLC) hydrolyzes phosphatidylinositol-4,5-bisphosphate to produce DAG, which is in turn phosphorylated by diacylglycerol kinase (DGK) to yield PA [16]。

Previous studies have demonstrated that PA interacts with ABI1 phosphatase 2C, thus promoting ABA signaling inA. thaliana[10]。PLDαgene expression, protein levels, and enzymatic activity all increase inA. thalianawhen the leaves are treated with ethylene [8]。InB. napusroots, exogenous application of 24-epibrassinolide (EBL) affects PA synthesis through the PLC/DGK pathway under optimal salinity conditions [17]。Furthermore,DGKgenes are rapidly but transiently regulated in different plant tissues in response to beneficial elements and other ions, including silver (Ag), aluminum (Al), arsenic (As), cadmium (Cd), chromium (Cr), mercury (Hg), and sodium (Na) [18]。Notably,DGKgenes may be candidates for influencing seed oil content inB. napus最近的未发表的数据显示,我们的洛杉矶boratory. Additionally, diacylglycerol transferase (DGAT) catalyzes DAG to generate triacylglycerol (TAG), which is a major component of vegetable oils in oilseed crops [19]。Therefore, the relative flux of TAG synthesis from de-novo-synthesized or PC-derived DAG can greatly affect the final seed oil content. TheDGKgene family has been widely characterized in the context of plant stress tolerance, including tomato (Solanum lycopersicum,SlDGK) [20],Arabidopsis thaliana(AtDGK) [21], maize (Zea mays,ZmDGK) [22], rice (Oryza sativa,OsDGK) [23], apple (马吕斯有明显,MdDGK) [4], and soybean (Glycine max,GmDGK) [24]。However, the genomic analysis of theDGKgene family inBrassicaspecies has not been reported.

Brassica napus(AACC, 2n = 38), a typical allotetraploid of theBrassicagenus, is an important oil crop planted worldwide, which originated from the hybridization and polyploidization ofB. oleracea(CC, 2n = 18) andB. rapa(AA, 2n = 20) [25].Whole-genome sequences ofB. rapa,B. oleracea,andB. napushave recently been assembled [25,26,27], providing valuable resources for studying theDGKgene family inBrassicaspecies.In this study, we identifiedDGKsin theBrassicaspecies and investigated their gene structures, conserved domains, protein properties, evolution, and cis-acting elements through bioinformatics analysis. We also evaluated the expression patterns ofBnaDGKsin different tissues and their responses to hormone treatment and metal ion induction. Additionally, we studiedBnaDGKexpression patterns in twoB. napuscultivars with different seed oil contents. This genome-wide identification ofDGKgene family members in threeBrassicaspecies provides strong evidence of functional homologies among theseDGKgenes inBrassicaspecies.

Identification and chromosomal distribution ofDGKgenes inB. napus,B. rapa, andB. oleracea

We identified 21, 10, and 11DGKgenes inB. napus,B. rapa, andB. oleracea, respectively.B. napushad the same number ofDGKgenes asB. rapaandB. oleraceacombined, indicating thatDGKgenes may not have experienced a gene-loss event during polyploidization. In addition, theB. napus DGKproteins (BnaDGKs) were between 450 and 720 amino acids (aa) in length, corresponding to coding sequence (CDS) lengths of 1353 to 2163 bp; their molecular weights ranged from 49.93 kDa to 79.30 kDa; and their pIs from 5.88 to 8.93 (Additional file8: Table S1).

The number of exons in theDGKgenes varied from 7 to 15, with the maximum inBnaDGK6–2(Additional file8: Table S1). Analysis of the physical and chemical characteristics of 21DGKproteins inB. rapaandB. oleracearevealed that, upon calculation, the average values were approximately equal inB. napusand its diploid progenitors.

In plants, previous studies were conducted on the subcellular localization ofDGKsto the nucleus, plasma membrane, cytoskeleton, and chloroplast [28,29]。In present study, the subcellular location of allDGKswas predicted using the PSORT website (Additional file8: Table S1). Results showed that most members of cluster 1 were located in the endoplasmic reticulum, exceptBnaDGK1–2,BnaDGK2–1,BolDGK1–2, andBolDGK2–2in the nucleus, in accordance with thatAtDGK1andAtDGK2位于内质网的膜amino-terminal hydrophobic segments that are sufficient to target proteins to the cell membrane to play their role in many signal transduction processes [30]; the members of cluster 2 were located in different parts, for example, most ofDGK3andDGK7were located in the peroxisome andDGK4swere localized in the chloroplast, whileBnaDGK3–3andBolDGK3–2occurred in the mitochondrion; whileDGK5andDGK6in cluster 3 were predicted to be located in the peroxisome and cytoplasm, respectively (Additional file8: Table S1). These subcellular localization results were in accordance with the published report [24] and suggested thatDGKshave been widely conserved in the same clusters. Furthermore, the diverse localization of plantDGKsimplies that they might be actively involved in different cellular processes during development and abiotic stress.

Additionally, we investigated the chromosomal localization ofDGKgenes inBrassicaspecies based on their physical position in the GFF3 files. In total, 20 of the 21 full-lengthDGKsofB. napusmapped to the assembled 11 chromosomes (9 in An genomes and 11 in Cn genomes), while only one gene was assigned to Ann (random A-genome chromosome). InB. oleracea, 7 of the 11 full-lengthDGKswere distributed on 5 chromosomes, with the remaining 4 genes randomly located on 4 scaffolds. InB. rapa, all full-lengthDGKswere positioned on 6 of the 10 assembled chromosomes (Additional file1: Fig. S1). Notably, chromosomes An03, An09, Cn03, Ar03, and Ar09 contained the maximum number of genes, i.e., three, whereas chromosomes Cn02, Cn07, Cn09, Co03, and Co07 each contained two genes. This result suggests thatDGKshave uneven distributions in different species. Statistical analysis showed that seven gene pairs have maintained their relative positions between the An subgenome ofB. napusand the Ar subgenome ofB. rapa, and six gene pairs have maintained their relative positions between the Cn subgenome ofB. napusand the Co genome ofB. oleracea(Additional file1: Fig. S1).

Multiple sequence alignments and sequence characterization of theDGKgenes

To explore the structure of DGK protein sequences, we performed multiple sequence alignments using ClustalX and visualized the results with Jalview and DOG 2.0. The sequences of two diacylglycerol/phorbol ester (DAG/PE)-binding domains (C1 domain, PF00130; Additional file2: Fig. S2) indicate that the first and second DAG/PE-binding domains harbor the sequences HX₁₄CX₂CX_19 ~ 22CX₂CX₄HX₂CX₇C and HX₁₈CX₂CX₁₆CX₂CX₄HX₂CX₁₁C, respectively. The two extremely conserved C₆/H₂cores were also observed in cluster 1DGKsof other plants, such as soybean [24], apple [4], and maize [22]。Moreover, the alignment revealed that, except for the orthologs ofAtDGK2, almost allDGKspossess a presumed ATP-binding motif with a GXGXXG consensus sequence in their catalytic domain (DGKc), where G represents glycine and X represents any other amino acid (Additional file3: Fig. S3). Furthermore, in the homologs ofAtDGK2, glycine (G) is replaced by alanine (A), but this change does not render the enzyme inactive [21]。Interestingly,BnaDGK2–1, with the shortest protein sequence in cluster 1, lacked the C₆/H₂cores, upstream basic region, and an extended cysteine-rich (extCRD-like) domain; it seems probable thatBnaDGK2–1lost these domains during evolution.

c的功能域的原理图luster 1 DGK proteins drawn using DOG 2.0 showed that compared to the other two clusters, two DAG/PE binding domains (C1), the upstream basic region, and a conserved 15-amino-acid extension were specific to cluster 1 (Additional file4: Fig. S4). In the upstream basic region, YT and VP residues remained as front and back boundaries, respectively, except in homologs ofAtDGK1, in which VP was replaced with TP. Furthermore, the second DAG/PE-binding domain included a preserved 15-aa extension, an extCRD-like domain; whereas the combination of the extCRD-like domain and DGKc domain is indispensable to the function of the functionally active DGKs [31]。Overall, theDGKsof cluster 1 inA. thaliana, rice, apple, soybean, andB. napus, as well as its two diploid progenitors have a universal framework: (YT-upstream basic region-VP/TP) – (3 aa) – (DAG/PE-binding domain 1) – (12 aa) – (DAG/PE-binding domain 1/extCRD-like domain) – (~ 130 aa) – (DGKc/DGKa domain). These results suggest that the structure ofDGKs在集群1非常保守。

Phylogenetic tree and syntenic analysis of theDGKsinB. napus,B. rapa, andB. oleracea

To explore the evolutionary relatedness ofDGKsbetweenA. thalianaandBrassica, we performed the multiple sequence alignment ofDGKsusing MUSCLE. Subsequently, we constructed a phylogenetic tree based on the DGK protein sequences, including 7, 21, 10, and 11 protein sequences fromA. thaliana,B. napus,B. rapa, andB. oleracea, respectively (Fig.1). Generally, the gene number in theB. rapa, andB. oleraceagenome was notably less than three times theA. thalianagene number because some genes might be lost during polyploidy speciation. For example, theA. thalianagenesAtDGK4,AtDGK6, andAtDGK7have at most only two homologs in each of the threeBrassicaspecies. By contrast, the remainingAtDGKorthologs possess at least three homologs inB. napus. Additionally, the phylogenetic tree is clearly divided into three clusters, designated clusters 1, 2, and 3, in accordance with previous studies inA. thaliana[21]。Among them, theDGKgenes inB. napuswere closely associated with their corresponding genes inB. oleraceaandB. rapain each clade, suggesting that these genes might have undergone whole-genome duplication events from diploid parental species (B. oleraceaandB. rapa) to allotetraploidB. napus, which had similar functions.

We also performed collinearity analysis ofDGKsinA. thalianaand the threeBrassicaspecies to explore the evolutionary relationship ofDGKs(Fig.2). A total of 22 collinear gene pairs betweenBrassicaandA. thaliana被识别,包括9 AT&Bna基因对,7AT&Bra gene pairs, and 6 AT&Bol gene pairs (Additional file8: Table S2). Furthermore, these 21DGKsinB. napuswere identified as three homologous gene pairs on homologous chromosomes, while three pairs of genes were not located on homologous chromosomes. For example, the gene pairBnaDGK2–3&BnaDGK2–1was distributed on the An06 and Cn03 chromosomes, respectively, possibly as a result of chromosomal rearrangement and segmental duplication (Additional file8: Table S2). Our results suggested that theDGKsmight undergo gene duplication or loss from the diploid parental species (B. oleraceaandB. rapa) to the allotetraploidB. napus. Overall, these findings provide the clues for investigating the expansion mechanisms and functional characteristics ofDGKfamily genes inBrassicaspecies.

The expansion ofDGKgenes inBrassicaspecies

We also studied theDGKgene duplication types using MCScanX to determine the expansion patterns of theDGKgene family inBrassicaspecies. A total of 101,040 genes annotated in theB. napusgenome [25] were examined in this study, and 21,680 genes (21.46%) were dispersed and 74,035 genes (73.27%) might have undergone segmental duplication. Among them, all identifiedDGKgenes belong to segmental duplication except for one dispersed geneBnaDGK4–1(Additional file8: Table S1). Moreover, allDGKgenes ofB. oleraceaandB. rapawere derived from segmental duplication. Our findings showed that segmental duplication appears to play an important role inDGKgene expansion inBrassicaspecies.

The allopolyploidB. napuswas formed by hybridization ofB. rapaandB. oleracea, and its estimated formation time was approximately ~ 7500 years ago [25]。To explore the selective pressure on theDGKsafter duplication events, we calculated the nonsynonymous (K_a) and synonymous (K_s) substitution rates and theK_a/K_sratio for the 69 identified syntenic gene pairs in the threeBrassicaspecies andA. thaliana(Additional file8: Table S2).K_a/K_s = 1 signifies that genes have experienced neutral selection, whereasK_a/K_s > 1 orK_a/K_s < 1 indicate that genes have experienced positive or negative selection, respectively. As a result, theK_a/K_svalues between theBrassicaspecies andA. thalianaranged from 0.08 to 0.27. TheK_a/K_svalues for most duplicatedDGKgene pairs inB. napusand each of its diploid progenitors (B. rapaandB. oleracea) were < 1, except for one gene pair (BnaDGK7–2&BolDGK7–1) with a value > 1. Furthermore, for three gene pairs (BnaDGK1–1&BolDGK1–1,BnaDGK2–1&BolDGK2–2,BnaDGK7–1&BraDGK7–1), there were noK_a/K_svalues, because the two genes of each pair had exactly the same CDS. These results showed that mostDGKgenes have undergone purified selection, whereas the gene pairBnaDGK7–2&BolDGK7–1has been positively selected (Additional file8: Table S2).

In addition, the divergence time of the duplicated genes was estimated by calculatingK_svalues. Our estimation showed that the divergence time ranged from 12.11 to 16.75 million years ago (MYA) and averaged 14.43 MYA between theBrassicaspecies andA. thaliana(Additional file8: Table S2). This result indicated that the divergence ofBrassica DGKgenes occurred at ~ 14.43 MYA, in accordance with the WGT event betweenBrassicaspecies andA. thalianathought to have occurred approximately 9–15 MYA [32]。For the 6 paralogous gene pairs inB. napus, we usedK_svalues to estimate the time of the whole-genome duplication event. TheK_svalues ranged from 0.07 to 0.17, with an average of 0.12. The corresponding derived divergence time varied from 2.48 to 5.83 MYA, with an average of 3.84 MYA, which is considerably lower than the average value inB. napus, indicating that the divergence of theDGKgenes inB. napusoccurred well after theBrassicaWGT event (Additional file8: Table S2). This may be because the formation of theB. napusspecies and theBrassicaWGT event had little effect on the divergence of the syntenicDGKgenes inB. napus.

Gene structure, domain, and conserved motif analysis

We studied the intron–exon structures of theDGKgenes to determine the structural diversity ofDGKgenes in different clusters and to explore whether the gene structure changed during polyploidization (Fig.3). EachDGKcluster possesses a different gene structure. The structure ofDGKgenes from cluster 1 was significantly more conserved than the other two clusters, with an average of 7 exons. Of 16DGKgenes in cluster 1, 13 gene pairs had an identical gene structure with the correspondingA. thalianahomologs, with 7 exons, whereas 3 genes had 9 exons. By comparison, the number of exons varied from 9 to 12 and from 11 to 15 in clusters 2 and 3, respectively. In cluster 2, 12 out of 19DGKgenes had 10 exons. In cluster 3, out of 14DGKgenes, 6 genes were composed of 11 exons and 5 genes had 12 exons. In addition, we performed a comparative analysis of 17 gene pairs with the closest genetic distance in the phylogenetic tree, which we considered might have a direct evolutionary relationship. Of these 17 gene pairs, 16 gene pairs had an identical gene structure. Therefore, theDGKgenes shared a high similarity in exon number during polyploidization.

Furthermore, we analyzed the domain location of DGK proteins to identify the changes in the domain’s position between different clusters or different species (Fig.3). Although, all the DGK proteins assessed contained a catalytic domain (DGKc, PF00781) and an accessory domain (DGKa, PF00609), only cluster 1 had two C1_1 (PF00130) domains. Furthermore,BnaDGK2–1had no C1_1 domain and had the shortest sequence in cluster 1. The C1 (or DAG/PE binding) domain binds an important secondary messenger DAG, as well as the analogous phorbol esters (PE) [33]。Additionally, the positions and lengths of functional domains belonging to the same clusters were generally identical, although they varied between clusters. For instance, most DGKc domains of cluster 1 started from the 261st aa, while the DGKc domains of cluster 2 began at the 82nd aa and the DGKc domains of cluster 3 began at the 35th aa. The location of the DGKa domain was likewise different between the three clusters.

Next, we used the MEME online tool to search for 10 conserved motifs in DGK proteins (Fig.3and Additional file5: Fig. S5). We detected 10 motifs, among which 7 (motifs 1, 2, 3, 4, 5, 7, and 8) occurred in all DGK proteins. In general, DGK proteins in the same cluster displayed parallel motif components. Notably, motif 9 existed only in cluster 1 and constituted the second C1_1 domain ofDGKs,而主题10 existed in 28 DGK proteins, including 19DGKsin cluster 2 as well as theDGKshomologous toAtDGK1. Taken together,DGKgenes were highly conserved at the protein level amongBrassicaspecies andA. thaliana.

Expression profiles of DGK genes in different tissues of Brassica species

As we know, the gene functions were predicted by the expression patterns during plant development. To explore the expression patterns ofBnaDGKs在12迪斯汀,我们分析了他们的表达水平ct tissues at different developmental stages, including radicle, hypocotyl, cotyledon, root, stem, young leaf, mature leaf, petal, pistil, stamen, seed, and seed pericarp. All samples of 12 tissues at different developmental stages and time points inB. napusare detailed in Additional file8: Table S4. Based on the transcriptome sequencing datasets fromB. napusZS11 (the BioProject ID PRJNA358784), we found that 7BnaDGKswere highly expressed in distinct tissues and 7BnaDGKswere expressed in specific tissues, whereas 3BnaDGKshad low expression levels (Fig.4a). The remaining 4BnaDGKs(BnaDGK4andBnaDGK6) exhibited almost no expression in most tissues, except thatBnaDGK4was detected in Sta_i andBnaDGK6was detected in Se_10d (Fig.4a). Several genes were expressed in all tissues except LeY_f, Sta_i, and seed, includingBnaDGK1–1,BnaDGK1–2,BnaDGK1–3,BnaDGK3–3,BnaDGK3–4,BnaDGK5–1, andBnaDGK5–2. Notably,BnaDGK2–2andBnaDGK2–3had high expression in Sta_i (Fig.4a). Our results showed thatDGKshad relative expression levels in young tissues, in accordance with the previous results [18]。In addition, some genes displayed the tissue-specific expression profiles, for example,BnaDGK1–4in stem,BnaDGK7–1in leaves and silique pericarps, andBnaDGK7–2in silique pericarps, respectively. Overall,BnaDGKshad expression detected at different levels in various tissues.

To further study the function ofDGKs, we then analyzed the expression levels ofDGKsin six tissues of parental linesB. rapaandB. oleracea(Fig.4b, and c). We also observed that the orthologs members ofDGKsexhibited different tissue-specific expression patterns. For example, the orthologs ofDGK4andDGK6inB. rapaandB. oleraceadisplayed the same expression pattern as those ofB. napus, which had low or almost no expression in all tissues (Fig.4). This result suggests that they might be pseudogenes or expressed only at specific developmental stages or under special conditions. However, homologs ofDGK1,DGK3, andDGK5possessed different expression patterns among them. For example,BnaDGK5–1andBnaDGK5–2showed markedly higher expression thanBnaDGK5–3andBnaDGK5–4, as didBolDGK5–1andBraDGK5–1compared toBolDGK5–2andBraDGK5–2, respectively. The divergent expression patterns can be explained by the pseudogenization and functional divergence.

Cis-acting elements in the promoters ofDGKgenes

To explore the regulatory mechanism ofDGKgenes, we carried out an analysis of transcription cis-regulatory elements in the 2000-bp regions upstream ofDGKgene transcription start sites in threeBrassicaspecies. We identified and counted cis-acting elements associated with plant development and growth, abiotic and biotic stresses, phytohormone responses and light responsiveness in the promoters of allDGKgenes (Fig.5). Of 13 cis-acting elements related to plant development and growth, the as-1 motif (present in 83.3% of theBrassica DGKpromoters) is involved in the root-specific expression, while the GCN4/AACA motif (28.6%) is involved in endosperm expression (Additional file8: Table S3).

We also identified four types of abiotic stress elements in theDGKpromoters: defense and stress responsiveness (MYB/MYC/TC-rich repeats), drought responsiveness (MBS and DRE core), low-temperature responsiveness (LTR), and heat stress responsiveness (STRE) (Fig.5). Apart from this, two cis-elements related to wounding and pathogen response (W box and WUN motif) were also widely distributed in theDGKpromoters (Additional file8: Table S3). For phytohormone-response-related cis-acting regulatory elements, allDGKpromoters possessed abscisic-acid-responsive element (ABRE). Approximately 83.3 and 71.4% of theDGKpromoters contained methyljasmonate (Meja)-responsive elements (CGTCA motif and TGACG motif) and gibberellin-responsive elements (GARE-motif, TATC-box, and P-box), respectively, and approximately 66.7% possessed ethylene-responsive element (ERE). Auxin-response elements (AuxRR core and TGA elements) and the salicylic-acid-response element (TCA element) were also present in certainDGKgene promoters (Additional file8: Table S3). Notably,BraDGK1–1andBraDGK4–1had 11 ABREs and 9 EREs in its promoters, respectively. Moreover, light-responsive regulatory elements consisting of 24 types of different elements were predicted in theDGKpromoters (Additional file8: Table S3). These results confirm that theDGKgenes play a major role in stress resistance and hormone signaling pathways.

Expression analysis ofDGKgenes in cultivars with different oil contents

DGAT catalyzes the transfer of an acyl chain from a coenzyme A ester to the sn-3 position of sn-1, 2-diacylglycerol to form triacylglycerol [19], while DGKs can catalyze the conversion of DAG to PA. Meanwhile, DGAT and DGK could also compete for the substrate DAG, which is used for the synthesis of TAG and PA, respectively. However, TAG is the main lipid storage form in plants. Therefore, we analyzed the expression patterns of allBnaDGKs在样本两个品种与不同石油contents using qRT-PCR to explore the relationship betweenBnaDGKsand seed oil content. The results showed that mostBnaDGKgenes were differentially expressed in developmental seed and silique pericarp. Except forBnaDGK6–1had the highest relative expression in 20 DAF seeds, with a change of up to 33 fold, but it had almost no expression in silique pericarp (Fig.6). In addition, the relative expression of 8 of 21BnaDGKgenes in the seeds,BnaDGK2–2,BnaDGK2–3,BnaDGK3–1,BnaDGK3–2,BnaDGK3–3,BnaDGK3–4,BnaDGK5–2, andBnaDGK7–1, gradually decreased over time (Fig.6and Additional file6: Figure S6). The results suggest thatDGKgene expression functions in the early stages of seed development. Importantly, we found that the relative expression levels ofBnaDGK2–2,BnaDGK7–1, andBnaDGK7–2had the higher expression levels in development silique pericarp and seeds of high-oil cultivar than that in the low-oil cultivar, butBnaDGK5–3andBnaDGK5–4had conversely expression patterns (Fig.6), indicating that they might be involved in the accumulation of seed oil contents. Therefore, theseDGKscould be selected as excellent candidates for further functional characterization and application in rapeseed breeding programs.

Expression analysis ofBnaDGKsunder hormone treatment and heavy metal stress

Based on their expression levels as determined by RNA-seq and the closest genetic distance, five gene pairs were selected and their expression patterns were validated under hormone treatments (ABA, BR, and GA) with qRT-PCR. In general, most ofDGKsshowed a similar expression pattern during the hormone stresses that they showed the sharply up-regulation at 1 h and decreased to the relatively low expression between 3 h and 6 h, then their expression levels were gradually increased after 6 h under ABA and BR treatment (Fig.7), such asBnaDGK2–2,BnaDGK2–3, andBnaDGK3–3. These results suggested that the gene expression levels ofBnaDGKswere obviously induced by exogenous hormone and were affected by their own regulation.

For ABA treatment, expression levels ofBnaDGK1–1,BnaDGK2–2,BnaDGK2–3,BnaDGK3–3, andBnaDGK5–2increased to varying degrees at 1 h, by approximately 2.3-, 5.1-, 4.0-, 2.9-, and 1.5-fold, respectively (Fig.7). Among these,BnaDGK1–1andBnaDGK5–2possess most ABA response cis-elements, whereasBnaDGK2–2the middle andBnaDGK3–3the least (Fig.5). Furthermore,BnaDGK3–4andBnaDGK5–1expression profiles were significantly downregulated between 3 h and 12 h of ABA exposure, and that ofBnaDGK7–2was initially downregulated and then gradually increased to peak at 12 and 24 h (Fig.7). These results suggest that cluster 1 (BnaDGK1andBnaDGK2)BnaDGKsmay be responsive to ABA. We only found that the expression pattern ofBnaDGK1under the ABA treatment was perhaps correlated with the cis-element accumulation (Fig.5). In terms of BR induction, mostBnaDGKgenes displayed a similar expression pattern. As BR treatment continued,BnaDGKexpression levels decreased between 1 h and 6 h, and then increased from 6 h to 24 h. In particular, we observed that the expression levels ofBnaDGK2–2,BnaDGK2–3, andBnaDGK3–3obviously increased at 1 h, whereasBnaDGK7–1andBnaDGK7–2exhibited stronger increase at 24 h. The expression levels of mostBnaDGKswere higher at 24 h than at 0 h, indicating that the BR response mechanism is complex and the role ofDGKgenes is minimal and delayed (Fig.7). The GA treatment results revealed that except forBnaDGK1–3,BnaDGK3–4, andBnaDGK5–1, the remainingBnaDGKswere upregulated with various degrees, withBnaDGK5–2andBnaDGK7–2exhibiting 2.8-, and 4.0-fold upregulation, respectively. After the GA application,BnaDGK1–3expression gradually increased, peaking at 24 h. In addition,BnaDGK3–4andBnaDGK5–1were downregulated by GA treatment. We found in Fig.7thatBnaDGK5–1lacked GA cis-elements (Fig.5). Overall,BnaDGK2–2andBnaDGK3–3were upregulated, whileBnaDGK3–4andBnaDGK5–2were downregulated, in response to the various hormone treatments, implying that someBnaDGKgenes respond to hormone induction (Fig.7).

A previous study reported thatDGKgenes in land plants can respond to known beneficial elements as well as to other metal and metalloid ions [18]。We conductedBnaDGKexpression profiling based on RNA-seq data fromB. napuscultivars subjected to heavy metal induction (Additional file7: Fig. S7). Upregulation and downregulation were defined by a log₂ratio. The result showed that severalBnaDGKswere differentially regulated by As³⁺and Cd²⁺. Thus, induction with As³⁺resulted in lower expression ofBnaDGK1–4andBnaDGK3–2, but higher expression ofBnaDGK5–4,BnaDGK7–1, andBnaDGK7–2. Application of 30 mg/L Cd²⁺led to the upregulation ofBnaDGK1–1,BnaDGK2–2,BnaDGK7–2, andBnaDGK5–4, but downregulation ofBnaDGK1–4. Notably,BnaDGK1–4showed some degree of downregulation in response to all metal ion treatments in all cultivars. Our results suggested that someBnaDGKsmight take part in response to heavy metal stress and their expression might be induced by heavy metal treatments inB. napus.

Plants can respond rapidly to environmental stimuli and protect themselves from diverse stressors by activating a series of signal transduction cascades. PA, an emerging signal molecule, mediates signal pathways related to environmental stimuli, and is produced mainly by the PLD and coupled PLC/DGK routes in eukaryotic cells [2]。Moreover, DGK activity has been investigated in several plant species. However, few studies have examined the transcript levels and functions ofDGKgenes. In this work, we identified 21BnaDGKgenes inB. napus, 10BraDGKgenes inB. rapa,and 11BolDGKgenes inB. oleracea. Given that theA. thalianagenome encodes 7DGKgenes, it is clear that a WGT event has occurred inBrassicaspecies since they diverged fromA. thalianaand thatDGKgenes have been lost since the WGT event [32]。A report that 35% of genes were lost via a deletion mechanism after theArabidopsisandBrassicalineages diverged explains why at least 21DGKsare not observed in eitherB. rapaorB. oleracea[34]。We identified only 7 collinear gene pairs betweenA. thalianaandB. rapa, and 6 collinear gene pairs betweenA. thalianaandB. oleracea, perhaps as a result of gene identification or undergo the gene loss during polyploidy speciation. The chromosomal distribution analysis showed that mostDGKgenes in the Ar subgenome ofB. rapaand the Co subgenome ofB. oleraceamaintained their positions relative to the An and Cn subgenomes ofB. napus, respectively. Our synteny analysis betweenB. napusand its diploid progenitors indicated that mostDGKgenes inBrassicaare located in the syntenic regions, with the A subgenome sharing 16 gene pairs and the C subgenome sharing 15 gene pairs.

Polyploidy and WGT are prevalent in the evolutionary history of various species [35]。Most duplicated genes have arisen through whole-genome duplication, and these are often lost or nonfunctional inBrassicaspecies [36]。探索的功能和特点duplicate genes is important for understanding plant evolution. Our analysis of theK_a/K_svalues for 69 duplicate genes showed that except forBnaDGK7–2andBolDGK7–1, all theirK_a/K_svalues were < 1. This implies that these duplicated genes have undergone strong positive selection, which is consistent with previous research indicating that surviving duplicate genes have undergone strong purifying selection [36,37]。These results indicate that theDGKgene family has been well conserved during evolution over time.

Mammalian DGK enzymes are classified into five groups according to sequence homology, whereas plantDGKshave been divided into three clusters based on analyses in other plant species [38]。The catalytic region of all plantDGKsreported so far consists of a catalytic domain (DGKc, PF00781) followed by an accessory domain (DGKa, PF00609) near the C-terminus, where it may contribute to the normal function of DGKc. The cluster 1DGKsare more complex, with two copies of the DAG-binding domain (C1, PF00130), an upstream basic region, and an extended CRD in their N-termini. The DAG/PE-binding domain features the C₆/H₂cores. By contrast, cluster 2DGKsharbor only the DGKa and DGKc domains, whereas cluster 3DGKsmay display an extCRD-like domain generated by alternative splicing [2,39]。Moreover, since a GXGXXG consensus sequence in their catalytic domain (DGKc) is the ATP-binding motif,DGKscan use ATP as an energy generator to catalyze the conversion of DAG to PA. Interestingly,DGK1p, a novelDGKgene in the yeastSaccharomyces cerevisiae, utilizes CTP, rather than the standard ATP, as phosphate donor in forming phosphatidate [40]。The sequence ofS. cerevisiae DGK1p不像的DGKsfrom other species and contains a CTP transferase domain essential for the protein’sDGKactivity [40]。The combination ofDGK1pand phosphatidate phosphatase encoded byPAH1can regulate the levels of DAG and PA for phospholipid synthesis [41]。

Eukaryotic cells contain two categories ofDGKs. On the one hand, inactiveDGKs, stimulated by inositol phosphate metabolic signaling, can catalyze DAG to generate phosphoglycerols. On the other hand, activeDGKswith greater cellular activity can phosphorylate DAG to yield PA [42]。In eukaryotic species, PA is usually produced from the PLD pathway when plants are subjected to drought [43], oxidative stress, or physical damage [16], whereas PA accumulation via the PLC/DGK pathway is promoted by pathogens and xylanase [16]。Additionally, cold and salt stresses can lead to fast, transient production of PA through both pathways [44]。InA. thaliana, transcript expression profiling revealed thatAtDGK1andAtDGK2are implicated in the plant cold response [45]。Furthermore,AtDGK2expression can be transiently induced by the wounding signal [46]。A recent study showed thatAtDGK4is highly expressed in pollen, and a homozygousAtDGK4mutation affects not only the male germline but also the vegetative tissue [47]。过度的大米DGKgeneOsBIDK1can enhance disease resistance in transgenicN. tabacum[48], and suppressingOsDGKgene expression results in a distinct depletion of transcripts of the defense geneOsNPR1and the salt-responsive geneOsCIPK15[23]。In addition, transcript levels ofMdDGK4andMdDGK8in apple were induced by salt and drought stress, respectively [4]。In soybean,GmDGK10transcripts showed dramatic upregulation in response to PEG stress in root tissue, andGmDGK8andGmDGK9were significantly upregulated in the presence of NaCl [24]。Additionally, previous studies have indicated that PA is involved via the PLC/DGK route in various hormone pathways, such as the ethylene [8], ABA [10], and BR [17] pathways. Furthermore, theDGKgene response to beneficial elements and other ions has been demonstrated [18]。

Rapeseed is one of the most widely grown oil crops worldwide, and seed oil content is a crucial agronomic trait in rapeseed breeding. Therefore, exploring ways to increase the oil production ofB. napusis of great agricultural and economic significance. DAG can be catalyzed by DGAT to form TAG, a major contributor to vegetable oils in oilseed crops. Thus we infer that theDGKgenes may be candidate genes for improving seed oil content inB. napus. In this study, we investigated the effects ofBnaDGKgenes on seed oil content using qRT-PCR analysis. SeveralBnaDGKgenes were highly expressed in the sample from the high-oil-content cultivar, whereasBnaDGK5–3andBnaDGK5–4had distinctly high expression levels in the samples from the low-oil-content cultivar, indicating that differentially expressedBnaDGKgene family members should be further functional identified and applied in rapeseed breeding programs.

Nonetheless, transcription level and functional analysis of the response ofDGK激素海藻糖酶基因家族成员tments and metal stress are fragmentary. In terms of hormone treatments, PA interacts with ABI1 phosphatase 2C, thus promoting ABA signaling inA. thaliana[10]。InB. napusroot, exogenous application of EBL resulted in PA synthesis through the PLC/DGK pathway under optimal salinity conditions [17]。In our study, cluster 1BnaDGKgenes (BnaDGK1andBnaDGK2) may respond to ABA treatment.BnaDGK2–2,BnaDGK3–3,BnaDGK3–4, andBnaDGK5–2are potential candidate genes in response to hormone treatments, asBnaDGK2–2andBnaDGK3–3were upregulated, andBnaDGK3–4andBnaDGK5–2were downregulated in response to three types of hormone treatments.

Meanwhile, heavy metals, including copper, manganese, cobalt, zinc, and chromium could be toxic to plants, disrupting plant metabolic processes, inhibiting plant growth, or causing plant death [49,50]。Of these, excessive Cd can produce toxic effects in plants by reducing chlorophyll contents, inhibiting root growth, and perturbing water homeostasis [51]。However, root length, total biomass dry weight, and plant height ofLonicera japonicaplants increase in the presence of 0.5–5.0 mg/L Cd [52]。The expression levels of allAtDGKgenes exceptAtDGK6were repressed by Cd, whereasOsDGK8expression was slightly upregulated [18]。In our study,BnaDGK1–1,BnaDGK2–2, andBnaDGK5–4were upregulated, whileBnaDGK1–4was downregulated inB. napusexposed to 30 mg/L Cd²⁺. Our results thus contribute further evidence that Cd²⁺influencesDGKgene activities.

一项研究提供了直接证据As uptake by roots via high-affinityPHT1transporters inA. thalianaand rice [53]。Furthermore, As stimulates the uptake of Pi, providing some growth benefits at the modest concentrations [54]。As (in the form of NaHAsO4) upregulatesOsDGK2and downregulatesOsDGK3aandOsDGK8expression in rice [18]。In our study,BnaDGK1–4andBnaDGK3–2were downregulated, whereasBnaDGK5–4,BnaDGK7–1, andBnaDGK7–2were upregulated, in response to As³⁺treatment. These observations suggest a possible role ofDGKgenes in As³⁺absorption. The cis-acting elements were also predicted from the promoter regions ofDGKgenes, and a large number of abiotic and biotic stress elements in these promoters were identified, including MYB/MYC/TC-rich repeats, DRE, LTR, and the W box and WUN motif. Apart from these elements, mostDGKpromotors analyzed contain ABREs, Meja-responsive elements (CGTCA motif and TGACG motif), and gibberellin-responsive elements (GARE-motif, TATC-box, and P-box), as well as EREs.

Together, these results confirm thatDGKgenes play a major role in stress resistance and hormone signaling pathways. Our findings provide useful information on the evolutionary aspects of theDGKgene family inBrassicagenome and assist to reveal the biological functions ofDGKsin response to both oil metabolism and abiotic stress.

In this study, we provided a complete overview of theDGKgene family inB. napusand its diploid progenitors, and analyzed its evolutionary patterns, including gene duplication history and positive selection level. We identified 21, 10, and 11DGKsinB. napusand its diploid progenitorsB. rapaandB. oleracea, respectively. We also investigated their chromosome distribution, gene structure, domain, and motif pattern, as well as the presence of cis-regulatory elements. The expression patterns ofDGKgene family members in diverse tissues at different developmental stages showed that except for those homologous toAtDGK4andAtDGK6, theDGKgenes expressed in most tissues. Further, the expression patterns ofDGKgenes under hormone treatment and metal ion induction indicated that someDGKgenes respond to these treatments. Finally, severalBnaDGKsexpression profiles in cultivars with different seed oil contents were not completely consistent. Together, these results indicate thatDGKgenes have various roles in plant growth and development, hormone response, metal ion stress response, and the control of seed oil content. This work improves our understanding of the evolution of theDGKgene family and provides a reference for future studies.

Identification and nomenclature ofDGKsinB. napus,B. rapaandB. oleracea

Genomic sequences, coding sequences, and protein sequences ofA. thaliana,B. napus,B. rapa, andB. oleraceawere downloaded from the TAIR database (http://www.arabidopsis.org) and the Brassica Database (BRAD,http://brassicadb.org/brad). SevenAtDGKsequences fromA. thalianawere used as queries to identify the candidateDGKgenes inBrassicaspecies via a BLASTp search with a thresholde-value of 1e-20 [55]。Furthermore, the Pfam database (http://Pfam.sanger.ac.uk/) and the conserved domain database (CDD) of the US National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/Structure/cdd) were used to examine two conserved domains (DGKc, PF00781 and DGKa, PF00609) of the putativeDGKgenes. In this study, only proteins containing the complete DGKa and DGKc domains were considered asDGKgenes.

For nomenclature, a species prefix, such as ‘Bna’ forB. napus, was used, followed byDGKand two numbers. The first number refers to the homologousArabidopsisgene, and the second represents the degree of homology to the correspondingArabidopsisgene, where ‘-1’ represents the highest homology, followed by ‘-2’, and so on: for exampleBnaDGK1–1.

Multiple sequence alignments and protein properties of allDGKsinB. napus,B. rapa,B. oleracea, andA. thaliana

The multiple sequence alignments of theDGKproteins were established with the ClustalX [56] program, using the default parameter mode. Then, Jalview 2.11.0 [57] and DOG 2.0 [58] software were used to visualize the multiple sequence alignment. In addition, molecular weight (kDa) and isoelectric point (pI) values of eachDGKprotein sequence were predicted using the ExPASy server (http://expasy.org).

Chromosomal distribution, phylogenetic tree, gene duplication, and syntenic analysis of theDGKsinB. napus,B. rapaandB. oleracea

The chromosome location information for theDGKgenes ofB. napusand its two diploid progenitors was collected from the BRAD database, and the chromosome distribution was plotted with MapChart 2.0 [59]。To comprehensively understand the evolutionary relationship of theDGKgene family members, a phylogenetic tree was constructed with MEGA 7.0.26 using the neighbor-joining (NJ) method [60] among four species (A. thaliana,B. napus,B. rapa, andB. oleracea). A series of parameters were as follows: the Poisson model, pairwise deletion, and conserved sequences with coverage of 100% and 1000 bootstrap replicates. Finally, the phylogenetic tree was visualized using the EvolView website (http://www.evolgenius.info/evolview).

To identify forms of gene duplication ofBnaDGKs, only a total of 101,040B. napusannotated gene sequences were aligned using BLASTp, with ane-value of 1e-10. Then, the MCScanX software [61] with default parameters was used to classify the duplication patterns of theDGKgenes, and TBtools software [62] was used to calculate the synonymous (K_s) and nonsynonymous (K_a) mutation rates of the duplicatedDGKgene pairs. Divergence time was inferred using the formula T = K_s/2R, whereK_srepresents the synonymous substitutions per site and R is the rate of divergence. For dicotyledonous plants, the assumption R is 1.5 × 10^− 8synonymous substitutions per site per year [63]。For the syntenic analysis ofDGKs, we used Multiple collinear scanning toolkits (MCScanX) to produce the collinearity blocks across the whole genome with the default parameters [61]。Then, the collinearity gene pairs of theDGKfamily were visualized by Advanced Circos programs of TBtools [64]。

Protein domain, motif, and gene structure analysis of allDGKsinB. napus,B. rapa,B. oleracea, andA. thaliana

To identify conserved domains ofDGKs, the NCBI conserved domain search (www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) was used and the results were confirmed by conducting a Pfam database search. A search for conserved motifs ofDGKproteins was performed with MEME 5.0.1 (http://meme-suite.org/) using default settings, except that the number of motifs was set to 10. Subsequently, the Gene Structure Display Server (GSDS 2.0,http://gsds.cbi.pku.edu.cn/) was used to conduct an exon–intron structure analysis ofDGKgenes. An integrated schematic including the phylogenetic tree, gene structure, motif, and conserved domain was visualized using EvolView.

Promoter elements analysis and subcellular localization of allDGKsinB. napus,B. rapa,B. oleracea, andA. thaliana

A 2000 bp region upstream of the translation start sites of eachDGKgene was acquired from Brassica database (BRAD) as a promoter sequence, and the cis-acting elements were analyzed using the PlantCARE website (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). The subcellular localization was predicted using the PSORT online website (http://psort1.hgc.jp/form.html).

Expression analysis of allDGKsin different tissues ofB. napus,B. rapa, andB. oleracea

First, total RNA was isolated using the EZ-10 DNAaway RNA Mini-prep Kit (Sangon Biotech Co., Ltd., Shanghai, China) for the rapeseed cultivar Zhongshuang 11 (ZS11) in distinct tissues at different developmental stages. The cDNA library construction and RNA sequencing were performed as described previously [65] at Novogene Bioinformatics Technology (Beijing, China) and were deposited in the BioProject database (BioProject ID PRJNA358784). Then the relative expression values in the bars were calculated based on FPKM values (fragments per kilobase of exon model per million) using Cufflinks with default parameters [66], and expression patterns of the candidateBnaDGKswere analyzed at 12 tissues (radicle, hypocotyl, cotyledon, root, stem, young leaf, mature leaf, petal, pistil, stamen, seed, and seed pericarp; Additional file8: Table S4 and Table S5). Furthermore, six tissues (callus, flower, silique, stem, leaf, and root) RNA-seq data of cultivatedB. oleraceaandB. rapawere obtained from the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/browse/;accession numbers GSE42891 and GSE43245). AllBnaDGKs,BraDGKs, andBolDGKsexpression levels were also estimated using FPKM (Additional file8: Table S5 and Table S6).

In addition, transcriptome sequencing datasets of five rapeseed cultivars under heavy metal stress were generated and analyzed as described above, which were treated cultivars P070, P085, and P087 with 15 mg/L As³⁺and cultivars P085, P155, and P163 with 30 mg/L Cd²⁺, respectively (Additional file8: Table S7). Gene expression levels ofBnaDGKswere estimated using FPKM values. Heatmaps of allDGKswere drafted using TBtools [64], which was normalized by Log₂(FPKM values).

Plant materials, hormone treatments, and heavy metal stress

TwoB. napuscultivars were used in this study, P594 (high oil content, 43%) and P134 (low oil content, 34%), both provided by the Chongqing Rapeseed Engineering Technology Research Center. Both cultivars were grown under natural conditions in Chongqing, China, and inflorescences were bagged before blossoming to prevent pollen contamination. Finally, developing seeds and silique pericarp of two cultivars were collected at 20, 30, and 40 days after flowering (DAF) and preserved at − 80 °C for further analysis.

For hormone treatment, ZS11 seeds were sown and grown in an artificial climatic chamber at 25 °C with a 16 /8 h (day/night) photoperiod. Two-week-old ZS11 seedlings were treated with 50 μM ABA, GA, or BR for 0, 1, 3, 6, 12, and 24 h, respectively, and their leaves were immediately collected. All materials were frozen in liquid nitrogen and stored at − 80 °C until use.

Five rapeseed cultivars (P070, P085, P087, P155, and P163) with extreme phenotypes for heavy metal response (strong of weak resistance) were selected from two hundred rapeseed accessions under different concentrations of two heavy metals (As³⁺and Cd²⁺). Subsequently, we treated rapeseed cultivars (P070, P085, and P087) with 15 mg/L As³⁺and rapeseed cultivars (P085, P155, and P163) with 30 mg/L Cd²⁺as the optimal concentration in this study. In brief, the filled seeds for each accession were sown on the filter papers in the petri dishes (9 × 9 cm), which were treated with distilled water (Control) and 15 mg/L As³⁺and 30 mg/L Cd²⁺(Treatment), respectively. The growth conditions were kept at 25 °C with long-day conditions (16 h light/8 h dark, 5000 Lux). After 7 days, whole plants were frozen in liquid nitrogen and stored at − 80 °C until use.

RNA extraction, complementary DNA synthesis, and quantitative real-time PCR analysis ofBnaDGKgenes

Total RNA was extracted with the EZ-10 DNAaway RNA Mini-prep Kit (Sangon Biotech Co., Ltd., Shanghai, China) according to the manufacturer’s instructions. Subsequently, the gel electrophoresis and a NanoDrop 2000 spectrophotometer were used for detecting the quality and concentration of each RNA sample, and the qualified RNA were used for further analysis. The first-strand cDNA was synthesized from 1 μg RNA by reverse transcription PCR (RT-PCR) with a Perfect Real-Time Synthesis Kit (TaKaRa Biotechnology, Dalian, China). The diluted cDNA after reverse transcription was used as the template for real-time quantitative RT-PCR. Real-time quantitative PCR was conducted with a BIO-RAD CFX96 qPCR instrument and SYBR II (TaKaRa). Each 20 μl PCR mixture that included 10 μl of SYBR® Premix Ex Taq™ II, 2 μl (100 ng) of template cDNA, and 0.4 μM of each PCR primer. The RT-PCR protocol was set to 95 °C for 30 s and 35 cycles of 95 °C for 5 s, followed by 56–60 °C (depending on the primers used) for 30 s. All samples were amplified in three biological replicates and three technical replicates. The 2^−ΔΔCtvalue was used to measure the relative expression levels ofBnDGKsunder hormone treatment [67]; and the 2^−ΔCtvalue was used to measure the relative expression levels ofBnDGKsbetween low and high oil content rapeseed [68]。B. napus ACTIN7(BnaACTIN7, GenBank Accession No. AF024716) was used as the housekeeping gene. All qRT-PCR primers are listed in Additional file8: Table S8.

Statistical analysis

The expression levels were defined as mean ± standard error of mean (SEM), and three biological replicates were performed in each experiment. Data was subjected to a one-way ANOVA using SPSS 15.0 (SPSS, Inc., Chicago, Ill) to define the significance differences between the mean values. Differences withpvalues of ≤0.05 and ≤ 0.01 were considered significant and highly significant, respectively.

RNA-seq ofB. napusvariety Zhongshuang 11 (ZS11) in distinct tissues at different developmental stages are available in the NCBI Sequence Read Archive (SRA) database under the accession number PRJNA358784. RNA-seq data for expression profiles ofB. rapaandB. oleraceain different tissues were acquired from NCBI gene expression omnibus (GEO) database (GSE43245 and GSE42891). All other datasets supporting the results of this article are included within the article and its Additional files.

DGK:

Diacylglycerol kinase

DAG:

Diacylglycerol

PA:

Phosphatidic acid

DGKc:

Diacylglycerol kinase catalytic

DGKa:

Diacylglycerol kinase accessory

WGT:

Whole-genome triplication

qRT-PCR:

Quantitative real-time polymerase chain reaction

ABA:

Abscisic acid

BR:

Brassinolide

GA:

Gibberellic acid

PLC:

Phospholipase C

PLD:

Phospholipase D

PC:

Phosphatidylcholine

EBL:

24-epibrassinolide

As:

Arsenic

Cd:

Cadmium

DGAT:

Diacylglycerol transferase

TAG:

Triacylglycerol

DAG/PE:

Diacylglycerol/phorbol ester

Ka:

Nonsynonymous

Ks:

Synonymous

MYA:

Million years ago

ABRE:

Abscisic-acid-responsive element

ERE:

Ethylene-responsive element

DAF:

Days after flowering

NJ:

Neighbor-joining

FPKM:

Framents per kilobase of exon per million mapped framents

Not applicable.

This research was funded by National Key R&D Program of China (2018YFD0100500), the National Natural Science Foundation of China (31830067), Fundamental Research Funds for the Central Universities (XDJK2020B030, XDJK2019C099), and the “111” Project (B12006), Modern Agro-industry Technology Research System (CARS-13).

Affiliations

1. Chongqing Engineering Research Center for Rapeseed, College of Agronomy and Biotechnology, Southwest University, Chongqing, 400716, China

   Fang Tang, Zhongchun Xiao, Fujun Sun, Shulin Shen, Si Chen, Rui Chen, Meichen Zhu, Qianwei Zhang, Hai Du, Kun Lu, Jiana Li & Cunmin Qu

2. Academy of Agricultural Sciences, Southwest University, Chongqing, 400715, China

   Fang Tang, Zhongchun Xiao, Fujun Sun, Shulin Shen, Si Chen, Rui Chen, Meichen Zhu, Qianwei Zhang, Hai Du, Kun Lu, Jiana Li & Cunmin Qu

Contributions

CMQ and JNL conceived and guided the project. FT and ZCX conducted the experiments. FJS, SLS, and SC collected and analyzed the data; RC, MCZ, and QWZ carried out the experiments and performed software; HD and KL analyzed the transcriptome data; FT, KL, and CMQ wrote the manuscript; JNL and CMQ reviewed the manuscript. All authors have read and approved the final manuscript.

Corresponding authors

Correspondence toJiana LiorCunmin Qu.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visithttp://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions