Comparative Analysis of Phospholipase D (PLD) Gene Family in Camelina sativa and Brassica napus and Its Responses in Camelina Seedlings under Salt Stress

1. Introduction

Plants exploit lipid signaling to regulate a variety of cellular processes in response to external stimuli [1]. In this regard, phospholipids which are the main components of biological membranes can play as intracellular signaling molecules [2,3]. As part of maintaining membrane stability and homeostasis as well as biotic and abiotic stress alleviation mechanisms, these molecules are released directly and transiently by the activity of phospholipases or lipid kinases [4,5]. Mediating different types of catalytic reactions, plant phospholipases are divided into A, C, and D classes. In both prokaryotes and eukaryotes, phospholipases D (PLDs) are the main enzymes that break the terminal phosphodiester bonds of phospholipids in plasmalemma and organellar membranes yielding phosphatidic acid (PA) and free soluble head groups (e.g., ethanolamine or choline) [6,7]. It has been shown that PLDs can modulate various process in plants at cell scale including programmed cell dead [8], interaction with cytoskeleton [9], pollen tube growth [10], vesicular transport [11], and stomata closure [12] or at whole plant level including seed germination [13], wounding [14], and hair root growth and patterning [15]. There are three PLDs subfamilies including SP-PLD, PX/PH-PLD and C2-PLD, in which SP-PLD subfamily is classified into α, β, γ, δ, ε, and ζ classes on the basis of their physicochemical catalytic, and structural features and sequence homology [2,6]. While each PLD appears to have a different physiological function, all classes have a phospholipid-binding region at their N-terminus, and two highly conserved C-terminus domains that are catalytically involved in lipase activity through the interactions of two catalytic HxKxxxxD (HKD) motifs [16,17].

PLD enzyme is found in both plant and animal cells, and its activity is regulated by various signaling pathways and cellular factors. Several PLD cDNAs have been cloned from higher plant species since the cloning of the first eukaryotic PLD cDNA from castor bean (Ricinus sativus L.) [18]. The members of the PLD gene family can now be studied and discussed at the genomic level in many plant species thanks to advances in genome sequencing technologies. Du et al. (2013) studied phospholipase D gene family in three Rosaceae species and identified three pairs of segmental PLD gene duplications in Malus × domestica and Prunus mume and four pairs in Fragaria vesca L. [19]. Out of 40 identified PLD genes, 33 genes represented conserved synteny and the authors suggested that PLD genes in Rosacea family could have a common ancestor and purifying selection could be the primary driving force behind PLD gene evolution. Using genome-wide analysis, 17 PLDs were discovered in rice across four chromosomes. In addition to the C2-PLDs and PXPH-PLDs, the SP-PLDs has also bene identified as a novel PLD subfamily in rice (which contains a signal peptide sequence rather than the PXPH or C2 domains) based on the protein domain structural analyses [13]. Liu et al. (2010) identified 18 and 11 PLD members in grape and poplar, respectively. PLD genes were clustered into φ, ζ, ε, δ, β/γ, and α subgroups where they have evolved from four preliminary ancestors as a results of gene duplications [17]. They also found that most of PLD gene family members in grape (10 out of 11) and poplar (13 out of 17) showed closer phylogenetically relationship together than with rice and Arabidopsis PLD gene members. Genome-wide in silico analysis of PLD genes from Corchorus olitorius and C. capsularis resulted in identification of 20 different motifs where β1 and γ1, and all members of δ1 classes in both Corchorus species represented the highest number of motifs [1]. Tissue-specific analysis of δ2, α2, and β1 subclasses revealed their roles involved in drought and salinity stresses. Genome-wide analysis of PLD gene family members in tree cotton (Gossypium arboretum) revealed 19 non-redundant genes with similar exon-intron architecture and motifs with highly conserved structures. There are 59 PLD gene members in alfalfa (Medicago sativa L.) grouped into ε, δ, γ, β, α, and ζ subclasses [3]. It has been demonstrated that cis-acting elements are different in PLD members within common subfamilies, but they represent similar physicochemical characteristics, sequence feature and domains. The PLD protein family has, however, not been bioinformatically analyzed in important crops within Brassicaceae family, including rapeseed (Barssica rapa L.) and camelina (Camelina sativa L.) and to date, no genome-wide comparative analysis of PLD gene family has been carried out in these oilseed crops. In the present study, we identified and compared PLD gene members from the whole genome of Brassica rapa and Camelina sativa. The available bioinformatics tools were applied to analyze various characteristics of PLD genes and their expression levels were analyzed in response to salinity stress. Our findings provide a better understanding of the evolution and structure of sequence and function of PLDs, which will be useful in future functional studies.

