Development of a PacBio Long-Read Sequencing Assay for High Throughput Detection of Fungicide Resistance in Zymoseptoria tritici

Fungicide resistance has become a challenging problem in management of Septoria tritici blotch (STB), caused by Zymoseptoria tritici, the most destructive disease of winter wheat throughout western and northern Europe. To ensure the continued effectiveness of those fungicides currently used, it is essential to monitor the development and spread of such resistance in field populations of the pathogen. Since resistance to the key families of fungicides used for STB control (demethyalation inhibitors or azoles, succinate dehydrogenase inhibitors or SDHIs and Quinone outside Inhibitors or QoIs) is conferred through target-site mutations, the potential exists to monitor resistance through the molecular detection of alterations in the target site genes. As more efficient fungicides were developed and applied, the pathogen has continuously adapted through accumulating multiple target-site alterations. In order to accurately monitor these changes in field populations, it is therefore becoming increasingly important to completely sequence the targeted genes. Here we report the development of a PacBio assay that facilitates the multiplex amplification and long-read sequencing of the target gene(s) for the azole (CYP51), SDHI (Sdh B, C, and D), and QoI (cytochrome b) fungicides. The assay was developed and optimised using three Irish Z. tritici collections established in spring 2017, which capture the range of fungicide resistance present in modern European populations of Z. tritici. The sequences obtained through the PacBio assay were validated using traditional Sanger sequencing and in vitro sensitivity screenings. To further exploit the long-read and high throughput potential of PacBio sequencing, an additional nine housekeeping genes (act, BTUB, cal, cyp, EF1, GAPDH, hsp80-1, PKC, TFC1) were sequenced and used to provide comprehensive Z. tritici strain genotyping.


