Genome-Wide Identification and Characterization of the RWP-RK Proteins in Zanthoxylum armatum

1. Introduction

Gene transcription is activated or repressed by transcription factors, which bind to specific cis-elements of the promoter [1]. The RWP-RK protein family, with a conserved RWP-RK domain, is a class of plant-specific transcription factors that is expressed in green algae, mosses, and vascular plants [2]. It is divided into two subfamilies: NLPs (NIN-like proteins) and RKDs (RWP-RK domain proteins) [3]. Unlike RKDs, all members of the NLP subfamily possess the PB1 (Phox and Bem 1) domain at their C-termini [4].

Two subfamilies of RWP-RK have been to participate in the modulation of various processes in plants. Members of the RKD subfamily are primarily involved in embryonic development [5]. For example, AtRKD4 deficiency may lead to inhibition of zygote elongation and disrupt early cell division patterns [6]. Overexpression of AtRKD1 and AtRKD2 was found to induce the formation of egg-like structures [7]. Members of the NLP subfamily are also linked to the nitrogen response [8]. The AtNLP7 receptor directly binds nitrate in intracellular environments [9]. Loss of NLP2 might trigger a decreased in nitrogen fixation and nitrogen content in plants [10]. Moreover, evidence suggests that RWP-RK proteins modulate plant responses to abiotic stress, such as heat stress [11].

Zanthoxylum armatum DC., also known as ‘Tengjiao’ or ‘Qinghuajiao’, belongs to the rutaceae family and is widely distributed in southwest China [12,13]. The fruit of Z. armatum can be utilized to prepare food spice, medicine, and oil, highlighting its significant economic value [14]. In addition, Z. armatum, like citrus, exhibit apomixis traits [15]. CitRWP, a gene that encodes a protein with a RWP-PK domain, was associated with apomixis in citrus [16]. To date, analysis of the RWP-RK family has only been conducted for a few species such as soybean, tea, and Brassica napus [5,17,18]. However, and no family analysis of RWP-RK proteins in Z. armatum has been reported. Recently, the genome of the Z. armatum was published allowing researchers to analyze various gene families [15,19].

In this study, 36 ZaRWP-RK transcription factors were identified from the Z. armatum genome. Chromosomal localization, phylogenetic relationships, physicochemical properties, gene structure, and cis-element of the ZaRWP-RK transcription factors were explored. Furthermore, the expression levels of ZaRWP-RK in different tissues and protein-protein interaction predictions were conducted which showed that some ZaRWP-RK transcription factors may regulate apomixis process. Therefore, this study provides a theoretical basis for further investigations on RWP-RK genes and their regulatory mechanisms in the apomixis process in Z. armatum.

2. Materials and Methods

2.1. Identification and Physicochemical Properties of ZaRWP-RK Genes

