Comparative Genomic Insights into Pathogenicity of Diaporthe amygdali and D. eres Associated with Pear Canker in Korea

1. Introduction

Pears (Pyrus pyrifolia) are widely consumed fruits in East Asia and have gained increasing global attention due to their distinctive texture and nutritional value [1,2,3]. This species contains a wide range of phytochemicals, including flavonoids, phenolic acids such as chlorogenic and caffeic acids, arbutin, and tannin-related compounds. These bioactive constituents have been reported to exhibit antioxidant, anti-inflammatory, and metabolic regulatory activities. Given their substantial nutritional and physiological importance, they have been increasingly exposed to diverse abiotic stresses, including high temperatures, heat waves, unpredictable cold spells, drought, heavy rainfall, and frost damage, as global climate change drives a rise in extreme weather events. These environmental disturbances negatively affect physiological homeostasis, photosynthetic efficiency, and they also increase the infection pressure of pathogenic fungi and other pests [4,5,6,7].

Among the diseases affecting stressed trees, tree canker is one of the most prevalent and destructive [8,9,10]. Tree canker refers to a group of chronic woody diseases characterized by localized necrosis of bark, cambium, and underlying woody tissues on branches and trunks [11,12]. These necrotic lesions can expand over time, disrupting vascular function and impairing the transport of water and nutrients [13]. As cankers enlarge or coalesce, they cause branch dieback, structural weakening of the tree, and a progressive decline in overall tree vigor [11,12].

Pear tree canker diseases are associated with a diverse assemblage of fungal pathogens, most of which belong to the phylum Ascomycota. Several genera, Botryosphaeria [14], Neofusicoccum [15], Valsa [16], Diplodia [17] and Diaporthe [18,19] have been repeatedly isolated from cankered pear tissues. These fungi often act as opportunistic pathogens, colonizing host tissues asymptomatically before switching to a pathogenic lifestyle under favorable environmental conditions or host stress. Among these pathogens, the genus Diaporthe has emerged as a particularly important group associated with canker and dieback diseases in a wide range of woody plants [20,21,22]. Diaporthe species exhibit remarkable ecological and pathogenic diversity, functioning as pathogens, endophytes, or saprophytes depending on host and environmental context [23,24,25]. Genomic studies have suggested that Diaporthe pathogenicity is mediated by a combination of secreted effector proteins, cell wall–degrading enzymes, and diverse secondary metabolites, which collectively facilitate host colonization and tissue necrosis [22,26,27]. Pathogenicity and aggressiveness can vary substantially among Diaporthe species and even among isolates within the same species, indicating that host association alone may not adequately predict disease severity [25,28].

To date, genomic resources for Diaporthe species infecting pear remain extremely limited, and no comparative genomic study has directly addressed the genetic basis underlying pathogenicity differences among Diaporthe species isolated from canker-symptomatic pear trees. In this study, canker-symptomatic pear branches were collected from a pear orchard in Korea, leading to the isolation of D. amygdali and D. eres. To the best of our knowledge, this represents the first report of both D. amygdali and D. eres associated with pear canker in Korea. The pathogenicity of these isolates was subsequently evaluated to confirm their role as causal agents of pear canker. To investigate the genomic basis underlying their pathogenicity, we generated a high-quality genome assembly of the pear-derived D. amygdali strain and compared it with that of D. eres, which exhibited contrasting pathogenicity. Through genome-scale comparative analyses of genome architecture, secreted proteins, candidate effectors, carbohydrate-active enzymes, and secondary metabolite biosynthetic gene clusters, we identified key genomic features potentially associated with virulence differences. This study provides valuable genomic resources, reports new fungal pathogens of pear in Korea, and offers comparative insights into the genetic determinants of virulence in Diaporthe species, establishing a foundation for future functional and epidemiological investigations.

2. Materials and Methods

2.1. Sample Collection and Fungal Isolation

In 2020, severe canker symptoms were pear trees at a farm in Hyeonmae-ri, Seoun-myeon, Anseong-si, Korea. The symptoms were predominantly observed on trunks and branches, particularly around branch nodes, where brown to black discoloration and cracking were evident (Figure 1). To identify the causal agent of pear canker, five branch samples showing symptoms were collected from five individual trees. For fungal isolation, the sampled symptomatic wood tissues were cut into small pieces, surface-sterilized, and placed on potato dextrose agar (PDA). The plates were incubated at 25 °C for five days, and emerging fungal colonies were transferred to fresh PDA plates to obtain pure cultures.

2.2. Morphological Characterization

The morphological characteristics of the pure fungal isolates were initially examined from cultures grown on potato dextrose agar (PDA). For species identification of D. eres and D. amygdali, further examination was performed using cultures grown on malt extract agar (MEA), Czapek yeast extract agar (CYA), and oatmeal agar (OA). Colony morphology was recorded for each medium, whereas the morphology and dimensions of conidiomata, conidiophores, and conidia were examined using an optical microscope (BX53, OLYMPUS, Tokyo, Japan). Light microscopic images were taken for morphological documentation. Fifty measurements were obtained for each of the conidiomata, conidiophores, and conidia.

2.3. DNA Extraction, PCR Amplification, and Phylogenetic Analysis