INTRODUCTION
The application of fungicides for the control of foliar diseases has become an integral component of winter wheat production throughout Northern Europe. Amongst the targeted-diseases, septoria tritici blotch (STB) is the most commonly encountered and often the most economically destructive (Jørgensen et al., 2014;O'Driscoll et al., 2014). Unfortunately, the casual pathogen Zymoseptoria tritici has demonstrated a propensity to adapt and overcome the control provided by these fungicides. This is due to the combination of the high plasticity of the pathogen's genome and the increasingly specific nature of fungicides used for its control (Lucas et al., 2015;McDonald and Mundt, 2016). Over the past two decades, European Z. tritici populations have developed varying degrees of resistance to the majority of fungicides available for use in cereals, including those belonging to the Demethylation Inhibitors (DMIs) of which the azoles are the most abundant and relevant, the Quinone outside Inhibitors (QoIs) and Succinate Dehydrogenase Inhibitors (SDHIs) (Torriani et al., 2009;Leroux and Walker, 2011;Dooley et al., 2016a). The development and spread of such resistances have led to varying efficacy levels provided by the different fungicides (Torriani et al., 2009;Jørgensen et al., 2018Jørgensen et al., , 2021Joynt et al., 2019). For each fungicide group, multiple mechanisms, including target-site alterations, target-site overexpression and increased efflux activity, have been associated with resistance development (Leroux and Walker, 2011;Cools et al., 2012;Cools and Fraaije, 2013;Omrane et al., 2015;Kildea et al., 2019b).
Target-site alterations are the most common and widespread means through which resistance has developed in Z. tritici to the azole, QoI and SDHI fungicides. Nevertheless, in each instance, the nature and impact each have on their respective fungicide group can differ. In the case of QoIs, a single alteration in the pathogen's cytochrome b (G143A) confers complete resistance to all commercially available QoI fungicides (Torriani et al., 2009). On the other hand, resistance to azoles has developed through stepwise accumulation of alterations in multiple positions of their target-site 14α-demethylase (CYP51), with varying effects on the sensitivity to the different azoles (Leroux and Walker, 2011;Cools and Fraaije, 2013;Huf et al., 2018). For the SDHIs, multiple individual alterations of the succinate dehydrogenase (Sdh) subunits B, C, and D have been detected in field strains with contrasting impacts on SDHI sensitivity (Dooley et al., 2016a;Rehfus et al., 2018). Understanding the dynamics of how these alterations are selected in Z. tritici populations is a key component in developing strategies that limit their impacts, both from the aspect of control and further spread. Moreover, it is also highly relevant to be able to rapidly detect the emergence of new alterations or combinations of those already in populations.
Traditionally, the detection of fungicide resistance has relied upon in vitro sensitivity analysis. These assays involve measuring fungal growth on/in different media (solid or liquid) amended with increasing concentrations of the fungicide of interest and subsequently determining a given concentration required to inhibit growth by a desired amount (e.g., EC50 or effective concentration of a chemical required to reduce growth by 50%), which can then be used to compare between strains or populations. Although both laborious and economically cumbersome, such phenotypic screenings constitute the foundations of resistance monitoring. As increasing knowledge has become available on the development of fungicide resistance, so did the capacity to supplement these monitoring with detailed molecular studies that assess resistance based on changes at the genetic level. Molecular assays have been applied for the detection of fungicide resistance in Z. tritici from simple PCR assays, designed to amplify only in the presence of specific alterations and that can be readily applied in the most basic of laboratories (Siah et al., 2010) to more complex high throughput sequencing assays that generate millions of sequence reads per reaction and are often only available in specialised laboratories or institutes (Pieczul and Wąsowska, 2017). The strengths or weaknesses of each assay depend greatly on the study objectives and desired information.
In most instances, these assays are restricted in their capacity to detect changes in the DNA across short sequences of DNA. This hinders their ability to capture the full information on the haplotype of a strain, which might consist of a combination of multiple alterations which together confer increasing resistance in response to the commercialisation of increasingly active molecules (Cools and Fraaije, 2013). Full haplotype determination is increasingly required for resistance monitoring, as highlighted by the complexity of accumulation of CYP51 alterations (Leroux and Walker, 2011;Kildea et al., 2019b;Garnault et al., 2021). Whilst it is possible to augment these assays to capture multiple alterations (e.g., multiplexing of different dual-labelled probes to capture multiple alterations in a single qPCR assay as described by Hellin et al., 2020) or indeed entire genes greater than 1 kb (e.g., Frey et al., 2018), these are often either limited in their capacity to capture the diversity present or are highly complex and as such difficult to readily implement. The PacBio longread sequencing methods based on Single-Molecule Real-Time (SMRT) sequencing technology provides a potential opportunity to overcome both these limitations. SMRT sequencing is capable of generating continuous long-reads with an average read-length of 10 kb or more, and thus can be applied to directly examine the entire haplotypes of the fungicide resistance genes. In addition, an error correction method (circular consensus sequencing; CCS) has been established and repeated sequencing of a singlemolecule DNA template can result in extremely high accuracy (Wenger et al., 2019).
Here we report the development of a PacBio long-read sequencing assay built on multiplex amplification of target genes to enable the simultaneous detection of alterations across CYP51, Sdh B, C, and D subunits and cytochrome b, which, as described above, play a role in conferring resistance in Z. tritici to the azole, SDHI and QoI fungicides, respectively. The assay has been optimised to be completed in two multiplex PCRs to ensure high throughput, with an associated pipeline to identify non-synonymous mutations. It can handle 96 samples, with increased scalability possible.
To further exploit the high read accuracy associated with SMRT sequencing and the depth of sequencing reads generated in a single run, the multiplex amplification and sequencing of nine housekeeping genes (act, BTUB, cal, cyp, EF1, GAPDH, hsp80-1, PKC, TFC1), alongside those associated with fungicide resistance, have been included to provide comprehensive Z. tritici strain genotyping. The ability of the assay to be used for strain genotyping and population diversity analysis was validated using a panel of Z. tritici isolates collected from Irish fields in 2017, for which data from Sanger sequencing, PCR-RFLP, SSR, and microtitre plate assays were collected.

Fungal Isolates and DNA Extraction
In spring 2017, winter wheat leaves carrying STB lesions (final leaf 6) were sampled (a single infected leaf collected every 10 m following a W shape across each field) from commercial fields at three different locations in Ireland prior to fungicide applications ( Table 1). Fungal isolations (single pycnidia per infected leaf) were performed as described by Dooley et al. (2016a), after which single spore isolates were generated and subsequently stored at −80 • C in 30% glycerol (vv −1 ). DNA was extracted separately at both Teagasc and SLU from the 4-day old single spore isolates cultured on PDA using GeneElute TM (Sigma-Aldrich) and E.Z.N.A. SP Plant DNA Mini Kit (Omega Bio-Tek, Doraville, GA, United States), respectively, in accordance with the manufacturers' instructions. The collection has previously been described by Kildea et al. (2019a).

Development of the PacBio Sequencing Assay
The sequencing assay followed a four-step workflow (Figure 1): (i) amplification of the target genes in individual Z. tritici isolates using multiplex PCRs including indexing to allow pooling of isolates, (ii) PacBio long-read sequencing, (iii) de-multiplexing of sequencing reads and assigning consensus sequences to isolates, and (iv) alignment of sequencing reads to reference sequences and mutation/SNP scoring using Geneious Prime software.

Multiplex Amplification of Fungicide Resistance and Housekeeping Genes
Primers used for the amplification of the azole (CYP51), SDHI (Sdh subunits B, C, and D), and QoI (cytochrome b) target genes in Z. tritici and the nine housekeeping genes obtained from relevant literature or designed in the present study ( Table 2). The lengths of amplicons generated for the housekeeping genes were 594-750 bp and 528-2,192 bp for the fungicide target genes. Universal heel sequences, allowing the addition of specific tags (indexing primers), were added to each forward (ctctctatgggcagtc) and reverse primer (ctcgtgtctccgact; Table 2). PCR conditions used for amplification were modified from Chaudhary et al. (2020), originally based on a protocol from Nguyen-Dumont et al. (2013). Due to the varying amplicon sizes, two separate multiplex PCRs were used; PCR "A" for long amplicons (1,200-2,200 bp) and PCR "B" for short amplicons (500-1,000 bp), with extension times and primer concentrations, adapted to amplify all DNA fragments in satisfactory amounts.
Amplification conditions for PCR "A" (long amplicons): The PCR mixture (10 µl) consisted of 1X Phusion HF PCR buffer, 8.3 nM of each CYP51 primer, 6.7 nM of each SdhB primer, 0.2 mM of each dNTP, 0.06 U µl −1 Phusion Hot Start II High-Fidelity DNA Polymerase and 5 ng of genomic template DNA. PCR was performed using a 2720 Thermal cycler (AB Applied Biosystems, Life Technologies Corporation, Carlsbad, CA, United States). An initial denaturation step at 98 • C for 1 min was followed by 10 cycles of denaturation at 98 • C for 30 s, annealing at 61 • C for 45 s and extension at 72 • C for 2 min 30 s, followed by another 15 cycles of denaturation at 98 • C for 30 s, annealing at 59 • C for 30 s and extension at 72 • C for 2 min 30 s, and thereafter a final extension step at 72 • C for 5 min. A second PCR reaction was performed directly afterwards in the same reaction tubes, after adding unique forward and reverse tags (2 µM each) (Supplementary Table 1) to each of the 96 reaction tubes. Thermal cycling conditions consisted of an initial denaturation step at 98 • C for 1 min was followed by 3 cycles of denaturation at 98 • C for 30 s, annealing at 50 • C for 30 s and extension at 72 • C for 2 min 30 s, followed by another 10 cycles of denaturation at 98 • C for 30 s, annealing at 57 • C for 30 s and extension at 72 • C for 2 min 30 s, and thereafter a final extension step at 72 • C for 10 min.
All PCR products from the 96 wells were pooled (separately for PCR "A" and "B") and then purified with AMPure XP beads (Beckman Coulter Genomics, Danvers, MA, United States) using a ratio of 1:0.6 (PCR product:ampure solution) to remove DNA fragments shorter than 300 bp. The PCR products were then further purified with E.Z.N.A. Cycle Pure Kit (Omega Bio-Tek, Doraville, GA, United States). The DNA quality and amplicon sizes were analysed using an Agilent 2100 Bioanalyser (Agilent Technologies, Santa Clara, CA, United States). PCR products of long amplicons (PCR "A") and short amplicons (PCR "B") were pooled in equal volumes before PacBio sequencing.

PacBio Sequencing
Sequencing was performed at SciLifeLab in Uppsala (NGI, Sweden) with the PacBio Sequel platform (Pacific Biosciences, Menlo Park, CA, United States), which utilises the SMRT technology and circular consensus sequencing (CCS) to produce long reads (Eid et al., 2009). The pooled PCR products from the two multiplex PCRs (containing 14 different amplicons from 96 isolates) were sequenced in one SMRT cell. a Number of isolates included in each assay (out of those included in the PacBio sequencing). b PacBio sequencing of all 14 target genes. c Sanger sequencing of the CYP51, and SdhB, C, and D subunits. The number of isolates sequenced varies for the different target genes. d PCR-RFLP used for the detection of the cytochrome b substitution G143A as reported by Kildea et al. (2019a). e Sensitivity as determined using in vitro microtitre plate assays as described by Dooley et al. (2016b). The number of isolates tested varies for the different fungicides (epoxiconazole, metconazole, prothioconazole-desthio, tebuconazole, benzovindiflupyr, bixafen, boscalid, fluxapyroxad, isopyrazam and penthiopyrad).
FIGURE 1 | Illustration of the multiplex PCR principal and workflow applied in the PacBio long-read sequencing assay development. In the first multiplex PCR reaction, a number of gene-specific primers (gsp) with heel sequences are used. In the second multiplex PCR reaction, unique tags (one for each reaction tube) are added and used as indexing primers. For simplicity, only one multiplex PCR is shown in the figure, however, amplifications were divided into two PCRs ("A" and "B") in the current study.

De-Multiplexing
The multiplexed PacBio sequences were de-multiplexed using an in-house developed script ("detag.py") 1 implemented in Python 2.7 2 , previously described by Chaudhary et al. (2020), to sort reads based on their sequence tag and assign them to their sample of origin. Reads of too low quality (described in detail below), failing to match the heel or tag sequencing, or chimeric reads, were discarded. The script was run with the following arguments: unpairedfastq (unpaired fastq files); h5score 0.9 (proportion of matching 5 heel bases for acceptance); h3score 0.9 (proportion of matching 3 heel bases for acceptance); filtering 0:400:2200:30:20 (full filtering, done before the program searches for tags in the sequences : minimum length : maximum length : minimum average per-base phred-score a read needs : minimum phred-score a base needs). Input files with tags and corresponding sample ID (tags file), heel sequences (heel file), and primer sequences (gsp file; for statistics) were supplied to the script. In the processing, reads were organised into new fastq files (labelled with sample ID) with heel and tag sequences removed.

Alignments of Sequence Reads
The bioinformatic software Geneious Prime 2019.0.4 3 was used to align sequencing reads and identify single nucleotide polymorphisms (SNPs) and mutations. The output fastq files from the de-multiplexing script were used as input in Geneious Prime. Alignments were performed with the function "Map to reference." First, all reads were aligned separately, within sequence lists (i.e., within the sample; option "Assemble each sequence list separately"), to reference sequences of the 14 target genes (Supplementary Table 3). All consensus files were then aligned to the 14 reference sequences, resulting in a list for each target gene, with aligned consensus sequences from all samples. For genes with several exons (CYP51, SdhB, SdhC, and SdhD) the nucleotide sequences were aligned to each exon's reference nucleotide sequences. To find nonsynoymous mutations in the fungicide resistance genes, the consensus nucleotide sequences were translated to amino acid sequences using the function "Translate." All variable amino acid positions throughout the genes were identified by visual examination. The data were copied to an Excel spreadsheet and subsequently compared to those identified through Sanger sequencing. The nucleotide error rate of PacBio sequencing reads were estimated after quality filtering by the de-multiplexing script. For each targeted gene, the nucleotide error rate was calculated among all reads for ten randomly selected isolates, using the function "Find Variations/SNPs." The number of nucleotides examined ranged from 44,000 to 1,009,000 per gene, and in total 8,932,000 nucleotides.
The nucleotide consensus reads from the alignment in Geneious Prime were used to locate SNP markers in the housekeeping genes. Variable nucleotide positions in the Sanger Sequencing of CYP51, SdhB, C, D, and Cytochrome b Genes Amplification, Sanger sequencing and analysis of the CYP51, SdhB, C, and D subunits of the 96 isolates of the three collections were performed as previously described by Dooley et al. (2016b) and Kildea et al. (2019a), respectively. Due to the low amplification of the SdhC subunit in some isolates, the sequencing assay described by Rehfus et al. (2018) was also used. Sequencing was performed by LGC Genomics (Berlin, Germany). Based on the combination of alterations in CYP51 gene, isolates were assigned to a CYP51 haplotype following the nomenclature of Huf et al. (2018), with those not previously reported assigned a new haplotype accordingly. The presence of amino acid substitution G143A in the cytochrome b, known to confer QoI resistance, was determined using a PCR-RFLP assay as previously reported by Kildea et al. (2019a). Sanger re-sequencing of ambiguous samples were performed by Macrogen Europe (Amsterdam, Netherlands) using primers pairs CYP51-F1/CYP51-R1, CYP51-F3/CYP51-R3, and CYP51-F4/CYP51-R4 for gene CYP51 (Estep et al., 2015), and primer pairs according to Table 2 for the SdhB, SdhC, and SdhD genes.

Fungicide Sensitivity Assessment
All in vitro fungicide sensitivity assessments were determined using a microtitre plate assay described by Dooley et al. (2016b), with minor modifications depending on the fungicide. In total, the sensitivities of all 96 strains were determined to four azole fungicides (epoxiconazole, metconazole, prothioconazole-desthio, and tebuconazole) and six SDHI fungicides (benzovindiflupyr, bixafen, boscalid, fluxapyroxad, isopyrazam, and penthiopyrad). In all assays, technical grade fungicides were used and purchased from Sigma-Aldrich (St. Louis, MO, United States), except for benzovindiflupyr (kindly provided by Syngenta). Fungicides were initially dissolved in DMSO and subsequently aliquot to their test concentrations in potato dextrose broth. Fungicides were adjusted to test concentrations of 0, 0.04, 0.123, 0.33, 1.1, 3.3, 10, and 30 mg l −1 , with the exception of prothioconazole-desthio which concentration was adjusted to give test concentrations of 0, 0.01, 0.04, 0.123, 0.33, 1.1, 3.3, 10 mg l −1 . Microtitre plates were incubated in the darkness at 18 • C for 7 days, after which Z. tritici growth was measured as light absorbance using a Synergy-HT plate reader and Gen5 microplate software (BioTek Instruments). Sensitivity was expressed as the concentration of fungicide reducing fungal growth by 50% (EC 50 ), which was calculated by fitting a logistic curve to the absorbance data generated for each individual isolate using XL FIT (IDBS Inc.). The sensitivity of the collections to the QoI fungicide pyraclostrobin has been previously reported by Kildea et al. (2019a).
Analysis of variance was performed using Tukey's studentised range test to determine differences between the three locations regarding the sensitivity of their Z. tritici population to each of the 10 fungicides. The analysis was conducted in R 4.0.2 using the packages "car, " "agricolae, " and "ggplot2."

Population Genetic Analysis
The population structure of the Z. tritici samples was individually determined based on both SNP marker data from the PacBio sequenced housekeeping genes and microsatellite data. To calculate the level of differentiation among populations, analysis of molecular variance (AMOVA) was done with the Microsoft EXCEL add-in GenAlEx version 6.503 Smouse, 2006, 2012).