The data of Z. armatum genome was obtained from figshare (https://figshare.com/articles/dataset/Genome_Data_of_Chinese_pepper/20217635) as described by Hu et al. [19]. A hidden markov model (HMM) file of the RWP-RK domain (PF02042) was downloaded from the Pfam database (http://pfam-legacy.xfam.org/), and the candidate ZaRWP-RK proteins were identified using HMMER 3.0 software. All candidate ZaRWP-RK proteins were examined by NCBI CDD (https://www.ncbi.nlm.nih.gov/cdd/) and Pfam (http://pfam.xfam.org/). The number of amino acids, molecular weight (Mw), and isoelectric point value (pI) of ZaRWP-RK proteins were assessed and screened using the ExPasy tools (http://web.expasy.org/protparam/). The subcellular localization of ZaRWP-RK proteins was predicted using WoLF PSORT (https://wolfpsort.hgc.jp/).

2.2. Chromosomal Location and Gene Duplication Events

The chromosomal location of ZaRWP-RK genes was identified on the Z. armatum genome and visualized with TBtools-II software [20]. The gene duplication events of ZaRWP-RK genes were analyzed using Multiple Collinearity Scan toolkit (MCScanX) [21]. The synteny analysis of the ZaRWP-RK genes between Z. armatum and other species (Arabidopsis thaliana, Oryza sativa, and Citrus grandis) was visualized using TBtools-II software [20].

2.3. Phylogenetic analysis of ZaRWP-RK genes

Alignment of RWP-RK protein sequences from Z. armatum, Arabidopsis thaliana, Oryza sativa, and Citrus grandis was achieved using the MAFFT multiple-alignment software [22]. Phylogenetic trees were constructed using the Maximum likelihood (ML) method in raxmlGUI 2.0 with 1000 bootstrap replicates [23]. They were also visualized with the online tool EvolView (https://www.evolgenius.info/evolview-v2/).

2.4. Conserved Structural Domains, Conserved Motifs, and Cis-element Analyses of ZaRWP-RK genes

Conserved structural domains of ZaRWP-RK proteins were identified using the online tool CD-Search (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) [24]. Subsequently, conserved motifs were determined using the online program MEME (https://meme-suite.org/meme/tools/meme). Promoter regions 2000 bp upstream of the start codon (ATG) were subjected to cis-elements analyses using PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/). Next, the conserved structural domains, conserved motifs, and cis-element analyses of ZaRWP-RK genes were visualized using TBtools-II software [20].

2.5. Plant Materials, RNA Extraction and qRT-PCR Analysis

Different tissues (female flowers, young fruits, stems and leaves) of Z. armatum were collected from Juzi town in Leshan of Sichuan Province, China (29.32°N, 103.43°E). Total RNA samples were extracted using the FastPure ® Plant Total RNA Isolation Kit (Polysaccharides and Polyphenolics-rich) (Vazyme, Nanjing, China) according to the manufacturer’s instructions. 1 µg of total RNA was then used to synthesize cDNA using ABScript III RT Master Mix for qPCR with gDNA Remover (ABclonal, RK20428). The cDNA was subjected to qRT-PCR on a QX96M quantitative PCR system (JLM, Chengdu, China). The reaction protocol was as follows: one cycle at 95◦C for 30 s; 40 cycles at 95◦C for 10 s, 60◦C for 30 s, and melting curve generated automatically. All data presented was from at least three repeats. The results were analyzed using IBM SPSS Statistics 27 software. The primers for qPCR were designed using primer premier 5.0 (see Supplementary Table S2).

2.6. Protein-Protein Interaction Network Prediction of ZaRWP-RK Proteins

The ZaRWP-RK proteins from Citrus sinensis that were homologous to the ZaRWP-RK proteins in Z. armatum were submitted to the STRING 12 (https://cn.string-db.org/) to construct a protein-protein interaction network.

3. Results

3.1. Genome-Wide Identification and Physicochemical Properties of ZaRWP-RK Genes

Based on the RWP-RK domain (PF02042) in the Pfam database, 36 ZaRWP-RK genes were identified from the Z. armatum genome. The ZaRWP-RK genes were then named from ZaRWP-RK1 to ZaRWP-RK36, and two ZaRWP-RKs (ZaRWP-RK35 and ZaRWP-RK36) were not annotated to chromosomes based on the order of their chromosomal locations (Figure 1). Details of the physical and chemical properties of the ZaRWP-RK genes are presented in Table S1. The coding sequences (CDS) of ZaRWP-RK genes ranged from 438 bp (ZaRWP-RK12) to 3024 bp (ZaRWP-RK29). Similarly, ZaRWP-RK genes encode amino acids ranging from 145 aa (ZaRWP-RK12) to 1007 aa (ZaRWP-RK29), with molecular weights between 17.16 and 110.13 kDa. The isoelectric point (pI) distribution of ZaRWP-RK proteins varied from 4.56 (ZaRWP-RK17) to 9.57 (ZaRWP-RK12). Most of the ZaRWP-RK members were predicted to localize in the nucleus (32), and a few in the cytoskeleton (1), plasma membrane (1), cytosol (1), and chloroplast (1).

Next, we analyzed the duplication events of ZaRWP-RK genes (Figure 2). Ten pairs of segmental duplications were identified, while no tandem duplicated genes were found, indicating that segmental duplication likely had a substantial impact on the evolutionary trajectory of the ZaRWP-RK family. To gain further insights into the evolutionary relationships among the 36 ZaRWP-RK genes between Z. armatum and other plants, three collinearity maps were constructed using the RWP-RK members of Z. armatum, Arabidopsis thaliana, Oryza sativa, and Citrus grandis (Figure 3). The number of orthologous gene pairs between Z. armatum and A. thaliana, Oryza sativa, and Citrus grandis was 48, 11, and 35, respectively. Among the three species, ZaRWP-RK genes are most closely related to AtRWP-RK genes.

3.2. Phylogenetic Relationship, Conserved Structural Domains and Motifs Analysis of ZaRWP-RK proteins

To explore the phylogeny of the RWP-RK proteins, a phylogenetic tree was constructed based on the amino acid sequences from 36 Z. armatum ZaRWP-RK proteins, 14 Arabidopsis AtRWP-RK proteins, 13 rice OsRWP-RK proteins, and 8 pummelo CgRWP-RK proteins using raxmlGUI 2.0 (Figure 4). The ZaRWP-RK proteins were clustered into two subfamilies (NLP and RKD) using the results of A. thaliana as a reference [25]. Among them, NLP subfamily contained 21 members (ZaNLP1 - ZaNLP21), and RKD subfamily contained 15 members (ZaRKD1 - ZaRKD15) (Figure 4, Table S1).

To determine the relationships among ZaRWP-RK family members, the phylogeny, conserved structural domains and motif information were analyzed (Figure 5). Analysis of the protein conserved motifs led to identification of 10 motifs in ZaRWP-RK proteins, ranging from 21 aa to 50 aa in length (Figure 5B,D). Motifs 2-6 and motifs 8-10 were only found in the ZaNLP subfamily members. It appeared that motif 1 and motif 7 were the most important motifs in ZaRWP-RK family members given that it was identified in almost all proteins. In addition, five conserved structural domains were found in ZaRWP-RK family members. RWP-RK domains were present in almost all members of the ZaRWP-RK family, with the exception of ZaNLP5. The PB1 domain (including PB1_NLP and PB1 superfamily) was only identified in ZaNLP proteins, suggesting that the PB1 domain at ZaNLP proteins was highly conserved.

3.3. Cis-Element Analysis of ZaRWP-RK Promoters

To investigate the regulatory mechanisms of ZaRWP-RK genes, promoter regions 2000 bp upstream of the start codon (ATG) were subjected to cis-elements analyses (Figure 6, Table S2). Based on the results of cis-acting regulatory element analysis, 21 types of hormone-responsive elements were identified and classified into three categories, including abiotic stress, response to phytohormone, and plant growth and development. Among these elements, most of them were phytohormone-responsive elements. Research has demonstrated a close association of abscisic acid (ABA), auxin, and gibberellin with embryonic development. In terms of hormone response, ABRE (58) elements were the most prevalent and were identified as ABA response elements. Additionally, numerous auxin-responsive elements (14 TGA-element and 4 AuxRR-core) and gibberellin-responsive elements (12 P-box, 8 GARE-motif, and 7 TATC-box) were found in some ZaRWP-RK genes. Moreover, among the plant growth and development-related elements, unnamed__1 elements were the most abundant (42, 60k protein binding site) followed by CCAAT box (16, MYBHv1 binding site). Interestingly, some cis-elements were specific to certain ZaRWP-RK genes, such as MSA-like (ZaRKD2), DRE (ZaNLP10 and ZaRKD10), HD-Zip 3 (ZaNLP13) and Box III (ZaRKD1, ZaRKD11, and ZaRKD12). These results suggested that the ZaRWP-RK genes might be regulated by hormones and some transcription factors, thereby modulating embryonic development.

3.4. Tissue-Specific Expression of ZaRWP-RK Genes

To investigate the tissue-specific patterns of ZaRWP-RK genes, their expression profiles was explored in four different tissues of Z. armatum (female flower, young fruit, stem, and leaf) using qRT-PCR (Figure 7). It was observed that most ZaRWP-RK genes were highly expressed in young fruits, suggesting that these genes participate in the regulation of apomixis. Among them, ZaNLP6, ZaNLP10, ZaNLP18, and ZaRKD3 were significantly upregulated in the young fruit, whereas ZaRKD5 was highly expressed in female flower and young fruit, indicating their potential role in apomixis. In contrast, the ZaNLP7 expression was significantly lower in the female flower and young fruit than in the stems and leaf, suggesting it might regulate apomixis in Z. armatum plants and perform negative regulatory function.

3.5. Correlation Analysis of Apomixis-Related ZaRWP-RK Genes

Previous research has demonstrated that citrus apomixis may be caused by the expression of a gene known as CitRWP (Cg4g018970), which encodes a protein with an RWP-PK domain [16]. Similarly, OsRKD3 can regulate somatic embryogenesis in black rice [26]. According to the results of phylogenetic tree analysis, four genes (ZaRKD2, ZaRKD8, ZaRKD9, and ZaRKD15) in Z. armatum were most closely related to CitRWP, and three genes (ZaRKD10, ZaRKD13, and ZaRKD14) were most closely related to OsRKD3 (Figure 4). Zanthoxylum is a close relative of citrus in the Rutaceae family [27]. To further investigate the functions of these genes, the STRING database was employed to construct a protein-protein interaction network of these ZaRWP-RK proteins based on their CsRWP-RK homologs (Figure 8). Cs_ont_4g023450 (homolog of ZaRKD2, ZaRKD10 and ZaRKD14) and Cs_ont_5g040080 (homolog of ZaRKD8 and ZaRKD9) were directly connected to ATG5, ATG12, AP2A1, AP2S1, CHC, FIS1. Moreover, Cs_ont_5g040080 (homolog of ZaRKD8 and ZaRKD9) directly interacted with Alba3.

4. Discussion

The RWP-RK transcription factor is widely expressed in many plants and regulates various physiological processes such as nitrogen response, gametophyte development, and abiotic stress regulation [28]. In many plants, including Arabidopsis, rice, elephant grass and soybean, members of the RWP-RK transcription factor family have been reported [5,25,29,30]. However, no study has analyzed this family in Z. armatum. In this study, we performed the genome-wide identification of the RWP-RK family members in Z. armatum and analyzed their expression during reproductive processes, based on the recently published genome information of Z. armatum [19].

Using the conserved motif (RWPXRK) of RWP-RK transcription factor, 36 members of the RWP-RK family in Chili pepper genome were identified, which surpassed that reported in other plants, for example, 14 in Arabidopsis, 13 in rice, and 8 in pummelo (Figure 1 and Figure 4). Z. armatum, having undergone a whole-genome duplication, has higher number of chromosomes (2n = 4x = 132) [15,19]. This indicates that the ZaRWP-RK genes have extensively expanded throughout the evolutionary process, probably due to the diversity of their functions. Gene duplication is a significant driving force for the diversification of gene functions and species evolution [31]. In eukaryotes, there are two main types of duplication events: tandem and segmental [32]. Previous studies have shown that the NLP subfamily in Arabidopsis evolved through segmental duplication rather than tandem amplification [25]. Similarly, we found that the ZaRWP-RK family in Z. armatum evolved primarily through segmental duplication (Figure 2). Moreover, we observed that the number of orthologous gene pairs between Z. armatum and dicots (Arabidopsis and pummelo) exceeded that between Z. armatum and monocots (rice), indicating that after the divergence of monocots and dicots, the ZaRWP-RK transcription factors have undergone extensive evolution and duplication (Figure 3).

The RWP-RK transcription factor family is generally divided into two subfamilies: the RKD subfamily and the NLP subfamily [25]. The available studies indicate that in Arabidopsis, the RKD subfamily has 5 members and the NLP subfamily has 7 members [29,33]. In rice, there are 7 members in the RKD subfamily and 6 members in the NLP subfamily [25]. In this study, we analyzed 36 ZaRWP-RK transcription factors from the Z. armatum genome. Through a phylogenetic analysis, we classified 15 of these members into the RKD subfamily and 21 members into the NLP subfamily (Figure 4, Table S1). Notably, the transcriptional activity, DNA-binding activity, and protein-protein interactions of transcription factors are often controlled by their motifs [34]. Among the 10 conserved motifs identified in ZaRWP-RK family, motifs 1 and 7 were expressed in all ZaRWP-RK members, indicating that they might constitute the most conserved part of ZaRWP-RK (Figure 5B). Motifs 2-6 and 8-10 were only identified in members of the NLP subfamily, suggesting that these motifs might plat important role in unique biological functions of the NLP subfamily. Analysis of the conserved domains of ZaRWP-RK showed that the PB1 domain was present in all NLP subfamily members, which is the main distinguishing point between RKD and NLP subfamily (Figure 5C). Interestingly, our analyzes showed that ZaNLP2 contained a conserved GAF_2 domain. A recent study reported that nitrate-triggered ROS signaling and the detection of nitrate deficiency are dependent upon the GAF domain of AtNLP7 [35]. Therefore, we speculated that ZaNLP2 might be involved in the regulation of nitrogen response. In addition, the PoIC superfamily domain was detected in ZaRKD10. PoIC is a major DNA polymerase involved in genome duplication, suggesting that ZaRKD10 might participate in rapid cell division stage, such as embryonic development [36].

The activation of gene expression is regulated by Cis-elements present in the promoter region [37]. Therefore, the physiological functions of cis-elements on the promoter can be used to predict the biological function of the genes [38]. In ZaRWP-RK family, 21 types of cis-elements were identified, categorized into abiotic stress (6), response to phytohormone (9), and plant growth and development (6) (Figure 6, Table S2). Most of them were hormone response elements, with the ABRE (ABA response element) accounting for the largest proportion. ABA are key hormones that control plant growth, embryogenesis, seed physiology reactions, and stress tolerance [39]. Pretreatment with ABA in embryo maturation medium increased cassava somatic embryo conversion into plants [40]. In addition, a significant number of auxin response elements and gibberellin response elements were identified. Auxin and gibberellin have been shown to regulate plant embryo development [41,42]. Furthermore, numerous elements related to plant growth and development were identified in the promoters of ZaRWP-RK members, suggesting that ZaRWP-RK transcription factors participate in embryogenesis. Consistent with this view, numerous studies have shown that ZaRWP-RK modulates plant reproduction processes [7]. Among the five RKD members in Arabidopsis, AtRKD1 and AtRKD4 are highly expressed during early embryo development, and the AtRKD4 mutant leads to abnormal embryo development [6,29]. The MpRKD controls the development of gametophytes and contributes to the formation of gemma cups [43]. Therefore, our results indicate that ZaRWP-RK transcription factors may participate in plant reproductive processes by responding to and regulating hormones such as ABA, auxin, and gibberellin.

Apomixis, found in over 300 genera and 40 families throughout the plant kingdom, displays a wide taxonomic distribution [44]. In this mechanism, the fertilization step is skipped, allowing cells from nucellus or integument tissues to develop into somatic embryos. Previous studies have highlighted the crucial role of RWP-RK transcription factors in both somatic embryogenesis and apomixis processes [15,26,45]. Therefore, investigating the expression pattern of the ZaRWP-RK transcription factors during the apomictic reproduction process in Z. armatum is of great importance. In this study, we found that the expression levels of most ZaRWP-RK genes increase from flowers to young fruits significantly (Figure 7). Among them, the expression levels of ZaNLP6, ZaNLP10, ZaNLP18, and ZaRKD3 and ZaRKD5 show the most significant differences, indicating their potential functions in apomixis. In contrast, ZaNLP7 expression was significantly lower in the female flower and young fruit than in the stems and leaf, suggesting it might negatively regulate apomixis.

Autophagy proteins are required for multiple functions during embryogenesis [46]. Our predictions based on the protein-protein interaction network suggested that Cs_ont_4g023450 (homolog of ZaRKD2, ZaRKD10 and ZaRKD14) and Cs_ont_5g040080 (homolog of ZaRKD8 and ZaRKD9) may interact with ATG5 and ATG12, suggesting that ZaRKDs may be involved in autophagy-regulated embryonic development (Figure 8). In addition, our results showed that Cs_ont_4g023450 and Cs_ont_5g040080 may interact with AP2A1, AP2S1, and CHC1 (Figure 8). The adaptor protein complex 2 (AP-2) proteins participate in various biological processes, such as plant growth and development, and stresses response [47]. AP2σ and AP1/2β modulate embryogenesis and plant growth in Arabidopsis [48,49]. AP-2 has been shown to interact with CHC and regulate flower organ development. Therefore, we speculate that ZaRKDs may interact with AP2A1, AP2S1, and CHC1, thereby influencing the reproductive process in Z. armatum. Moreover, we also found that Cs_ont_4g023450 and Cs_ont_5g040080 interacted with MIS1 and Alba3. The mitochondrial fission complex has been implicated in Arabidopsis heat tolerance [50]. ALBA proteins modulate various stress responses [51]. Recent studies also have shown that overexpression of one RWP-RK gene in pearl millet can significantly enhance plant heat tolerance [11]. Therefore, it can be hypothesized that ZaRKDs may regulate heat tolerance in Z. armatum.

5. Conclusions

In this study, we identified 36 ZaRWP-RK transcription factors in the genome of Z. armatum. Most ZaRWP-RK genes were distributed on 24 chromosomes expect ZaRWP-RK35 and ZaRWP-RK36. The ZaRWP-RK genes were into two subfamilies: NLP and RKD, through phylogenetic analysis. The conserved structural domain and motif analysis revealed the functional similarity and specificity of these ZaRWP-RK genes. Meanwhile, cis-element analysis showed that the ZaRWP-RK genes had multiple functions in plant reproductive processes. Analysis of gene expression in different tissues and protein-protein interaction network further highlighted the potential regulatory role of ZaRWP-RK genes in the apomixis of Z. armatum. In summary, through genome-wide identification and characterization of the RWP-RK transcription factor in Z. armatum, our results not only improve our understanding of the ZaRWP-RK family but also offer theoretical support for further studies into the regulatory mechanism of apomixis in Z. armatum.