Genome-wide identification and expression patterns of the aspartic protease gene family in Epimedium pubescens

Abstract

Background

Aspartic proteases (APs), proteolytic enzymes involved in protein maturation, degradation, and signaling, are found in various organisms. The plant Epimedium pubescens is known for its pharmacologically active flavonoids and its use in traditional Chinese medicine. Despite this, to date, the AP gene family in this species has not been functionally analyzed. This study aimed to uncover the roles of AP genes in E. pubescens (EpAPs), focusing in particular on their involvement in light stress responses.

Results

Genome-wide analysis identified and characterized a total of 103 EpAPs, which were categorized into four phylogenetic groups and revealed conserved motifs crucial for their catalytic function. Structural analysis highlighted the diversity of intron-exon arrangements and the predominant role of tandem duplication in gene expansion. Promoter analysis showed an enrichment of light-responsive elements, indicating potential involvement in light stress responses. Tissue-specific expression patterns revealed specialized roles in various organs, whereas several EpAPs exhibited stage-specific expression during the formation of abscission zones. The analysis of protein–protein interactions identified links to reproductive development, programmed cell death, and stress responses. Under light stress, selected AP genes exhibited dynamic changes in expression, with some showing transient upregulation or recovery phases, which suggests their involvement in short-term adaptation or sustained light stress responses.

Conclusions

This study provides the first comprehensive analysis of AP genes in E. pubescens, highlighting their potential roles in development and stress adaptation. The presence of light-responsive elements and changes in expression under light stress suggest that AP genes may serve as key regulators of environmental responses in this species. Further validation studies could inform strategies to improve light stress resistance in shade-adapted plants.

Introduction

Aspartic proteases (APs; EC 3.4.23) are proteolytic enzymes found in various organisms, including bacteria, fungi, animals, and plants. Their catalytic activity relies on highly conserved Asp–Thr/Ser–Gly (DTG/DSG) motifs, which enable their functions in protein maturation, degradation, and signaling [1, 2, 3, 4]. Within the MEROPS protease database, APs are classified into 17 families based on their structural characteristics and evolutionary origins and are further grouped into five following clans: AA, AB, AC, AD, and AF [5]. Most plant APs belong to the A1 family of Clan AA [4]. The APs of this family, commonly referred to as pepsin-like APs, exhibit catalytic activity under acidic conditions and are inhibited by gastric protease inhibitor A [2, 6]. However, this is not the case for all plant APs; for example, the Arabidopsis CDR1 protein displays peak activity at higher pH values ranging from 6.0 to 6.5 [7].

In 2004, 59 A1 members were identified in Arabidopsis thaliana and classified into groups A1-1 through A1-5 [8]. In 2005, another study identified 51 A1 members and divided them into three categories, namely typical APs, nucellin-like APs, and atypical APs [9]. Subsequently, in 2008, 69 AP genes were identified in Arabidopsis, 63 of which belonged to the AP family and 6 were classified as non-peptidases [10]. Phylogenetic and structural analyses by Duan et al.. further classified 70 APs of the A1 family in A. thaliana into four distinct groups (Groups I–IV), with Group IV representing a conserved group in Spermatophyta that lacks the classical DTGS/DSGT motif [11]. Moreover, studies have identified 104 A1 APs in rice (Oryza sativa) [12], 50 in grape (Vitis vinifera) [13], 67 in black cottonwood (Populus trichocarpa) [14], and 129 in Moso bamboo (Phyllostachys edulis) [15], highlighting the expansion and functional diversity of AP genes in plants.