Mycelia of D. eres and D. amygdali were cultured on PDA plates for five days at 25°C. Genomic DNA was extracted using the Doyle’s method [29]. DNA quality and integrity were evaluated by 1% agarose gel electrophoresis, and DNA concentration was measured with a Qubit 4 Fluorometer (Thermo Fisher Scientific, USA). Only samples with an A260/A280 ratio between 1.8 and 2.0 were used for sequencing. Partial regions of the internal transcribed spacer (ITS), translation elongation factor 1-α (TEF), calmodulin (CAL), and β-tubulin (Tub2) genes were amplified by PCR using the primer pairs ITS1/ITS4 [30], TEF1/TEF728 [31], Cl1/Cl2a [32] and Bt2a/Bt2b [33], respectively, and subsequently sequenced. PCR amplification was performed in a thermal cycler with an initial denaturation at 95 °C for 5 min, followed by 35 cycles of denaturation at 95 °C for 1 min, annealing at 54 °C for 45 s for ITS, TEF, and Tub2, or at 52 °C for CAL, and extension at 72 °C for 1 min, followed by a final extension at 72 °C for 5 min. The amplified PCR products were purified and subjected to sequencing. The nucleotide sequences obtained for ITS, TEF, CAL, and Tub2 were deposited in GenBank under accession numbers PZ867154, PX508498, PX508497, and PX508496, respectively. Individual gene sequences were aligned using MAFFT v. 7 and subsequently concatenated for phylogenetic analyses [34]. Maximum likelihood (ML) analysis was performed using RAxML under the GTR+G model with 1,000 bootstrap replicates [35,36]. Bayesian inference (BI) was conducted using MrBayes v. 3.2.7 under the GTR+G model [37]. The resulting phylogenetic trees were used to determine the phylogenetic placement of strain DUCC20226 among related Diaporthe species. Reference sequences of Diaporthe species were downloaded from the GenBank Database in NCBI (https://www.ncbi.nlm.nih.gov/genbank/).

2.4. Pathogenicity Test

Pathogenicity was evaluated using two-year-old pear (cultivar Shingo) seedlings. Three seedlings were used for each fungal strain and the control treatment. A 1 × 1 cm wound was artificially made on the stem of each seedling using a sterile drill bit, and a 1 × 1 cm PDA plug containing actively growing mycelium of D. amygdali DUCC20226 and D. eres DUCC15256 was placed onto the wound. The inoculated sites were sealed with parafilm to maintain moisture and facilitate infection. Seedlings inoculated with sterile PDA plugs served as controls. Seedlings were maintained in a greenhouse at 24 °C and 50% relative humidity and examined for symptom development two weeks after inoculation. To fulfill Koch’s postulates, fungal isolates were re-isolated from symptomatic inoculation sites and compared with the original strain based on morphological and molecular characteristics.

2.5. Genome Sequencing

Genomic DNA (gDNA) was extracted for genome sequencing. For PacBio sequencing, a SMRTbell library was prepared and sequenced using Sequel II chemistry. For Illumina sequencing, paired-end libraries of 2 × 150 bp were constructed from randomly fragmented gDNA and sequenced on the NovaSeq 6000 platform after adapter ligation. In addition, total RNA was extracted and sequenced to provide transcript evidence for genome annotation. Whole-genome and RNA sequencing was performed by Macrogen Inc. (Seoul, Republic of Korea) using a hybrid approach combining long-read PacBio Sequel II (Pacific Biosciences, USA) and short-read Illumina NovaSeq 6000 (Illumina, USA) technologies.

2.6. De novo Genome Assembly and Evaluation

De novo genome assembly was carried out using FALCON-integrate v2.1.4 [38]. The draft assembly was first polished using Arrow to refine the long-read consensus and then further corrected with Pilon v1.21 [39] using high-quality Illumina short reads to reduce base-level errors and fill small gaps. K-mer analysis was performed with KMC v3.2.1 [40] to estimate genome size, heterozygosity, and repeat content from Illumina short reads. A k-mer size of 21 was applied, and the distribution of k-mer frequencies was used to infer the estimated genome length and proportion of repetitive sequences. Genome completeness was further evaluated with BUSCO v5.8.0 [41] using the sordariomycetes_odb10 lineage dataset. BUSCO identified single-copy and duplicated orthologs based on hidden Markov model profiles to estimate the proportion of conserved genes present in each assembly. The percentages of complete, fragmented, and missing BUSCOs were calculated to determine overall assembly quality.

2.7. Mitochondrial Genome Assembly, Annotation and Phylogenetic Analysis

The raw short reads data generated by HiSeqX-nano sequencer (Illumina) were subjected to de novo assembly done by SOAPdenovo GPL-3.0 [42]. MITOS web server [43] was used to annotate its mitogenome and protein coding genes, tRNA genes and rRNA genes were predicted by MitoZ [44]. The mitochondrial genome was visualized using OGDRAW [45]. The mitochondrial phylogenetic positions of D. amygdali DUCC20226 and D. eres DUCC15256 within Ascomycota were inferred using concatenated sequences of 15 core mitochondrial protein-coding genes from D. amygdali DUCC20226, D. eres DUCC15256, and 21 phylogenetically closely related species. The sequences were aligned using MAFFT v7 [34] Maximum likelihood analysis was performed using RAxML under the GTR+G model with 1,000 bootstrap replicates [35,36], and Bayesian inference was conducted using MrBayes v3.2.7 under the GTR+G model [37].

2.8. Gene Prediction and Functional Annotation

Prior to gene prediction, repetitive elements were identified using RepeatMasker v4.1.5 [46] to soft-mask the assembled genome, preventing false gene predictions within repetitive regions. Gene prediction was performed using BRAKER2 v2.1.6 [47], which combines ab initio prediction with RNA-seq evidence to improve gene model accuracy. RNA-seq reads were aligned to the assembled genome and used as extrinsic evidence during training and prediction. Functional annotation was conducted using BLAST+ v2.7.1 [48] with an E-value threshold of 1e−3 against the Swiss-Prot (release 13 October 2024) and RefSeq (release 13 October 2024) databases. Additional functional and domain annotation was performed using InterProScan v5.59-91.0 [49] and EggNOG v5.0.2 [50] databases to assign Gene Ontology (GO) terms, protein domains, and orthologous group information.

2.9. Prediction of Carbohydrate-Active Enzymes and Secreted Proteins

To further characterize the functional potential of the genomes, carbohydrate-active enzymes (CAZymes) and secreted proteins were systematically analyzed. CAZyme annotation was performed using dbCAN2 meta server (release v14) [51], which integrates multiple approaches including HMMER [52] and DIAMOND [53] for accurate identification of glycoside hydrolases (GH), glycosyltransferases (GT), carbohydrate esterases (CE), polysaccharide lyases (PL), auxiliary activity (AA) enzymes, and carbohydrate-binding modules (CBM). Only genes supported by at least two of the three prediction tools were retained as confident CAZyme candidates. Secretome prediction was carried out using a stepwise computational pipeline combining signal peptide prediction, transmembrane topology analysis, and effector protein identification. Signal peptides were identified using SignalP v6.0 [54], while transmembrane helices were predicted with DeepTMHMM v1.0 [55] to exclude membrane-anchored proteins. Candidate secreted proteins were further screened with EffectorP v3.0 [56] to predict potential effector proteins involved in plant–pathogen interactions. Orthologous genes between D. amygdali and D. eres were identified using OrthoFinder v2.5.5 with default settings [57]. Protein sequences from both genomes were compared in an all-versus-all manner to infer orthogroups and single-copy orthologs.

2.10. Prediction of Biosynthetic Gene Cluster and Network Analysis of BGCs

Prediction of biosynthetic gene clusters (BGCs) was performed using the fungal version of antiSMASH (v8.0) to identify putative secondary metabolite biosynthetic regions in the genome sequences [58]. The analysis was conducted with default parameters, including the detection of core biosynthetic domains and cluster types such as polyketide synthase (PKS), nonribosomal peptide synthetase (NRPS), terpenes, hybrid clusters and others. To investigate the relationships among the predicted BGCs, pairwise similarity and clustering analyses were conducted using BiG-SCAPE (v1.1.9) [59]. This program calculates similarity scores based on domain composition, gene order, and sequence identity to group BGCs into gene cluster families (GCFs). The resulting network was visualized to compare BGC diversity and to identify conserved or species-specific clusters across D. amygdali and D. eres genomes.

3. Results

3.1. Fungal Isolation

A total of 38 fungal strains were recovered from symptomatic tissues. Among them, five strains belonged to the genus Diaporthe, comprising four D. eres and one D. amygdali. The Diaporthe strains were individually distributed across five branches, each from a different pear tree. Isolates DUCC15256 and DUCC20226 were designated as representative strains of D. eres and D. amygdali, respectively, and were compared in subsequent analyses.

3.2. Morphological Characterization of D. eres and D. amygdali

Diaporthe amygdali DUCC20226 formed white to beige, flat colonies with undulate margins on PDA and MEA (Figure 2A-B). On CYA, colonies were flat and radially arranged, with a dark gray center and white to beige margins (Figure 2C). On OA, colonies were circular, radial, flat, and beige in color (Figure 2D). Descriptions of the microstructures of pycnidial conidiomata, conidiophores, and alpha conidia grown on PDA are provided (Figure 2E-H). The strain produced dark brown, aggregated, and immersed pycnidial conidiomata measuring 250–415 μm (mean = 350 μm, n = 50) in diameter. Conidiophores were phialidic and hyaline, measuring 8.5–17.5 × 2.3–4.1 μm (mean = 13.7 × 3.4 μm, n = 50). Alpha conidia were hyaline and ellipsoid, measuring 4.7–9.0 × 2.5–3.9 μm (mean = 7.5 × 3.2 μm, n = 50). Beta conidia were not observed. These morphological characteristics align with known descriptions of Diaporthe spp. associated with canker diseases of woody hosts [18,60].

Diaporthe eres DUCC15256 grew radially and rapidly on OA, PDA, CYA, and MEA at 25 °C, reaching 67–69, 79–80, 80–82, and 90 mm in diameter after 7 days, respectively (Figure 3A-D). Colonies were circular, white, and slightly fluffy, with sparse aerial hyphae. Black pycnidia developed after 26 days on OA and 25 days on PDA, while olive-colored spots appeared on PDA after 17 days. Hyphae in the central colony area showed slight to severe degeneration on CYA and MEA, respectively. Descriptions of the microstructures of pycnidial conidiomata, conidiophores, and alpha conidia grown on PDA are provided (Figure 3E-L). The produced dark brown, aggregated, and immersed pycnidial conidiomata measuring 185–350 μm (mean = 260 μm, n = 50). Hyphae were hyaline, branched, septate, and thin-walled, 2–5 μm wide. Conidiophores were septate, smooth, densely packed, and mostly bottle-shaped or cylindrical, measuring 13–24 × 1.4–3.5 μm (mean = 18.7 × 2.9 μm, n = 50). Alpha conidia were aseptate, smooth, thin-walled, cylindrical with one pointed end, and contained two guttules, measuring 5.6–9.4 × 2.0–3.2 μm (mean = 8.2 × 2.6 μm, n = 50). Beta conidia were aseptate, smooth, thin-walled, filiform, and measured 21.8–34 × 1.2–2.5 μm (mean = 29.8 × 2.1 μm, n = 50). The sexual state was not observed (Fig. 3). These cultural and conidial characters agreed with descriptions of Diaporthe eres [61].

3.3. Molecular Identification and Phylogenetic Analysis

The obtained sequences comprised 583 bp for ITS, 682 bp for TEF, 785 bp for Tub2, and 805 bp for CAL. BLAST searches of the sequences showed that the TEF sequence of DUCC20226 shared 97.94% sequence identity with D. amygdali strain BR1 (GenBank accession: OQ238794), while the ITS, CAL, and Tub2 sequences showed 100% identity with D. amygdali strain CBS 126679, 98.88% identity with D. amygdali strain DAL-3 (GenBank accession: MT006694), and 99.43% identity with strain DAL-4 (GenBank accession: MT006469), respectively. Following multiple-sequence alignment, 463 bp of ITS, 513 bp of CAL, 362 bp of TEF, and 437 bp of Tub2 were retained for phylogenetic analyses. The CAL alignment contained the highest number of variable sites (383), followed by Tub2 (283), TEF (281), and ITS (83). The concatenated alignment consisted of 1,775 bp, of which 1,030 were variable sites. The concatenated phylogenetic analysis placed D. amygdali DUCC20226 within a well-supported D. amygdali clade together with strains CBS 111811, CBS 136969, CGMCC 3.17087, CGMCC 3.17092, and CGMCC 3.17088. Within this clade, DUCC20226 formed a sister group with D. amygdali CGMCC 3.17088 with strong bootstrap support of 98%. The D. amygdali clade was distinct from closely related species, including D. chinensis and D. decedens, supporting the identification of DUCC20226 as D. amygdali (Figure 4).

3.4. Pathogenicity

The pathogenicity test demonstrated that D. amygdali DUCC20226 was pathogenic to pear seedlings. Two weeks after inoculation, black discoloration and water-soaked lesions developed around the inoculation sites, accompanied by collapse and necrosis of the internal stem tissues. The lesions extended above and below the inoculation sites. Seedlings inoculated with D. eres DUCC15256 exhibited similar symptoms, including black discoloration and necrotic lesions around the inoculation sites; however, the lesions were restricted to a narrower area and were less severe than those caused by D. amygdali DUCC20226, indicating weak pathogenicity. No symptoms were observed in the control seedlings (Figure 5). The inoculated fungi were successfully re-isolated from the symptomatic tissues and identified as D. amygdali DUCC20226 and D. eres DUCC15256, respectively, based on their morphological and molecular characteristics. These results fulfilled Koch’s postulates and confirmed the pathogenicity of D. amygdali DUCC20226 and the weak pathogenicity of D. eres DUCC15256 on pear seedlings.

3.5. Genome Assembly and Genomic Characteristics

The assembly of D. eres DUCC15256 resulted in 64 contigs with a total length of 61.16 Mb and an N50 of 4.08 Mb (GenBank accession: JCDZBA000000000), whereas D. amygdali DUCC20226 comprised 15 contigs totaling 49.82 Mb with an N50 of 6.15 Mb (GenBank accession: GCA_032433605). The GC content was 50.3% in D. eres and 52.4% in D. amygdali. The largest contigs reached 6.62 Mb and 10.15 Mb in the respective assemblies. Illumina reads were mapped to the assembled genomes, achieving mean sequencing depths of 149.1× for D. eres and 146.8× for D. amygdali, supporting the high quality of the genome assemblies (Table 1).

K-mer analysis estimated genome sizes of 61.16 Mb for D. eres and 49.82 Mb for D. amygdali, with heterozygosity levels of 0.159 and 0.089, respectively. BUSCO evaluation indicated high completeness in both genomes, reaching 98.7% for D. eres and 98.9% for D. amygdali (Table 1).

The genome of D. eres contained 32,075 tandem repeats with a total length of 2,323,842 bp. It possessed 260 tRNA genes that together comprised 23,442 bp. The genome included 47 SINEs measuring 3,674 bp and 299 LINEs measuring 22,761 bp. Four LTR elements were identified with a combined length of 302 bp, and 59 DNA transposons accounted for 5,573 bp. Six unclassified repeat elements were detected with a total of 468 bp. It also contained 106 small RNA elements measuring 26,705 bp, two satellite sequences totaling 165 bp, 13,949 simple repeats measuring 604,368 bp, and 1,318 low-complexity regions totaling 64,175 bp.

The genome of D. amygdali contained 8,374 tandem repeats with a total length of 465,175 bp. It had 184 tRNA genes comprising 17,234 bp. The genome included 46 SINEs measuring 3,171 bp and 257 LINEs measuring 19,891 bp. Six LTR elements were identified with a combined length of 468 bp, and 77 DNA transposons accounted for 6,358 bp. Seven unclassified repeat elements were detected with a total of 565 bp. The genome also contained 71 small RNA elements measuring 6,881 bp, five satellite sequences totaling 371 bp, 11,714 simple repeats measuring 479,384 bp, and 957 low-complexity regions comprising 45,313 bp (Table 2).

3.6. Structural Features of the Mitochondrial Genome

The length of the D. amygdali DUCC20226 mitochondrial genome (GenBank accession: OQ349375) is 49,922 bp. The mitochondrial genome of D. amygdali encodes 15 protein-coding genes, two rRNA genes, and 18 tRNA genes. The protein-coding genes include the NADH dehydrogenase subunits (nad1, nad2, nad3, nad4, nad4L, nad5, and nad6), cytochrome b (cob), the cytochrome c oxidase subunits (cox1, cox2, and cox3), the ATP synthase subunits (atp6, atp8, and atp9), and the ribosomal protein (rps3). The rRNA region comprises rrnL and rrnS, and the genome also contains 18 tRNA genes (trnW-CCA, trnR-ACG, trnK-UUU, trnL1-UAG, trnY-GUA, trnF-GAA, trnV-UAC, trnE-UUC, trnA-UGC, trnS1-GCU, trnS2-UGA, trnI-GAU, trnM-CAU, trnG-GCC, trnD-GUC, trnN-GUU, trnP-UGG, and trnT-UGU) (Figure 6). The protein-coding region of the mitochondrial genome is 17,293 bp, accounting for 34.64% of the total. The GC content of the entire mitochondrial genome is 32.17%.

The length of the D. eres mitochondrial genome is 93,346 bp (GenBank accession: OQ349376). The mitochondrial genome of D. eres encodes 15 protein-coding genes, two rRNA genes, and 22 tRNA genes. The protein-coding genes include the NADH dehydrogenase subunits (nad1, nad2, nad3, nad4, nad4L, nad5, and nad6), cytochrome b (cob), the cytochrome c oxidase subunits (cox1, cox2, and cox3), the ATP synthase subunits (atp6, atp8, and atp9), and the ribosomal protein (rps3). The rRNA region comprises rrnL and rrnS, and the mitochondrial genome also contains 22 tRNA genes (trnY-GTA, trnN-GTT, trnW-TCA, trnK-TTT, trnD-GTC, trnS1-GCT, trnG-TCC, trnV-TAC, trnI-GAT, trnS2-TGA, trnP-TGG, trnT-TGT, trnE-TTC, trnM-CAT, trnL2-TAA, trnA-TGC, trnF-GAA, trnL1-TAG, trnQ-TTG, trnR-TCT, trnC-GCA, and trnH-GTG) (Figure 7). The protein-coding region of the mitochondrial genome is 14,640 bp, accounting for 15.68% of the total. The GC content of the entire mitochondrial genome is 31.91%.

3.7. Phylogenetic Analysis Based on Mitochondria Genome Sequence

Phylogenetic analysis based on the combined sequences of 15 conserved mitochondrial protein-coding genes showed a clear and well-supported tree. Based on the phylogenetic tree, the family Diaporthaceae, including D. amygdali DUCC20226, was clearly separated from Nectriaceae, Ophiostomataceae, and Cryphonectriaceae, each forming well-supported distinct lineages. In contrast, Diaporthaceae was not clearly resolved from Cytosporaceae, which was represented solely by Valsa mali in this analysis. Within Diaporthaceae, D. amygdali DUCC20226 clustered with D. phaseolorum ZM33-4, D. eres ZM79-3, and D. nobilis. However, it formed an independent monophyletic lineage within this clade (Figure 8).

3.8. Gene Prediction and Analysis of Putative Secreted Proteins

Gene prediction identified 16,506 and 14,825 protein-coding genes in D. eres and D. amygdali, respectively. Among them, 1,950 proteins in D. eres and 1,656 in D. amygdali were predicted to contain N-terminal signal peptides, indicating their potential for secretion via the classical secretory pathway. A total of 1,635 and 1,357 proteins were estimated to have a high probability of being secreted after excluding those with transmembrane domains or GPI anchors. Effector prediction further identified 494 candidate effectors in D. eres and 442 in D. amygdali.

3.9. Orthologous Gene Composition and Candidates of Species-Specific Effector

Orthologous gene clustering revealed a high degree of genomic conservation between D. eres and D. amygdali. A total of 31,331 protein-coding genes were identified across the two genomes, of which 29,641 genes (94.6%) were assigned to 11,949 orthogroups, indicating that most genes are shared between the two species. Only 1,690 genes (5.4%) were unassigned, suggesting limited gene novelty or lineage-specific expansion. Among the identified orthogroups, 11,722 (98.1%) were present in both species, and 9,402 were classified as single-copy orthogroups. A total of 227 species-specific orthogroups containing 711 genes (2.3% of all genes) were detected. Orthology-based clustering revealed 575 protein-coding genes that were specific to D. eres and 184 that were unique to D. amygdali (Table 3).

In D. eres, SignalP and TMHMM detected 42 proteins containing N-terminal signal peptides among the 575 species-specific proteins. EffectorP predicted 10 of these 42 proteins with a high probability of being secreted as putative effectors. BLASTP annotation predicted one effector as a pectate lyase, while the remaining nine proteins were considered hypothetical.

In D. amygdali, SignalP and TMHMM detected 11 proteins among the 184 species-specific proteins. EffectorP predicted four of these 11 proteins as candidate effectors. BLASTP annotation showed that two corresponded to putative pathogen effector proteins categorized as candidate necrosis-inducing effectors, whereas the other two remained unannotated. The lengths of the two proteins annotated as putative necrosis-inducing factors were 202 aa and 191 aa, and their structural features suggested that they likely function as apoplastic effectors. In contrast, the two unannotated proteins, which were 77 aa and 161 aa in length, were predicted to be cytoplasmic effectors. A search of these candidate effector sequences using PHI-BLAST showed no hits for all four proteins. RNA sequencing of total RNA extracted from D. amygdali grown on PDA revealed no detectable reads for these four species-specific effector genes.

3.10. CAZyme Composition

A total of 961 carbohydrate-active enzyme (CAZyme) genes were identified in D. eres, while 863 were detected in D. amygdali. Both species showed similar overall CAZyme distributions. In D. eres, the CAZyme repertoire consisted of 437 GHs, 254 AAs, 141 GTs, 62 CEs, 40 PLs, and 27 CBMs. In D. amygdali, 405 GHs, 226 AAs, 124 GTs, 56 CEs, 39 PLs, and 13 CBMs were identified.

Among the GH families detected in both species, GH43 was the most abundant, with 38 copies in D. eres and 33 in D. amygdali, associated with β-xylosidase and α-L-arabinofuranosidase activities for hemicellulose degradation. GH5 followed, represented by 28 copies in D. eres and 24 in D. amygdali, containing endo-β-1,4-glucanase activity involved in cellulose breakdown. GH3 was similarly abundant in both species, with 25 and 26 copies, respectively, encoding β-glucosidases that act on cellulose-derived sugars. GH18 contained 22 and 14 copies (chitinase activity), while GH16 was represented by 21 copies in each species (β-1,3-glucanase activity). Beyond these dominant groups, additional GH families such as GH28, GH13, and various mannosidases, xylanases, rhamnosidases, and glucanases were also present, indicating a broad enzymatic capacity for degrading diverse plant cell wall components.

PL1 was the most abundant polysaccharide lyase family in both species, with 19 copies in D. eres and 18 in D. amygdali, associated with pectate lyase activity involved in pectin degradation. PL4 and PL3 followed, represented by 8 and 7 copies in D. eres, and 7 and 8 copies in D. amygdali, respectively (pectate and rhamnogalacturonan lyase activities). Lower-frequency families such as PL9, PL26, PL11, and PL42 were detected in both species, whereas PL27 was found exclusively in D. eres.

3.11. Prediction of Biosynthetic Gene Cluster and Network Analysis of BGCs

AntiSMASH analysis of the whole genomes of D. amygdali and D. eres predicted 87 and 116 biosynthetic gene clusters (BGCs) for secondary metabolites, respectively. Both species contained clusters corresponding to PKS, NRPS, terpene, indole, and NRPS-like types, with PKS-related clusters representing the largest proportion.

In D. amygdali, most of the predicted BGCs were classified as type I PKS (T1PKS), with several showing similarity to known polyketide metabolites such as alternariol, alternapyrone, depudecin, zearalenone, and scytalone. The NRPS and NRPS–PKS hybrid clusters were predicted to be associated with peptide metabolites such as fusarin C, neosartorin, and BAB/BAA, while the terpene clusters were linked to compounds including PR-toxin, phyllocladan-16α-ol, and copalyl diphosphate. In addition, clusters related to cytokinin, isocyanide, and siderophore biosynthesis were detected. MIBiG comparison revealed that many of these clusters shared genetic organization with known pathogenic polyketide and terpene biosynthetic pathways.

In D. eres, 29 additional BGCs were predicted compared to D. amygdali, with particularly high diversity among PKS and terpene clusters. The predicted major metabolites included alternapyrone, fusarin C, nectriapyrone, zearalenone, and depudecin-like compounds, as well as terpene derivatives such as clavaric acid, koraiol, and squalestatin S1. The NRPS–PKS hybrid clusters showed combinations similar to those responsible for the biosynthesis of fusarin C, wortmanamide A/B, and AKML A–D. In addition, unique clusters belonging to siderophore, RiPP, and isocyanide classes were identified.

Network analysis using BiG-SCAPE revealed that only 16 BGCs shared similar cluster or domain architectures between D. eres and D. amygdali, suggesting potential conservation of specific secondary metabolite pathways (2 terpene, 7 NRPS, 2 NRPS–T1PKS hybrid, 3 T1PKS, 1 T3PKS, and 1 indole-derived cluster). In contrast, the remaining 100 BGCs in D. eres and 71 BGCs in D. amygdali were classified into distinct groups, indicating that the two species likely synthesize largely divergent sets of secondary metabolites.

4. Discussion

This study provides the first report of Diaporthe amygdali and D. eres as fungal species associated with pear canker on Pyrus pyrifolia in Korea, based on morphological, molecular, and pathogenicity evidence. In addition, this study reports high-quality genome assemblies of D. eres and D. amygdali isolated from canker-symptomatic pear trees and aims to elucidate the molecular genetic basis underlying the marked differences in pathogenicity between the two species through comparative genomic analyses, despite their high degree of genomic conservation.

Several pathogens have been reported to cause important diseases of Pyrus pyrifolia in Korea. Erwinia pyrifoliae is known to cause shoot blight [62], whereas Erwinia amylovora causes fire blight [63]. Botryosphaeria dothidea has been reported as a causal agent of black rot and canker [65], Nectria cinnabarina causes Nectria canker [64], while Diaporthe ambigua is associated with dieback of pear trees [64]. In addition, Colletotrichum acutatum has been reported as a causal agent of anthracnose [65], and Venturia nashicola [66] and Gymnosporangium asiaticum [67] are important fungal pathogens associated with diseases of pear leaves and fruits.

In contrast, D. amygdali has not previously been reported as a pathogen on P. pyrifolia in Korea. However, D. amygdali has previously been reported as an important plant pathogen of several woody fruit crops worldwide. It is well known as a causal agent of twig canker and shoot blight of almond (Prunus dulcis) [21] and peach (Prunus persica) [68]. It has also been reported to cause twig canker and shoot blight of nectarine [20]. Therefore, the present study expands the known diversity of fungi associated with pear diseases in Korea by providing the first evidence of D. amygdali and D. eres as pathogenic fungi associated with pear canker. The pathogenicity assays, together with the successful re-isolation and molecular identification of the inoculated fungi, support their pathogenicity on pear and their association with pear canker.

Interestingly, D. amygdali exhibited stronger pathogenicity toward pear than D. eres under the tested conditions. Such variation in pathogenicity against seedlings of five Asian pear cultivars among isolates of the same or closely related 19 Diaporthe taxa has also been documented in previous studies [18]. But this previous work did not include ‘Singo’ cultivar which is most popularly cultivated in Korea and fairly in Japan. This study extends the pathogenicity of D. eres to Singo’ cultivar. In addition, Guo isolated 20 Diaporthe species from shoot canker pear. But there was no D. amygdali [18]. To our knowledge, this is the first time to isolate D. amygdali from canker of P. pyrifolia.

Additionally, Dissanayake and Guo reported substantial variation in aggressiveness among isolates within the D. eres complex, with different phylogenetic lineages producing markedly different lesion sizes on peach shoots and pear seedlings [18,25,28]. This suggests that pathogenicity in Diaporthe may depend on lineage- or isolate-specific genetic traits rather than species identity alone. Furthermore, the cross-family infectivity of D. amygdali from Rosaceae stone fruits (Prunus) to pome fruits (Pyrus) underscores its host adaptability and potential threat to fruit orchards.

D. eres exhibited a relatively larger genome size, a higher proportion of repetitive elements, and a greater number of predicted genes compared with D. amygdali. However, under the conditions tested in this study, D. eres showed weaker pathogenicity toward pear than D. amygdali, suggesting that a larger genome does not necessarily confer greater pathogenicity. Although genome enlargement has been proposed to facilitate diversification and adaptive evolution of pathogenicity-related genes [69], our results instead revealed an opposite trend, in which the more virulent species possessed a smaller genome with fewer repetitive elements. This pattern is consistent with previous comparative genomic studies indicating that enhanced pathogenicity can be achieved not through genome expansion but rather via selective retention of specific functional gene sets and reduction of nonessential genes [70]. In some pathogenic fungi, genome reduction or the retention of a relatively compact genome has been proposed as an adaptive feature associated with specialized pathogenic lifestyles. For example, Geosmithia morbida, a host-specific tree pathogen, possesses a smaller genome than closely related nonpathogenic species, and genome reduction was proposed as a potential adaptive force during the evolution of its specialized lifestyle [71]. Similarly, the smaller genome and lower repetitive content observed in D. amygdali, despite its high pathogenicity toward pear, may indicate that pathogenicity in this lineage is not necessarily associated with genome expansion.

The comparative analysis between Diaporthe amygdali DUCC20226 and D. eres highlighted a conserved core set of mitochondrial genes alongside striking variations in total genome size and non-coding proportions. Similar patterns of structural variation and core gene conservation have been widely observed across the genus Diaporthe [72]. In our study, D. amygdali exhibited a substantially smaller genome size (49,922 bp) compared to D. eres (93,346 bp), while retaining a higher proportion of protein-coding region (34.64% vs. 15.68%). Xie similarly demonstrated that Diaporthe mitogenomes display exceptional size plasticity, ranging from 53,646 bp in D. longicolla to 108,865 bp in D. sojae, while maintaining an identical complement of 15 core protein-coding genes (PCGs), two rRNAs, and a set of tRNAs [72]. Their findings revealed that intron gain/loss events and variations in intergenic or intronic regions—rather than expansion of core PCGs—drive the broad spectrum of mitogenome sizes within Diaporthe. The lower proportion of coding sequence in D. eres relative to D. amygdali aligns with this phenomenon, indicating that the larger size of D. eres is primarily driven by intron accumulation or intergenic sequence expansion. Furthermore, we observed a discrepancy in tRNA gene numbers (18 in D. amygdali vs. 22 in D. eres). While standard fungal mitogenomes typically harbor a relatively conserved tRNA set to support mitochondrial translation, variations in tRNA copy number and unidentified open reading frames (uORFs) are frequently reported during interspecific fungal comparative genomics [72]. Despite these genomic rearrangements and intron dynamics, the core functional components required for oxidative phosphorylation remain strictly conserved under purifying selection across Diaporthales [72]. Taken together, our results corroborate that Diaporthe mitochondrial genomes preserve high functional conservation in core energetic machinery while undergoing extensive dynamic expansion and contraction in non-coding architectures.

Secretome and effector analyses revealed that although both species harbor a large repertoire of secreted proteins and candidate effectors, they display striking differences in pathogenic outcomes. This finding indicates that virulence is not determined by the sheer size of the effector repertoire, but rather by the functional contribution of specific effector combinations. This observation is in agreement with experimental evidence demonstrating that host-specific functional effectors play a decisive role in pathogenicity [73,74]. Notably, among the species-specific effectors identified in D. amygdali, two were annotated by BLAST-based analysis as putative necrosis-inducing factors, suggesting a potential direct involvement in host tissue necrosis. Necrosis-inducing elements are known to play a critical role in triggering disease development during the early stages of infection [75,76,77], and thus may represent candidate virulence determinants contributing to the enhanced pathogenicity of D. amygdali. Furthermore, the lack of detectable expression under PDA culture conditions suggests that these genes may be subject to condition-dependent regulation and may be preferentially expressed during host infection (in planta). This possibility is consistent with the host- and infection-dependent expression patterns reported for effector genes in various plant-pathogenic fungi [78,79], emphasizing the necessity of transcriptomic profiling during host tissue colonization to validate their functional roles.

CAZyme analysis showed that D. eres possesses a greater number of carbohydrate-active enzymes, including plant cell wall–degrading enzymes, than D. amygdali, yet exhibits lower pathogenicity, as reported in this study. Comparative genomic analyses across diverse fungi have similarly demonstrated that the total number of hydrolytic enzymes does not consistently explain virulence levels; instead, pathogenicity is influenced by lifestyle, enzyme composition specificity, and regulatory patterns of functional gene groups [80]. These findings suggest that pathogenicity is determined not by the overall abundance of CAZymes, but rather by the expression of specific glycoside hydrolase and polysaccharide lyase families, as well as their functional interactions with effectors and toxins [81,82].

Secondary metabolite analysis based on antiSMASH predicted broadly similar classes of polyketides and other metabolites in both species. However, BiG-SCAPE analysis revealed that only 16 biosynthetic gene clusters (BGCs) shared similar cluster architectures between the two species. This indicates that despite the presence of homologous BGCs, differences in gene order and domain composition may substantially limit functional equivalence between clusters [83]. Moreover, even when identical BGC types are present, variation in enzymatic modifications and tailoring reactions during downstream biosynthetic pathways can result in pronounced differences in the toxicity, host reactivity, and functional roles of the final metabolites during infection [84,85]. Such metabolic divergence therefore represents a plausible contributor to the observed differences in pathogenicity between the two species.

Taken together, despite the high level of genomic conservation between D. eres and D. amygdali, the observed difference in pathogenicity was not explained by overall genome size, mitochondrial genome size, gene number, or CAZyme abundance alone. Instead, species-specific effector candidates (particularly putative necrosis-inducing factors), divergent secondary metabolite BGC repertoires, and differences in the composition of specific CAZyme families may contribute to the contrasting pathogenic phenotypes of the two species. These findings suggest that qualitative differences in a subset of functional genes and biosynthetic pathways may be more relevant to pathogenicity than overall genomic size or gene content. However, the genomic comparisons in this study were based on a single isolate representing each species, and the predicted functions of the candidate effectors and secondary metabolite BGCs were not experimentally validated. Therefore, these genomic features should be considered potential contributors to pathogenicity rather than direct determinants of virulence. Further studies involving multiple isolates, infection-stage dual RNA-seq profiling, functional validation via gene knockout of candidate virulence genes, and experimental characterization of secondary metabolites will be necessary to clarify their direct roles in the pathogenicity of D. amygdali.

5. Conclusions

This study provides the first report of Diaporthe amygdali and D. eres as fungal pathogens causing pear canker on Pyrus pyrifolia in Korea, expanding the known etiology of pear diseases. Despite high genomic conservation and conserved mitochondrial core gene composition, the two species exhibited contrasting pathogenic phenotypes, with D. amygdali causing significantly more severe necrosis than D. eres. Comparative genomic analyses revealed that higher virulence in D. amygdali is not driven by overall genome size, gene count, or CAZyme abundance. Instead, specific genomic features—such as candidate species-specific necrosis-inducing effectors, distinct carbohydrate-active enzyme family profiles, and divergent secondary metabolite biosynthetic gene clusters—likely contribute to its enhanced pathogenicity. Overall, these findings provide high-quality genomic resources for pear-associated Diaporthe species and offer valuable insights into the genetic determinants of virulence, serving as a foundation for future functional studies and disease management strategies.