Genomewide analysis of the Class III peroxidase gene family in apple (Malus domestica)

Gene identification

In the present study, the hidden Markov model (HMM) profile of PRX (PF00141) (Mistry et al., 2021) was downloaded and queried against the apple GDDH13 genome database (https://iris.angers.inra.fr/gddh13/index.html) (Daccord et al., 2017) using HMMER 3.0 software. Default parameters were used with a cutoff value set to 0.001 to identify PRX genes. To verify peroxidase-specific domains, candidate sequences were analyzed using three databases: (1) NCBI’s Conserved Domain Database (CDD) with an E-value threshold of 0.01, (2) the Pfam database (release 35.0), and (3) the SMART database (version 9.0). Sequences lacking a peroxidase-specific domain were excluded from further analysis. Candidate MdPRXs were cross-referenced with RNA-Seq and NCBI EST databases, removing those without detectable expression. The chromosomal locations of all identified PRX members were then mapped using MapChart 2.30 for further analysis.

Phylogenetic analysis

The full length sequences of the characterized PRX proteins were aligned using ClustalW with default parameters. Based on these alignments, phylogenetic trees were constructed in MEGA 7.0, employing the JTT model, pairwise deletion, and 1,000 bootstrap replications. Additionally, for the composite phylogenetic tree, the core peroxidase domain sequences of PRX proteins from Arabidopsis and Malus were aligned, and a neighbor-joining (NJ) tree was constructed using ClustalW in the same manner.

Sequence analysis, gene structure and motif analysis of MdPRXs

The sequence lengths, molecular weights, and isoelectric points of the identified PRX proteins were calculated using TBtools (Chen et al., 2023). Subcellular localization analyses were performed using WoLF PSORT II (https://www.genscript.com/wolf-psort.html). The exon-intron organization of Malus PRX genes was determined by comparing the predicted coding sequences with their corresponding full-length sequences using the Gene Structure Display Server (GSDS: http://gsds.cbi.pku.edu.cn) (Guo et al., 2007). Conserved motifs in MdPRX proteins were identified using MEME Suite version 4.12.0 (Bailey et al., 2009), with the analysis parameters set to detect a maximum of 10 distinct motifs.

Gene duplication

Gene duplication events were analyzed using the Multiple Collinearity Scan toolkit (MCScanX) with default parameter settings (Wang et al., 2012). A Circos plot (Kärkönen et al., 2009) was generated using TBtools software (Chen et al., 2020) to visualize the chromosomal distribution and relationships of segmental duplication gene pairs.

Transcriptome data

MdPRX expression was profiled using RNA-Seq data (GSE253335/276181/274104, PRJNA308148/826123/801073) across rootstocks, tissues, and developmental stages. The detailed RNA-Seq data information was listed in Table 1.

All samples were collected from Qingdao Academy of Agricultural Sciences (36°24′N, 120°58′E) under standard cultivation (fertilization: N in February, nutrient solution in May; pest control: imidacloprid March–May, pyrethroids May–July).

All RNA-seq datasets contain three biological replicates per sample, except for GSE253335, which has two. Samples were sequenced on an Illumina platform with 150 bp paired-end reads. Raw reads were processed with fastp software, to generate clean reads. The clean reads were mapped to the apple genome using HISAT2. The RNA-Seq data were normalized by Fragments Per Kilobase of Exon Per Million Mapped Fragments (FPKM). Differential expression analysis was performed with DESeq2.

RNA extraction and gene expression analysis

50 mg of samples were macro-dissected and immediately flash-frozen in liquid nitrogen, and stored at −80 °C. Total RNA was extracted using the Biospin Plant Total RNA Extraction Kit (Bioer Technology, Hangzhou, China) according to the manufacturer’s instructions. RNA integrity was checked by gel electrophoresis, and concentration was measured with a QuickDrop™ spectrophotometer (Molecular Devices). First-strand cDNA was synthesized using the HiScript III RT SuperMix + gDNA wiper (Vazyme), following the manufacturer’s recommendations. qRT-PCR was performed using the Roche LightCycler^® 480 system with SYBR Green I, with apple 18S rRNA as the internal reference. Primers were designed using NCBI Primer-BLAST, targeting the CDS of genes in Malus domestica (amplicons: 100–300 bp; Tm: 60 °C; GC: 40–60%; length: 18–25 bp). Intron-spanning was disabled as gDNA wiper effectively removes gDNA contamination during reverse transcription. Primers with off-target risks were excluded after genome-wide screening. Primers (Table S1) were HPLC-purified (Sangon Biotech). PCR conditions: 95 °C (30 s), then 40 cycles of 95 °C (10 s) and 60 °C (30 s). Reactions were run in triplicates, and gene expression was analyzed by 2^−ΔΔCT. Additional data (primer details, PCR specificity) are in File S1. Results (mean ± SD, n = 3) were analyzed with SPSS 18.0 and graphed using ORIGIN 8.0.

Protein structure prediction and comparison

Protein structure prediction of MdPRXs was performed using the Swiss-Model server (https://swissmodel.expasy.org/). The highest-quality models, as determined by GMQE (Global Model Quality Estimation) scores and confirmed as plant class III peroxidases, were downloaded in PDB format and analyzed using PyMOL (version 2.4.0; Schrödinger, LLC) for structural visualization and alignment (Schrödinger, 2021).

Identification of the PRX proteins in apple

The initial HMMER search identified 155 candidate gene models matching the PRX HMM profile (PF00141). PRX domain sequences were extracted from these candidates and aligned using ClustalW, followed by construction of an apple-specific PRX HMM profile using hmmbuild function in HMMER 3.0. A subsequent HMMER search with an E-value cutoff of 0.001 yielded 138 non-redundant candidates. Thirty-nine sequences lacking conserved peroxidase domains in NCBI CDD, Pfam, or SMART databases were excluded. A total of 99 PRX family genes were identified. These MdPRX genes were mapped to the Malus chromosomes based on their physical locations and were renamed MdPRX1 to MdPRX99 according to their order on the 17 apple chromosomes. During the subsequent gene structure and SMART analysis, we identified that MD11G1060900, (renamed as MdPRX61) possessed an unusually long sequence with 19 exons but lacked a signal peptide domain. Furthermore, no peroxidase domain was detected within the first 15 exons. IGV analysis revealed that only the first 15 exons were transcribed, whereas exons 16–19 exhibited no expression in the specific RNA-Seq data, suggesting that MD11G1060900 may consist of two independent transcripts. SoftBerry analysis predicted a novel transcript consisting of the last four exons, which contains both a signal peptide domain and a secreted peroxidase domain. This transcript was manually annotated as MdPRX61 using IGV-GSAman software (Fig. S1).

Physical distribution of MdPRX genes

The original gene IDs and their corresponding renamed IDs, along with their positions on the chromosomes, are available in Table 2. The FASTA information for all MdPRX coding sequences (CDS) is provided in File S2. The lengths of the identified CDS sequences ranged from 489 bp (MdPRX48) to 1,932 bp (MdPRX80), with an average length of 1,000.6 bp. A maximum of 15 MdPRX members were distributed on chromosome 3, followed by chromosomes 13 and 15, with 10 and nine members, respectively. Chromosome 14 contains the fewest MdPRX genes, with only two members, while the other chromosomes harbor between four and eight MdPRX genes (Fig. 1).

Phylogenetic relationships of MdPRX proteins in apple and Arabidopsis

To investigate the evolutionary relationships, we used the PRX domain amino acid sequences of 99 MdPRXs and 73 AtPRXs from Arabidopsis thaliana to construct a neighbor-joining (N-J) tree using MEGA 7.0. These PRX proteins were clustered into six groups. According to the phylogenetic distances, Group 1 was further subdivided into four subgroups (A, B, C, and D), and it represented the largest group, comprising 73 PRX members (46 from apple and 27 from Arabidopsis). Group 4 ranked as the second-largest cluster, containing 35 PRX proteins (18 from apple and 17 from Arabidopsis). Group 4 consisted of 26 PRX proteins, with 11 originating from Arabidopsis and 15 from apple. In contrast, Group 2, 3, and 6 were the smallest groups, each containing only 13 PRX members (Fig. 2).

Absolutely purifying selection for MdPRX genes

To elucidate the evolutionary expansion of MdPRX gene family, McScanX was employed to analyze duplication events of MdPRX genes. The results showed that 56 MdPRX genes were involved in duplication events, resulting in 38 pairs of segmental duplications among MdPRX members. Additionally, seven pairs of segmental duplications were identified between MdPRX members and other regions of the apple genome. No tandem repeat events were detected, suggesting that segmental duplication plays a dominant role in the expansion of the MdPRX gene family. A Circos map illustrating these duplication events was generated based on collinearity analysis. In apple, Chr 03 and Chr 11, Chr 13 and Chr 16, and Chr 05 and Chr 10 are copies of the same ancestral chromosome (Velasco et al., 2010). In this study, we observed strong collinearity of MdPRX members between Chr 13 and Chr 16, Chr 03 and Chr 11, and Chr 05 and Chr 10, whereas the MdPRX members on Chr 08 exhibited no collinearity with other chromosomes (Fig. 3).

The selection pressure among different types of duplication was also inferred by calculating the rates of synonymous substitution (Ks) and non-synonymous substitution (Ka). During evolution, genes typically undergo purifying selection (Ka/Ks < 1), positive selection (Ka/Ks > 1), or neutral selection (Ka/Ks = 1) (Khan et al., 2019). In this study, most of the Ka/Ks values for the MdPRX collinear pairs were less than 0.5, excerpt for two pairs (MdPRX27 vs MD12G1121300 and MdPRX77 vs MdPRX91) showed values around 0.5–0.6 (Fig. 4, File S3). The divergence time of each pair ranged from 3.61 million years ago (Mya) to 116.28 Mya. Therefore, we propose that the MdPRX gene family has undergone strong purifying selection at a slow evolutionary rate.

The first class III peroxidase appeared approximately 450 Mya, coinciding with the emergence of terrestrial plants. Since then, the number of gene copies has increased significantly. This expansion appears to be correlated with the evolution of plant architecture and complexity, as well as the diversification of biotopes and pathogens (Passardi et al., 2005). The number of MdPRXs (99) exceeds that of Arabidopsis thaliana (73), grapevine (47), and pear (94). This difference stems from to apple’s more complex plant architecture in contrast to the herbaceous nature of Arabidopsis thaliana and the vine-like structure of grapevine. Notably, both pear and apple evolved from a common ancestor that existed around 5–20 Mya. However, apple possesses a greater number of PRX genes than pear, which likely contributes to its superior adaptability—particularly in terms of cold and drought tolerance—as well as its broader ecological adaptability. This phenomenon of gene family expansion due to enhanced adaptability is also observed in the WRKY gene family, with apple possessing 118 members compared to pear’s 103 (Huang et al., 2015; Lui et al., 2017).

Gene characteristics and cellular sublocallization of MdPRXs

Gene characteristics, including protein sequence length, molecular weight (MW), isoelectric point (pI), and subcellular localization, were analyzed (File S4). Protein lengths ranged from 162 aa (MdPRX22) to 643 aa (MdPRX80), with an average length of 332.5 aa. The MW ranged from 17,702.1 Da (MdPRX22) to 70,893.1 Da (MdPRX80), with an average of 36,268.7 Da. The isoelectric point (pI) ranged from 4.03 (MdPRX12) to 9.63 (MdPRX41), with an average of 7.40. Most MdPRXs (84 out of 99) exhibited hydrophobicity. According to TargetP (https://services.healthtech.dtu.dk/services/TargetP-2.0/), 88 out of 99 MdPRX members possess a signal peptide sequence, whereas 11 MdPRXs lack signal peptide sequences (File S5). According to Wolf PSORTII predictions, two MdPRXs (MdPRX29 and MdPRX48) of them were predicted to localize to the cytoplasm, whereas the other nine MdPRXs were predicted to be chloroplast-localized. 33 MdPRXs were predicted to localize to extracellular regions, while 50 (including the nine members without signal peptides) and one member were predicted to localize to chloroplasts and mitochondria, respectively. Notably, 10 MdPRXs were predicted to localize to vacuoles, however, all of these vacuole-localized PRXs possess signal peptides. The remaining MdPRXs were predicted to localize to the endoplasmic reticulum (ER), plasma membrane, and mitochondria (File S4).

Gene structure, motif and conserved domain analysis of MdPRXs

The 10 conserved motifs (motifs 1-10) of  MdPRX were identified using the MEME program. Notably, Motif 1, Motif 10, Motif 2, and Motif 9 were the most common among MdPRXs (Fig. 5B). LOGOs for these motifs were also generated from MEME (Fig. S2). The highest consensus sequence (41) was observed in Motif 1, while the lowest (11) was recorded in Motif 9. All MdPRX members, except for MdPRX11, MdPRX19, MdPRX22, MdPRX48, MdPRX67, MdPRX87 and MdPRX80, contain all ten conserved motifs arranged in the following order: motif 1, motif 10, motif 2, motif 9, motif 5, motif 8, motif 3, motif 4, motif 7, and motif 6. Notably, MdPRX80 featured two sets of motifs (Motifs 1 to 5 and Motifs 7 to 10), but contained only a single copy of Motif 6 (Fig. 5A).

Gene structure analysis revealed that MdPRX members had between one and four exons, with most containing four. Diverging from high structural conservation observed in protein architecture across the whole gene family, MdPRX gene structures exhibited partial phylogenetic coherence. Paralogous clusters within subclades maintained conserved exon-intron architectures, for example:six genes (MdPRX50/96/72/91/36/37) which clustered within the same subclade, all contained only one exon (Fig. 5C).

Batch SMART analysis revealed that all MdPRXs, contain a single Pfam peroxidase domain, except for MdPRX80, which contains two distinct Pfam peroxidase domains (Table 2). In all cases, the peroxidase domains encompass all exons (Fig. 5C). The majority of MdPRXs contain a signal peptide domain at the N-terminus (Fig. 5C; File S6).

Cis-regulatory element analysis of MdPRXs

The analysis of cis-regulatory elements within the 2,000 bp upstream regions of the sequences is presented. MdPRXs clustered within the same phylogenetic clades exhibit distinct cis-element profiles in their upstream regulatory regions, suggesting divergent regulatory mechanisms and transcriptional responsiveness (Fig. 6).

Among the identified cis-elements of MdPRXs, light-responsive elements constituted 61% (Table S2), suggesting a potential role for MdPRXs in photoprotection through ROS scavenging. O₂^•− could be produced during photosynthesis in the chloroplasts, but their accumulation can lead to oxidative stress, damaging various cellular components (Asada, 1999). Once O₂^•− is converted to H₂O₂ by SOD, CATs and PRXs further degrade H₂O₂ into water and oxygen, thereby completing the detoxification process. This step is essential because H₂O₂, although less reactive than O₂^•−, can still cause oxidative damage if not removed (Zeng et al., 2017).

Approximately 13% of identified cis-elements were methyl MeJA-responsive, 5% drought-responsive, and 4% ABA-responsive elements. These findings align with previous studies on the enzymatic and molecular functions of PRXs, providing additional evidence for their dual roles in biotic and abiotic stress responses (New, Barsky & Uhde-Stone, 2023). Additionally, 4% and 2% of the identified cis-elements are associated with gibberellin- and auxin-responsive elements, respectively, suggesting that these PRXs may be involved in growth regulation. Elements specific to meristem, endosperm, and seed tissue expression, as well as those involved in circadian control, were also identified within the promoters of MdPRXs.

Spatio-temporal expression patterns of MdPRXs

The average FPKM value from three replicates for each MdPRX was obtained from transcriptomic data of axillary buds, leaves, stems (young shoots), and flowers of Malus spectabilis ‘Bly114’ grafted onto M. robusta, as well as from M. robusta root tips (GSE274104 and PRJNA801073). A total of 68–89 members were expressed in these tissues, with a minimum of 68 in buds and a maximum of 87 in roots, indicating a preferential accumulation of MdPRXs in roots. More than half (55 out of 99) of MdPRXs exhibited constitutive expression patterns, while some members showed strict tissue-specific expression patterns. For example, MdPRX62, MdPRX63, MdPRX22, and MdPRX30 were exclusively expressed in roots. Others exhibited significantly elevated expression levels (log2 fold change >2) in particular tissues relative to others, as evidenced by their distinct red coloration in the heatmap (data were log2-transformed and subjected to row normalization). For example, MdPRX 38/83/95/99 exhibited predominant expression in stem tissues, whereas MdPRX4 and MdPRX97 displayed leaf-preferential expression patterns (Fig. 7A).

Interestingly, the expression of MdPRX15 was not detected in any of the tissues analyzed in this study. However, a BLAST search of its CDS against the NCBI EST database revealed that MdPRX15 is expressed in xylem tissue of Royal Gala plants subjected to 5 °C for 24 h at a specific developmental stage (EST: lIBEST_024520).

The widespread distribution of MdPRX s across various non-stressed apple tissues and organs (with the exception of MdPRX15) strongly supports their multifunctional roles in developmental regulation, particularly through tissue-specific specialization of individual members. This conclusion is consistent with and substantiated by previous studies in model systems, including Arabidopsis leaf, stem, and root development, as well as transgenic analyses in both tobacco and apple systems. In Arabidopsis thaliana, a set of PRXs has been implicated in regulating leaf cell expansion and determining final organ size by modulating ROS homeostasis in the apoplas (Lu et al., 2014). In stems, peroxidases contribute to lignin deposition. They play a role in xylem differentiation and vascular tissue development (Hoffmann et al., 2020). Arabidopsis PRX8 transgenic apple showed an increase in the number of xylem vessels in the stem (Vicuna, 2005). Arabidopsis PRX01, PRX44, and PRX73 participate in root hair formation and root tissue differentiation (Marzol et al., 2022). Based on microarray analysis, AtPrx13, AtPrx30 and AtPrx55 potentially involved in Arabidopsis flower and fruit development (Cosio & Dunand, 2010).

The expression patterns of MdPRXs during axillary bud sprouting

Quiescent multibranching mutant (MB) axillary buds, minimally branching buds, and elongation zones of new branches were designated as MB1, MB2, and MB3, respectively (Fig. S3). Meanwhile, WT (Malus spectabilis ‘Bly114’) axillary buds lacking branching potential were designated as WB. RNA-Seq analysis was performed on MB1, MB2, and WT axillary buds (PRJNA801073). As illustrated in Fig. 7B, Seventy-five MdPRXs were expressed in these three types of buds and clustered into five clades. MdPRXs in Clade 1 were downregulated in MB2, while those in Clade 2 were upregulated in MB2, indicating their involvement in the axillary bud-breaking process. MdPRXs belonging to Clade 3 exhibited higher expression in WB but lower expression in MB1 and MB2, suggesting their role in maintaining the stemness of WB.

Traditionally, PRXs are considered H₂O₂ consumers, transferring electrons from other substrates to H₂O₂ and reducing it to water. The function of Clade 3 MdPRXs aligns with this traditional role. However, Class III PRXs participate in three distinct cycles—peroxidative, oxidative, and hydroxylic (Abbas et al., 2018). While they reduce H₂O₂ to water in the peroxidative cycle, they also produce hydroxyl radicals (OH^•) in the hydroxylic cycle and superoxide anions (O₂^•−) in the oxidative cycle (Passardi et al., 2005; Shigeto & Tsutsumi, 2016). Thus, PRXs function as bifunctional reductases, both consuming H₂O₂ and generating ROS species such as OH^• and O₂^•−. In this study, Clade 4 and 5 PRXs are upregulated in MB2 and/or MB1 (Fig. 7B). Buds must maintain high O₂^•− levels and low H₂O₂ levels to preserve stemness. However, to break the bud dormancy, H₂O₂ should accumulated in the peripheral (elongation) zone to inhibit WUS (WUSCHEL), a key regulator of stem cell fate in plant meristems, thus differentiation take place (Zeng et al., 2017). Given that elevated H₂O₂ concentrations in the peripheral zone are associated with differentiation, these upregulated MdPRXs may contribute to H₂O₂ production rather than its elimination.

The expression profile of MdPRX87, a Clade 3 member, was analyzed via qRT-PCR, revealing a gradual decrease across the SCN region in MB1 (quiescent lateral bud), MB2 (sprouting axillary bud), and the elongation zone in MB3 (developed lateral shoot). The highest expression level was observed in WT axillary buds, which are incapable of branching (Fig. 7C). This suggests that MdPRX87 may contribute to H₂O₂ elimination in the QC region of axillary buds. As its expression declines, H₂O₂ concentration increases in the peripheral zone of MB buds, promoting lateral branch differentiation. In contrast, the high expression level of MdPRX87 in WT axillary buds may suppress differentiation of the lateral branches.

The expression patterns of MdPRXs in dwarfing and vigorous rootstocks across different growth stages

A total of 78 MdPRXs were expressed in phloem tissues. In RNA-seq dataset GSE253335, 45, 50, 69, 45, and 54 MdPRXs were detected in the actively growing phloem of A1d, Budagowski 9 (B9), M. sylvestris, M. hupehensis (PYTC), and M9, respectively. Meanwhile, 32 and 49 MdPRXs were expressed in the bud-breaking phloem of A1d (A1d0) and M. hupehensis (PYTC0) (RNA-seq dataset GSE276181), respectively. Overall, more MdPRXs were expressed during active growth than at the bud-breaking stage. Remarkably, five MdPRX genes—MdPRX1, MdPRX28, MdPRX50, MdPRX57, and MdPRX61—exhibited constitutive high-level expression patterns across all examined growth stages in both vigorous and dwarfing rootstock genotypes. In contrast, MdPRX79, MdPRX96, MdPRX58, MdPRX39, and MdPRX53 were also ubiquitously expressed but at lower levels. MdPRX5, MdPRX35, MdPRX71, and MdPRX7 were highly expressed in the actively growing phloem of both vigorous and dwarfing rootstocks; however, they exhibited little to no expression at the bud-breaking stage (Fig. 7D).

MdPRX99 is highly expressed in the dwarfing rootstock M26, in contrast, MdPRX59 was exclusively expressed in vigorous rootstocks (M. hupehensis and M. sylvestris) across both bud-breaking and actively growing stages. Further RT-PCR analysis of MdPRX59 in the phloem of four dwarfing rootstocks (A1d, B9, M9, QZ1) and four vigorous rootstocks (M. hupehensis, M. sylvestris, M. sieversii, M. baccata) at the active growth stage (Fig. 7E) confirmed that MdPRX59 was expressed exclusively in vigorous rootstocks. These findings posit MdPRX59 expression deficiency as both a predictive biomarker for dwarfing rootstock selection and a candidate breeding target. Mechanistically, targeted downregulation of MdPRX59—achievable through CRISPR/Cas9-mediated knockout or RNAi silencing—represents a strategy for engineering rootstocks with elevated dwarfing ability, warranting prioritized functional validation in apple breeding programs. Concurrently, complementary investigations should evaluate whether MdPRX99 over-expression conversely enhances dwarfing capacity.

Based on total FPKM values, the previous proposal that higher PRX activity correlates with dwarfing ability was supported by comparisons between WT M. hupehensis (PYTC) and its GA-insensitive dwarf mutant (A1d), as well as between M9 and the two types of vigorous rootstocks. However, this correlation did not hold when comparing B9 with vigorous rootstocks. B9 exhibited the lowest total PRX transcription level, however it is known for its strong dwarfing ability among several rootstocks (e.g., M.9T 984, M.9T 337, Jork 9, Mark 9, B9, M.9 EMLA, Pajam 1, Pajam 2, and Supporter 4) (Gjamovski & Kiprijanovski, 2011). In addition, there was little difference between total PRX transcription of the vigorous M. sylvestris and the dwarfing rootstock A1d.

Previous studies have demonstrated that PRX activity progressively increases during plant development, correlating with the accumulation of lignin content. In accorded with this proposal, we found that the total of PRX transcripts of A1d at the bud-breaking stage was much lower than that of at actively growing stage. However, in the vigorous rootstock M. hupehensis, higher total PRX expression was observed at bud-breaking stage than that of at actively growing stage.

These discrepancies suggest that PRXs may have bifunctional or antagonistic roles in regulating plant growth vigor, an assertion supported by existing studies. Different PRX members exhibit opposing functions in cell wall modification, as well as in auxin (IAA) metabolism and signaling (Chen & Schopfer, 1999; Francoz et al., 2015; Kawano, 2003; Schopfer et al., 2002). Some PRXs oxidize aromatic cell wall compounds in the presence of H₂O₂, thereby stiffening the cell wall and causing growth retardation (Francoz et al., 2015; Passardi, Penel & Dunand, 2004). In contrast, other PRX members may promote cell division and growth by reducing local H₂O₂ concentrations or generating OH^• that break covalent bonds in cell wall polymers (Schopfer, 2001). In relation to IAA metabolism and signaling, PRX-mediated signaling primarily facilitates the degradation (oxidation) of both IAA and its downstream signaling molecule, H₂O₂, through the conventional PRX cycle. Meanwhile, extracellular IAA activates a plasma membrane-bound NADH oxidase, stimulating the apoplastic production of superoxide anions (O₂^•−), a portion of which is subsequently converted into H₂O₂. The generated H₂O₂ then serves as a substrate for PRX enzymes bound to cell wall polymers, leading to the production of hydroxyl radicals (OH^•). These OH^• radicals cleave the backbone of cell wall polysaccharides, thereby mediating IAA-induced cell elongation and ultimately promoting plant growth (Chen & Schopfer, 1999; Kawano, 2003; Schopfer et al., 2002).

Rootstock-Scion interaction on MdPRX expression characteristics

The expression profiles of MdPRXs in Fuji6 scions grafted onto vigorous rootstock WT M. hupehensis (PF) and its GA-insensitive dwarfing mutant A1d (AF) were analyzed using RNA-Seq data (PRJNA826123). The results indicated that rootstocks influence the expression of MdPRXs in scions. Specifically, MdPRX26, MdPRX9, MdPRX5, MdPRX4, MdPRX95, and MdPRX38 were significantly downregulated in AF compared to PF, whereas MdPRX35 exhibited significant upregulation in AF. Additionally, the total MdPRX transcription level, calculated as the sum of FPKM values, was lower in AF than in PF (Fig. 7F).

Notably, our previous study documented a statistically significant reduction (p < 0.05, Student’s t-test) in internode length between AF and primary PF (Shi et al., 2023). This observation aligns with reports indicating that Norway spruce plants overexpressing the peroxidase-like gene spi2 exhibited higher guaiacol peroxidase activity concomitant with reduced height and cell length (Elfstrand et al., 2001). Despite these precedents, the correlation between altered MdPRX expression and scion performance still needs to be validated through transgenic experiments.

It was also observed that the root tips of M. robusta rootstocks grafted with scions from WT M. spectabilis and a more-branching (MB) mutant exhibited distinct MdPRX expression patterns (Fig. 7G). Additionally, the primary root length, diameter, and root weight of rootstocks grafted with MB scions were significantly reduced compared to those grafted with WT scions (Li et al., 2016). Thirty-one PRX members were significantly upregulated, while six were downregulated in the root tips grafted with MB scions compared to those grafted with WT scions. This suggests that the interaction between scions and rootstocks influences MdPRX expression patterns, which in turn affects root growth and development. The result is consistent with findings that overexpression of a tobacco anionic peroxidase gene led to a 50% reduction in root weight in transgenic tobacco plants (Lagrimini et al., 1997). It also consistent with the study proposed that the accumulation of peroxidase-related ROS in trichoblasts correlates with root hair initiation in barley (Kwasniewski et al., 2013). These results suggest that scions may regulate MdPRX expression in the RAM, thereby influencing the differentiation of rootstock root tips.

MdPRX protein structures

The three-dimensional structures of six representative MdPRX proteins were predicted using homology modeling based on high-quality templates retrieved from the Protein Data Bank (PDB) (Berman et al., 2000) or from the AlphaFold Protein Structure Database (https://alphafold.ebi.ac.uk; DeepMind and EMBL-EBI, CC-BY 4.0) (Jumper et al., 2021). The AlphaFold templates (AFDB IDs M5XQS7, M5XC68, A0A2P6RAP4, and A0A834ZTL9, corresponding to MdPRX76, MdPRX90, MdPRX27, and MdPRX99, respectively) and the PDB templates (PDB IDs 4cuo and 1qgj, corresponding to MdPRX59 and MdPRX2, respectively) were utilized, with all models meeting a stringent Global Model Quality Estimation (GMQE) threshold of ≥0.85 to ensure structural reliability. The results revealed that MdPRXs exhibit highly similar structural architectures. As illustrated in Fig. 8A, six MdPRXs homologous structures were aligned with near-perfect precision. These proteins are characterized by two predominantly α-helical domains, along with a smaller third domain, the “β-domain”. Collectively, these three structural domains constitute a conserved pocket-like architecture at their interface. The heme ligand is anchored to the PRXs via a network of both polar and nonpolar interactions within the pocket-like interdomain configuration. This arrangement not only stabilizes the heme but also creates an accessible environment for substrates to enter and engage with the heme, thereby, facilitating their catalytic function. Notably, the size of these pockets varies among different MdPRX members, with some displaying larger pockets compared to others. Specifically, MdPRX2, MdPRX76, and MdPRX99 exhibit larger pockets, as shown in Fig. 8B, while MdPRX27, MdPRX59, and MdPRX90 feature smaller pockets, as depicted in Fig. 8C. Studies suggest that PRX-oxidized compounds (e.g., lignin precursors) require heme-edge accessibility to undergo catalytic transformation (Kwon, Moody & Raven, 2015). This suggests that variations in pocket size among different PRXs may influence their substrate accessibility and catalytic efficiency, thereby modulating their biological roles. Our results are align with the proposal derived from AlphaFold modeling (New, Barsky & Uhde-Stone, 2023) that PRXs capable of producing ROS exhibit distinct interdomain and protein–heme interactions, potentially functioning as ‘gatekeepers’ by preventing larger substrates from accessing the heme.

In the context of axillary bud sprouting, MdPRX59 /27/90, which are expressed at higher levels in the sprouting MB2 compared to the quiescent MB1 and WB, feature a smaller heme-containing interdomain. This structural characteristic may enhance their capacity for ROS production, leading to an increase in local H₂O₂ concentration within the peripheral zone of the SCN. Such a shift could disrupt the delicate balance of high O₂^•− and low H₂O₂ levels in the SCN, thereby promoting differentiation and initiating the sprouting process. Conversely, MdPRX2/76, which have larger heme-substrate pockets, may act as ROS scavengers, thereby reducing the concentration of H₂O₂ in the SCN. Notably, these genes demonstrated the highest expression levels in WB, while displaying relatively lower expression in both MB1 and MB2, implying their potential role in maintaining low H₂O₂ concentrations in the dormant WB by consumping it.

Concerning the regulation of rootstock growth vigor, MdPRX99, which exhibited the highest expression in the dwarfing rootstock M26, is classified as a large-pocket PRX. It may function as a conventional PRX, oxidizing larger substrates such as IAA and lignin precursors, thereby contributing to growth inhibition. In contrast, MdPRX59, which is specific to vigorous rootstocks, has been identified as a smaller-pocket PRX. It may play a role in the production of ROS, specifically OH^• and O₂^•−, thereby promoting growth by modulating cell wall elasticity.

Purifying selection versus cis-regulatory divergence drive functional diversification

Phylogenetic analysis of MdPRXs coupled with McScanX-based genomic collinearity assessment revealed distinct evolutionary trajectories of gene duplication events. Notably, most segmentally duplicated gene pairs were clustered within the same phylogenetic subclades and exhibited minimal genetic divergence, suggesting functional conservation during evolution. Only a limited subset of duplication pairs exhibited sequence divergence after duplication. Subsequent evolutionary rate analysis (Ka/Ks < 1) further confirmed that MdPRXs evolved under strong purifying selection. These findings were supported by conserved structural patterns in gene motif architecture, where paralogs across divergent phylogenetic groups retained nearly identical motif organization and composition profiles, indicating strong protein level sequence conservation. Beside the segmental duplication gene pairs, several groups of consecutive MdPRX members were observed in the same branches of phylogenetic tree, such as MdPRX11-15 in Group 1a and MdPRX45-47 in Group 2 formed tight phylogenetic clusters (Fig. 3). These consecutive MdPRXs were not classified as tandem duplication events due to McScanX’s inherent filtering parameters. McScanX employs a consolidation protocol for collinear gene identification wherein consecutive BLASTP matches sharing a common anchor gene are algorithmically processed. Specifically, if paired homologs reside within a genomic interval of fewer than five intervening genes, these matches are computationally merged, retaining only the representative gene pair exhibiting the most significant sequence similarity (i.e., the lowest BLASTP E-value). This filtering mechanism prevents overcounting of tightly linked homologs (Wang et al., 2012). These local tandem arrays are evolved from very early duplications in the phylogenies, for which mechanisms are difficult to infer (Cannon et al., 2004). Along with segmental duplicates, these several sets of local tandem arrays contribute to the highly homologous sequences of MdPRX gene family.

However, despite their high degree of protein sequence homology, 3D structural predictions revealed substantial variations in the dimensions of the heme-substrate binding interdomain among PRX members, even within the same phylogenetic clades. For instance, although MdPRX90 and MdPRX2 cluster within the same phylogenetic group (Fig. 2, group 1), they exhibit distinct structural characteristics, with MdPRX90 representing a small-pocket PRX variant and MdPRX2 belonging to the large-pocket PRX category. These observations suggest that even minor variations in amino acid sequences can dramatically alter the architecture of the heme-substrate binding domains. Our observations are consistent with the findings of (New, Barsky & Uhde-Stone, 2023), who identified three key amino acid positions (Alpha1, Alpha2, and Beta Buttons) that significantly influence the size of the heme-substrate binding pocket. Specifically, PRXs with smaller heme-substrate pockets tend to feature arginine (R) or lysine (K) at the Alpha1 Button, serine (S) at the Alpha2 Button, and arginine (R) or glutamine (Q) at the Beta Button. In line with this research, we observed a strong correlation between the size of the heme-substrate binding domain and the expression patterns of MdPRXs in WB, MB1, and MB2 (Fig. 7B). MdPRXs belonging to clade 2 of Fig. 7B, such as MdPRX2 and MdPRX76, are exhibit larger heme-substrate binding pockets, whereas those from clade 5, including MdPRX59, MdPRX27, and MdPRX90, display smaller pockets. This clustering pattern further supports the functional relevance of these structural variations in MdPRXs.

Additionally, MdPRX promoters showed cis-regulatory divergence, with segmental duplicates and some tandem arrays having distinct cis-element combinations and motif arrangements. Subcellular localization divergence was further observed between segmentally duplicated pairs and localized tandem arrays, likely resulting from evolutionary diversification of sorting signal peptides. This spatial differentiation implies post-duplication functional specialization, with paralogs gaining distinct sorting signals. As consequence, despite phylogenetic clustering, most MdPRX paralogs exhibited divergent spatiotemporal expression patterns, as evidenced by incongruent profiles in expression heatmaps (Fig. 7). The observed dichotomy suggests that while MdPRXs retained ancestral peroxidase activity (as indicated by phylogenetic conservation), they evolving novel regulatory logics. They have acquired specialized regulatory codes enabling tissue- and stage-specific responses to developmental cues. Thus, the highly homologous and conserved MdPRX proteins response to specific developmental signals in different organs/tissues at different stages. For example, MdPRX13 and MdPRX59 are a pair of segmental duplication genes, with only MdPRX59 exhibiting exclusive expression in vigorous rootstocks. Given the observed cis-regulatory diversification among MdPRX paralogs, future transgenic validation studies should incorporate both functional characterization of promoter architectures and coding sequence analyses. This dual-focused approach is essential for untangling the relative contributions of protein sequence conservation and regulatory element innovation in shaping the spatiotemporal functional diversification of the MdPRX family.

MdPRXs integrate hormone signals and regulate plant development through ROS homeostasis

Since indoleacetic acid (IAA) is one of the substrates of PRX, the relationship between IAA and PRXs has been extensively investigated. However, in this study, we revealed that MdPRXs may also play a significant role in homeostasis, transportation and signal transduction of strigolactone (SL), ABA, and cytokinin (CK), through the processes of reactive oxygen species (ROS) production and scavenging.

The MB apple mutant used in this study is a SL-insensitive mutant. Compared to WT, the IAA content and transcription of the IAA transporter PIN were significantly lower, while the ABA and CK contents were higher in the axillary buds of the MB mutant (Tan et al., 2019). Extensive research has established intricate connections between plant hormone signaling and ROS metabolism. Specifically, apoplastic ROS can be induced and regulated by ABA (Li et al., 2022), while auxins can increase O₂^•− production and simultaneously induce the expression of reductases such as PRXs, SODs, and CATs, thereby maintaining ROS homeostasis (Mir et al., 2020). This ROS equilibrium reciprocally regulates polar auxin transport dynamics (Cséplő et al., 2021). Additionally, exogenous SL (GR24) treatment has been shown to reduce H₂O₂ levels in cucumber seedlings (Zhou et al., 2022).

In the tomato system, two NADPH oxidase isoforms, RBOH1 and WFI1, serve as the primary enzymatic sources of apoplastic H₂O₂ production. The accumulation of H₂O₂ in the apoplast functions as a critical secondary messenger that activates auxin biosynthesis pathways in the shoot apical meristem, subsequently downregulates CK biosynthesis in stem tissues (Chen et al., 2016). The consequent hormonal imbalance, characterized by elevated auxin-to-cytokinin ratios, effectively suppresses axillary bud development, thereby maintaining apical dominance. Given the dual enzymatic functionality of MdPRXs-functioning either as H₂O₂-scavenging enzymes that antagonize NADPH oxidase activity or as H₂O₂-generating enzymes analogous to NADPH oxidases, as previously discussed, we propose that the tissue-specific expression patterns of MdPRXs may mediate the observed physiological changes in WB tissues. Specifically, these expression patterns likely contribute to the downregulation of PIN-formed auxin transporters and the concomitant elevation of cytokinin levels, potentially through modulation of ROS homeostasis and subsequent hormonal signaling pathways. Furthermore, the resulting hormonal imbalance, involving dynamic changes in IAA, CK, and SL levels, appears to establish a feedback regulatory mechanism that modulates MdPRX expression.

Supporting this hypothesis, cis-element analysis revealed the presence of ABA-, IAA-, and SA-responsive elements in the upstream promoter sequences of MdPRXs. These findings strongly suggest the existence of a feedback regulation network between hormonal signaling pathways and peroxidase-mediated ROS homeostasis. Within this regulatory circuitry, PRXs function as molecular hubs that process hormonal signals to finely tune ROS dynamics, thereby regulating lateral shoot differentiation and growth vigor (Fig. 9). According to our proposed regulatory model (Fig. 9), distinct MdPRX members exhibit specific responsiveness to different plant hormones, enabling precise modulation of ROS (H₂O₂ and O₂^•−) concentrations within OC and elongation zones of axillary buds. This ROS-mediated regulation subsequently influences multiple physiological processes, including auxin biosynthesis, polar auxin transport through PIN proteins, and cytokinin metabolism. The resulting hormonal crosstalk ultimately determines axillary bud developmental fate as well as growth vigor. However, further studies are needed to validate this hypothesis. For instance, over-expressing specific MdPRXs that are up-regulated in sprouting MB2 buds such as MdPRX90/59/27 could provide critical insights into their functional roles and confirm the proposed regulatory mechanism.