Plant APs play multifaceted roles in plant growth, development, and stress responses. In vegetative tissues, they are crucial for protein processing and turnover. For instance, phytepsin, a vacuolar AP in barley, participates in the autolysis of tracheary elements and sieve cells [16, 17]. During leaf senescence in tobacco, CND41, an AP found in the chloroplast, facilitates the turnover of ribulose-1,5-bisphosphate carboxylase/oxygenase (Rubisco) [18]. During seed germination, ASPARTIC PROTEASE IN GUARD CELL 1 (ASPG1) in A. thaliana breaks down storage proteins, assisting in seedling development [19]. APs also play critical roles in plant reproduction by regulating programmed cell death (PCD) in tapetal cells and participating in pollen development. In rice, OsAP25 and OsAP37, which are regulated by the transcription factor ETERNAL TAPETUM1 (EAT1), control PCD in tapetal cells, a function that is essential for pollen viability [20], whereas OsAP65 contributes to the biosynthesis of key compounds for pollen germination and tube growth [21]. In A. thaliana, the AP known as PROMOTION OF CELL SURVIVAL 1 (PCS1) regulates cell fate during embryonic development and reproduction [22], while A36 and A39 may be involved in gametogenesis and pollen guidance [23]. Studies have suggested that plant APs also play important roles in defense against biotic and abiotic stresses. APCB1 is a nucellin-like AP that induces Bcl-2-ASSOCIATED ATHANOGENE 6 (BAG6) to initiate autophagy and defense mechanisms [24]. An extracellular AP in A. thaliana and rice known as Constitutive disease resistance 1 (CDR1) enhances immunity through peptide signaling pathways [25, 26]. Furthermore, overexpression of VqAP13 in grape has been shown to improve resistance to powdery mildew by modulating the salicylic acid and methyl jasmonate signaling pathways [27]. In A. thaliana, ASPG1, which is induced by abscisic acid (ABA), plays a role in enhancing drought resistance through ABA-dependent signaling pathways [28]. Additionally, the overexpression of APA1 in A. thaliana confers greater tolerance to drought stress than that exhibited by wild-type plants [29].

Despite significant progress being made in understanding APs in model species, their roles in non-model plants remain underexplored. Epimedium pubescens, a medicinal herb belonging to the Berberidaceae family, is rich in pharmacologically active flavonoids such as icariin and epimedin A, B, and C [30, 31]. This species is widely used in traditional Chinese medicine to treat kidney-yang deficiency and rheumatism [32]. Although it is a shade-loving plant, recent studies have demonstrated that controlled light enhancement can significantly increase its flavonoid content. A study on E. sagittatum has shown that total flavonoid content peaked at 40–100 µmol·m⁻²·s⁻¹, was lowest at 20–35 µmol·m⁻²·s⁻¹, and declined above 160 µmol·m⁻²·s⁻¹ due to inhibited synthesis [33]. In E. pseudowushanense, flavonoid content has also been shown to vary with light levels, with the highest epimedin C yield being recorded at 54.6 ± 2.5 µmol·m⁻²·s⁻¹ [34]. In contrast to these two species, E. pubescens exhibits higher tolerance to strong light. A specific study has revealed that as light intensity increased from L1 (36.4 ± 5.0 µmol·m⁻²·s⁻¹) to L4 (200.2 ± 5.0 µmol·m⁻²·s⁻¹), both branch and leaf numbers increased significantly. Moreover, the content of Epimedium-derived flavonoid glycosides also increased with light intensity, peaking under L3 (145.6 ± 5.0 µmol·m⁻²·s⁻¹), suggesting that this level is optimal for maximizing the accumulation of medicinal compounds [35]. Some studies have found that certain APs, such as NANA [36] and APP-A1 [37], play important roles in photosynthesis. Therefore, the present study focused on the AP genes in E. pubescens (EpAPs) to explore their roles in light stress responses. Specifically, EpAPs were identified and their phylogeny, structure, conserved motifs, duplication, and expression patterns were analyzed. Notably, promoter analysis revealed a substantial number of light-responsive elements, prompting us to examine their expression under light stress via real-time quantitative PCR (qPCR) to better understand their potential functions. The findings obtained contribute to the identification of key genes that could be used to enhance light stress resistance in E. pubescens.

Results

Genome-wide identification and characterization of AP genes in E. pubescens

A total of 109 and 143 candidate AP proteins in E. pubescens were initially identified using HMMER v3.4 and BLASTP searches against AP proteins in A. thaliana, respectively. After Batch-CDD domain verification, 103 EpAPs were confirmed (Table S1) and further validated using the SMART and InterProScan tools (Figure S1, Table S2). These genes were subsequently renamed EpAP1–EpAP103 based on their chromosomal scaffold positions.

The predicted EpAP proteins exhibited diverse characteristics. Their amino acid lengths ranged from 211 to 795, with molecular weights ranging from 22,810.02 to 86,282.81 Da. The isoelectric points (pIs) varied between 4.46 and 9.63. The instability index varied from 24.82 to 53.48, with 63 proteins being classified as stable (instability index values below 40), while the rest exhibited slightly higher values for this index. Overall, 63 proteins exhibited hydrophilicity, as indicated by negative GRAVY scores, whereas the remaining 40 were classified as hydrophobic.

Subcellular localization predictions using WoLF PSORT indicated that the EpAP proteins were likely localized primarily in the chloroplast (40), followed by extracellular space (25), and vacuoles (17), while DeepLoc predictions revealed a predominant extracellular localization (73), followed by localization to cell membrane (9) and lysosome/vacuole (11) (Table S3). These differences are likely due to variations in the training datasets and algorithms of the two prediction tools, affecting in particular proteins with dual localization potential, such as those predicted to localize both in the chloroplast and extracellular space.

