A well-annotated Arabidopsis, rice sPLA₂-α and β protein sequences were retrieved from Phytozome v13 (https://phytozome.jgi.doe.gov/pz/portal.html) and compiled as initial query sequences. These AtsPLA₂ query sequences were used to perform BlastP search against thirty-four plant genomes, including Chlorophyta, Bryophyta, Pteridophyta, Gymnosperm, and Angiosperm species in Phytozome v 13 (https://phytozome.jgi.doe.gov/pz/portal.html), ONEKP database (One Thousand Plant Transcriptomes Initiative, 2019) and NCBI database. Moreover, sPLA₂ sequences of several Nicotiana species (N. tabacum TN90, N. sylvestris, and N. tomentosiformis) were searched in the SolGenomics database (Mueller et al., 2005). In subsequent rounds of blast searches, PLA2 sequences from bryophytes, and charophyte and chlorophyte algae were also used as queries. Full-length protein and nucleotide sequences were retrieved and manually checked for the presence of an N-terminal secretory peptide, highly conserved PA2c domain in the ScanProsite (https://prosite.expasy.org/scanprosite/), and Interproscan 5 (https://www.ebi.ac.uk/interpro). Further, partial and truncated sequences that suggested incomplete gene predictions were curated using the SoftBerry FGENESH+ gene prediction algorithm (https://www.softberry.com); short and dubious sequences were removed from the dataset. All sPLA₂ and PLA₂-like sequences were compiled in the table with their genomic information, including gene id, protein length, and chromosomal location ( Supplementary Table 1 ).

Protein sequences were aligned using the MAFFT E-INS-I algorithm (Katoh and Standley, 2013) in Jalview (Waterhouse et al., 2009). All sequences were manually checked for gaps and non-conserved regions, which were eliminated from the alignment and exported in FASTA file format. The resulting alignment was used to build a sequence logo of the Ca²⁺ binding motif and catalytic motif (https://weblogo.berkeley.edu/logo.cgi) (Crooks et al., 2004). The evolutionary relationship between plant sPLA₂ and PLA₂-like was inferred using the IQ-TREE algorithm with Maximum likelihood (ML) supported by the ultrafast bootstrap method (1000 replicates). Model finder was performed and the Whelan and Goldman model with Invariable and gamma (WAG+I+G4) was selected as a best model based on the Bayesian information criterion (BIC) score (Minh et al., 2020). The Interactive Tree Of Life (iTOL v5) (https://itol.embl.de/) online tool was used to display and annotate the sPLA₂ phylogenetic tree. Likewise, species trees were constructed in NCBI taxonomy and edited in the iTOL server.

Tissue-specific expression of the sPLA₂ and PLA₂-like members have been searched in various vegetative (leaf, stem, root, and seeds) and reproductive (flower, anther, pollen, pollen tube, carpels, pistil, ovary, ovule, and egg cells) tissues of Arabidopsis, Amborella, tomato, grape, rice, and maize using the CoNekT database (https://conekt.sbs.ntu.edu.sg/) (Proost and Mutwil, 2018). Gene expression was represented in transcripts per kilobase million (TPM)-based normalization because it can be used for both gene count comparisons within a sample or between samples of the same sample group (Abrams et al., 2019). The expression values were analyzed in the CIMminer one matrix server (discover.nci.nih.gov/cimminer).

Total RNA was isolated from tobacco leaves, roots, buds, flowers, imbibed pollen, germinating pollen grains and growing pollen tubes using Qiagen RNAeasy Kit, and Turbo DNA-free Kit (Applied Biosystems, Waltham, MA, USA) was used for DNA removal. cDNA synthesis was carried out using Transcriptor High Fidelity cDNA Synthesis Kit (Roche, Penzberg, Germany) with anchored-oligo (DT)₁₈ primer according to manufacturer’s instructions. Semi-quantitative RT-PCR was performed with NtPLA₂ gene-specific oligonucleotides 1-6 ( Supplementary Table 2 ) designed to span an intron in the corresponding genomic DNA sequence. Actin7 (Bosch et al., 2005) was used as load control. Amplification conditions were 94°C for 30 sec, 55°C for 30 sec, 68°C for 30 sec and final extension 68°C for 10 min for 28 or 34 cycles.

The promoter regions of sPLA₂ and PLA₂-like genes (~ 2kbp upstream of the start codon) were retrieved from the Phytozome v13 database. Promoter sequences were then analyzed via the PlantCARE server with default parameters (http://bioinformatics.psb.ugent.be/webtools/plantcare/html). The obtained cis-acting elements data were processed using the CIMminer server.

Homology models of N-terminal truncated peptides signals of Amborella trichopoda (ATR0564G112-α; 129 aa, ATR0789G151-β; 118 aa), Arabidopsis thaliana (AT2G06925-α; 128 aa, AT2G19690-β; 119 aa), and Nicotiana tabacum TN90 (Gene_60450(SS4740)-α; 127 aa, Gene_37244 (SS1768)-β; 117 aa) were obtained through AlphaFold server and also predicted using Robetta server (https://robetta.bakerlab.org/), and RaptorX server (http://raptorx.uchicago.edu/). Moreover, homology models of C-terminal regions (containing PA2c domain) of PLA2-like sequences were predicted for Amborella trichopoda (AmTr _scaffold00063; 168 aa), Arabidopsis thaliana (AT4G29070- 145 aa), and Nicotiana tobacumTN90 (mRNA_158062; 149 aa). For Robetta and RaptorX, five models of each protein were generated, and the model qualities were compared with the AlphaFold predictions using the SWISS-refinement server (https://swissmodel.expasy.org/assess). MolProbity Score, clash score, Ramachandran scores, and QMEAN were considered to select the best model among all generated predictions. Further, electrostatic potentials for the best models were calculated using APBS web-server (Adaptive Poisson-Boltzmann Solver, https://www.poissonboltzmann.org/) and plotted on the predicted protein structures (Unni et al., 2011). The conserved residues and non-conserved residues were analyzed using the Consurf web-server (https://consurf.tau.ac.il/) for the best models. All the generated models were visualized and edited using the CHIMERA (https://www.cgl.ucsf.edu/chimera/).

To get a deeper insight into the evolution of the plant sPLA2 family, we have selected thirty-four plant genomes representing a diverse and balanced sample of Viridiplantae, including algae, bryophytes, lycophytes, pteridophytes, gymnosperms, and angiosperms. We searched within this genome sample for the presence of the sPLA₂ homologs, using known dicot and monocot sPLA₂ as input and employing various homology-based searches like Blast and HMMER ( Figure 1A ). We ultimately identified 113 non-redundant canonical plant sPLA₂ genes based on the following criteria: the presence of N-terminal signal peptide and PA2c domain with Ca²⁺ binding motif and catalytic HD dyad site. Surprisingly, we also uncovered 37 additional genes coding for proteins that lack the N-terminal signal peptide but contain a well-conserved PA2c domain, including the calcium-binding motif and the catalytic dyad ( Supplementary Figure 1 , Supplementary Tables 1 , 3 ). The data showed a ubiquitous distribution of both sPLA₂ and PLA₂-like genes in plants, from unicellular algae to multicellular flowering plants ( Figure 1A ). Remarkably, the two groups exhibit a clear difference in evolutionary dynamics. PLA₂-like orthologs were found as single-copy genes in most diploid species and never exceeded two copies per genome. Two PLA2-like genes were found in either polyploid species (N. tabacum) or species that underwent relatively recent whole genome multiplications (B. rapa, G. max, P. patens). This suggests that PLA2-like gene number is under the purifying selection to remain singleton. On the other hand, sPLA₂ genes display much wider genomic plasticity, showing gene numbers ranging from one (chlorophyte and streptophyte algae) to twelve (gymnosperm T. plicata). This higher evolutionary dynamics of sPLA₂ is also apparent for individual plant clades like ferns (two to nine genes) and gymnosperms (four to twelve genes). Within angiosperms, Basal lineages (represented by A. trichopoda) seem to contain two sPLA₂ genes, while distinct dynamics can be found between dicots (two to six genes) and monocots (four to five genes). Despite the general increase in sPLA₂ gene complexity during plant evolution, reductive events can also be observed. In bryophytes, moss P. patens and liverwort M. polymorpha have two sPLA₂ isoforms, whereas only one sPLA₂ member was observed in hornwort A. angustus. Similarly, just one sPLA₂ was retained in the lycophyte S. moellendorffii ( Figure 1A ).

To get more insight into the sequence characteristics of sPLA₂ versus PLA₂-like subfamilies, we performed multiple sequence alignments of the two groups and analyzed their overall protein sizes, motif occurrences, and the distribution of conserved residues. All sPLA₂ sequences have N-terminal signal peptides and range between 90-191 residues, while PLA₂-like members range between about 143-320 residues, have extended N- and C-terminal regions, and lack any recognizable signal peptide sequence.

The sequence analysis revealed other notable differences among sPLA₂ and PLA₂-like members at amino acid levels. Significantly, all sPLA₂ members possess twelve cysteine residues forming six disulphide bridges and providing extra structural stability  (Mariani and Fidelio, 2019). In contrast, PLA₂-like members have only six C residues, possibly forming up to three disulphide bridges ( Supplementary Table 3 ). Moreover, although the highly conserved Ca²⁺ binding loop is present in all members within the PA2c domain, PLA₂-like has a subtle variation in Ca²⁺ binding motif with the insertion of an extra residue (YGHYCGxxxxxGK vs YGKYCGxxxxGC) ( Figure 1B ). Another cysteine loss occurred near the catalytic site where sPLA₂-α and β show invariable DxCCxxHDxC motif, whereas PLA₂-like has DxCCxxHDxG.

Collectively, our data show a widespread and dynamic occurrence of canonical sPLA₂ genes and document a highly conserved, evolutionary-constrained subfamily of non-characterized PLA₂-like genes.

To shed light on the evolutionary relationships of sPLA₂ and PLA₂-like families, a curated dataset of 113 sPLA₂ and 38 PLA₂-like genes was compiled from a diverse set of plants with sequenced genome drafts ( Figure 2 ). It included evolutionary important chlorophytic and streptophytic algal species such as Chlamydomonas reinhardtii, Volvox carteri, Penium margaritaceum, Klebsormidium nitens, Mesotaenium endlicherianum, and Chara braunii. Moreover, we have included several bryophytes (moss Physcomitrium patens, liverwort Marchantia polymorpha, and hornwort Anthoceros angustus), lycophyte (Selaginella moellendorffii), pteridophytes (Azolla filiculoides, Salvinia cucullata, and Ceratopteris richardii), gymnosperms (Thuja plicata and Ginkgo biloba), Basal angiosperm (Amborella trichopoda), and multiple dicots, and monocots species ( Figure 2 ). The sequences were aligned using the MAFFT E-INS-I algorithm, unalignable regions were removed, and a maximum likelihood (ML) tree was constructed using the IQ-tree software with ultra-fast bootstrap support ( Figure 2 ).

Despite a limited evolutionary signal (due to relatively short sequence lengths), the resulting tree showed clearly that plant PA2c domain-containing genes separated into three well-supported clades. All sPLA₂ sequences were clustered into two groups, referred to as sPLA2-α​ clade and β​ clade and corresponding to XI-A (sPLA₂-β​) and XI-B (sPLA₂-α​) classification of sPLA₂ genes (Mansfeld, 2009; Mariani and Fidelio, 2019). Significantly, PLA₂-like members formed a third independent clade, possibly separating sPLA₂-α and sPLA₂-β ( Figure 2 )​​. It should be noted that the two canonical clades seem to possess several specific enzymological characteristics, which were best explored in Arabidopsis and elm. Most importantly, sPLA₂-α members are show optimal activation by millimolar calcium, while sPLA₂-β show maximal activity already at micromolar calcium levels. Compared to sPLA₂ members, newly identified PLA₂-like members are large in length (>240 aa), high molecular weight, disordered N-terminal region (~100 aa) without signal peptide, and residual variation in Ca²⁺ binding loop. Moreover, the presence of clear sPLA₂-α and sPLA₂-β orthologs in gymnosperms and ferns clearly shows that the diversification and stable retention of an α- and β​-clade occurred relatively early in land plant evolution. Interestingly, the contrasting sPLA₂ distribution in algae and bryophytes (all chlorophyte and streptophyte algae retain 1 sPLA₂ isoform clustering in the β​-clade, while mosses, liverworts and hornworts possess only α​-clade members), suggests that the two clades may have emerged very early in plant sPLA₂ evolution, and the single clades were lost in distinct lineages. On the whole, the evolution of the β​-clade is clearly under the evolutionary constraint, both in gene number (except for Brassicaceae) and mutation rate ( Figure 2 ). On the other hand, rapidly evolving α​-clade is highly expanded in several lineages, including pteridophytes (C. richardii), gymnosperms (T. plicata and G. biloba), and angiosperms (E. grandis, P. trichocarpa, G. max, and M. truncatula).

As mentioned above, PLA₂-like clade, which consists of primarily single-copy genes, branched-off from sPLA₂ at the earliest stage of green plant evolution, although its presence in green algae was detected universally ( Supplementary Table 1 ). Generally, PLA₂-like phylogeny roughly follows the species evolution with a mutational rate similar to sPLA₂-α, suggesting a gradual evolution constrained at the gene number level ( Figure 2 ).

To get a comprehensive evolutionary insight into sPLA2 expression, publicly available tissue-specific RNA-Seq data of representative members of eudicots (A. thaliana, S. lycopersicum, and V. vinifera), monocots (O. sativa and Z. mays), and Basal angiosperms (A. trichopoda) have been extracted and arranged clade-specific on the phylogenetic tree ( Figure 3A ). In total, we analyzed the expression of eleven α​-clade members, eight β​-clade members, and six PLA₂-like members in four vegetative tissues (leaves, stem, root, and seeds) and ten reproductive tissues (flower, anther, pollen, pollen tube, carpels, pistil, ovary, ovules, and egg cells). The data strongly suggest the existence of conserved clade-specific expression patterns across angiosperms. Members of the sPLA₂-α clade in all six analyzed species exhibited significant expression in the vegetative tissues, including leaf, stem, root, and seeds. Notably, while most dicot, monocot and Amborella sPLA₂-α members are expressed in the various maternal reproductive tissues, such as flowers, anther, carpels, ovules, and egg cells, they are all conspicuously absent from pollen and pollen tube.

On the other hand, sPLA₂-β members across angiosperms showed significant (often the strongest) expression in male gametophytic tissues such as pollen, and pollen tube (and also anther), in addition to variable sporophytic expression. This is well illustrated on the three β-clade and one α-clade Arabidopsis members, where our global expression data analysis is corroborated by earlier RT-PCR and promoter studies  (Bahn et al., 2003; Wang et al., 2008). Thus, AtsPLA₂-α​ and β​ show similar expression in all tissues except male gametophyte, where AtsPLA₂-β​ is significantly expressed. Conversely, AtsPLA₂-γ (and to a lesser extent also AtsPLA₂-δ) displayed significant expression predominantly in the male gametophyte tissues, including anther, pollen and pollen tube, suggesting that these isoforms might play an essential role in plant reproduction (Bahn et al., 2003; Lee et al., 2005; Wang et al., 2008). The general expression profile of angiosperm PLA₂-like then suggests a ubiquitous expression with the highest values typically occurring in the sporophyte and species-diversified expression in the male and female gametophytes ( Figure 3A ).

To experimentally corroborate these in silico studies, we performed a semiquantitative RT-PCR analysis of tobacco PA2c homologs (please note that the sequences of tobacco sPLA₂ and PLA₂-like homoeologous gene pairs are nearly identical (> 98% nucleotide identity) in tobacco amphidiploid genome and were considered as single cDNAs in the RT-PCR analysis). Analysis of leaf, root, bud, flower, dry pollen, germinating pollen and growing pollen tubes showed that tobacco PA2c genes show expression patterns strongly supporting the RNAseq data on other species, particularly the absence of AtsPLA₂-α in pollen and pollen tubes, where PLA₂-β is strongly present ( Figure 3B ). Similarly, tobacco PLA₂-like was expressed predominantly in the sporophyte but could be clearly detectable also in the male gametophyte.

Taken together, our comprehensive examination of available transcriptomic data of six angiosperm species and RT-PCR analysis of sPLA₂ and PLA₂-like genes in tobacco showed evolutionarily conserved expression patterns for the three clades, highlighting their overlapping expression in most sporophytic tissues, and suggesting a significant role for the sPLA₂-β clade in pollen.

Having established conserved expression patterns for the three clades of the angiosperm PLA₂ family, we next sought to understand their transcription regulatory networks. Therefore, we predicted cis-elements for the promoter sequences of sPLA₂ and PLA₂-like genes from Arabidopsis, wine, rice, and Amborella ( Figure 4 ). Various categories of cis-elements were found, including common elements TATA and CAAT box and tissue-specific, light, stress, and hormone-responsive elements.

While no obvious pattern of evolutionarily conserved subset of cis-acting motifs could be attributed to individual sPLA₂ or PLA₂-like clades, the analysis strongly suggested that specific stress- and hormone-responsive elements are the most significantly overrepresented category in angiosperm sPLA₂ and PLA₂-like genes. These include abiotic stress-responsive elements, like dehydration-responsive element MBS (TAACTG), Myb (CAACTG), MYB (CAACAG), MYB recognition sites (CAACAG), MYB-like sequence (TAACCA) and MYC (CATTTG).

In addition to the abiotic stress-related elements, motifs implicated in biotic interactions, such as wound-responsive elements W-box (TTGACC), WRE3 (CCACCT), and as-1, were detected in most promoter regions for distinct species. Among those, the WRE3 motif is the most abundant and reported in almost all sPLA₂-α and β gene members except Amborella and grape genes. W-boxes were found to interact with transcription factors belonging to the WRKY family, regulate defence-related genes, and play a vital role in biotic and biotic stress, senescence and seed dormancy.

Several hormone-responsive elements were documented in the promoter regions, among which abscisic acid-responsive element ABRE (TACGTG) and methyl jasmonate-responsive cis-elements (CGTCA and TGACG) are the most abundant and conserved, particularly within the sPLA₂-α clade. On the other hand, ethylene-responsive element ERE (ATTTCAAA), gibberellic acid-responsive elements (GARE: AAACAGA, PA box, and TATC box), auxin responsiveness core element (auxRR: GGTCCAT, TGA elements), show much lower abundance and are scattered among species ( Figure 4 ).

Lastly, many light-responsive cis-elements have been documented in the sPLA₂ and PLA₂-like promoter regions, including I-box, G box, GATA-motif, GA motif, TCT motif, Box 4, and MRE (MYB binding light responsive elements). The most abundant and conserved motif is G-box (CACGTG), which is involved in the light, abscisic acid, methyl-jasmonate, ethylene and anaerobiosis responses (Sibéril et al., 2001).

Collectively, the abundantly predicted cis-acting elements corroborate earlier functions experimentally attributed to selected sPLA₂ members and suggest a conserved transcriptional control. Indeed, sPLA₂ transcription was activated in response to blue light (Seo et al., 2008), auxin (Scherer, 2002; Lee et al., 2010), abiotic stresses  (Chapman, 1998), wound stress and pathogen elicitors  (Creelmen and Mullet, 1997; Lee et al., 1997; Laxalt and Munnik, 2002; Ellinger et al., 2010).

Despite a wealth of knowledge about PLA₂ physiological function, there is limited structural data on the plant sPLA₂-α​, β​, and PLA₂-like members, thus impeding our understanding of sPLA₂ substrate preferences, interfacial recognition surface (IRS), and catalytic sites. To get better insight into the structural features of sPLA₂ and PLA₂-like proteins, we thoroughly analyzed structural models predicted for representative members of sPLA₂-α​, β, ​and PLA₂-like from Arabidopsis, tobacco, Amborella ( Figure 5 ; Supplementary Figure 2 , 3 ). Since the current databases of automatically-generated structural models include the signal peptide, we generated de novo models of the processed sPLA₂ forms. Moreover, since the N-terminal extensions in PLA₂-like sequences show the characteristics of an intrinsically disordered region, we analyzed only the C-terminal portions. In a pilot analysis, we generated the models using three top-ranked algorithms that are not constrained by existing experimental structures (AlphaFold2_mmseq2 via the ColabFold infrastructure, RoseTTAFold, and RaptorX). We assessed them for their folding accuracy based on the Molprobity score, Ramachandran score and Q means plus Z-score criteria (Wang et al., 2016; Waterhouse et al., 2018; Jumper et al., 2021; Baek et al., 2021). Almost universally, AlphaFold-generated models displayed the best criteria and were therefore selected for further analyses ( Supplementary Table 4 ).

Comparative analysis of the final validated models revealed that all display a rather tightly-packed globular structure corresponding to the general sPLA₂ fold ( Figure 5A ; Supplementary Figure 2 , 3 ). Two main structural regions are present across the analyzed species in all sPLA₂ and PLA₂-like models. The N-terminal part contains mainly a loop region, including a conserved Ca²⁺-binding loop, and is preceded by two antiparallel beta strands. The C-terminal segment is then represented by three antiparallel α-helices, of which the two first are also present in other secreted PLA₂s and contain the conserved catalytic histidine and calcium-coordinating aspartate residues (Guy et al., 2009). Similar structural features were documented in PLA₂-like sequences, including the N-terminal loop region with Ca²⁺ binding motif and three antiparallel helices at the C-terminal region with the extended region. Notably, in canonical sPLA₂s, the calcium-binding loop and the N-terminal and C-terminal parts are held together by six disulfide cysteine bridges that stabilize the overall structure. However, not all cysteines are conserved in PLA₂-like sequences, and the structural models suggest that only two cysteine bridges may be retained in PLA₂-like structures, giving them higher conformational flexibility.

Since sPLA₂ is an interfacial enzyme interacting with the phospholipid bilayer, we next sought to see the electrostatic potential distribution in all analyzed models that could indicate distinct membrane-binding patterns for separate clades. Our data show that sPLA₂-β members in all analyzed species contain positively-charged pockets, which may be involved in anionic lipid binding. On the other hand, all PLA₂-like models lack this feature and produce mostly neutral surfaces. The highest diversity was found within sPLA₂-α​ clade, where significant charge differences can be found between Arabidopsis, Nicotiana and Amborella members ( Figure 5B ; Supplementary Figure 2 , 3 ). The low evolutionary conservation of the charge distribution was also corroborated by the visualization of evolutionarily conserved surface residues, suggesting that among sPLA₂-α​, β​, and PLA₂-like members, there is rather low conservation beyond the Ca²⁺ binding motif and the catalytic region ( Figure 5C ).

Despite the critical importance of PLA₂ documented in higher plants and the availability of detailed molecular information for phospholipases from non-plant organisms, evolution-based knowledge of phospholipase A2 from plants is still meagre. To fill this gap, we have performed a comparative phylogenomic, sequence and structural analysis of the plant’s secretory PLA₂ family, along with previously uncharacterized PLA₂-like members. This is the first systematic analysis of the plant family of PA2c domain-containing proteins, covering chlorophyte and streptophyte algal species, bryophytes, lycophyte, pteridophytes, gymnosperms, Basal angiosperm, dicots and monocots. However, the detailed evolutionary history and clade separation of the plant PA2c domain superfamily is murky. Our data demonstrated that one genuine sPLA₂ can already be documented in multiple algae species (chlorophytes, charophytes and zygnematophytes), and all canonical algal sPLA₂ fall into the XI-A (sPLA2-β) clade. On the other hand, all canonical sPLA₂ genes from three main bryophyte groups (hornworts, liverworts, and mosses) belong to the XIB (sPLA₂-α) clade, and the simultaneous presence of both clades can be first documented in ferns. Three evolutionary scenarios can explain this diverse distribution: (i) sPLA₂-β is the ancestral plant form of secretory phospholipase A₂, which through gene duplication and subsequent diversification, gave rise to the sPLA₂-α clade during the plant colonization of earth but was secondarily lost in bryophytes and lycophytes; (ii) both clades were present already in the common ancestor of chlorophyte and streptophyte algae and were sometimes lost in particular groups; (iii) individual sPLA₂ clades were subject to the independent horizontal gene transfer events during the evolution of separate Viridiplantae groups. These scenarios are not necessarily mutually exclusive, especially considering the widespread occurrence of horizontal gene transfer in green algae (Ma et al., 2022).

One notable element of plant sPLA₂ evolution is the diversity in gene copy numbers for the two clades. The terminal duplication and diversification events of sPLA₂ may be attributed to whole-genome duplication (eg. in some bryophyte species), polyploidization or amphidiploidization events (eg. N. tabacum or Z. mays). In addition, in some pteridophyte and gymnosperm species, sPLA₂ underwent massive multiplication. Several hypotheses explain plant gene family expansion besides the whole-genome duplications and hybridization events. Gene families can expand through either segmental, tandem, or retro-transposition (RT) mechanism, but segmental and tandem duplication events were more predominant than RT (Kondrashov, 2012; Panchy et al., 2016). These new paralogs may perform an existing gene function (sub-functionalization) or acquire a novel role (neo-functionalization)  (Panchy et al., 2016). These gene duplication mechanisms may often co-occur, as evidenced in the Eucalyptus grandis, which genome (n=11) has been shaped by lineage-specific genome duplication events and a high rate of tandem gene duplication  (Myburg et al., 2014), leading to the highest number of sPLA₂ members found in the analyzed angiosperms. As noted above, most of these within-family multiplication events occurred in the sPLA₂-α clade, which is evolving considerably faster than the sPLA₂-β clade. The notable exception is Arabidopsis (and the whole Brassicaceae clade, Supplementary Table 5 , Supplementary Figure 5 ), where the expansion occurred predominantly in the sPLA₂-β clade.

Our analysis of the sPLA2 expression revealed that while both -α​ and -β clades are present in various sporophytic tissues, sPLA2-α​ genes seem to be absent in the male gametophyte in dicots, monocots and Basal angiosperms ( Figure 3A ). This contrasts with the expression of the angiosperm sPLA₂-β clade members, expressed either exclusively in the pollen and pollen tubes (Arabidopsis sPLA2-γ and -δ), or showing the highest expression in the male gametophyte (sPLA₂-β from Arabidopsis, tomato, monocots, and Amborella). Significantly, our high-throughput data analysis is corroborated by the RT-PCR analysis of tobacco sPLA₂ members ( Figure 3B ) and by multiple reports on Arabidopsis sPLA₂ members (Lee et al., 2005; Kim et al., 2011). Although the split between the two sPLA2 clades occurred before the emergence of sexual reproduction via pollen, the presence of only sPLA2-β ortholog already in Amborella pollen and pollen tube suggests that the β-clade is the major sPLA2 in all angiosperms. Indeed, Kim et al. (2011) showed that when all three Arabidopsis PLA2-β clade members were suppressed by RNA interference, pollen development and germination were severely affected. Moreover, lysophosphatidylethanolamine, the product of sPLA2 activity, likely plays a vital role in pollen germination and pollen tube growth (Kim et al., 2011). Significant changes in lysophospholipid levels, including plasma membrane lysophosphatidylethanolamine and lysophosphatidylcholine, were recently described during tobacco pollen germination (Serrano et al., 2022). These findings further corroborate the conserved role of sPLA2 in the male gametophyte. Interestingly, compared to other beta-clade members, Arabidopsis sPLA2-γ shows the highest expression in the growing pollen tubes, suggesting that the three different β-clade members in Brassicaceae may have distinct roles during pollen development, germination, and tube growth.

A feature reportedly distinguishing plant sPLA2-α​ and -β clades is their distinct localization, suggesting that sPLA2-β may act primarily inside the endomembrane system and are not secreted to the apoplast (Fujikawa et al., 2012);. Three Arabidopsis sPLA2-β clade paralogs (β, γ, and δ) localized mainly to the ER and/or Golgi (Seo et al., 2008; Kim et al., 2011) while sPLA₂-α localizes either to the apoplast or to the nucleus, depending on the plant status (Froidure et al., 2010; Jung et al., 2012). However, a canonical C-terminal ER-retention signal (KxEL) can be found only in 14 angiosperm sPLA₂-β members and is missing even in Arabidopsis sPLA₂-γ, and -δ. Therefore, the subcellular localization of sPLA₂s, particularly the β-clade, might be variable and needs to be determined in other angiosperms.

While annotating the plant’s sPLA₂ sequences, we came across a subfamily of proteins with unusually long sequences (~250 aa), that possess a conserved PA2c domain containing both Ca²⁺ binding motif and HD catalytic dyad but lacking a signal peptide, which we termed PLA₂-like. While the presence of the PLA₂-like subfamily was briefly noticed before (Gupta et al., 2017), its phylogenomic distribution and sequence-structural properties were not explored. Our genome-wide analysis demonstrated that at least one PLA₂-like member is consistently present from charophyte algae to higher plants, highlighting an ancient origin. It should be noted that PLA₂-like genes can also be found in chlorophyte algae, but their distribution is patchy, suggesting either frequent losses or horizontal gene transfer in the plant group. Phylogenetically, PLA₂-like sequences clustered into one clade, representing the independent evolutionary history of PLA2-like members and separating XI-A and XI-B sPLA₂s. As noted above, the low gene-copy number of PLA₂-like genes is under strong selection, effectively keeping PLA₂-like as a single-copy gene in most species. While the putative enzymatic activity of PLA₂-like remains obscure, the single-copy status is typically linked with genes often involved in essential metabolic processes and/or the formation of macromolecular complexes. Interestingly, the N-termini of PLA₂-like proteins show intrinsically disordered region features, which function in protein-, DNA-, or RNA- binding (Han et al., 2014).

Our analysis of cis-acting elements in promoters of evolutionarily-distinct angiosperms species suggested that the transcription of both sPLA₂ and PLA₂-like genes is regulated by various biotic, abiotic, and developmental stimuli ( Figure 4 ). We have reported different cis-elements (W-box, WRE2, and WUN motif) involved in biotic stress responses, and found several methyl-jasmonate responsive elements in the sPLA₂ promoters, further suggesting their active involvement during plant reactions to biotic stresses. The transcriptional activation of sPLA₂ genes after a pathogen attack was shown in diverse species such as Arabidopsis (Froidure et al., 2010) and grapevine (Laureano et al., 2018), supporting our in-silico analysis. Analogously, sPLA₂ promoters have multiple abiotic stress-responsive, abscisic acid-responsive and ethylene-responsive elements in the promoter region, suggesting that sPLA₂ may have some regulatory roles during abiotic stresses. Verlotta et al., 2013 characterized durum wheat sPLA₂s and assessed their involvement in drought stress (Verlotta et al., 2013). Similarly, Rice sPLA₂ also showed a differential expression pattern in response to abiotic stress (Singh et al., 2012).

Besides the involvement in the stress responses, several case studies from distinct species corroborate the conserved transcriptional regulation of sPLA₂ after developmental cues ( Figure 4 ). Arabidopsis sPLA₂-β is upregulated by auxin and in the curved regions of the peduncle, which were undergoing the gravitropic response (Lee et al., 2003). In Citrus sinensis, the sPLA2-α and β genes expression and enzyme activity in leaves and fruits exhibited diurnal rhythmic changes and light regulation, which suggested that diurnal fluctuations in lipophilic second messengers are involved in the regulation of physiological functions (Liao and Burns, 2010).

Very few high-throughput studies described the transcriptional regulation of the Arabidopsis PLA₂-like gene. Interestingly it was demonstrated among genes with expression changes between pollen germination and tube growth  (Wang et al., 2008), and listed among genes involved in Arabidopsis acyl lipid metabolism  (Beisson et al., 2003). In the present study, we have predicted that PLA2-like members exhibited ubiquitous expression throughout the vegetative and reproductive tissues. In addition, it has light, stress, and hormonal-responsive elements in the promoter region that may be involved in diverse regulatory mechanisms and cellular functions.

The comparative analysis of structural models for diverse plant sPLA₂ and PLA₂-like members strongly suggested that there are only subtle differences in the catalytic sites among typical members of sPLA₂-α and -β clades. All tested models displayed characteristic sPLA₂ features, shared with the experimentally-determined structures of non-plant sPLA₂s from groups I, II, X and also rice XI-B sPLA₂- namely three α-helices, a short two-stranded β-sheet, and a conserved calcium-binding loop  (Dijkstra et al., 1981; Holland et al., 1990; White et al., 1990; Jabeen et al., 2006; Guy et al., 2009). On the other hand, the total conformation (due to the flexibility of N-terminal and C-terminal regions) and charge distribution differ substantially between the two clades or between different phylogenetic groups. For example, the electrostatic charge differences would affect the membrane-binding and protein-binding properties and thus may impact the sPLA₂ catalytic functions. Notably, surface charge differences affecting the interactions with negatively-charged phospholipids were also described for non-plant sPLA₂ members from groups IIA and X (Quach et al., 2014).

Importantly, template-free modelling of the C-terminal half of PLA₂-like proteins suggested a fold similar to canonical sPLA₂ structures, including a calcium-binding loop and catalytic dyad. However, the striking difference between the sPLA₂ and PLA₂-like is the loss of the disulfide bridges, where only four out of twelve Cys residues - forming two disulfide bridges - are conserved in PLA₂-like sequences. On the other hand, the two retained disulfide bridges (C140-C167, C166-C192) - which connect helix 1 with the calcium loop and the helices 1 and 2 together - structurally correspond to those that are vital for the PLA₂ catalytic activity of porcine PLA₂ (C29-C45, C44-C105, Zhu et al., 1995). Therefore, despite the lower number of disulfide bridges retained in PLA₂-like structures, PLA₂-like members may still possess an enzyme catalytic activity towards lipidic substrates analogous to canonical sPLA₂ (Mansfeld, 2009). In addition to the phylogenomic survey of canonical sPLA₂ genes, our study may also serve as a call for further enzymological and physiological characterization of the enigmatic PLA₂-like gene subfamily.