PacBio Sequencing
A single PacBio sequencing run resulted in a total of 375,804 reads with Phred quality score above Q20 (99% base call accuracy), covering a total of 396,117 Mbp. After processing and checking the quality of the reads using the de-multiplexing script, 88,481 reads (23.5%) were accepted. Reads not accepted were discarded due to low minimum quality (41%), missing primer(s) (17%), missing tag(s) (17.2%), chimeric tags (1.3%), too short length (0.2%), and truncation (<0.1%). Across all isolates, the average number of reads per targeted gene was between 5.6 and 171.6, with a range between 0 and 556 reads (Table 3). A high overall similarity was found among nucleotides in repeated reads of the same gene for the ten randomly selected isolates, with an overall error rate of 2.21 × 10 −4 , ranging from 0.95 × 10 −4 to 5.02 × 10 −4 for the different genes, i.e., 1-5 errors per 10,000 nucleotides ( Table 3).

Cytochrome b Sequence and Sensitivity to QoI Fungicides
Using PacBio sequencing, a 521 bp fragment encoding amino acids 96-268 of the cytochrome b was successfully amplified and sequenced with a mean read depth of 171.6 reads in each of the 96 isolates ( Table 3). The alteration G143A was detected in 92 of the 96 isolates sequenced, which is identical to the detection obtained with the PCR-RFLP assay. No further alterations related to QoI resistance were identified in these isolates or in those without G143A. The sensitivity of 85 of these isolates to the QoI pyraclostrobin was determined with the four isolates with G143 deemed sensitive (EC50 < 0.04 mg l −1 ) and the remaining 81 isolates with G143A ranging from 0.37 to 4.30 mg l −1 as previously described by Kildea et al. (2019a).

