During the 2022 growing period, S. pistaciarum isolates were obtained in June from pistachio leaves showing leaf spot symptoms (cv. Aeginis) across three distinct sites located in the areas of Amfikleia, Fthiotida (38.629963 N, 22.681427E; code FTH), Kiparissia, Larissa (39.5138684 N, 22.5733432E; code LAR) and Kitros, Pieria (40.3785945 N, 22.607257E; code PIE). The trees were over 20 years old and trained under the standard open vase system to an average height of 10 m. A standard fungicide program based on the fungicides boscalid + pyraclostrobin (Signum 26,7/6,7 WG; BASF), fluxapyroxad (Sercadis 30 SC; BASF), pyrimethanil (Scala 40 SC; BASF), dodine (Syllit 544 SC; UPL), and copper (Cuprofix Ultra 40 WG; ADAMA) was applied on all orchards in spring. Each diseased leaf sample was transferred to the lab in a separate moist polyethylene bag to prevent cross-contamination and stored at room temperature for 72 h to induce the formation of cirrhi in the pycnidia. From each sample a single-spore isolate was obtained by slight touching a flamed wire loop onto a freshly formatted cirrhus of conidia from one spot per leaf randomly picked up. A sparse pycnidiospore suspension in 0.5 mL distilled water was prepared and spread onto Potato Dextrose Agar (PDA) in Petri dishes amended with 100 mg/L of streptomycin sulfate. After 48 h of incubation at 24°C in the dark, individual germinated single-spore conidia were transferred in glass test tubes with potato dextrose agar and stored until use (Figure 1). Pycnidiospores were hyaline, curved, with obtuse ends, having 1 to 5 septa (Chitzandis, 1956) typically of S. pistaciarum. In total, 27 isolates were obtained (with nine isolates representing each region; Table 1) and they were maintained in PDA tubes at 4°C. Because of the very slow growth of the hyphae, in order to obtain larger colonies in plates for further analysis each representative plate was seeded with 0.4 mL of a conidial suspension from each isolate (Dhingra and Sinclair, 1985). After 4 weeks of incubation in the dark, fungal DNA was extracted from conidiomata (Figure 1) of all isolates using the Quick-DNA™ Fungal/Bacterial Miniprep Kit (Zymo Research).

Fungal DNA from all isolates was randomly fragmented into 350 bp inserts, genomic libraries were constructed using the Novogene NGS DNA Library Prep Set (Cat No.PT004) and sequenced by a PE150 sequencing strategy. Illumina reads were trimmed by cutadapt v3.7 (-a AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA-A AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG --minimum-length = 50 -- max-n = 3 -n2 -q30; Martin, 2011). All isolates were assembled by SPAdes v3.15.4 (--only_assembler –cov-cutoff auto; Bankevich et al., 2012). Assembly quality and gene set completeness was predicted via BUSCO v5.4.6 (Simao et al., 2015). Initial genome assembly metrics were assessed for contiguity and BUSCO completeness (Supplementary Data Sheet 1), from which the “SPF8” isolate was selected as a reference isolate. Reads for isolate ‘SPF8’ were merged and gap-filled with BBMerge (Bushnell et al., 2017) then de novo assembled with SPAdes as above. Genome assemblies of all isolates were further scaffolded were possible using pairwise comparisons vs. other isolates with ragtag v2.1.0 (scaffold: –remove-small -f 75 -r, merge; Alonge et al., 2022). Mitochondrial genomic DNA (mtDNA) was assembled via MitoZ v3.6 (Meng et al., 2019), and assembled mtDNA contigs were appended to the genome assemblies, while previously-assembled sequences matching mtDNA contigs were removed. Contamination checks for each assembly was performed via BLAST v2.12.0 (blastn -max_target_seqs 5 -evalue 1e-100 -perc_identity 95; Altschul et al., 1990) to the UNITE database (Abarenkov et al., 2010). Overall G:C content and AT-rich compartments were predicted with OcculterCut v1.1 using the reference isolate SPF8 genome assembly (Testa et al., 2016).