Phylogenetic analysis

A total of 173 AP proteins, 70 from A. thaliana and 103 from E. pubescens, were used to construct a phylogenetic tree, where they clustered into four distinct groups (Groups I–IV) (Fig. 1). Group III, with 77 members, represented the largest group, potentially reflecting gene expansion or functional diversification. In contrast, Groups I, II, and IV contained 2, 12, and 12 members, respectively.

Fig. 1

Phylogenetic tree based on AP protein sequences in E. pubescens and A.thaliana

Analysis of gene structure, conserved domains, and motifs

The gene structure and motif composition of the examined EpAPs varied significantly among groups (Fig. 2). Group I contained motifs 10, 1, 7, 6, 3, and 2, while Group II exhibited motifs 10, 1, 6, 4, 3, and 2, with motifs 4, 7, or 9 being present in most members. Group III was highly diverse and the only one exhibiting motif 8. Group IV predominantly contained motifs 6, 4, 3, and 2. Motifs 1 and 6, each 15 amino acids in length, were the shortest and most conserved (Fig. 3). Motif 1 was shown to encode the AP active site DTGS, a signature feature of APs that is essential for enzymatic activity, and was identified in 93 EpAPs. Motif 3, comprising 29 amino acids, was one of the most conserved, identified in 101 EpAPs. Another highly conserved motif, also present in 101 genes, was motif 2, which likely contributes to protein stability and enzymatic integrity. In contrast, motif 5 was the least conserved, being found in only 32 EpAPs and predominantly within Group III, highlighting a potential role in subgroup-specific functions.

Fig. 2

Phylogenetic tree, gene structure, protein conserved motifs, and conserved domain analysis of the EpAP gene family. (A) Phylogenetic tree constructed by IQtree in TBtools with 1000 bootstrap number. (B) Analysis of motifs in the proteins encoded by EpAP genes. (C) Conserved domains of EpAPs predicted by Batch CD-Search in NCBI. (D) Gene structure of EpAP genes

Fig. 3

Motif expression pattern in EpAP proteins. (A) Conserved motif logos by TBtools. (B) Length and distribution of motif

All EpAP protein sequences contained the AP domain (PLN03146), although some were not displayed due to structural overlap (Fig. 2C). According to the results obtained via InterProScan, all 103 EpAP proteins possessed a Peptidase Family A1 domain. Additionally, Group I members EpAP103 and EpAP90 contained a Saposin-like domain. SMART domain analysis revealed that certain EpAP proteins (e.g., EpAP12, EpAP28, EpAP35, and EpAP48) lacked an N-terminal signal peptide. Furthermore, other members, including EpAP16, EpAP24, EpAP35, and EpAP40, contained transmembrane domains at either the N- or C-terminus.

The analysis of gene structure revealed intron numbers ranging from 1 to 15. A total of 48 EpAPs (46.6%) contained a single exon and no introns, whereas 32 had two exons and one intron, and 23 contained more than two introns (Table S4). Untranslated regions (UTRs) were absent in 78.6% of the genes, whereas 4.9% had one UTR, and 16.5% had more than one UTRs. The absence of UTRs in most EpAPs suggests that post-transcriptional regulation may play a significant role in controlling their expression, which is consistent with the stress-responsive roles of these genes.

Chromosomal locations, genomic duplications, and Ka/ks ratios

The chromosomal distribution of EpAPs revealed a non-uniform pattern across six chromosomes (Fig. 4A), with chromosome 4 (Chr4) harboring the highest number of genes (30 in total; EpAP50–EpAP79), predominantly clustered in the upper and middle regions. In contrast, Chr5 exhibited the lowest number, with only six evenly distributed genes (EpAP80–EpAP85). Chr1 and Chr2 each contained 21 genes, wherein they were clustered in the upper regions of Chr1 (EpAP1–EpAP21) and central region of Chr2 (EpAP22–EpAP42). Chr3 contained seven genes (EpAP43–EpAP49) scattered on the chromosomes, whereas Chr6 contained 18 genes (EpAP86–EpAP103), mainly clustered at its terminal region. The abundance of tandem repeats within the examined EpAPs contributed to their diversity.

Fig. 4