2. Materials and Methods

2.1. Identification and Characterization of PLD Family Members

In order to identify the members of the PLD gene family, the known sequence of these genes in Arabidopsis was run as queries against the genome of Brassica napus and Camelina sativa, using BLASTp tool of Ensembl Plants database [20]. The peptide sequences of the PLDs in Brassica napus (BnPLDs) and Camelina sativa (CsPLDs) were checked by CCD database [21] for confirmation. In addition, PLDs were identified in Glycine max and Oryza sativa similar to the previous method. In this study, the physicochemical properties, including number of amino acids, isoelectric point (pI), grand average of hydropathicity (GRAVY), molecular weight (MW), and instability index were enumerated for the candidate PLDs, using the ProtParam tools [22].

2.2. Phylogenetic Analysis of PLDs and Prediction of Phosphorylation Sites

To investigate the evolution of the PLD gene family, peptide sequences of CsPLDs and BnPLDs with the PLD characterized in Oryza sativa, Glycine max, and Arabidopsis thaliana were aligned using online tool, the Clustal Omega [23]. In the next step, the aligned sequences were introduced to IQ-TREE web server [24] to construct a phylogeny tree based on the Maximum likelihood (ML) method with 1000 bootstrap. Finally, the phylogeny tree of PLDs was visualized by Interactive Tree Of Life (iTOL) v6.3 [25]. In addition, the strong potential phosphorylation sites (potential value ≥ 0.90) of CsPLDs and BnPLDs were identified using the NetPhos 3.1 server [26].

2.3. Duplication Analysis of CsPLD and BnPLD Genes

In the present study, the coding sequence (CDS) of pairs of CsPLD and BnPLD genes were aligned and of pairs with an identity ≥0.85 were selected as a duplicated genes [27,28]. In addition, the synonymous (Ks) and non-synonymous (Ka) indexes were calculated for duplicated pairs using TBtools software [29]. Besides, the time of divergence for each duplicated pairs was estimated using the following equation, T = (Ks/2λ) × 10−6. (λ = 6.5 × 10−9).

2.4. Conserved Motifs and Promoter Analysis

To recognize the conserved protein motifs, CsPLDs and BnPLDs were analyzed using the MEME (version 5.0.5) [30], based on the default settings. Moreover, an upstream region (1500 bp) of the start codon in each CsPLD and BnPLD was screened, as a promoter site, to identify the putative cis- regulatory elements using the PlantCARE database [31].

2.5. Interaction Networks and Structure Analysis of CsPLD and BnPLD

In the current study, the interaction networks for CsPLDs and BnPLDs were constructed using the STRING database (v11.5) [32]. To better understand direct and indirect interactions, the first and second shell were adjusted to ≤ 20 and ≤ 5, respectively. Finally, the interaction networks were visualized using Cytoscape v3.8.2 [33]. In order to draw the three-dimensional structure, five CsPLDs and BnPLD proteins (one candidate from each subfamily) were selected and their structure of these proteins were predicted using the Phyre2 server [34].

2.6. Expression Analysis of CsPLD and BnPLD in Different Tissues

