Genome-wide identification and molecular evolution analysis of BPA genes in green plants

Abstract

Reactive oxygen species (ROS) signaling plays a central role in plant immune response. BPAs, referred to as binding partner 1 of accelerated cell death 11 (ACD11) (BPA1) and BPA1-like proteins, regulate ROS-mediated defense responses in Arabidopsis thaliana. However, their distribution and evolutionary characteristics in the plant lineage remain unexplored. In this study, we demonstrated that most BPA genes form a plant-specific family with expansion events observed. We found that BPA and ACD11 genes co-exist in all land plants, suggesting that this immune regulatory module may originate at the early stage of land plant emergence and contribute to their colonization. Angiosperm BPAs can be classified into four distinct groups (I-IV) in our analysis. Domain organization and motif composition are highly conserved within each group but divergent across different groups. In certain species, BPAs undergo complex alternative splicing, suggesting their regulatory and functional divergence. The protein-protein interaction network we constructed predicted additional acting partners of BPAs. The yeast two-hybrid assay revealed 15 BPA interaction pairs forming homo- or hetero-dimers. Taken together, our results provide the first synopsis of BPA evolutionary pattern and adaptations to green plant colonization.

Background

Host plants and microbial pathogens are engaged in a constant evolutionary arms race. To counteract pathogen invasion, plants have evolved a two-tier immunity system (Dangl et al. 2013). In the first tier, pattern recognition receptors (PRRs) located in the plasma membrane recognize conserved microbe-associated molecular patterns (MAMPs), such as bacterial flagellins and oomycete elicitins (Mukhtar et al. 2016). The recognition leads to MAMP-triggered immunity (MTI), which is sufficient to halt most pathogens. Some highly adapted pathogens secrete effector proteins to interfere with MTI. These effectors can be directly or indirectly recognized by intracellular nucleotide-binding domain leucine-rich repeat containing (NLR) receptors, which constitute the second tier of plant defense known as effector-triggered immunity (ETI). ETI is a robust response that often includes programmed cell death (PCD) and primers systemic acquired resistance (SAR) (Cui et al. 2015). Despite their significant differences in the activation mechanisms, MTI and ETI, however, share some vital signaling pathways including the burst of reactive oxygen species (ROS) (Torres et al. 2006).

Initially recognized as toxic by-products of aerobic metabolism, ROS are now considered as a major class of signaling molecules in plant immune response (Baxter et al. 2014; Mittler 2017; Waszczak et al. 2018). The balance between ROS producing and scavenging determines the two-faced roles of ROS as either suppressing or promoting pathogen infection (Waszczak et al. 2018). Intracellular ROS are generated primarily in chloroplasts, mitochondria and peroxisomes/glyoxysomes, whereas plasma membrane localized NADPH oxidases, amine oxidases and cell wall peroxidases are responsible for the generation of apoplastic ROS (Mignolet-Spruyt et al. 2016). Major ROS-scavenging enzymes include mitochondrial oxidase (AOX), catalase (CAT), copper/zinc superoxide dismutase 2 (CSD2) and ascorbate peroxidase (APX1) (Mittler et al. 2004). Besides ROS-producing and scavenging enzymes, several protein regulators essential for maintaining ROS homeostasis have been identified in Arabidopsis thaliana (Qi et al. 2017). BOTRYTIS-INDUCED KINASE 1 (BIK1), which belongs to the receptor-like cytoplasmic kinase (RLCK) family, interacts with respiratory burst oxidase homolog protein D (RBOHD) to enhance extracellular ROS production (Li et al. 2014). The calcium-dependent protein kinase CPK28 suppresses RBOHD-mediated ROS production by interacting with BIK1 to facilitate its turnover (Monaghan et al. 2014). Interestingly, CPK28 also modulates vegetative stage transition via tissue-specific balancing of jasmonic acid (JA) and gibberellic acid (GA) (Matschi et al. 2015). We recently reported that the Arabidopsis binding partner 1 of accelerated cell death 11 (ACD11) (BPA1) and BPA1-like homologs are novel regulators of ROS accumulation and cell death under biotic stresses (Li et al. 2019).