Repetitive sequences were predicted using Dfam TE Tools 1.88 (Lerat et al., 2016). Gene prediction was performed across all isolates in 2 rounds. In the first round, funannotate v1.8.15 (predict; max_intronlen 1000; Palmer and Stajich, 2020) was used to predict genes for all isolates. In order to provide further evidence for the accurate support of gene loci and their exon boundaries, transcriptomic data was obtained for the isolate SPF8. Thus, total mRNA was extracted in triplicate from the mycelial phase of this isolate grown on PDA medium from a fresh two-weeks-old culture. The Quick-RNA™ Fungal/Bacterial Miniprep kit (Zymo Research) was employed, and RNAseq reads were generated through Illumina platform (PE150 Novaseq 6000) using the Novogene NGS RNA Library Prep Set (PT042). RNAseq alignment vs. the reference isolate was performed with HiSAT2 v2.2.1 (max-intronlen 5000; dta; Kim et al., 2019), converted to GFF via Stringtie v2.2.1 (Levy Karin et al., 2020) and provided as input to funnannotate. Funannotate-predicted proteomes were clustered into orthogroups (including singletons) with ProteinOrtho v6.3.1 (selfblast; singles; Lechner et al., 2011). A representative pan-genome proteome dataset was selected from the longest member of each orthogroup, and used as input to MetaEuk Release 6-a5d39d9 (--easy-predict –max-intron 500; Levy Karin et al., 2020) for a second round of gene prediction in all isolates, and orthogroups with MetaEuk matches were retained as a representative ‘pan-genome’ proteome set. Functional annotations were predicted across the pan-genome proteome dataset with InterProScan v5.63-95.0 (Jones et al., 2014). Secretion and effector-like properties were predicted using Predector v1.2.7 (Jones et al., 2021). Trophic niche (i.e., biotrophy, necrotrophy, etc.) was predicted via CATAStrophy v0.1.0 (Hane et al., 2020; using HMMER 3.3) vs. dbCAN v10. Mating type genes were screened vs. a representative dataset (Wilken et al., 2017) using MetaEuk (as above). Secondary metabolite synthesis gene clusters were predicted in the reference isolate assembly (SPF8) with antiSMASH v6.1.1 (Blin et al., 2021). Fungicide-resistance mutations were predicted with the fungicide-resistance allele screening tool (FRAST) (Oliver et al., 2024).

Gene ortholog group presence-absence variation (PAV) across the S. pistaciarum pan-genome was determined based on matches to a pan-genome-wide representative proteome dataset using MetaEuk (see above). Genome sequencing reads were aligned to the SPF8 assembly with BWA v0.7.17-r1198-dirty (Li and Durbin, 2009). Genome sequencing reads were processed fastqToSam to generate raw ubam inputs, and passed through MarkIlluminaAdapters, MarkDuplicates, and combined with BWA alignmed bam data with MergeBamAlignment (-CREATE_INDEX true -ADD_MATE_CIGAR true) to generate alignmed uban inputs for variant calling. Variant calling was performed with GATK v4.2.6.1 (HaplotypeCaller -ERC GVCF –minimum-mapping-quality 20 –min-base-quality-score 20 -G StandardAnnotation -G AS_StandardAnnotation -G StandardHCAnnotation; McKenna et al., 2010). Bi-allellic single nucleotide polymorphisms (SNPs) were filtered and randomly selected for 1 per 5 kb with BCFTools v1.15 (+prune -w 5,000 bp -n1 -N rand; Danecek et al., 2021), and used to generate an unrooted distance-based tree via IQTree v2.2.2.7 (1,000 iterations, -bb 1000 -alrt 1000; Minh et al., 2020). The tree was visualized alongside selected and summarized predicted protein function data (mating types, effector-like proteins) with iTOL (mid-point root; Letunic and Bork, 2021).

