The protein sequences of Arabidopsis thaliana PAO proteins were obtained from The TAIR database (https://www.arabidopsis.org/). The protein sequence FASTA files of eight Pyrus species including Cuiguan, Shanxi Duli, Zhongai 1, Nijisseiki, Yunhong No.1, d’Anjou, Bartlett v2.0, Dangshansuli’ v.1.1 were used as subject sequences to run the command-line tool, BLAST+. A BLASTp search was conducted against the eight Pyrus proteomes using AtPAO protein sequences as the queries. The hits from this search were then refined, including the removal of isoforms and duplicates.

The Pyrus PAO candidate sequences were examined for domains through the NCBI Conserved Domain Database (CDD; https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) (Marchler-Bauer et al., 2011) and InterPro (https://www.ebi.ac.uk/interpro/) (Hunter et al., 2009) to identify the definitive protein family sequences. Physicochemical properties such as molecular weight (MW), isoelectric point (pI), aliphatic index (AI), instability index (II), and the grand average of hydropathicity (GRAVY) were forecasted using the ExPASy ProtParam tool (https://web.expasy.org/protparam/) (Gasteiger et al., 2005). The subcellular localization of these PAO proteins was forecasted using the WoLF PSORT tool (https://wolfpsort.hgc.jp/) (Horton et al., 2007).

Phylogenetic analysis was carried out to assess the intra-specific evolutionary connections among Pyrus PAOs. A multiple sequence alignment of 11 Cuiguan, 8 Shanxi Duli, 10 Zhongai 1, 10 Nijisseiki, 10 Yunhong No.1, 9 d'Anjou, 5 Bartlett v2.0, 8 Dangshansuli’ v.1.1, 4 A. thaliana, 5 O. sativa, and 4 Z. mays protein sequences was performed using ClustalW (Thompson et al., 2003). The phylogenetic tree was constructed using the IQTREE Web Server (http://iqtree.cibiv.univie.ac.at/) (Trifinopoulos et al., 2016) with the maximum likelihood (ML) method and bootstrap replicates of 1000. The iTOL: Interactive Tree of Life (https://itol.embl.de/) (Letunic and Bork, 2007) was used for editing and visualization.

To identify the conserved common motifs in all eight Pyrus PAO sequences, the MEME tool (https://meme-suite.org/meme/) (Bailey et al., 2009) was utilized. The number of conserved motifs was set to 20 for each sequence. Gene structures were constructed using CDS and genomic sequences through the Gene Structure Display Server (GSDS; https://gsds.gao-lab.org/) (Hu et al., 2015). The identified motifs and gene structures were visualized using TBtools-II v2.067 (C. Chen et al., 2018).

The chromosomal positions of each Pyrus PAOs were extracted from GFF/GFF3 files and mapped to chromosomes using TBtools-II v2.067: Gene Location Visualize advanced (C. Chen et al., 2018). By assessing whether the shorter gene’s length covered 70% of the longer gene and if the alignment similarity between the two genes was 70% or higher, instances of PAOs gene duplication were identified (ul Qamar et al., 2023). The duplication pattern, whether segmental or tandem, was also examined. To predict the selection pressure on the duplicated genes, Ka/Ks values were calculated using DnaSP v.6 software (Le et al., 2011; Rozas et al., 2017). Based on whether the Ka/Ks ratio was greater than, equal to, or less than one, purifying, neutral, or positive selection was inferred (Zia et al., 2022). Additionally, the divergence time for the duplicated gene pairs was estimated using the formula “t = Ks/2λ×10^–6″ in million years (Mya), where the λ value is 1.5 × 10⁻⁸ for dicots (Fatima et al., 2023; Sadaqat et al., 2023).

The amino acid sequences of the Pyrus PAOs proteins were analyzed for protein-protein interactions (PPIs) using the STRING database (Mering et al., 2003). The top ten interactions were selected for prediction, with an interaction threshold set at 0.4. The PPI network was then visualized using Cytoscape software (Shannon et al., 2003). Additionally, the PANNZER database (http://ekhidna2.biocenter.helsinki.fi/sanspanz/) (Törönen et al., 2018) was utilized to predict the GO enrichment analysis components, including biological processes (BPs), cellular components (CCs), and molecular functions (MFs).

To predict the cis-regulatory elements, the 2 kb sequences upstream of the translation start site of Pyrus PAOs genes were obtained and analyzed using the PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Rombauts et al., 1999).

To understand the expression patterns of Pyrus Dangshansuli PAOs, transcriptomic RNA-seq data from different developmental stages of pear fruit (BioProject: PRJNA309745), under drought stress (BioProject: PRJNA655255), and Fruit hardening disease (BioProject: PRJNA763913) were retrieved from the SRA-NCBI database (https://www.ncbi.nlm.nih.gov/sra). The genome and annotation (GFF) files of Dangshansuli were downloaded from “The pear genomics database” (PGDB; http://pyrusgdb.sdau.edu.cn/) (S. Chen et al., 2023). The quality of reads was assessed using the FastQC tool (Wingett and Andrews, 2018). HISAT2 (Kim et al., 2019) was used to build the indexes of the Pyrus genome, and the high-quality paired-clean reads were then mapped onto the indexed genome. The abundance estimation of gene family members was carried out using StringTie (Pertea et al., 2016). Finally, a heatmap was generated using the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values (Zameer et al., 2021).

To ensure proper functionality, proteins require a three-dimensional (3D) structure. The 3D structures of eight Dangshansuli PAOs (DaPAO1-DaPAO8) were predicted using AlphaFold2 (https://alphafold.ebi.ac.uk/) (Jumper et al., 2021). The accuracy of these predicted structures was assessed through the SAVES server (https://saves.mbi.ucla.edu/) (Zameer et al., 2022) and MolProbity (http://molprobity.biochem.duke.edu/) (Davis et al., 2007). Finally, the UCSF ChimeraX software (Goddard et al., 2018) was utilized to visualize these structures (Neupane et al., 2018; Zia et al., 2024).

A total of 11 genes from the Cuiguan genome (CuPAO), eight from Shanxi Duli (ShPAO), 10 from Zhongai1 (ZhPAO), 10 from Nijisseiki (NiPAO), 10 from Yunhong No.1 (YuPAO), nine from d’Anjou (AnPAO), five from Bartlett v2.0 (BrPAO), and eight from Dangshansuli’ v.1.1 genome (DaPAO) were identified. All of these members were confirmed to have the Rieske [2Fe-2S] iron-sulphur domain superfamily and Pheophorbide a oxygenase (PAO) domain. The protein names of each member were assigned based on their position on chromosomes (Table 1).

The physical and chemical properties of all Pyrus PAO proteins were computed and compared. The protein length ranged from 174–1628 aa and molecular weight ranged from 19.99 to 18.50 kDa. The isoelectric point showed that most of the proteins are basic in nature. The instability index denoted that most of these proteins are unstable. Aliphatic index showed that proteins are good thermostable. GRAVY values indicates that the proteins have hydrophilic behaviour. Subcellular localization analysis revealed that most proteins were found in the chloroplast, with a few in the nucleus, plasma membrane, mitochondria (Table 1).

To examine the evolutionary connections among the PAO members from eight Pyrus genomes, a phylogenetic tree was created using 84 amino acid sequences from 11 species. All PAO proteins were grouped into four sub-family clusters: Group A, B, C, and D. The Group D cluster was the largest, with 24 members. After that the Group A was the second largest with 23 members. Group C and B has the members 22 and 15 respectively (Figure 1).

To understand the evolutionary patterns among Pyrus PAOs, their conserved motifs and gene structures were studied. The gene structure was observed to be highly conserved within members of the same subfamily. The number of exons ranged from 3–18. Group A and B has more numbers of exon as compared to Group C and D (Figure 2A).

All members of each subfamily shared highly conserved motifs. Members of Group A, B, C, and D had exactly the same motif pattern (conserved motifs 2, 3, 4, 5, 8, 11, 12, 14, 16). Motif seven was only present in Group D and the motifs 6, 9, 10, 13, 15, 17, and 18 were present in all members except few members of Group D (Figure 2B). This high conservation of motifs suggests no major differences in the structure and functions of Pyrus PAOs.

To assess the gene distribution pattern of Pyrus PAOs across the 17 chromosomes of each Pyrus genome, their chromosomal gene localization was determined. This analysis revealed an uneven distribution of genes across chromosomes. In the Cuiguan genome, 11 CuPAOs were localized on five out of seventeen chromosomes (Chr3, 8, 11, 13, and 15). The remaining chromosomes did not contain any CuPAOs genes. Similarly, in the Shanxi Duli genome, eight ShPAOs were scattered across five out of seventeen chromosomes. In the Zhongai1 genome, 10 ZhPAOs were distributed across four out of seventeen chromosomes. In the Nijisseiki genome, 10 NiPAOs were distributed across four chromosomes. In the Yunhong No.1 genome, 10 YuPAOs were unevenly distributed across four chromosomes. In the d’Anjou genome, nine AnPAOs were distributed across five chromosomes. In the Bartlett v2.0 genome, five BaPAOs were distributed across four chromosomes. Finally, in the Dangshansuli’ v.1.1 genome, eight DaPAOs members were distributed across five out of seventeen chromosomes (Figure 3).

Gene duplication events were investigated within each Pyrus PAO gene family member. In the Cuiguan genome, eight pairs of genes were found to be duplicated, in which two pairs duplicated through tandem duplication and six pairs duplicated through segmental duplication. The Shanxi Duli genome contained two segmentally duplicated pairs of ShPAOs. In the Zhongai1 genome, three pairs of ZhPAOs were duplicated, with two pairs being tandemly duplicated and the one pair being segmentally duplicated. The Nijisseiki genome had two pairs of segmentally duplicated NiPAOs. In the Yunhong No.1 genome, three pairs of YuPAOs were found, all originating from segmental duplication. The d'Anjou genome exhibited four pairs of duplicated genes, with two pairs being tandemly duplicated and the other two pairs segmentally duplicated. The Bartlett genome contained only one pair of segmentally duplicated genes. The Dangshansuli genome had four pairs in which three pairs are segmentally duplicated and one gene pair is tandemly duplicated (Table 2).

To investigate the evolutionary pressures acting on the duplicated Pyrus PAOs genes, the Ka, Ks, and Ka/Ks ratios were computed for all para-homologous gene pairs. Across the eight genomes, the Ka/Ks ratio varied from 0.48 to 1.05, indicating a mix of positive and negative selection events. The divergence time of the 27 duplicated gene pairs of Pyrus PAOs ranged from 0.5 to 91.51 million years ago (MYA) (Table 2).

A protein-protein interaction (PPI) network of the Pyrus PAOs was constructed to explore the functional diversity among its members. Among the five Dangshansuli members, DaPAO2, DaPAO5, DaPAO6, DaPAO7, and DaPAO8, interactions with several other proteins were observed. The majority of these interactions were identified with PORA (degree 9) and NYC1 (degree 9) proteins, suggesting the potential roles of these Pyrus members in skotomorphogenesis, photomorphogenesis and throughout the plant life under specific light conditions (Figure 4A).

GO enrichment analysis was conducted to determine the molecular functions of Pyrus PAOs in a dynamic context. This analysis classified the genes into three main categories: biological processes (BPs), cellular components (CCs), and molecular functions (MFs). The primary BPs identified included the protein targeting to chloroplast, chlorophyll biosynthetic process, and flower and fruit Development. The genes were found to be present in CCs such as the chloroplast, plastid envelope, and thylakoid membrane. The MFs associated with these genes included 2 Iron, two sulfur cluster binding, metal ion binding, and pheophorbide A oxygenase activity (Figure 4B, Supplementary Material S1).

To gain a deeper understanding of the diverse stress responses exhibited by Pyrus PAOs, the cis-regulatory elements in their promoter sequences were analyzed. Across all genomes, cis-elements associated with stress responses, including light, hormones, and development, were found in abundance. Specifically, G-box, GT1-motif, and GATA-motif (cis-element Box 4) were identified as being involved in the regulation of light-stress. Hormone responsiveness was linked to five cis-elements: P-box, TGA-element, ABRE, CGTCA-motif, and TCA-element. Additionally, the GC-motif, LTR, TC-rich repeats, and MBS were identified as the four cis-elements associated with stress responsiveness. For developmental processes, five elements were involved: CAT-box, MBSI, circadian, HD-Zip 1, and o2-site. The presence of these elements in Pyrus PAO suggests their involvement in hormone, stress, and development-related responses (Figure 5, Supplementary Material S2).

Transcriptome expression data were utilized to assess the expression levels of eight DaPAOs under fruit hardening disease to evaluate the expression of DaPAOs in diseased conditions. DaPAO4 and DaPAO5 were highly expressed in all conditions. DaPAO1 was downregulated in disease conditions as compared to normal condition (Figure 6A). Under drought stress, different expression levels were observed, with DaPAO2 being upregulated. On the other hand, the DaPAOs have no change in expression during normal and drought condition (Figure 6B). Transcriptome expression data was also utilized to assess the expression levels of eight DaPAOs across various tissues, including fruit, leaves, petal, sepal, ovary, stem, and bud. DaPAO1 was highly expressed in ovary. Three genes named as DaPAO2, 3, and 8 were highly expressed in stem. Two genes named as DaPAO5 and DaPAO6 were highly expressed in most of the tissues (Figure 6C). DaPAO1, DaPAO3, DaPAO4, and DaPAO5 showed fluctuated and high expression patterns which make these genes potential candidates for future research.

To delve deeper into the structural and functional diversity, the 3D structures of the eight DaPAO proteins were modeled. These proteins exhibited a high degree of conservation in their 3D structures, with similar patterns of helices and turns. Additionally, all eight proteins shared a similar helical structure at both the C and N termini. This structural conservation implies that the DaPAO proteins may have similar functions. DaPAO7 has a slightly changed structure as compared to the other DaPAOs (Figure 7).