Succinate Dehydrogenase Subunit B, C, and D Sequences and Sensitivity to the SDHI Fungicides
A 1,217 bp sequence covering the entire SdhB subunit was successfully obtained for each of the 96 isolates, with a mean read depth of 418.1 reads/isolate ( Table 3). Only three non-synonymous mutations were detected amongst the isolates, with the alterations B-N225T identified in two isolates (M16 and C21), and B-R265L and B-R265Q in a single isolate each (M6 and KB31, respectively) (Figure 2A). All three alterations were also identified using Sanger sequencing (a total of 79 isolates compared with both PacBio and Sanger sequencing). Sensitivity for one of the isolates with B-N225T was obtained and resulted in a 3.0-9.8-fold reduction in its sensitivity to the six SDHIs tested compared to those without any Sdh subunit alteration. The isolate with B-R265L exhibited a 4.7-19.8-fold reduction in sensitivity (Figure 3). The sensitivity of the isolate with B-R265Q was not determined.
For the SdhC subunit, a 903 bp sequence was successfully sequenced in 94 of the 96 isolates tested with a mean read depth of 7.6 reads/isolate. For two isolates (M21 and KB31) no reads were produced. Only a single read was produced for a further seven isolates (M37, C3, C5, C14, C27, KB19, and KB36) and as such, these were not included in the analysis. Amongst the remaining 87 isolates, 10 non-synonymous mutations were identified (Figure 2A). All mutations were confirmed using Sanger sequencing with exception of C-R22H (isolate KB39), for which Sanger sequencing was not performed (a total of 68 isolates compared with both PacBio and Sanger sequencing). As anticipated, only specific alterations, C-T79N, C-N86S, C-V150L & C-V166M, and C-W155R, impacted sensitivity to the different SDHIs, with 3.3-10.5-fold reductions in sensitivity depending on the alteration and SDHI fungicide (Figure 3).
A 857 bp sequence capturing the entire SdhD subunit was successfully amplified and sequenced in 92 of the 96 isolates, with a mean read depth of 12.4 reads/isolate. For four isolates (C16, C33, KB5, and M21) no reads were produced. Amongst the 92 isolates successfully sequenced, two nonsynonymous mutations were identified (Figure 2A), and their presence was confirmed by Sanger sequencing (a total of 75 isolates compared with both PacBio and Sanger sequencing). In two of the isolates (M23 and M24), this led to the alteration D-T45A, which conferred a 5.2-12.9-fold reduction in sensitivity to the different SDHIs (Figure 3). The alteration D-I95T, found in a single isolate (KB1), did not result in any sensitivity change.
Differences in the frequencies of the various alterations were detected amongst the three sites, with C-T79N most prevalent in Knockbeg and Meath, whilst absent in Cork ( Figure 2B). These differences reflected the differing SDHI sensitivities between the three sites (Supplementary Image 1).