ACD11 encodes a sphingosine transfer protein which belongs to the glycolipid transfer protein (GLTP) superfamily. Its null mutant acd11 exhibits accelerated PCD and constitutive immune response activation phenotypes in the absence of pathogen attack (Brodersen et al. 2002; Braun et al. 2011). BPA1 was initially reported to interact with ACD11 in yeast two-hybrid (Y2H) screen and co-immunoprecipitation assay (Petersen et al. 2009). Likewise, the six BPA1-like homologs in Arabidopsis, namely BPL1–6, all interact with ACD11 (Li et al. 2019). All seven BPAs negatively regulate plant resistance to Phytophthora capsici in a functionally redundant manner (Li et al. 2019). Furthermore, BPA1 and BPL2/3/4 can stabilize ACD11 to suppress ROS production and cell death (Li et al. 2019). Besides BPAs, the alternative splicing isoform of a Golgi-located E3 ligase, XBAT35.2, also interacts with ACD11 and promotes its 26S proteasome-dependent turnover (Liu et al. 2017).

Besides being partners of ACD11 in regulating ROS production and cell death, BPAs may interact with additional unknown proteins to modulate plant immunity and other biological processes as well. Their evolutionary pattern across green plants is also to be explored. The increasing availability of sequenced genomes enables us to perform a genome-wide analysis of BPA gene repertoires across the tree of life. Here we reported key evolutionary features detected in the BPA gene family including phylogeny, conserved domains and motifs, and alternative splicing events. We also predicted additional interacting proteins and the regulatory network of BPAs. In particular, we illustrated the interaction map of all Arabidopsis BPAs. Taken together, our results revealed the evolutionary pattern of BPAs and provided clues for further investigation of their functions, interacting partners and regulatory mechanisms.

Results

Most BPAs form a plant-specific gene family

Due to their novelty, BPA genes have been previously identified only in Arabidopsis. To characterize BPAs across the tree of life, seven known Arabidopsis homologs (Li et al. 2019) were used as queries to perform BLASTP and PSI-BLAST searches against the National Center of Biotechnological Information (NCBI) non-redundant protein database with an e-value cutoff of 1e-5. The hit sequences were filtered by the presence of an RNA recognition motif (RRM_1) domain. A total of 1868 candidate BPA-encoding genes were detected exclusively in two eukaryotic lineages, with 1541 homologs from green plants and 327 homologs from fungi (Fig. 1 and Additional file 1: Table S1). Specifically, all fungal candidates belong to the previously reported Vip1 gene family (Rhind et al. 2011).

Fig. 1