Distribution of 103 EpAP genes on chromosomes and collinearity analysis. (A) EpAP genes on chromosome scaffolds visualized by TBtools v2.110. The central gray lines represent the collinearity of the E. pubescens genes, whereas the cyan line indicate the collinearity of the EpAP gene. a, chromosome scaffolds; b, gene density; c, GC ratio; d, GC skew; e, N ratio; f, localization of EpAPs. (B) Synteny analysis of APs between E. pubescens and A. thaliana and Coptis chinensis

Our analysis identified one pair of collinear EpAPs (EpAP1 and EpAP26; Ka/Ks = 0.179) (Fig. 4A, denoted by the green line) and 17 pairs of tandemly repeated EpAPs with Ka/Ks ratios ranging from 0.10 to 1.12 (Table S5). The majority of gene pairs (16/17) exhibited Ka/Ks ratios < 1, indicating that these genes have primarily undergone purifying selection during evolution. Notably, the gene pair EpAP70/EpAP71 (Ebr02G018630.1/Ebr02G018640.1) had a Ka/Ks ratio slightly greater than 1 (1.125), suggesting that this duplication event may have experienced a brief period of positive selection or functional divergence. The low Ks values (0.045–1.710) reflected recent tandem repeat events that are potentially related to recent genomic expansion events in E. pubescens. In the comparison of collinear genes between E. pubescens and A. thaliana, all gene pairs showed Ka/Ks ratios significantly lower than 0.5 (0.104–0.419), indicating that AP genes exhibit a high degree of functional conservation in core eudicots (Fig. 4B). In particular, the Ka/Ks ratio between EpAP103 and AT1G11910.2 (APA1) was only 0.105. A total of 20 pairs of collinear AP genes were identified between E. pubescens and Coptis chinensis (Table S5). The Ka/Ks ratios of all gene pairs were significantly lower than 1 (range: 0.076–0.381), indicating that these genes were subjected to strong purifying selection during evolution. Among them, EpAP37 (Ka/Ks = 0.077) and EpAP38 (Ka/Ks = 0.076) showed extremely low ratios, suggesting that they may have assumed highly conserved physiological functions after the two species diverged.

Identification of Cis-elements

Promoter analysis of the EpAPs revealed diverse cis-elements associated with environmental and hormonal responsiveness within the 2000-bp upstream regions of their transcription start sites (Fig. 5; Table S6). Notably, all 103 EpAPs contained light-responsive elements in their promoter sequences. Other commonly identified elements included the anaerobic response element (present in 97 genes), the low-temperature response element (68 genes), and TC-rich repeats associated with defense and stress reactions (41 genes). In addition, elements related to hormonal regulation were also prevalent, with 79 genes containing methyl jasmonate-responsive motifs such as CGTCA and TGACG. These findings suggest that EpAP genes may play roles in the responses to environmental stressors and regulation of hormonal signaling pathways.

Fig. 5

Cis-element analysis of EpAP genes. This analysis was conducted on the promoter regions of 2000 bp upstream of 103 EpAP genes to predict the types and numbers of cis-elements

Expression patterns of EpAP genes

The expression profiles of EpAPs across various tissues and during the formation of abscission zones indicated their involvement in distinct biological functions. Transcriptome data revealed that 12 EpAPs showed no detectable expression in any tissue, while other genes exhibited tissue-specific expression patterns or high expression levels (Fig. 6). For instance, EpAP54, EpAP87, and EpAP95 were predominantly expressed in the flowers, suggesting roles in reproductive development. EpAP40 and EpAP90 showed high expression levels in fruit, suggesting functions in fruit maturation or stress adaptation. Leaf-specific expression was observed for EpAP2, EpAP5, and EpAP27, while EpAP30, EpAP47, and EpAP43 were highly expressed in the roots, which suggests a potential role in nutrient uptake or environmental responses. Additionally, distinct expression patterns were observed in the three stages of abscission zone formation during fruit development (Fig. 7). Genes such as EpAP17, EpAP45, EpAP49, and EpAP59 were highly expressed prior to abscission, which is indicative of their involvement in the early development of abscission zones. In contrast, EpAP15, EpAP46, EpAP68, and EpAP73 were predominantly expressed during abscission, indicating a role in the active separation process. Following fruit detachment, genes like EpAP16, EpAP70, EpAP78, and EpAP93 showed elevated expression, which highlights their roles in cellular reorganization or recovery post-abscission.

Fig. 6

Expression patterns of EpAPs in different tissues of E. pubescens. The data were obtained from the transcriptome of fully developed leaves, roots, shoots, leaves, flowers, and fruits [56 ]. The color scale illustrates log2 expression values, ranging from blue for low expression levels to red for high expression levels

Fig. 7