The reference isolate (SPF8) assembled into a total of 43.1 Mb. All isolate assemblies exhibited ~70-120X sequencing depth (ex. SPP26), an average N50 of 80.4 kb, an average contig number of 9,084, and an average of 42% and 20% gene-coding and repetitive regions, respectively (Table 2) with LTR retrotransposons being the most common repeat type (Table 2; Supplementary Data Sheet 1). The genome metrics above were comparable to pan-genome datasets from the sister species Zymoseptoria tritici, which had an average core chromosome coverage of ~30X, an average N50 of 84.7 kb and average contig number of 1,088 (Feurtey et al., 2023). CATAStrophy (Hane et al., 2020) predictions for S. pistaciarum indicated non-haustorial hemibiotrophy (Table 3), which may correspond with prior reports of a 2 week latent phase (Crous et al., 2013). OcculterCut analysis of G:C content of the reference isolate assembly (Figure 2B) also indicated a bi-modal G:C distribution typical of hemibiotrophic fungi (Testa et al., 2016), with the AT-rich peak at approximately 38.6% G:C containing over a third (38.3%) of the genome length. Blastn vs. UNITE (Abarenkov et al., 2010) indicated >99.25% identity to the ITS region for all isolates to Mycosphaerella pistaciarum, however isolate 6 was also contaminated with plant DNA, and isolate 26 was contaminated with a Basidiomycete matching Tomentella fuscograulosa [NUC: UDB028526] (Supplementary Data Sheet 2). Despite varying levels of contamination in these two assemblies, the pan-genome survey approach of this study (Figures 2C,D) was capable of focusing on genome features relevant to the S. pistaciarum population as a whole, as genes specific only to isolates 6 and 26 were not the subject of further study. Assessment of the reference assembly via BUSCO indicated 98.7% completeness relative to the capnodiales_odb10 dataset (2020-08-05, genomes:13, BUSCOs:3578), with 3738/3786 BUSCOs detected completely (Supplementary Data Sheet 3).

Based on analysis of SNP-based distance trees the two regions—Pieria (north) and Larissa (central)—were distinct, whereas Fthiotida (south) was similar to both Pieria and Larissa regions (Figures 2A,C). While the geographic range and number of isolates sampled in this study is limited, based on the available data, we speculate that various events of long-distance dispersal (Golan and Pringle, 2017) caused by anthropogenic influences, spore dispersal, and infected nursery stock may be the cause of this clustering of isolates from the Fthiotida region with a relatively higher genetic divergence than their expected geographic distance. Such a clustering, where the populations were not strictly clustered based on their geographical origins, has also been reported in Zymoseptoria tritici populations (Mekonnen et al., 2020; Chedli et al., 2022). It may also indicate that these isolates belong to a population that might has experienced range expansions across the other two regions, which has resulted in admixture complemented by gene flow (Rogers and Rogers, 1999), gradually reducing the isolation by distance (Mills et al., 2007) of Fthiotida isolates. On the other hand, the extensive clustering into geographically-restricted lineages of isolates from Larissa and Pieria regions may indicate that these populations experienced a geographic isolation leading to local adaptation (Bazzicalupo, 2022). Future expansion of pan-genomic resources with a broader sampling to include surrounding regions or a globally-representative set of isolates of S. pistaciarum, might be able to confirm if isolates of the southern Fthiotida region represent a “wilder” population, relative those of the genetically narrower Pieria and Larissa regions. At present it is not clear whether pathways of movement through neighboring regions that were not sampled may have contributed to the transportation history of these isolates. Assemblies for all isolates possessed two versions of MAT1-1 type mating type genes, homologous to MAT-1-5/COX13 [ABS19615] and MAT1-1-6/APN2 [ABS19616] (Wilken et al., 2017; Supplementary Data Sheet 4). No MAT1-2 homolog was detected, suggesting limited potential for meiotic recombination across these three regions.