CYP51 Sequence and Sensitivity to the Azole Fungicides
A 2,194 bp sequence capturing the entire CYP51 gene was successfully amplified and sequenced in 92 of the 96 isolates, with a mean read depth of 42 reads/isolate ( Table 3). No reads were produced for three isolates (M32, C3, and C5), and the only read for one isolate (M17) was lost due to failure to map to the reference CYP51 gene. For an additional two isolates (M30 and C26), only a single read was produced, and whilst each mapped to the reference, they were not included in further analysis. Comparison of the translated sequences generated using PacBio sequencing with those generated using Sanger (n = 69) showed that the amino acid sequences were identical.
Amongst the 90 isolates successfully sequenced using PacBio, 11 non-synonymous mutations and a 6-bp deletions were detected, resulting in 22 different CYP51 amino acid haplotypes ( Figure 2C). All amino acid alterations detected have previously been reported in Z. tritici (Cools and Fraaije, 2013). The frequencies of both CYP51 amino acid alterations and haplotypes differed between the three sites (Figures 2B,C), reflected in differing levels of sensitivity between the sites to the four azoles tested (Supplementary Image 2). The sensitivity profile of the isolates to the different azoles varied depending on the specific alteration present at CYP51 amino acid positions 136 (Figure 4A), 381 (Figure 4B), and 524 ( Figure 4C). The sensitivity profiles to combinations of alterations at the three amino acid positions are shown in Figure 4D.