Expression patterns of EpAPs in abscission zone of E. pubescens. EpA, pre-abscission stage (before abscission zone formation); EpB, abscission stage (abscission zone formation); and EpC, post-abscission stage (after abscission zone formation and pod detachment). The color scale illustrates log2 expression values, ranging from blue for low expression levels to red for high expression levels

Expression analysis of EpAPs at different periods after light stress

Given the abundance of light-responsive cis-acting elements in the promoter regions of EpAPs, it was hypothesized that these proteases might be involved in the plant’s response to light stress. To explore gene expression patterns under light stress, real-time quantitative PCR (qPCR) was performed for 12 specific EpAPs (EpAP9, EpAP11, EpAP19, EpAP20, EpAP37, EpAP38, EpAP42, EpAP68, EpAP78, EpAP85, EpAP88, and EpAP90). These genes were selected based on the results of subcellular localization predictions, which suggested their localization to chloroplast, making them relevant candidates for studying responses to light stress. The qPCR results are shown in Fig. 8. EpAP42 showed no response to light stress, with its expression remaining stable throughout the entire period. EpAP9 exhibited an increased expression at L3 and L6, followed by a decrease back to the initial level. The other genes showed a significant reduction in expression at L3. Among them, EpAP11, EpAP19, EpAP20, EpAP68, and EpAP78 displayed increased expression levels after 9 h of light stress (L9), but these decreased thereafter. However, only EpAP11 and EpAP19 exhibited expression levels at L9 that surpassed their pre-stress levels at L0. The expression levels of the other genes, while higher than those at L3, remained below or at the same level as before the treatment. Both EpAP85 and EpAP90 showed sustained suppression of expression under light stress.

Fig. 8

qPCR result of 12 selected EpAP genes. Error bars indicate Standard Error (SE) from three technical replicates. Asterisks or n.s. indicate significant differences between different time point, *p < 0.05, **p < 0,01, **** P < 0.0001, ns, no significant difference, by One-way ANOVAs

Prediction of Protein–Protein interactions (PPIs)

The observed PPIs highlighted the potential roles of EpAP proteins in critical biological processes (Fig. 9). For instance, EpAP64 and EpAP102 were predicted to interact with Egg cell-secreted protein 1.5 (EC1.5), a protein that mediates gamete interactions during double fertilization, suggesting their potential involvement in reproductive success. Several EpAPs, including EpAP16, EpAP46, and EpAP92, were predicted to interact with RDUF1 and RDUF2 (E3 ubiquitin-protein ligases), which are implicated in protein degradation pathways essential for stress responses and development. Additionally, EpAP20, EpAP24, and other EpAP proteins were predicted to interact with cathepsin B-like protease 3 (CATHB3), a key regulator of PCD in plants. These findings collectively indicate that EpAP proteins have diverse functional roles in development, stress adaptation, and cell death regulation.

Fig. 9

Prediction of the protein–protein interaction network of EpAP proteins via STRING. The interaction score was set to medium confidence 0.400

Discussion

Plant APs are critical enzymes involved in numerous biological processes, such as development, pathogen defense, and stress responses [38, 39]. Despite extensive research being conducted on APs in model plants, genome-wide identification and functional analysis of the AP gene family in E. pubescens have remained unexplored. The present study comprehensively analyzed the AP genes in E. pubescens, providing insights into their structural diversity, evolutionary conservation, and potential functional roles.

Tandem duplications are the main driver of EpAP gene expansion

In this study, a total of 103 EpAPs were identified and classified into four phylogenetic clades. Notably, proteins EpAP90 and EpAP103 in Group I were shown to contain Saposin-like domains, which play multifaceted roles in plant APs, primarily in intracellular sorting, host defense, and membrane interactions [40, 41]. Additionally, consistent with previous findings in A. thaliana [9], Moso bamboo [15], and potato (Solanum tuberosum) [42], some EpAP proteins possessed signal peptides and transmembrane domains. Interestingly, except for EpAP24 in Group III, all EpAP proteins with transmembrane domains (EpAP16, EpAP35, EpAP40, EpAP79, and EpAP80) belonged to Group II. Furthermore, most EpAP proteins contained xylanase inhibitor domains (TAXI_C and TAXI_N), which overlap with the Asp domain, consistent with previous reports [15, 42].