In the current study, the available transcriptome data (RNA-seq information) related to different stages of growth and development for B. napus and C. sativa were used to screen the expression profile of the BnPLD and CsPLD genes. The Camelina eFP Browser (https://bar.utoronto.ca/efp_camelina/cgi-bin/efpWeb.cgi) was used for CsPLD genes and the BrassicaEDB was used for BnPLD genes to extract expression data in terms of FPKM (fragments per kilobase of exon model per million mapped reads) value. Finally, data were illustrated in heatmaps based on the log2 transformed method using the TBtools software.

2.7. Plant Materials and Treatment

In this study, expression patterns of candidate CsPLD genes were evaluated in camelina seedlings under salt stress. In the first step, sterilized seeds were planted in pots containing peat moss and perlite (2:1) and were kept under a temperature of 24±2 °C and 16 hours of light and pots were irrigated every three days. Then, after five weeks, seedlings were irrigated with salt (200 mM of NaCl) and was repeated after 24 hours. Finally, leaves were harvested at different time points (6, 24, and 72 h after treatment). The experiment was carried out in three biological replicates.

2.8. Expression Analysis Using qPCR

The total RNA samples from leaves were extracted by an RNX kit (Sinaclon, Iran) according to manufacturer protocols and the cDNA was produced by a reverse transcriptase kit (Roche, Germany), based on manufacturer protocols. In this study, six CsPLD genes were selected for qPCR analysis based on phylogenetic analysis, one member from delta, beta, zeta, and gamma subfamily and two members from alpha subfamily. Besides, actin-2 gene (Csa15g026420) was selected as a housekeeping gene for normalizing the raw expression data. The Primer3 online tool (version 4.1.0) [35] was applied to design the specific primers of candidate CsPLDs (Table S1). Maxima SYBR Green/ROX qPCR Master Mix kit (Thermo Fisher, France) was used to analyze the expression levels of the CsPLDs by the ABI Step One, according to manufacturer protocols. Finally, the raw data were analyzed by the delta delta ct method [36].

3. Results

3.1. Physiochemical Properties of PLDs

In this study, 33 CsPLDs and 41 BnPLDs were recognized in genome of Camelina sativa and Brassica napus, respectively (Table S2). PLDs were compared in five plant species that diverse physiochemical properties were observed between PLDs (Table 1). For instance, 33% of CsPLDs were predicted as stable protein, while 61% of BnPLDs were stable. All studied PLDs showed a negative GRAVY value, indicating that PLDs are probably hydrophilic proteins. Based on the exon number, less exon number was recognized in alpha subfamily members.

3.2. Phylogenetic Analysis and Prediction of Post-Translational Modifications of PLD Gene Family

Phylogenetic analysis of PLD gene family in five plant species, including C. sativa, Arabidopsis, B. napus, O. sativa, and G. max could separate the five PLD subfamilies, including gamma, delta, beta, alpha, and zeta (Figure 1). According to phylogeny, subfamilies gamma, delta, and beta are more closely related. In addition, the most of PLDs were in alpha subfamily group. CsPLDs showed more close relationships to Arabidopsis PLDs than oilseed crops. The potential phosphorylation modifications, as an important post-translational modification, were predicted in PLD proteins (Figure 1). Alpha subfamily members showed less potential sites than other subfamilies. It seems that alpha proteins are less influenced by post-translational modifications.

3.3. Conserved Motifs in PLDs

In this study to investigate the structure of PLDs, the conserved motifs of the CsPLDs and BnPLDs were identified. A total of 12 motifs were detected that motifs 1, 6, 10, 5, 4, 3, and 2 were observed in almost of PLDs (Figure 2). However, the members of zeta subfamily were different from other PLD subfamilies based on the distribution of conserved motifs. For instance, motifs 8, 9, 11, and 12 were observed in subfamilies alpha, beta, delta, and gamma. Structural differences between PLDs reveal the functional specificity between different subfamilies. On the other hand, the PLD members that have the same pattern of motif distribution probably had a similar evolutionary process and have a common ancestor.

3.4. Duplication Events in PLD Gene Family

The CsPLDs and BnPLD genes were unevenly located on chromosomes of C. sativa and B. napus (Figure 3 and Table S3). For the genome of C. sativa, the highest number of CsPLDs were distributed on chromosomes 1, 8, 13, 15, and 19 (Figure 3a). BnPLDs were mainly distributed on chromosomes A03, A05, A06, C03, and C07 (Figure 3b). However, the chromosomal location for 9 genes from B. napus is unknown, presented with An and Cn in the Figure 3b. Duplication events between PLDs revealed that many duplication events have occurred in PLD gene family from C. sativa and B. napus (Figure 3 and Table S3). Besides, the results disclosed that all CsPLD genes have been segmentally duplicated, while in BnPLDs, three tandem duplications were observed. According to ka/ks index, selective pressure of all studied PLDs was purifying. Divergence time of duplication events estimated ranging 1.38 to 8.21 and 1.56 to 219.36 MYA in CsPLDs and BnPLDs, respectively.

3.5. Promoter Analysis

The upstream of CsPLDs and BnPLDs were screened to identify the putative cis-regulatory elements having roles in stress and defense responsiveness, hormone responsiveness, protein binding site, and plant growth and light responsiveness (Figure 4a,b; Table S4). Cis-regulatory elements related to hormones were more observed in promoter region of CsPLDs, while in BnPLDs, cis-regulatory elements (about 44%) had roles in stress responses. The promoter elements related to ABA and auxin responsiveness were more distributed in CsPLDs, while elements linked with MeJA responsiveness were more observed in promoter region of BnPLDs. This is suggestive of the participation of CsPLDs in response under abiotic stress and BnPLDs are associated with defense responses against biotic stresses by inducing the MeJA mechanism.

3.6. Interaction Network of CsPLDs and BnPLDS

In order to predict the interactions of PLDs with other proteins, an interaction network for PLDs in camelina and B. napus was constructed based on the available data (Figure 5a,b). As presented in Figure a, CsPLDs showed interactions with non-specific phospholipase C (nsPLC), guanine nucleotide-binding protein (GNBP), 1-acyl-sn-glycerol-3-phosphate acyltransferase (AGPAT), diacylglycerol kinase (DAGK), and phosphatase 2C. Besides, BnPLDs interacted with nsPLC, phosphatidic acid phosphohydrolase 2 (PAH2), AGPAT, lipid phosphate phosphatase 2 (LPP2), lysophosphatidyl acyltransferase 2 (LPAT2), monogalactosyldiacylglycerol synthase type C (MGD3), phosphatidic acid phosphatase (PAP2), alpha/beta-Hydrolases, and phosphatidyl-N-methylethanolamine (Figure 3b). PLDs strongly showed interactions with proteins that are involved in phosphonate and phosphinate metabolism, ether lipid metabolism, phosphatidic acid metabolic, and inositol phosphate metabolism. These results disclose that PLDs can be involved in various pathways in response to upstream signals. However, interaction network of CsPLDs and BnPLDs was somewhat different.

3.7. 3D Structure and Pocket Sites of Candidate PLD Subfamilies

Three dimensional structure of each subfamily of PLD gene family was predicted based on homology. Results revealed that PLDs in C. sativa and B. napus have similar 3D structure, while they showed different binding points according to the location and type of sequence in the pocket region (Figure 6). Based on the abundance of amino acids in the pocket region, PRO, TRP, LEU, and GLU were more present in zeta subfamily of BnPLD, while TRP and GLU more observed in CsPLD zeta. Moreover, LEU, and GLU in BnPLD alpha and GLY in CsPLD alpha were more present in pocket site. Interestingly, ALA was frequently distributed in pocket site of beta subfamily of both BnPLD and CsPLD. In candidate gamma subfamily proteins, LYS, GLN, and ALA in BnPLD and SER, and ASN in CsPLD were more observed. LYS and THR were highly distributed in BnPLD delta and CsPLD delta, respectively.

3.8. Expression Profile of PLDs in Different Tissues and Organs

Expression levels of PLDs were investigated in different tissues and organs of C. sativa and B. napus (Figure 7 and Figure 8). In C. sativa, three CsPLD genes from alpha subfamily including Csa15g020390, Csa19g022540, and Csa01g018370, showed high expression in all tissues and organs (Figure 7). Besides, two genes, Csa13g049040 and Csa02g016220, from gamma subfamily were more expressed in senescing leaf tissues. In B. napus, five genes from alpha subfamily, including BnaC01g35830D, BnaA01g28530D, BnaA01g28540D, BnaC05g37540D, and BnaA05g23740D, showed high expression in all studied tissues and organs, except for anther tissues (Figure 8). A delta subfamily gene, BnaA08g04540D, and two beta subfamily genes, BnaC05g00230D and BnaA10g00160D, were more expressed in anther. In bud tissues of B. napus, two genes from alpha subfamily, BnaC03g69220D and BnaA05g35700D, were more induced.

3.9. Expression Profile of CsPLD in Response to Salt Stress

To better understand the potential roles of the PLDs in camelina, the expression profile of six candidate CsPLD genes were investigated in response to 200 mM of NaCl, as salinity stress. The results indicate that CsPLD candidate genes are involved in the response of Camelina seedlings to salt stress, although they had different expression patterns (Figure 9). PLD beta1 (Csa06g048050) and PLD delta (Csa10g006270) showed a similar expression pattern, so that they their expression levels were upregulated after 24 hours of salt stress and showed the highest level of expression after 72 hours of salt stress. The expression pattern of one PLD gamma gene (Csa13g049040) was significantly down-regulated after 6 hours and its expression level was not significantly induced by salinity. Besides, the expression level of Csa01g006520 (which is a PLD zeta) was reduced early stages of salt stress (after 6 and 24 hours), while its expression was increased after 72 hours. Two members of the alpha subfamily presented a different expression pattern after 6 hours of salt stress, alpha 4 (Csa03g061750) increased expression and alpha 1 (Csa15g020390) decreased expression, while their expression levels of both genes were increased after 72 hours.

4. Discussion

Oilseed plants such as rapeseed play an important role in the production of healthy oils in the food basket. Camelina, as an almost new agricultural plant, also has a good percentage of healthy oil and also has high resistance to adverse environmental conditions. Besides this, phospholipases have a key role in regulating cellular signaling pathways in response to stresses. PLD, as a key element linked in lipid mediated signaling, which involved in various physiological and developmental processes [5,37]. In this study, two oilseed plants were compared based on the PLD gene family to gain a better understanding of their structures and functions. The results showed that PLD family has 41 members in B. napus and 33 members in C. sativa, which were more members compared to other plants such as rice, 17 PLDs, [13] and Arabidopsis, 12 PLDs [38], 13 in maize [39] and 19 in potato [37]. Ploidy level and genome size can have an effect on the number of gene family members in plant species [40]. Evolutionary analysis was able to separate the subfamilies (gamma, delta, beta, alpha, and zeta) well. These findings suggest that PLDs in the same group share conserved sequences, disclosing that they can have similar functions. However, less relationship was observed between CsPLDs and BnPLDs. These results indicate that the evolutionary process of PLD gene family has been somewhat different in C. sativa and B. napus. Moreover, it seems that the ancestral gene of PLD family were formed before C. sativa and B. napus, and after the separation of the species, a different evolutionary process has been passed in each species [41]. In addition, the Ka/Ks value of all duplicated CsPLD and BnPLD pairs was estimated to be less than 1.0, revealing that the CsPLD and BnPLD family members survived under strong selection pressure of purification [42,43]. However, it was stated that purifying selection can generate pseudogenes and/or genes with conserved functions [44,45]. Besides, it seems that the expansion in CsPLDs happened much earlier than BnPLDs, and the evolution process of PLD was different in CsPLD and BnPLD family.

Signaling of phospholipases activates pathways related to growth and development and response to stresses [12,46]. Direct and indirect interactions can be effective in expanding the effect of phospholipases. By examining the interaction network, it was determined that the members of CsPLD and BnPLD family interact more with other elements of lipid-dependent cell signaling such as non-specific phospholipase, diacylglycerol kinase, 1-acyl-sn-glycerol-3-phosphate acyltransferase, phosphatidic acid phosphohydrolase, lipid phosphate phosphatase, lysophosphatidyl acyltransferase, and phosphatase. However, the interaction network of BnPLDs was estimated to be different from CsPLDs. It was previously reported that the interaction of PLD alpha1 with GPA is effective in increasing resistance to osmotic stress by affecting the opening and closing of stomata through ABA signaling pathway [47]. Moreover, it was stated that PLDs can participate in the response to abiotic and biotic stresses by cooperating with abscisic acid and jasmonic acid-dependent pathways [13,14,48]. The presence of cis-regulatory elements in the upstream region of PLD genes can strengthen this hypothesis, although more additional studies are needed. Post-translational modifications such as phosphorylation are effective in activating or deactivating target proteins [51]. The evaluation of post-translational changes in BnPLDs and CsPLDs showed that they have different potentials and thus can contribute in different ways. For example, less phosphorylation site was predicted in PLD alpha proteins, suggesting that this subfamily is less affected by upstream signaling elements.

The results of the expression profile of CsPLD and BnPLD genes in different tissues indicated that the members of this gene family are expressed in all tissues and developmental stages. A comparison between subfamilies revealed that the alpha subfamily has a greater range of activity. Also, the analysis of the structure of CsPLD and BNPLD genes disclosed that the members of alpha subfamily have fewer introns/exons. It has been reported that genes with fewer exons/introns are induced faster [45,52]. It seems that the smaller number of introns has an effect on the duration of the mRNA splicing process before leaving the nucleus, and as a result, they are translated earlier [43,53]. The quantitative findings under salt stress revealed that candidate CsPLDs were significantly induced. It was reported that PLDs are affected by abiotic stresses such as drought, high and low temperature, salt stress [37,46,54,55]. CsPLD gamma1, CsPLD zeta2, CsPLD alpha1, CsPLD alpha4 are recognized as early elements in response to salt stress. With the persistence of salinity stress, all studied CsPLD genes, except for CsPLD gamma1 gene, showed upregulation, suggesting that CsPLDs are probably involved in both the initial responses to activate the downstream pathways and in the resistance response to reduce the damages of the salinity stress.

5. Conclusions

Phospholipases D (PLD) is an enzyme that plays a crucial role in cell signaling, lipid metabolism, and membrane dynamics. In the current study, PLD gene family is characterized in Brassica napus and Camelina sativa. CsPLD and BnPLD have had different evolutionary processes, although both of them have been under almost similar evolutionary pressures, which has caused their structure to be conserved. Among the PLD subfamilies, members of the alpha subfamily have more obvious differences from the others, which can be considered more in future functional works. According to the expression data, PLDs were identified as widely expressed genes that can be involved in most of the developmental stages of plants. Our findings revealed that CsPLDs are involved in the response to high salinity stress and are probably active in a multi-branched interaction network in increasing resistance and reducing stress damages.