Population Structure Based on Housekeeping Gene Markers and SSR Markers
In total 105 polymorphic SNP markers were identified in the nine housekeeping genes, with 3-25 markers per gene (Supplementary Table 2). The housekeeping gene markers revealed 92 unique multilocus genotypes (MLGs) among the 96 isolates, with 4 MLG's found twice.
The 13 SSR markers could identify 1-9 alleles per marker, and totally 56 alleles in the entire dataset. The SSR markers distinguished 45 unique MLGs, of which 3 MLGs were found twice, 4 three times, 2 four times, 1 nine times, and 1 was observed fifteen times.
The AMOVA showed that 100% of the variation for the SNP marker data (housekeeping genes) and 97% of the variation for the SSR marker data were found within populations, and no significant differentiation among the three locations was indicated ( Table 4).

DISCUSSION
With increasing restrictions surrounding the registration and usage of pesticides in the European Union, it is essential to ensure those fungicides that continue to be available and those that are approved into the future remain effective (Jess et al., 2014). Understanding the dynamics of fungicides resistance in target pathogen populations is therefore a key pillar of integrated pest management programmes to control the ensuing diseases. The development of the long-read PacBio assay and associated downstream analysis for detecting target-site fungicide in Z. tritici, as detailed here, provide tools to aid such essential research.
To validate the assay, Z. tritici populations were collected from three different regions in Ireland collected in spring 2017. They represent a period when the Irish Z. tritici population was undergoing considerable changes in terms of fungicide resistance (Dooley et al., 2016a;Kildea et al., 2019b) and, therefore, ideal for testing the robustness of the assay. Based on the phenotypic sensitivity assays varying levels of resistance were present in the collections. Sanger sequencing and a PCR-RFLP assay confirmed the presence of multiple alterations previously identified as conferring azole, SDHI and QoI resistance in their respective target-sites. In all instances, the PacBio assay detected these alterations.
The housekeeping gene markers developed in this study showed high variability and were able to distinguish among multilocus genotypes, even to a greater extent than the SSR markers. The high genotypic diversity indicates that the Z. tritici populations are highly recombined as a result of regular cycles of sexual recombination (McDonald and Mundt, 2016). Neither the SSR markers nor the housekeeping gene markers revealed any population structure, i.e., no differentiation between populations and no subgrouping within populations, and the majority of genetic diversity was identified within the fields. The conformity of the two marker sets demonstrates that the housekeeping gene markers, like SSR markers, are neutral and unlinked to the genes involved in fungicide resistance. The lack of population structure among the three Irish fields agree with previous studies, where little geographic differentiation among Z. tritici populations was found within Ireland (Welch et al., 2018;Kildea et al., 2019a) as well as on a wider geographic scale (Zhan et al., 2003;McDonald and Mundt, 2016) and is suggested to be explained by high gene flow between populations combined with regular sexual recombination. The differences observed in patterns of fungicide sensitivity and haplotype distribution of fungicide resistance genes among locations, despite the high similarity for neutral markers, suggest that the rate of gene flow may not be sufficient to counteract gene frequency changes caused by ongoing local selection for fungicide resistance. The differences in sensitivity identified between the sites may reflect differences in local fungicide usage. However, it should be noted that the populations examined in this study were collected prior to fungicide application. Even though subtle differences in fungicide usage may occur between the geographic regions from which the samples were obtained, Irish winter wheat fungicide programmes have been broadly similar, relying on combinations of SDHIs and azoles, applied at final leaf three emergence and again at flag leaf emergence. Furthermore, in all instances, the populations were collected from fields that were previously planted to a non-cereal which suggests initial infections would primarily have resulted from airborne ascospore inoculum. As identified by Hellin et al. (2020) Z. tritici ascospores are abundant in Ireland in late autumn and early winter and the frequencies of alterations conferring fungicide resistance in these can vary  over this period. Further research is required to dissect these differences in fungicide sensitivity.
In the initial design of the assay, a single multiplex PCR assay was trialed; however, due to varying amplicon sizes, larger amplicons were under-represented once sequenced. As a compromise, two initial PCRs were performed with larger fragments (CYP51 and SdhB) amplified together. Nevertheless, as the SdhB amplicon is about half the size of the CYP51, its average read depth following sequencing was almost 10× greater than for CYP51 (Table 3). Similar differences were observed for the other amplicons highlighting that further optimisation of the assay should be possible, especially where individual specific targeted amplicon sequencing is required (e.g., CYP51 only).
The amplification efficiency of individual targeted genes in a multiplex PCR reaction may be impaired in an unpredictable way because of interactions and unequal competition among amplicons for reaction components. For example, shorter amplicons will amplify at a higher rate than longer ones; some primers might be constrained because of primer dimer formations, amplification of certain templates may be favoured due to the properties of the target, the target's flanking sequences, or differential accessibility of targets within the genome due to secondary structures (Elnifro et al., 2000). In the current assay, individual primer concentrations were adjusted based on empirical testing to promote disadvantaged amplicons. An increase of primer concentrations must, however, be done with caution. The approach of the multiplex PCRs is to keep the concentrations of primers low in the first PCR, ensuring they are depleted before the second PCR, where the added indexing tags will act as primers. If the primers from the first PCR still remain at this stage, they may produce new fragments without tags. Seventeen percent of the reads generated in the present study were rejected in the de-multiplexing step due to lack of tags, which may partly be a result of incomplete depletion of primers in the first PCR.
Depending upon user requirements, the number of targets can be restricted to those relating to fungicide resistance or other specific traits such as those used for genotyping. By reducing the number of targets, the read depth for those remaining will increase. This can be used either to increase the overall number of isolates included or to provide greater depth on those individuals assessed. As the assay is also designed with an initial multiplex PCR, reducing the number of targets in this PCR will also increase the read number of those amplicons under-represented due to the initial PCR process. This potential to increase read depth also opens the assay to be used as a possible quantitative detection method. Phelan et al. (2017) have previously demonstrated the potential of high throughput sequencing to be applied to fungicide resistance detection, with 454-sequencing used to detect the cytochrome b alteration G143A in the Irish field populations of the barley pathogen Rhyncosporium commune. As identified in the present study, multiple alterations across the entire CYP51 gene impact the sensitivity to the azole fungicides in a cumulative manner, and similar impacts from multiple Sdh subunit mutations were recently reported (FRAC, 2021). Being able to quantify the specific CYP51 or Sdh haplotypes directly from field samples will be invaluable in understanding the dynamics of fungicide resistance following varying treatments.
The fungicide resistant component of the assay described here is specific to target site fungicide resistance. Undoubtedly, this is the predominant form of resistance contributing to erosion in fungicide sensitivity of European Z. tritici populations and associated decrease in efficacy observed in the field (Jørgensen et al., 2021). Resistance mechanisms other than target-site alterations, such as target-site overexpression and/or overexpression of efflux pumps, have also been reported in European populations of Z. tritici and have the potential to contribute to resistance development (Cools et al., 2012;Omrane et al., 2017). In the case of CYP51 overexpression conferred by inserts in the promoter region, it should be possible to increase the size of the targeted CYP51 amplicon to capture some of these key inserts (e.g., 120 bp as described by Cools et al., 2012). Given the increased size of this amplicon, it might be required to amplify it separately as previously highlighted to ensure  its adequate representation in the final sequencing. Additional CYP51 promoter inserts, some of which potentially impacting CYP51 expressions, have also been identified (Kildea et al., 2019b;Huf et al., 2020). However, as some of these inserts can range in size up to 866 bp, careful amplicon design will be required to ensure the region is effectively captured and sequenced. Given the likely varying fragment sizes that would be amplified in a mixed DNA sample and the resulting bias that may occur, the potential applicability of the assay would be restricted to qualitative detection. This would also be relevant for the detection of inserts in the promoter region of the efflux pump Mfs1, which have been associated with increased efflux of a diverse range of fungicides (Omrane et al., 2017). Although initial setup costs for the PacBio sequencing assay (excluding the costs of the sequencing machines) are slightly more expensive due to the additional PCR step and adapters and indexes, when utilised to its full capacity to screen multiple targets in large collections of isolates, costs per strain can be expected to be considerably lower compared to traditional Sanger sequencing. Future research is ongoing to apply the assay to contemporary collections of Z. tritici collected throughout Europe to provide a comprehensive overview of the frequency and spread of fungicide resistance throughout the region.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA720680.