In E. pubescens, the expansion of the AP gene family is primarily driven by tandem duplication, which is particularly prominent in the members of Group III. Notably, the expansion patterns of AP gene families vary significantly among plant species. For example, in P. trichocarpa [14] and P. edulis [15], segmental duplication is the primary driver of AP gene family expansion, whereas in S. tuberosum, tandem duplication is more prevalent [42]. This difference may reflect distinct evolutionary strategies for environmental adaptation. As a perennial herbaceous plant adapted to shaded environments, E. pubescens may favor rapid, localized gene duplication through tandem repeats to enhance its adaptability. Tandem duplications often facilitate the rapid expansion of gene families within a short evolutionary timespan, forming a flexible and efficient regulatory network in response to external stimuli such as light stress and pathogen invasion [43, 44].

Additionally, the numbers of introns in the EpAPs examined exhibited a distinct pattern across different phylogenetic groups, with an average of 13 in Group I, 5.75 in Group II, 0.65 in Group III, and 1 in Group IV. This distribution pattern is generally consistent with APs in other plant species [11]. Notably, it has previously been shown that 46.6% (48/103) of the AP genes in E. pubescens are intronless, a proportion comparable to that observed in S. tuberosum (43.5%, 27/62) [42]. This simplified gene structure may confer an adaptive advantage to E. pubescens, allowing it to respond swiftly to environmental changes and maintain strong competitiveness under variable ecological conditions.

EpAPs are enriched in light-responsive elements

It has been shown that the promoters of all EpAPs are enriched with light-responsive elements, including the G-box, GT1-motif, and GATA-motif [45, 46, 47]. As a shade-loving plant, E. pubescens likely employs EpAPs to mediate responses to light stress. Studies on other species have demonstrated the critical roles of APs in chloroplast function and leaf senescence [18, 48, 49, 50]. For example, the chloroplast-localized protein CND41 in tobacco is linked to the degradation of Rubisco and exhibits a negative correlation with transcriptional activity within chloroplasts [18, 51]. In S. tuberosum, the canonical AP SPAP1 protein participates in ethylene-mediated leaf senescence, functioning in coordination with multiple members of the ethylene signaling pathway [52]. In addition, the atypical AP NANA, which is localized in the chloroplast in Arabidopsis, plays a crucial role in maintaining chloroplast homeostasis and starch metabolism. Its precise substrate and mechanisms of action remain unclear, but it has been shown that its loss reduces chlorophyll content, impairs electron transport, and alters sugar metabolism [36]. The presence of light-responsive elements in EpAP promoters suggests potential involvement in these pathways, but further functional validation studies are needed to elucidate specific regulatory networks.

Tissue-specific expression and abscission zone formation

The EpAPs examined in this study exhibited dynamic expression patterns across various plant tissues and participated in multiple biological processes. Tissue-specific analysis showed that EpAP102 and EpAP64 were predominantly expressed in the flowers and fruit, and their predicted interaction with EC1.5 suggests functions in reproductive development [53]. Similarly, the significant upregulation of EpAP15 and EpAP46 during the formation of fruit abscission zones indicates the involvement of these genes in the regulation of organ abscission. Such role has been suggested by some studies, although this remains a poorly documented subject. For instance, in a proteome analysis of ethylene-induced pedicel abscission in tomato, AP 1 (F2VJ74) showed significant changes in abundance during the formation of the pedicel abscission layer [54]. Additionally, a transcriptomic analysis of the sepal abscission zone (AZ-C) in ‘Hamlin’ oranges identified 16 differentially expressed APs by comparing fruit that detached from a Huanglongbing-infected tree upon shaking with those that remained attached to the infected tree [55]. Furthermore, in the present study, PPI analysis predicted interactions between several EpAP proteins and key regulatory proteins. For example, EpAP20/24/34 were predicted to interact with the PCD-related protease CATHB3 [53], potentially contributing to cell separation in the abscission zone via substrate hydrolysis. Additionally, the observed interactions with the E3 ubiquitin-protein ligases RDUF1 and DURF2 suggest a role in stress responses, including those to drought and salinity. Collectively, these interactions highlight the multifaceted roles of EpAPs in maintaining cellular homeostasis and facilitating stress adaptation. However, it is important to note that PPIs can vary across species, and interactions predicted based on model plants may not always directly apply to E. pubescens. The PPI results obtained in this study require experimental validation, for example via yeast two-hybrid assays, co-immunoprecipitation, or bimolecular fluorescence complementation, to confirm the validity of predictions.

EpAPs participate in light stress responses