AntiSMASH prediction of secondary metabolite synthesis gene clusters (SMCs) in the reference isolate (SPF8) assembly, indicated homologs of 4 T1PKS (polyketide synthase) regions potentially encoding toxins similar to cercosporin (Chen et al., 2007), fusarubin/oxyjavanicin (Chowdhury et al., 2017), aspyridone A (Bergmann et al., 2007), fumonisin (Kamle et al., 2019), as well as melanin which may have a role in strengthening the pathogen cell wall (Jacobson, 2000). Also predicted were 15 non-ribosomal peptide synthase (NRPS)-like clusters potentially encoding the cytotoxic cyclic peptide serinocyclin (Krasnoff et al., 2007) and other unknown metabolites, and 5 terpene sythase clusters (Supplementary Data Sheet 5). Accurate prediction of secondary metabolite biosynthesis products remains challenging and requires further validation, however the clusters for serinocyclin and melanin were highly conserved. Overall, the predicted SMC profile of S. pistaciarum was similar to that previously reported for the relatively closely-related species Zymoseptoria tritici (Hassani et al., 2022), however S. pistaciarum had more NRPS SMCs whereas Z. tritici had more PKS SMCs (Supplementary Data Sheet 5).

Although testing of differential fungicide efficacy and resistance across S. pistaciarum isolates was beyond the scope of this current study, it was possible to predict potential fungicide resistance mutations from pan-genomic data. Analysis of amino acid changes in the products of known fungicide resistance loci (Mair et al., 2016) revealed several potential resistance adaptations, inferred from reported resistance associated with equivalent mutations in other fungal species (Supplementary Data Sheet 6; Oliver et al., 2024). The CYP51A protein of all isolates contained amino acid residues corresponding to P216L, M220I, and H147Y of Aspergillus fumigatus CYP51A [NCBI: AF338659], which may confer azole resistance (Howard et al., 2009). Some isolates had mutations in CYP51B corresponding to deletions at Q287 and G412, and a mutation at S208T (relative to Zymoseptoria tritici [NCBI: AY253234]), that may be involved in DMI resistance (Stammler and Semar, 2011), CytB loci were not well represented across this pan-genome with data for only two isolates, however for these no intron/intein mutations were detected, and one deletion mutation corresponded to site D203 (relative to Z. tritici [NCBI: AY247413]), which in Plasmopora viticola, was involved in resistance to cyazofamid (Mounkoro et al., 2019). Beta-tubulin exhibited mutations corresponding to E198A and M257L (relative to Aspergillus nidulans [NCBI: M17519]), which may confer benzamidazole resistance (Leroux et al., 2002). A single isolate had a mutation in the OS1 protein corresponding to A350S of Botrytis cinerea [NCBI: AF435964] which may confer resistance to Fludioxonil (Ren et al., 2016). Overall, the mutations detected above and their corresponding resistance phenotypes inferred from related fungal species may reflect the broad application of boscalid + pyraclostrobin, fluxapyroxad, pyrimethanil, dodine and copper that was applied to all isolates in this study.

Comparative analysis of orthologous groups across the pan-genome indicated 21,174 ortholog groups (Supplementary Data Sheet 7), with isolates encoding an average of 11.2 K genes each, and 6,805 (32%) core groups containing a single member that was present in all isolates (Supplementary Data Sheet 8). Candidate secreted effector-like protein (CSEP) prediction across the representative pan-genome orthologous gene set, which in an effort to highlight high-priority effector candidates, was filtered for: predicted secretion (Predector—any method), Predector Score ≥ 2, and cysteine residues ≥ 2, which resulted in 59 CSEP-orthogroups (Figure 3) (Supplementary Data Sheet 8). Of these, notable CSEPs with conserved functional domains and/or with homology to confirmed effectors of other pathogen species (Supplementary Data Sheets 9, 10) are summarized below according to their speculative roles at various phases of S. pistaciarum infection. Overall, the predicted CSEP set indicated that S. pistaciarum may employ a combination of defensive effectors with roles in suppression of host defenses, and offensive effectors with a range of cytotoxic activities. Notably, some effector-like ortholog groups may have presented as divergent versions of the same protein (Table 4; Supplementary Data Sheet 8), suggesting region-specific adaptations may have occurred. There were 23/59 highly-conserved CSEP orthogroups, common to all 27 isolates, with another 5 similarly conserved sets of groups (WSC, PBP, DLH, ZtNip2 and Cutinase, Table 4) presenting as separate orthogroups that may be highly divergent versions of the same ortholog.