The distribution of BPA and ACD11 homologous genes in the tree of life. The left panel depicts eukaryotic phylogenetic relationships based on the Tree of Life phylogeny (http://tolweb.org/). Lineages with bold fonts indicate the presence of BPA or ACD11 homologs. The right panel shows the numbers of BPA or ACD11 homologs identified in each lineage

Being the only known partner of BPA (Petersen et al. 2009; Li et al. 2019), ACD11 has 630 candidate homologs in two eukaryotic lineages as revealed by our search using similar criteria. 410 and 220 putative ACD11-encoding genes were found in green plants and animals, respectively (Fig. 1 and Additional file 2: Table S2). Despite their wide distributions in the plant kingdom, both BPA and ACD11 genes are absent in Rhodophyta and glaucophytes.

Expansion of BPAs in land plants

A total of 160 BPA homologous sequences from 22 plant species were kept after manual curation (Fig. 2 and Additional file 3: Table S3). These BPAs are distributed in monocots (6 species: 61 sequences), dicots (12:88), basal angiosperms (1:4), bryophytes (1:5) and chlorophytes (2:2), with no homologs found in animals. Regarding protein sizes, most predicted BPAs are similar to their Arabidopsis homologs (Li et al. 2019) with an average length of 286 amino acids (Additional file 3: Table S3). BPA copy numbers vary across plant species, ranging from 0 in two chlorophytes (Volvox carteri and Chlamydomonas reinhardtii) to 18 in wheat (Triticum aestivum). Every land plant species examined has 4 or more BPA copies (Fig. 2 and Additional file 3: Table S3). In contrast, only a single copy of BPA can be detected in two chlorophytes and all fungi species examined. These results indicate that BPA gene duplication events likely occurred in land plants after their divergence from chlorophytes.

Fig. 2

The schematic of species phylogenetic relationships and distribution of BPA and ACD11 homologous genes in different species. The species tree was generated using the NCBI Taxonomy tool. Classified and total numbers of BPA and ACD11 homologs identified in each species were listed

Similarly, 35 ACD11 homogeneous sequences were retrieved from 27 species (Fig. 2), including animals (7 species: 7 sequences), monocots (6:8), dicots (12:17), basal angiosperms (1:2) and bryophyte (1:1). Interestingly, none of the fungi or chlorophytes examined harbors ACD11. Unlike BPAs, ACD11 gene duplication events can only be detected in six land plant species. The observation that BPA and ACD11 genes co-exist in all land plants we surveyed indicates the establishment of their interaction in the early stage of land plant emergence.

BPA genes exhibit early divergence in angiosperms

A maximum-likelihood phylogenetic tree was constructed based on 160 BPA genes from representative species. Fungal Vip1 genes were included as an outgroup. BPAs in green plant were clustered into three distinct clades, consistent with their classifications in angiosperm, bryophyte, or chlorophyte (Fig. 3a). Notably, angiosperm BPAs can be further split into four groups designated as BPA-I to -IV (Fig. 3a). BPAs from monocots, dicots and basal angiosperms can be found in all four groups, suggesting the existence of four ancestral BPA paralogs in the most recent common ancestor (MRCA) of angiosperms. The 5 BPAs in Amborella trichopoda form a distinct group, indicating the independent expansion of bryophyte BPAs after their divergence from angiosperms.

Fig. 3

The phylogenetic relationship of BPAs and their conserved motifs. a The maximum-likelihood phylogenetic tree of BPA sequences from representative species. Fungal sequences were used as the outgroup (green). Angiosperm BPAs were classified into four clades (I-IV). Red, purple, orange, aqua green or blue labels indicate proteins from monocots, dicots, basal angiosperms, bryophytes and chlorophytes, respectively. b Domain structures and conserved motifs in each clade/group. Motifs are distinguished by different colors as shown in the legend. Domain structures were predicted by scanning against the Pfam database (El-Gebali et al. 2019). Motifs were identified using MEME (Bailey et al. 2009)

After scanning green plant BPA proteins against the Pfam database, we found that they all contain an RRM_1 domain with exception of Bra004270.1 from Brassica rapa (Fig. 3b and Additional file 4: Table S4). Bra004270.1 harbors a DUF747 domain with unknown function (Li et al. 2011). In addition, 10 significantly-overrepresented (E value <1e-5) novel motifs of 11–50 residues were identified in BPAs using the MEME motif detection software (Fig. 3b). Motifs 1, 2, 4, 5 and 6 were present in angiosperms IV and bryophyte group. Motifs 1, 2 and 4 were present in all groups of green plants, whereas motif 10 was specific to fungi. Motifs 1, 2, 4 and 5 were present in chlorophytes group. Motif 1 corresponds to RRM_1 (RNA recognition motif) domain. It was recently reported that the domain-containing gene may play a key role in plant immunity (Zhai et al. 2019) while none of other motifs can be found in the Pfam database. They distribute unevenly in different clades and angiosperm groups with motif 10 being specific to fungi. The distribution pattern of these conserved motifs may reflect the functional divergence of BPA proteins during green plant evolution.

Alternative splicing may enhance the functional diversity of BPAs

Besides gene duplication, alternative splicing is another evolutionary mechanism that increases functional diversity (Krylov et al. 2003; Reddy et al. 2013), which may be critical for plant stress responses (Syed et al. 2012; Filichkin et al. 2015). The occurrence of alternative splicing in BPAs was inspected in our data set (Fig. 4a). Totally 60 alternative splicing events that lead to the peptide change from 39 BPA genes were detected in 13 land plants (Fig. 4b), ranging from 1 to 18 events in each species.

Fig. 4

Alternative splicing of the plant BPA genes. a Alternative splicing events identified in land plant BPAs. Alternative splicing events were mapped onto the BPA phylogenetic tree and labeled with their pattern numbers. ①, lacking peptide in the non-domain region; ②, having supplementary peptide in the non-domain region; ③, lacking peptide in the domain region; ④, having supplementary peptide in the domain region; ⑤, lacking domain region. b Summary of alternative splicing events in selected species. c Alternative splicing patterns detected in this study

The alternative splicing that leads to the peptide change of BPAs has five patterns (Fig. 4c). Among the 60 BPA isoforms, 17 proteins lack peptide in the non-domain region (Patterns No. 1), 20 proteins have supplementary peptide in the non-domain region (No. 2), 5 proteins lack peptide in the domain region (No. 3), 1 protein has supplementary peptide in the domain region (No. 4), and 17 proteins are truncated with domain being removed (No. 5). In particular, the RRM_1 domain is disrupted in the last three alternative splicing patterns. Patterns No. 1, 2 and 5 are widely distributed across all land plant clades and groups. Pattern No. 3 is present in all four angiosperm groups, whereas pattern No. 4 can only be found in angiosperm BPA-II (Fig. 4a). In addition, we detected 28 BPA genes undergoing alternative splicing events that change the UTR region, such as Arabidopsis BPL5 and BPL6. The alternative splicing events detected may increase the functional diversity of BPA isoforms.

The protein-protein interaction network of Arabidopsis BPAs

Besides ACD11 (Petersen et al. 2009; Li et al. 2019), BPAs may have additional interacting partners. The interaction network of seven Arabidopsis BPAs was generated using the STRING software (Fig. 5a). A total of 70 candidate proteins potentially interact with five BPAs (BPA1, BPL1, 3,5 and 6) with no partners being predicted for the two members left. GO annotation showed that the interacting candidates participate in a wide range of developmental, signaling and immune response processes (Fig. 5b). There are some interesting clues in the network. For example, BPA1 and BPL1 were predicted to interact with SUPPRESSOR OF NPR1–1 CONSTITUTIVE 4 (SNC4), which is an atypical receptor-like kinase essential for PTI response (Bi et al. 2010; Li et al. 2014). Nine GLYCEROPHOSPHODIESTER PHOSPHODIESTERASE (GDPD) family lipid metabolic proteins, including SHAVEN 3 (SHV3), SHV3-LIKE (SVL) 1–5 and GDPD4–6, are involved in cell wall organization and root hair morphogenesis (Hayashi et al. 2008; Cheng et al. 2011). They are all predicted interactors of BPA1 and BPL1. Furthermore, both BPA1 and BPL1 were also predicted to interact with METALLOTHIONEIN 1A (MT1A), which plays an important role in copper homeostasis and seed development (Benatti et al. 2014). BPL3’s interacting partners include MITOCHONDRIAL GRPE 1 (MGE1), a contributor to plant high-temperature adaptation (Chen et al. 2019). Overall, this predicted interaction network greatly expands our future research directions on BPA functions.

Fig. 5

The protein-protein interaction network of Arabidopsis BPAs. a The protein-protein interaction network of Arabidopsis BPAs. The network was predicted using STRING (https://string-db.org/). There are 77 protein nodes and 297 protein-protein association edges in the network. Green, red, and blue lines indicate interactions predicted from gene neighborhood, fusions, and co-occurrence, respectively. Light green, black, and dark blue lines indicate additional interactions inferred from text mining, co-expression, and protein homology, respectively. b GO annotation of predicted BPA interacting proteins

To further investigate the biological roles of these BPAs, we analyzed the expression of these genes using the Genevestigator database (Additional file 5: Figure S1). Generally, the BPAs show expression across all developmental stages except that BPL5 is not detectable (Additional file 5: Figure S1a). BPL4 shows the highest expression level in the first eight developmental stages, while BPA1 exhibits the highest expression level in the last two developmental stages. BPA1, BPL3 and BPL6 display up-regulation throughout the entire life at different developmental stages, while BPL1, BPL2 and BPL4 display down-regulation. Next, we investigated stress-responsiveness of BPAs to 13 different abiotic and biotic stress conditions (Additional file 5: Figure S1b). BPA1, BPL1 and BPL6 are up-regulated while BPL3 and BPL4 down-regulated under cold treatment. In response to Pseudomonas syringae infection, BPA1, BPL1, BPL3 and BPL4 are up-regulated while BPL2 down-regulated. Notably, among all BPAs, BPL1 was found to be up-regulated in response to all stresses except heat stress. Taken together, expressions of BPA genes are tightly regulated in different development stages and in response to biotic and abiotic stresses.

Interactions between Arabidopsis BPAs

The interaction networks predicted indirect interactions between BPA1, BPL1, 3 and 6 (Fig. 5a). Since some plant proteins with binding activity often form dimers (Feller et al. 2011), we performed a Y2H assay to test BPA interactions in vivo. The result revealed 15 pairs of BPA homo- and hetero-dimers. Taking the pair of BPA1 and BPL1 as an example, we conducted a Y2H assay using BPL1 as bait and BPA1 as individual preys. Our results showed that BPL1 interacts with BPA1 (Fig. 6a). To confirm this association, we cloned BPA1 into the bait vector pGBKT7, and BPL1 into the prey vector pGADT7 for reciprocal Y2H assay. The result clearly showed that BPA1 associates with BPL1 in yeast (Fig. 6a). As shown in Fig. 6a and b, yeast two-hybrid assay showed four BPAs (BPA1, BPL2, 3 and 6) could form homodimers (Fig. 6a). Our results also revealed that mutilple hetero-dimers were formed between BPAs, including five for BPA1, two for BPL1, three for BPL2, four for BPL3, four for BPL4, one for BPL5 and three for BPL6. These BPA dimers may play similar and/or diverse biological roles via different combinations.

Fig. 6

Physical interactions between Arabidopsis BPAs. a Physical interactions between Arabidopsis BPAs revealed by Y2H assay. Arabidopsis BPA genes were cloned into the bait vector pGBKT7 (BD) and the prey vector pGADT7 (AD). The combination of BD-53 and AD-T was used as the positive control. BD-Lam and AD-T combination was used as the negative control. Yeast transformants were grown on selective medium (SD) lacking Tryptophan (T) and Leucine (L) and selected on SD lacking Tryptophan (T), Leucine (L), Histidine (H) and Adenine (A). Plates were photographed 3 days after inoculation. b Summary of the interactions. “+” represents interaction. “-” represents no interaction

Discussion

BPAs was initially described as a group of genes encoding RRM_1 domain-containing proteins (Petersen et al. 2009). BPA1 and its homologs were recently found to regulate plant immunity and ROS via interacting with ACD11 (Li et al. 2019). In this study, we performed a genome-wide analysis on BPA gene repertoires in green plants to infer their evolutionary history and molecular interactions.

In the present study, BPA copy numbers were systemically surveyed in 24 representative plant species, ranging from chlorophytes to terrestrial higher plants. Our analysis indicates that BPA genes exist in all land plants as well as some chlorophytes. BPA-like sequences identified in fungi were previously reported as Vip1 homologs (Rhind et al. 2011). Totally, 160 BPA sequences were identified from the species examined. The number of BPA genes ranges from 0 in the two chlorophytes (Volvox carteri and Chlamydomonas reinhardtii) to 18 in wheat (Triticum aestivum). Surprisingly, only 0 or 1 BPA gene was detected in four chlorophytes and 4 to 18 genes were identified in land plants, strongly suggesting that BPA genes may have emerged before land plants and expanded during the evolution of land plant species. Unlike BPA genes, the number of ACD11 is relatively conservative in the representative plant species. Notably, The BPA-ACD11 pair can be found universally and exclusively in land plants, suggesting that this immune regulatory module may originate at the early stage of land plant emergence and contribute to their colonization.

In the phylogenetic tree, BPA genes cluster into angiosperm, bryophyte and chlorophyte clades, indicating that the evolution of BPAs is in accordance with their taxonomic classifications. Angiosperm BPAs can be further divided into four distinct phylogenetic groups with each group containing members across monocots, dicots and basal angiosperms. This finding indicates that angiosperm BPAs may originate from four ancestral genes in their MRCA. We also found that BPA gene duplication events in angiosperms and bryophytes are independent. BPA protein domain and motif organization patterns are highly conserved within groups/clades but more diversified across them, indicating the functional divergence of BPAs during land plant evolution. Motif composition in fungal BPAs is an outlier when compared with those of the plant groups, which implies that BPA-like proteins in fungi may also be functionally different from their plant counterparts. This claim is supported by the observations that no ACD11 homolog can be found in fungi and all fungal BPAs belong to the Vip1 family.

More than 60% of plant intron-containing genes may undergo alternative splicing (Barbazuk et al. 2008; Syed et al. 2012), and play important roles in modulating plant development, pathogen response and stress tolerance. Regarding BPA genes, 60 alternative splicing events that lead to the peptide change were identified from 24% (39 out of 161) plant members, with 23 events causing the disruption of RRM_1 domain. RRM domain containing protein has important roles in regulating plant defense (Zhai et al. 2019), and domain-disrupted isoforms often exhibit remarkably different functions when compared with their corresponding normal proteins (Finet et al. 2013). Therefore, the widespread occurrence of alternative splicing in plant BPAs may also increase their functional diversity.

In the protein-protein interaction network we established for Arabidopsis BPAs, 70 proteins are interacting candidates of BPA1, BPL1, 3, 5 and 6. Out of these candidates, ACD11 is a demonstrated interactor regulating ROS and cell death (Li et al. 2019). Other potential interactors, including SNC4, GDPDs, MT1A and MGE1, are involved in PTI response (Bi et al. 2010; Li et al. 2014), root hair development (Hayashi et al. 2008; Cheng et al. 2011), copper homeostasis (Benatti et al. 2014), and plant heat adaptation (Chen et al. 2019), respectively. The diverse roles of these BPA-interacting candidates provide new clues for exploring BPA functions beyond plant immunity regulator. We also utilized a Genevestigator analysis to gain insight into the expression profiles of the BPA genes. We found that most BPAs show abundant expression across all developmental stages, suggesting broad roles of BPAs in plant development. Indeed, we noticed that silencing of BPL4 alone or BPL1 and BPL4 together in a bpl2 background has a negative effect on plant growth, while silencing or deletion of single gene has no visible growth phenotypes, indicating of functional redundancy (Li et al. 2019). Moreover, the results of the Genevestigator analysis showed that most BPA genes were predicted to be regulated by various stresses. In addition, we found that BPA genes have very diverse expression patterns. For example, BPL1 is up-regulated when treated with stresses such as cold, UV-B, ozone, Botrytis cinerea and Phytophthora infestans, whilst BPL3 is down-regulated. These results indicate that BPA genes may play important roles in stress response.

Another possible way for BPAs to enhance and/or expand their functions is to form homo- and hetero- dimers within the family. Arabidopsis BPA1 and BPL4 are known to function redundantly in modulating immunity against P. capsici (Li et al. 2019). In our study, we found BPA1 and BPL4 form a heterodimer, which may explain their functional redundancy. BPL1/2/4 also functions redundantly in modulating plant immunity (Li et al. 2019). In our observation, only BPL2 and 4 can form a heterodimer, indicating the existence of additional mechanisms for BPL1/2/4 genetic and/or physical interactions. Overall, the 15 dimers identified in our assay imply the biochemical and genetic complexity of BPA interaction and function network.

Conclusions

In the study, we report the genome-wide analysis of BPA repertoires across the tree of life. According to our results, most BPAs are plant-specific and enriched in land plants. Their sequences exhibit multiple evolutionary features including early divergence, conserved domain/motif organization at the clade/group level, and complex alternative splicing patterns. In Arabidopsis, the predicted protein-protein interaction network for four BPAs and the multiple homo−/hetero- BPA dimers identified indicate their broader roles in plant development, immunity and abiotic stress response. Taken together, our findings for the first time reveal the evolutionary pattern and interaction map for BPAs, which provide clues for further investigation of their diverse functions.

Methods

Sequence retrieval and homolog identification of BPAs and ACD11

Several resources were utilized to build a broad-scale initial data set. Sequenced genomes and predicted proteomes of 47 species (Additional file 1: Table S1) were downloaded from Phytozome (version 12.1; http://www.phytozome.net), Ensembl Genomes (release 97; http://www.ensembl.org), Joint Genome Institute (JGI) (http://genome.jgi.doe.gov) or The Arabidopsis Information Resource (TAIR) (https://www.arabidopsis.org). The downloaded protein sequences were integrated into a local protein database for homolog identification. When alternative splicing isoforms were annotated at the same locus, the longest one was selected.

BPA and ACD11 homologs were identified in three steps. Firstly, the protein sequences of Arabidopsis BPAs (BPA1, BPL1–6) and ACD11 were employed as queries to perform BLASTP and PSI-BLAST searches against the NCBI non-redundant protein database (https://www.ncbi.nlm.nih.gov) with an e-value threshold of 1e-5. Then the same BLASTP searches were performed against our local protein database with identical settings. Finally, sequences obtained from both databases were verified using NCBI CDD (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi), SMART (http://smart.embl-heidelberg.de), and PFAM (http://pfam.xfam.org/search). Proteins harboring intact RRM_1 (Pfam accession no. PF00076) or Glycolipid transfer protein domain (GLTP, Pfam accession no. PF08718) were identified for subsequent analyses.

Evolution analysis

Alignments of full-length protein sequences were performed using MUSCLE v3.8.31 (Edgar 2004) with default setting. Maximum-likelihood phylogenetic trees were constructed using IQ-TREE v1.6.8 (Nguyen et al. 2015) with automatic selection of optimal model for protein substitution and rate heterogeneity. For tree construction, the SH-aLRT test and ultrafast bootstrapping (Hoang et al. 2018) were conducted with 1000 replicates. FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree) was used for tree visualization and editing.

Conserved motifs in BPA proteins were identified using MEME 5.0.5 (Bailey et al. 2009) with motif length setting at 6–100 amino acids and number limit ≤30.

Identification of alternative splicing events

Alternative splicing analyses were performed as previously described (Finet et al. 2013). Briefly, alternative splicing events and the sequences of multiple isoforms were obtained from Phytozome (version 12.1) or Ensembl Genomes (release 97). The isoform exhibiting similar gene structure to that of Arabidopsis was selected as the canonical pattern, which was used as a reference for other isoforms to determine the patterns of alternative splicing. Specifically, the alternative splicing events were classified into six patterns: ①, lacking peptide in the non-domain region; ②, having supplementary peptide in the non-domain region; ③, lacking peptide in the domain region; ④, having supplementary peptide in the domain region; ⑤, lacking domain region; ⑥, occurring in the UTR region.

Protein-protein interaction network construction

Protein-protein interaction network was constructed using the STRING database (http://string-db.org) with default setting. Sources in STRING include experimentally determined interactions, curated databases, and information of co-expression, fusion, text mining and co-occurrence (Szklarczyk et al. 2019).

Expression profile analysis

The expression profiles of BPA genes in different developmental stages, biotic and abiotic stress conditions were retrieved from the Genevestigator database (https://genevestigator.com/gv/). For the developmental stages, the raw expression values were log2 transformed. For biotic and abiotic treatments, expression was indicated as fold-change relative to a control treatment. Heatmap was generated using the HemI software (Deng et al. 2014).

Yeast two-hybrid (Y2H) assay

Y2H assay was done as previously described (Luban and Goff 1995). Briefly, BPA1, BPL1, 2, 3, 4, 5 and 6 coding regions were PCR-amplified by using A. thaliana cDNA as templates with the reported primers (Li et al. 2019). PCR was performed in a reaction volume of 50 μL containing 10 μL 5 × PsBuffer, 200 μM each of dNTPs, 0.2 μM primers, 1.25 U of PrimeStar polymerase and 50 ng template DNA. Then the corresponding PCR products were cloned into both pGBKT7-BD and pGADT7-AD vectors. The Y2H assay was performed using the Gold Yeast Two-Hybrid System (Clontech). Yeast cells were co-transformed with the indicated plasmid combinations. Transformed cells were selected using the synthetic dropout (SD/−Leu/−Trp) medium and transferred to the SD/−Leu/−Trp/−His/−Ade selective medium for growth analysis. The BD-53 and AD-T were also co-transformed as a positive control, while BD-Lam and AD-T were co-transformed as a negative control. All Y2H experiments were repeated three times independently.