Consistent with the presence of light-responsive cis-elements in their promoters, several EpAPs showed dynamic changes in expression in response to different periods of exposure to light stress. EpAP42 exhibited a stable expression pattern throughout the entire treatment period, indicating that this gene is not responsive to light stress. EpAP9 showed a transient increase in expression at the early stages of light stress (L3 and L6), followed by a return to the initial levels. This transient upregulation may reflect a short-term adaptation or regulatory response during initial acclimation to varying light conditions. Most of the tested genes, including EpAP11, EpAP19, EpAP20, EpAP68, and EpAP78, displayed a significant decrease in expression at L3, followed by a rebound at L9. Notably, only EpAP11 and EpAP19 surpassed their baseline expression levels at L0 during this recovery phase, which indicates their potential roles in sustained light stress adaptation or recovery processes. The partial recovery of expression levels observed for the other genes suggests that they may function during specific stages of the light stress response while remaining inactive for prolonged periods.

To further elucidate the roles of EpAPs, it is necessary to conduct functional validation studies using reverse genetics approaches, such as RNA interference or CRISPR/Cas9. Investigations into protein targets and interaction networks will also enhance our understanding of how these proteases contribute to light stress adaptation. Additionally, examining their expression under other abiotic stresses, such as drought or temperature extremes, could reveal shared or stress-specific regulatory pathways. Considering that E. pubescens is a shade-adapted plant, understanding how EpAPs respond to different light intensities could offer valuable insights for the optimization of its growth conditions. Targeted regulation of these genes to enhance light tolerance may assist in balancing biomass accumulation and the biosynthesis of secondary metabolites, such as flavonoids. Future studies could explore genetic or agronomic approaches, such as the overexpression of specific EpAPs or the adjustment of light exposure strategies, to improve both stress resilience and the yield of medicinal compounds.

Conclusion

This study represents the first comprehensive genome-wide analysis of the AP gene family in E. pubescens. It identified 103 EpAP genes and provided insights into their structural diversity, evolutionary conservation, and potential functional roles. These genes were shown to exhibit significant phylogenetic diversity as well as diverse structural domains and duplication patterns, which reflects their evolutionary adaptations to various physiological and environmental challenges. Notably, the presence of light-responsive cis-elements in EpAP promoters highlighted their potential involvement in light stress responses. Furthermore, tissue-specific expression patterns and predicted interactions with key regulatory proteins suggested roles in reproductive development, PCD, and abscission zone formation. Our findings establish a foundation for future functional studies aiming to explore the precise roles of EpAPs, potentially contributing to improving stress tolerance and reproductive efficiency in crops.

Methods

Identification and physicochemical characterization of epap proteins

The AP genes in E. pubescens were identified using a comprehensive approach combining HMMER v3.4 and BLASTP in order to ensure both analytical sensitivity and specificity. The hidden Markov model profile for the AP domain (PF00026) from the Pfam database was used as a query to search the E. pubescens genome (downloaded from the GigaDB database [56]) in HMMER v3.4, which effectively detects conserved domains even in divergent sequences. Protein sequences with an expected value (E) < 1E − 5 were retained. In addition, 70 AP protein sequences from A. thaliana [11] were used as queries in BLASTP (E < 1E − 5) to identify potential EpAPs based on sequence similarity. To improve accuracy, the sequences retrieved using both methods were merged, and duplicates were removed. NCBI CD-Search [57], SMART [58], and InterProScan [59] were employed to further validate the identified AP genes and minimize false positives, which allowed to confirm the presence of conserved domains and structural integrity. This multi-step approach ensured a robust and reliable identification of the EpAP gene family. The physicochemical properties of EpAP proteins were predicted in ProtParam [60]. Subcellular localization was predicted using the WoLF PSORT [61] and DeepLoc 2.0 [62] web tools. Additionally, TargetP 2.0 [63] and DeepLoc 2.0 were used to predict the presence and location of signal peptides in EpAP protein sequences.

Phylogenetic analysis

The 103 EpAPs identified in this study were aligned with 70 AtAPs in MUSCLE v5.0, and the data were subsequently refined using TrimAl. A maximum-likelihood phylogenetic tree was constructed using IQ-tree with 1,000 bootstrap replicates in TBtools v2.110 [64] and visualized via iTOL (https://itol.embl.de/) [65].

Analysis of gene structure, motifs, and domains

Exon-intron positions were extracted from the GFF annotation file obtained for the E. pubescens genome. The structural domains of EpAPs were analyzed using the NCBI Conserved Domain Database (https://www.ncbi.nlm.nih.gov/cdd/). Gene structures and protein domains were visualized using TBtools.

The conserved motifs in the EpAPs were predicted using the MEME online tool (http://meme-suite.org/tools/meme), with the predicted number of motifs set to 10. Gene and motif structures were mapped and visualized using TBtools.

Cis-elements in the promoter regions of EpAPs

To identify conserved cis-elements in the promoter regions of EpAPs, a screening was conducted of the sequences located 2000 bp upstream of the initiation codons, which represent the proximal promoter regions. Specifically, these sequences were analyzed using the PlantCARE database to predict promoter regions (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) [66]. Subsequently, the functional elements were categorized, and the results were visualized using TBtools.

Chromosomal locations, genomic duplications, and Ka/ks ratios

The chromosomal locations of EpAPs were determined using CDS coordinates from the E. pubescens genome. The genes were mapped onto their respective chromosomes based on their physical location and visualized using TBtools. Gene duplications within the EpAP gene family were illustrated using the Advanced Circos tool. One Step MCScanX was applied to identify homologous regions between E. pubescens and A. thaliana and between E. pubescens and C. chinensis. The Ka/Ks ratios were calculated to assess homologous relationships and estimating the divergence times of EpAPs. Homologous gene pairs were identified in BLAST, and their Ka/Ks ratios were computed using KaKs_Calculator. Based on the absolute molecular evolutionary rate (r) for Ranunculales plants (i.e., 6.98 × 10⁻⁹), the origin time of each duplication event was calculated using the formula T = Ks / (2r). The above analyses were all conducted using TBtools.

Expression patterns of EpAPs

To investigate the expression patterns of EpAPs, RNA-seq data (leaf: SRR15561380, SRR15562432, SRR15562431, SRR15562430, SRR15562429, SRR15562521, SRR15562520, SRR15562519, SRR15562518; Flower: SRR15569255, SRR15569254, SRR15569253, SRR15569252, SRR15569236, SRR15569235, SRR15569234, SRR15569232; fruit: SRR15569250, SRR15569249; root: SRR15561383, SRR15561382, SRR15561381; shoot: SRR15569248, SRR15569247) were obtained from the NCBI Sequence Read Archive database (Table S7) [56]. The transcriptome data related to fruit abscission in E. pubescens were obtained from the CNCB database (GSA number: CRA022394).

Predictive analysis of PPI networks

The EpAP protein sequences were analyzed on the STRING web platform (https://string-db.org/) using A. thaliana as a reference. The minimum required interaction score was set to 0.400, and the maximum number of interactors in the first shell was set to 20. Disconnected nodes were hidden, and the other parameters were maintained as default.

Plant materials and light stress treatments

In June 2024, 1-month-old seedlings of E. pubescens were obtained from the Xianling Medicinal Valley Yinyanghuo Cultivation Base (29°20 × 41.75´´ N, 103°33 × 39.60´´ E) in Shawan District, Leshan, Sichuan Province, and identified based on their morphological characteristics using the iPlant database (https://www.iplant.cn/info/Epimedium%20pubescens). The seedlings were then transferred to a greenhouse at the Southwest Characteristic Economic Plant Hybridization and Breeding Research Center of Leshan Normal University. Uniformly growing seedlings were selected and planted in pots measuring 11 cm in height, 11.5 cm in top diameter, and 8 cm in bottom diameter. The growth medium consisted of a 1:3 mixture of vermiculite and nutrient soil. Pre-treatment growth conditions were maintained at 22 °C under a light intensity of 25 µmol·m⁻²·s⁻¹. After an acclimation period of 3 weeks, the seedlings were exposed to light stress (330 µmol·m^–2·s^–1) for 12 h. Leaf samples were collected at 0 h (pre-treatment), 3 h, 6 h, 9 h, and 12 h post-treatment. These samples were subsequently flash-frozen in liquid nitrogen and stored at − 80 °C.

RNA extraction and qPCR

Total RNA was isolated using the FastPure^® Plant Total RNA Isolation Kit (Cat#RC401-01, Vazyme, Nanjing, China). cDNA was synthesized using the ABScript III RT Master Mix for qPCR with gDNA Remover (Cat#RK20428, ABclonal, Wuhan, China). qPCR amplification was conducted using the ChamQ Universal SYBR qPCR Master Mix kit (Cat#Q711-02, Vazyme, Nanjing, China), with Actin as the internal control. The qPCR primer sequences are listed in Table S8. The amplification protocol consisted of an initial denaturation step at 95 °C for 30 s, followed by 40 cycles at 95 °C for 10 s and at 60 °C for 30 s, with the melting curve being automatically generated. Relative gene expression levels were calculated using the 2^−ΔΔCt method [67], with three technical replicates per sample. Data analysis and visualization were performed using GraphPad 10. Statistical significance was assessed using one-way ANOVA followed by Dunnett’s multiple comparison test to determine differences between time points at specific light stress conditions. qPCR data were normalized based on the relative expression of EpAPs at time point L0.