A fungal avirulence factor encoded in a highly plastic genomic region triggers partial resistance to septoria tritici blotch.

Cultivar-strain specificity in the wheat-Zymoseptoria tritici pathosystem determines the infection outcome and is controlled by resistance genes on the host side, many of which have been identified. On the pathogen side, however, the molecular determinants of specificity remain largely unknown. We used genetic mapping, targeted gene disruption and allele swapping to characterise the recognition of the new avirulence factor Avr3D1. We then combined population genetic and comparative genomic analyses to characterise the evolutionary trajectory of Avr3D1. Avr3D1 is specifically recognised by wheat cultivars harbouring the Stb7 resistance gene, triggering a strong defence response without preventing pathogen infection and reproduction. Avr3D1 resides in a cluster of putative effector genes located in a genome region populated by independent transposable element insertions. The gene was present in all 132 investigated strains and is highly polymorphic, with 30 different protein variants identified. We demonstrated that specific amino acid substitutions in Avr3D1 led to evasion of recognition. These results demonstrate that quantitative resistance and gene-for-gene interactions are not mutually exclusive. Localising avirulence genes in highly plastic genomic regions probably facilitates accelerated evolution that enables escape from recognition by resistance proteins.

Introduction

Regardless of whether they are mutualistic or parasitic, colonising microbes evolve a high degree of specialisation to recognise and infect their hosts and overcome host‐inducible defences (van der Does & Rep, 2017). Host manipulation is frequently achieved by the secretion of effectors, which are often small secreted proteins (SSPs) that support growth and development of the microbe by conferring protection against host antimicrobial compounds or by altering host metabolism (Lo Presti et al., 2015). Although effectors are beneficial for host colonisation, some are specifically recognised by certain host genotypes, triggering an immune response (Jones & Dangl, 2006; Lo Presti et al., 2015). This interaction typically follows the gene‐for‐gene model, in which a so‐called resistance protein recognises an effector, which is then called an avirulence factor (Avr) (Flor, 1971; Jones & Dangl, 2006). A common assumption is that resistance/avirulence gene interactions confer complete resistance, whereas quantitative resistance, understood here as incomplete or partial resistance that allows some pathogen infection and reproduction, is based on different, race‐nonspecific and therefore avirulence‐independent mechanisms. This paradigm originated from work on biotrophic pathogens, where avirulence recognition often leads to complete immunity via induction of a hypersensitive response (Cook et al., 2015; Niks et al., 2015). However, it is often overlooked that gene‐for‐gene interactions could also lead to quantitative resistance, as suggested by several studies (Antonovics et al., 2011; Rietman et al., 2012; Chen et al., 2013). Recently, more refined concepts such as the ‘invasion model’ (Cook et al., 2015) or ‘effector‐triggered defence’ (Stotz et al., 2014) emphasised a broader perspective for the gene‐for‐gene model, in which resistance gene‐based effector recognition and quantitative resistance are not mutually exclusive (Niks et al., 2015). However, avirulence factors leading to quantitative resistance have only rarely been described (Schirawski et al., 2010; Rietman et al., 2012).

Host recognition of effectors exerts an evolutionary pressure that favours sequence modification, deletion or acquisition of new effectors to overcome the immune response. Thus, genes encoding effectors are among the most polymorphic found in pathogen genomes (Win et al., 2012). The mechanisms underlying effector diversification remain largely unexplored. Many pathogen genomes are compartmentalised into highly conserved or rapidly evolving regions, often described as the ‘two‐speed genome’ (Raffaele & Kamoun, 2012). Effector genes are frequently localised in the highly variable compartments, which are often rich in transposable elements (Ma et al., 2010; Soyer et al., 2014; Plissonneau et al., 2018). Transposable elements are thought to contribute to genome evolution and the diversification of effector genes (Raffaele & Kamoun, 2012). They translocate within a genome, causing gene disruption, duplication or deletion of genomic sequences. In addition, transposable elements contribute to variability by favouring nonhomologous recombination or through repeat‐induced point mutations (RIPs) (Möller & Stukenbrock, 2017; Seidl & Thomma, 2017). Pathogens carrying these highly plastic genome regions are thought to benefit from an increased versatility to adapt to different conditions or to an evolving host (Dong et al., 2015; Faino et al., 2016).

The most damaging pathogen of wheat in Europe is Zymoseptoria tritici, an ascomycete fungus that causes septoria tritici blotch (STB; Fones & Gurr, 2015). It is an apoplastic pathogen with a latent necrotrophic lifestyle (Sánchez‐Vallet et al., 2015). Fungal hyphae penetrate the stomata and colonise the apoplast during a long asymptomatic phase that lasts between 7 and 14 d, depending on the weather conditions, the host cultivar and the pathogen strain. This long latent period is followed by a rapid induction of necrosis that is accompanied by the development of asexual reproductive structures called pycnidia, which contain asexual spores that spread the disease during a growing season (Kema et al., 1996; Duncan & Howard, 2000). The genetic basis of Z. tritici virulence is poorly understood as a result of its largely quantitative nature (Hartmann et al., 2017; Stewart et al., 2018). Two highly conserved lysin motif (LysM) effectors, Mg1LysM and Mg3LysM, prevent fungal recognition and shield the fungal cell wall from degradation by host hydrolytic enzymes (Marshall et al., 2011). The other known effectors of Z. tritici, Zt80707, AvrStb6 and Zt_8_609, are rapidly evolving small secreted proteins (Poppe et al., 2015; Hartmann et al., 2017; Zhong et al., 2017; Kema et al., 2018). The latter two were identified because they are specifically recognised by certain wheat cultivars, and they were found to be located in transposable element‐rich genomic regions (Brading et al., 2002; Hartmann et al., 2017; Zhong et al., 2017). AvrStb6 is recognised by the resistance protein Stb6 in a gene‐for‐gene interaction that leads to a strong resistance response, completely blocking the progression of the infection (Kema et al., 2000, 2018; Brading et al., 2002; Saintenac et al., 2018). In addition to Stb6, 19 other race‐specific Stb resistance genes with large effects have been mapped, but their corresponding avirulence factors remain unknown (Brown et al., 2015). We hypothesised that one of these Stb genes might be responsible for the differences in resistance of cultivar Runal to two Swiss strains (3D1 and 3D7, Stewart et al., 2018). The more virulent strain 3D7 produced necrotic lesions faster than strain 3D1. The less virulent strain 3D1 was successful in producing pycnidia, but at a lower density and with a less uniform distribution across the leaf surface than 3D7 (Fig. 1a,b). A single, large‐effect quantitative trait locus (QTL) encoding differences in lesion size and pycnidia density between 3D1 and 3D7 was mapped to a region on chromosome 7 (Stewart et al., 2018). However, the genes responsible for the differences in virulence were not identified.

Figure 1

Z09_7_00581 encodes the avirulence factor Avr3D1. (a) The second true leaves of the wheat cultivar Runal spray‐infected with Zymoseptoria tritici strains 3D1 and 3D7 and harvested at different days post‐inoculation (dpi). (b) Fully necrotised leaves with pycnidia of cv Runal infected by 3D1 and 3D7 at 21 dpi. The leaves shown in (a) are from a different experiment than the leaves shown in (b). Experiments were repeated twice with similar results. (c) Percentage of leaf area covered by lesions (PLACL) produced by the wild‐type (wt), the Avr3D1 knockout (KO, ∆avr3D1) and the ectopic mutants expressing the Avr3D1 allele of either 3D1 (Avr3D1 ) or 3D7 (Avr3D1 ) in the knockout background. Left panel: mutants in the 3D1 background at 15 dpi. Right panel: mutants in the 3D7 background at 14 dpi. Red dots represent the median of at least 15 leaves (except for the mock treatment, for which at least eight leaves were used), error bars represent 95% confidence intervals of the medians and black dots represent individual data points. Asterisks indicate statistical differences between wild‐type and knockout (P < 0.01, Kolmogorov–Smirnov test). (d) Amino acid sequence alignment of Avr3D1 variants of 3D1 and 3D7. The signal peptide sequence is highlighted in yellow and sequence polymorphisms between both alleles are shown in red.

Here we aimed to broaden our knowledge of the genetic basis of host‐race specificity in Z. tritici. First, we showed that Avr3D1 is the gene responsible for the differences in quantitative virulence between 3D1 and 3D7. We then demonstrated that Avr3D1 is an avirulence factor whose recognition is host‐specific, but triggers an incomplete, quantitative resistance. We next studied the evolutionary trajectory of Avr3D1 by combining population genetic and comparative genomic analyses involving 132 Z. tritici strains originating from four field populations on three continents as well as 11 strains of the closest known relatives of Z. tritici. We found that Avr3D1 is a member of an effector gene cluster that is located in a highly dynamic genomic region containing many independent insertions involving different families of transposable elements. Because an intact and presumably functional version of Avr3D1 was found in all strains of Z. tritici and in its closest relatives, we conclude that Avr3D1 plays an important role in the life history of Z. tritici. Maintaining Avr3D1 in a highly plastic genomic region probably provides an advantage by accelerating evolution that enables an escape from recognition in wheat populations carrying the corresponding resistance gene.

Materials and Methods

Quantitative trait locus mapping

To generate a genetic map, we used the previously generated restriction site associated DNA sequencing (RADseq) data from the progeny of the cross between 3D7 and 3D1 (Lendenmann et al., 2014). Quality trimmed reads were aligned to the genome of 3D7 (Plissonneau et al., 2016) using bowtie2 with default parameters (Langmead & Salzberg, 2012). Single nucleotide polymorphisms (SNPs) were called in each progeny with the HaplotypeCaller tool from Gatk v.3.3 (McKenna et al., 2010) and further filtered for their quality using the following parameters: > QUAL 5000, QD > 5, MQ > 20, and ReadPosRankSum, MQRankSum and BaseQRankSum between −2 and 2. We constructed the linkage map using R/qtl v.1.40‐8 (Arends et al., 2010). We retained only progenies for which 45% of all SNPs were genotyped, and we then removed SNPs genotyped in < 70% of the progenies. Potential clones (i.e. progenies with > 90% shared SNPs) were excluded. We removed adjacent nonrecombining markers. QTL mapping was performed with the qtl package in R (R Core Team, 2013), similar to the procedure described by Lendenmann et al. (2014) using the pycnidia density dataset (Stewart et al., 2018).

Zymoseptoria tritici and bacterial strains

The Swiss strains ST99CH_3D1 (3D1) and ST99CH_3D7 (3D7, described by Linde et al., 2002) or mutant lines derived from them were used in this study. Standard conditions for Z. tritici cultivation consisted of yeast‐sucrose broth (YSB) medium (10 g l−1 yeast extract, 10 g l−1 sucrose, 50 μg ml−1 kanamycin sulfate) at 18°C or yeast‐malt‐sucrose (YMS) medium (4 g l−1 yeast extract, 4 g l−1 malt extract, 4 g l−1 sucrose, 12 g l−1 agar) at 18°C. For molecular cloning and plasmid propagation, Escherichia coli strains HST08 (Takara Bio, Shiga, Japan) or NEB 5‐alpha (New England Biolabs, Ipswich, MA, USA) were used. Agrobacterium tumefaciens‐mediated transformation was performed with A. tumefaciens strain AGL1. If not stated otherwise, E. coli and Agrobacterium lines were grown in Luria Bertani (LB) medium containing kanamycin sulfate (50 μg ml−1) at 37°C or in LB medium containing kanamycin sulfate (50 μg ml−1), carbenicillin (100 μg ml−1) and rifampicin at 28°C (50 μg ml−1), respectively.

Generation of plasmid constructs for targeted gene disruption and ectopic gene integration

All PCRs for cloning procedures were performed using NEB Phusion polymerase (New England Biolabs) with primers listed in Supporting Information Table S1. All DNA assembly steps were conducted with the In‐Fusion HD Cloning Kit (Takara Bio) following the manufacturer's instructions. To create constructs for targeted gene disruption, two flanking regions of at least 1 kb in size for homologous recombination were amplified from Z. tritici genomic DNA. The hygromycin resistance gene cassette, used as a selectable marker, was amplified from pES6 (E. H. Stukenbrock, Department of Environmental Genomics, Christian‐Albrechts‐University of Kiel and Max Planck Institute for Evolutionary Biology, Plön, Germany, unpublished). The three fragments were assembled into pES1 (E. H. Stukenbrock, unpublished) linearized with KpnI and PstI (New England Biolabs), resulting in pES1∆5813D1 and pES1∆5813D7. To create the construct for ectopic integration of Avr3D1 , a fragment containing Avr3D1 including the 1.3‐kb sequence upstream of the start codon and the 1‐kb sequence downstream of the stop codon was amplified and cloned into pCGEN (Motteram et al., 2011) that had been linearised with KpnI, resulting in pCGEN‐5813D1ect. To exchange the coding sequence (CDS) in pCGEN‐5813D1ect, we first digested it with XhoI (New England Biolabs) to linearise it and remove the Avr3D1 CDS. In this digestion, the promoter was partially removed from the vector. In a second step a fragment containing the CDS and intron 1 of Avr3D1 (amplified from 3D7 genomic DNA) and a fragment to reconstitute the promoter sequence of Avr3D1 (amplified from pCGEN‐5813D1ect) were assembled into the linearised pCGEN‐5813D1ect, resulting in pCGEN‐5813D7ect. Constructs were transformed into E. coli by heat shock transformation, mini‐prepped and verified by diagnostic digests and Sanger sequencing (MicroSynth, Balgach, Switzerland). Confirmed plasmids were transformed into A. tumefaciens cells by electroporation.

Agrobacterium tumefaciens‐mediated transformation of Z. tritici cells

A. tumefaciens‐mediated transformation (ATMT) of Z. tritici was performed according to Zwiers & De Waard (2001) with the following modifications: A. tumefaciens lines were grown as liquid cultures for c. 24 h. Cell concentrations were estimated by measuring the optical density (OD600) and the cultures were diluted to an OD600 of 0.15 in induction medium (pH 5.7, 50 μg ml−1 kanamycin sulfate, 100 μg ml−1 carbenicillin, 50 μg ml−1 rifampicin, 10 mM glucose, 200 μM acetosyringone). These cultures were incubated at 28°C until they reached an OD600 of 0.25–0.35 and 100 μl was mixed with 100 μl of Z. tritici cell suspensions (cells grown on YMS for 4–6 d and washed off with water) and plated on induction medium covered with Amersham Hybond‐N+ nylon membranes (GE Healthcare Life Sciences, Chicago, IL, USA) membranes. After 3 d of incubation at 18°C, the nylon membranes were placed on YMS medium containing cefotaxime (200 μg ml−1) and either hygromycin B (100 μg ml−1) or geneticin (150 μg ml−1), depending on the resistance cassette of the construct, and incubated at 18°C until colonies appeared. Colonies were streak‐plated on the same selective medium to isolate single colonies before the mutant lines were grown on YMS without selection. For knockout lines, disruption of the target genes was verified using a PCR‐based approach. We determined the copy number of the transgene by quantitative PCR (qPCR) on genomic DNA extracted with DNeasy Plant Mini Kit (Qiagen). The target used was the selection marker and the reference gene was TFIIIC1 (Mycgr3G110539, Table S1). Only single insertion lines were selected for further experiments.

Infection assays

Seeds from wheat (Triticum aestivum L.) cultivars Runal, Titlis, Drifter, Chinese Spring and Arina were purchased from DSP Ltd. (Delley, Switzerland). Seeds were sown in peat substrate Jiffy GO PP7 (Jiffy Products International, Moerdijk, the Netherlands) and grown for 17 d in a glasshouse at 18°C (day) and 15°C (night) with a 16‐h photoperiod and 70% humidity. For all infection experiments, square pots (11 × 11 × 12 cm; Lamprecht‐Verpackungen GmbH, Göttingen, Germany ) containing 16–18 seedlings or 2 × 3 pot arrays (7 × 7 cm and 200 ml each, Bachmann Plantec AG, Hochdorf, Switzerland) containing two seedlings per unit were used. The infection procedure for the two pot types was identical. Z. tritici inoculum was prepared as follows: 50 ml of YSB medium was inoculated in 100‐ml Erlenmeyer flasks from Z. tritici glycerol stocks stored at −80°C. After 4–6 d of incubation (18°C, shaking at 120 rpm), liquid cultures were filtered through sterile cheesecloth and pelleted (3273 , 15 min, 4°C). The supernatant was discarded and the cells were resuspended in sterile deionised water and stored on ice until infection (0–2 d). The concentrations of the spore suspensions were determined using KOVA Glasstic counting chambers (Hycor Biomedical, Inc., Garden Grove, CA, USA) and adjusted to 106 spores ml−1 in 0.1% (v/v) Tween 20. Spore viability and concentration was analysed by performing a developmental assay on YMS medium as described below. Plants were sprayed until run‐off with 15 ml spore suspension per pot/array. Square pots were placed in plastic bags (PE‐LD, 380 × 240 mm) to support the leaves and stems. Subsequently, they were placed in a second plastic bag (PE‐LD 650 × 400 mm, two pots each), which was sealed to keep humidity at 100%. Pot arrays were placed directly into the sealing bags. After 3 d, the sealing bags were trimmed to a height of c. 27 cm and then opened, in the case of the 2 × 3 pot arrays, or completely removed in the case of the square pots, leaving the supporting bags intact in the latter case. For symptom quantification, second or third leaves were mounted on paper sheets, scanned with a flatbed scanner (CanoScan LiDE 220) and analysed using automated image analysis (Stewart et al., 2016). Data analysis and plotting was performed using RStudio v.1.0.143. Confidence intervals of the medians were determined using the ‘boot’ package and Kolmogorov–Smirnov (KS) tests for statistical significance with the ‘matching’ package in RStudio.

RNA isolation and quantitative reverse transcription PCR

Second leaves from cv Runal were infected with 3D1 or 3D7, harvested and scanned as described. Immediately after scanning, the top 2 cm of the leaves was excised and discarded and the adjacent 8.5 cm sections were frozen in N2. Three biological replicates were harvested. Leaf tissue was homogenised using a Bead Ruptor with a cooling unit (Omni International, Kennesaw, GA, USA) and zirconium oxide beads (1.4 mm). RNA was isolated using the GENEzol reagent (Geneaid Biotech, Taipei, Taiwan) and purified with the RNeasy Mini kit (Qiagen) including an on‐column DNase treatment with the RNase‐Free DNase Set (Qiagen) according to the manufacturer's instructions. cDNA was produced with the RevertAid First Strand cDNA Synthesis Kit (Invitrogen), using up to 900 ng RNA (estimated with NanoDrop) per reaction. To determine expression of Avr3D1 relative to the 18S reference gene, quantitative reverse transcription PCR (qRT‐PCR) was performed with a LightCycler 480 (Roche) using white 384‐well plates. Each reaction consisted of 250 nM of each primer, template cDNA generated from 11 to 30 ng of RNA and 1× HOT FIREPol EvaGreen qPCR Mix Plus mastermix (Solis BioDyne, Tartu, Estonia) in a total volume of 10 μl. Amplification was performed with a 10‐min step of initial denaturation and enzyme activation and 40 cycles of 95°C (15 s) and 60°C (60 s). Each sample was run in technical triplicates. Relative expression was calculated with LightCycler 480 software using the ‘advanced relative quantification’ tool. The mean and confidence interval of the mean was calculated with RStudio v.1.0.143.

Stress and development assay

The obtained Z. tritici mutant lines were tested for an altered, plant‐unrelated phenotype under various conditions including stress by growing them on potato dextrose agar (PDA), YMS and YMS supplemented with H2O2 (2 mM for 3D1 lines and 1 mM for 3D7 lines) or 1 M NaCl at 18°C. All media contained kanamycin sulfate (50 μg ml−1). An additional stress condition consisted of growth at 28°C on YMS. Inoculum preparation and quantification were the same as for the infection assays. Drops (2.5 μl) of spore suspensions of 107, 106, 105 and 104 spores ml−1 were plated on the media described above. Plates were assessed after 6 d of upside‐down incubation. Mutant lines exhibiting abnormal development or growth deficiencies were excluded from further experiments.

Manual annotation of three small secreted proteins in the QTL for virulence

We used RNA sequencing (RNA‐seq) raw data of IPO323 infecting wheat seedlings (Rudd et al., 2015) to manually annotate the gene Zt09_7_00581. To annotate the other genes in the cluster, we used RNA‐seq raw data of 3D7 from two different experiments and at six different time points (Palma‐Guerrero et al., 2016). The data were previously deposited in NCBI with experiment numbers SRP061444 and ERP009837. RNA‐seq reads were analysed as described by Hartmann & Croll (2017). Possible reading frames were manually examined using Integrative Genomics Viewer (IGV; Broad Institute; Robinson et al., 2011). Signal peptides were predicted using Signal P 4.1 (CBS; Petersen et al., 2011).

Zymoseptoria tritici strain collections

We used 132 strains collected in four different countries (Switzerland, Israel, USA and Australia; Zhan et al., 2005). Whole‐genome Illumina sequencing data of the 132 strains were previously deposited on the NCBI Short Read Archive under BioProject ID numbers PRJNA178194 and PRJNA327615 (Torriani et al., 2011; Croll et al., 2013; Hartmann & Croll, 2017; Hartmann et al., 2017). We used complete genome assemblies of IPO323, ST99CH_3D1 (3D1), ST99CH_3D7 (3D7), ST99CH_1E4 (1E4) and ST99CH_1A5 (1A5) previously described by Goodwin et al. (2011) and Plissonneau et al. (2016, 2018). BLAST searches were performed using the blastn command of the ncbi‐blast‐2.2.30+ software (Camacho et al., 2009). Synteny of the QTL between IPO323, 3D1, 3D7, 1E4 and 1A5 was analysed using blastn and visualised using the R package genoPlotR v.0.8.4 (Guy et al., 2010). Homologues of Avr3D1 were identified by blastn using CLC Genomic Workbench 9 (Qiagen) in strains 3D1, 3D7, 1E4 and 1A5.

We searched for orthologues of Avr3D1 using the blast algorithm implemented in CLC Genomics Workbench 9 (Qiagen) in one strain of Zymoseptoria passerinii (NCBI genome accession no. AFIY01 (fungal strain SP63)), four strains of Z. ardabiliae (STIR04 1.1.1 (accession no. AFIU01), STIR04 1.1.2 (AFIV01), STIR04 3.13.1 (AFIW01), STIR04 3.3.2 (AFIX01)), one strain of Z. brevis Zb18110 (LAFY01) and five strains of the sister species Z. pseudotritici (STIR04 2.2.1 (AFIQ01), STIR04 3.11.1 (AFIO01), STIR04 4.3.1 (AFIR01), STIR04 5.3 (AFIS01), STIR04 5.9.1 (AFIT01)). The genomes were downloaded from NCBI under accession numbers PRJNA63035, PRJNA277173, PRJNA63037, PRJNA63039, PRJNA343335, PRJNA343334, PRJNA343333, PRJNA343332, PRJNA63049, PRJNA273516 and PRJNA46489.

Presence/absence polymorphism of transposable elements and annotation

Repetitive DNA was identified for the 132 strains. For 3D1, 3D7, 1E4 and 1A5, full genome annotations were already available (Plissonneau et al., 2016, 2018). We annotated and masked repetitive elements in the genomes of the remaining 128 strains using RepeatModeler v.1.0.8, as described before (Plissonneau et al., 2016) and we masked the genomes using RepeatMasker v.4.0.5 with the library previously obtained for Z. tritici strain IPO323 (Grandaubert et al., 2015) according to transposable element nomenclature defined by Wicker et al. (2007).

DNA and protein alignments and phylogenetic tree

DNA and protein sequence alignments of Avr3D1 and the other SSPs from different strains were obtained using CLC Genomics Workbench 9 (Qiagen). For the phylogenetic analysis, amino acid sequences of Avr3D1 were aligned using Muscle. The maximum likelihood phylogeny reconstruction was performed applying the WAG model, with the software Mega6 (Tamura et al., 2013).

Population genetic analysis

DnaSP v.5 (Librado & Rozas, 2009) was used to calculate summary statistics of population genetic parameters associated with Avr3D1. Sliding window analyses of π were conducted using DnaSP with a window length set to 20 bp and a step size of 5 bp. The haplotype alignment of the coding region was used to generate a parsimony haplotype network using the TCS method (Clement et al., 2000) as implemented in the PopArt package v.1.7 (Leigh & Bryant, 2015). TCS utilises statistical parsimony methods to infer unrooted cladograms based on Templeton's 95% parsimony connection limit. Mutational steps resulting in nonsynonymous changes were identified with DnaSP.

The degree of selection was estimated by comparing dN (the number of nonsynonymous changes per nonsynonymous site) with dS (the number of synonymous changes per synonymous site) for all pairwise sequence comparisons using DnaSP. A dN : dS ratio of 1 (ω = 1) indicates neutrality, while ω < 1 suggests purifying, and ω > 1 suggests diversifying selection. Because diversifying selection is unlikely to affect all nucleotides in a gene, ω averaged over all sites is rarely > 1. We focused on detecting positive selection that affects only specific codons in Avr3D1 by applying the maximum likelihood method CodeML implemented in the Paml software (phylogenetic analysis by maximum likelihood, Yang, 1997, 2007).

Results

Differences in virulence map to an effector gene cluster on chromosome 7

To identify the gene(s) responsible for the differences in virulence between 3D1 and 3D7, we generated a new linkage map based on the completely assembled genome of the parental strain 3D7 (Plissonneau et al., 2016). Mapping onto the new genome sequence provided twice as many SNP markers and enabled the identification of additional crossovers that allowed us to reduce the number of candidate genes in the previously identified virulence QTL on chromosome 7 (Stewart et al., 2018). The new map yielded a narrower QTL interval (logarithm of the odds, LOD = 41.5, P < 10−15) located within the original QTL interval. The 95% confidence interval for the new QTL in 3D7 spanned 75 kb and contained only four of the 35 genes identified in the original QTL, including Mycgr3T105313, Zt09_7_00581, Mycgr3T94659 (Zt09_7_00582) and the predicted SSP‐encoding gene QTL7_5. A manual RNA‐seq‐supported reannotation in 3D7 of the confidence interval revealed two additional genes predicted to encode SSPs, which were named SSP_3 and SSP_4 (Fig. S1; Table S2). Zt09_7_00581 was reannotated as also encoding a predicted SSP after identifying an upstream start codon (Fig. S1; Table S2). The four genes predicted to encode SSPs in 3D7 formed a cluster of putative effectors.

Avr3D1 recognition contributes to quantitative resistance

In contrast to SSP_3 and SSP_4, the genes Zt09_7_00581 and QTL7_5 are highly expressed during infection (Stewart et al., 2018; Fig. S1b). Therefore, we considered them as the best candidate genes to explain the virulence QTL and they were selected for functional validation. Knockout mutants in both parental strains were generated by targeted gene disruption and used for virulence assessments in cv Runal. Mutants in QTL7_5 in the 3D7 and 3D1 backgrounds (3D7∆qtl7_5 and 3D1∆qtl7_5) did not show an altered phenotype when they were scored for host damage (Fig. S2), suggesting that QTL7_5 is not involved in virulence on cv Runal. Similarly, the virulence phenotype of the Zt09_7_00581 mutant in the 3D7 background (3D7∆avr3D1) was unaltered compared to the wild‐type (Figs 1a, S3), but disrupting Zt09_7_00581 in 3D1 (3D1∆avr3D1) led to faster development of necrotic lesions and to the production of more pycnidia compared to the wild‐type 3D1 (Figs 1a, S3, S4). Phenotypic alterations of the knockout lines in 3D1 were specific to in planta conditions, as no developmental alterations were observed when the mutants were grown on solid media used for stress assays (Fig. S5). The fact that Zt09_7_00581 negatively affects virulence in 3D1 but not in 3D7 and that in vitro growth is unaffected by gene deletion suggests that this gene encodes an avirulence factor, so we renamed this gene Avr3D1. Although Avr3D1 hinders the progression of the infection by 3D1, the avirulent strain is able to infect and produce pycnidia. Thus, Avr3D1 triggers a quantitative resistance response.

To determine if 3D7 modulates the expression of Avr3D1 to escape recognition, we quantified expression levels during infection for both strains. The expression pattern of Avr3D1 in the virulent 3D7 strain was similar to 3D1, demonstrating that 3D7 is able to infect despite highly expressing Avr3D1. Avr3D1 expression was high during the entire asymptomatic phase, peaking before the switch to the necrotrophic phase but dropping rapidly after the first symptoms appeared (Fig. S6), indicating a role for this SSP in host colonisation, possibly during the asymptomatic phase, the switch to necrotrophy, or both.

Avr3D1 is recognised by different wheat cultivars harbouring Stb7

To determine if recognition of Avr3D13D1 is mediated by a specific resistance protein, a set of 16 additional wheat cultivars was assessed for resistance against 3D1 and 3D1∆avr3D1. Three (Estanzuela Federal, Kavkaz‐K4500 L.6.A.4 and TE‐9111) of 16 cultivars exhibited a significantly lower level of resistance against 3D1∆avr3D1 compared to 3D1 (Figs 2, S7), suggesting the presence of a host‐specific factor contributing to resistance against 3D1, possibly a resistance protein. In none of these three cultivars did the presence of Avr3D1 completely abolish lesion development and pycnidia production, demonstrating that the quantitative nature of Avr3D13D1‐induced resistance is a general phenomenon and not restricted to cv Runal. All three cultivars that exhibited Avr3D13D1‐induced resistance were reported to carry the resistance gene Stb7 (Brown et al., 2015) and are also likely to carry the linked resistance gene Stb12 (Chartrain et al., 2005), leading us to propose Stb7 and Stb12 as putative candidate resistance proteins recognising Avr3D13D1.

Figure 2

Avr3D1 is specifically recognised by four wheat varieties. Violin plots showing the percentage of leaf area covered by lesions (PLACL) produced by Zymoseptoria tritici wild‐type (wt) 3D1, the Avr3D1 knockout (3D1∆avr3D1) and the mock control in 17 wheat varieties. Harvesting time points varied because of cultivar‐specific infection dynamics. Red dots represent the median of at least 10 leaves (except for the mock treatments, for which at least four leaves were used), error bars represent 95% confidence intervals of the medians and black dots represent individual data points. Asterisks indicate statistical differences between wild‐type and knockout (P < 0.01, Kolmogorov–Smirnov test). Synth. CS, Synthetic Chinese Spring; ST6, Estanzuela Federal; M6, M6 synthetic (W‐7984); Kavkaz‐K4500 Kavkaz‐K4500 L.6.A.4; dpi, days post‐inoculation. This experiment was repeated with cultivars Runal, Kavkaz‐K4500 L.6.A.4, ST6, TE‐9111, Arina, Titlis, M6 and Bulgaria‐88 and similar results were obtained.

The effector cluster resides in a highly dynamic region of the genome

Effector genes are located in plastic, transposable element‐rich regions of the genome in many fungal pathogens (Soyer et al., 2014; Dong et al., 2015; Faino et al., 2016). We explored the plasticity of the genomic region harbouring the effector gene cluster in order to understand the evolution of Avr3D1. With this aim, we performed alignments of the QTL of the 3D7 genome to the genomes of 3D1, the reference strain IPO323 and Swiss strains 1E4 and 1A5. These alignments revealed the absence of SSP_3 and SSP_4 in 3D1, IPO323 and 1E4 and the absence of SSP_3 in 1A5 (Figs 3a, S8). To gain further insight into the plasticity of this effector cluster, we extended our analysis using Illumina genome sequences of 128 Z. tritici strains obtained from four different field populations located on three continents (Hartmann & Croll, 2017; Hartmann et al., 2017). SSP_3 and SSP_4 were absent in 65% and 42% of the strains, respectively, whereas Avr3D1 and QTL7_5 were present in all or 95% of the strains, respectively (Fig. 3b). The presence/absence polymorphisms exhibited by several SSP‐encoding genes in this cluster highlight the dynamic nature of the genomic region harbouring the virulence QTL.

Figure 3

A dynamic and effector‐rich region on chromosome 7 is associated with quantitative virulence. (a) Synteny plot comparing the quantitative trait locus (QTL) for virulence between Zymoseptoria tritici strains 3D7 and 3D1. The borders of the 95% confidence interval of the QTL in 3D7 are marked by black vertical lines. Genes are represented by red arrows and transposable elements are represented by blue blocks. Collinear sequences between the two strains are shown in different shades of brown indicating sequence identity. (b) Presence/absence polymorphisms of genes predicted to encode small secreted proteins (SSPs) and transposable elements (TEs) in different populations of Z. tritici and in four closely related species. The pie charts in red shades show the presence/absence polymorphisms of genes encoding SSPs. The blue pie charts display the presence/absence polymorphisms of transposable elements in Z. tritici strains up‐ and downstream of the gene Zt09_7_00581. The transposable elements were classified according to the three‐letter code described by Wicker et al. (2016): the first letter indicates the class (R, RNA class; D, DNA class); the second letter indicates the order (L, LTR; I, Line; T, TIR; Y, Crypton); and the third letter indicates the superfamily (C, Copia; G, Gypsy; L, L1, I, I; H, ; M, Mutator; T, Tc1‐Mariner; X, unknown). nd, Not detected; na, not analysed; n, number of strains; NoCat, no category; CH, Switzerland; ISR, Israel; AUS, Australia.

To investigate whether Avr3D1, SSP3, SSP4 and QTL7_5 originated after speciation, we analysed Z. tritici sister species to determine if they contained homologues of the genes. A homologue of Avr3D1 was identified in all examined strains of Zymoseptoria pseudotritici and Zymoseptoria ardabiliae, but not in Zymoseptoria brevis or Zymoseptoria passerinii, suggesting that Avr3D1 originated before Z. tritici speciation. Homologues of QTL7_5 and SSP_3 were found in only two of four strains of Z. ardabiliae but not in Z. pseudotritici (Fig. 3b). Homologues of SSP_4 were not identified in the other Zymoseptoria species, indicating that this gene may have originated after Z. tritici speciation.

We extended our investigation on the genomic plasticity of the effector gene cluster to consider the presence of repetitive elements and transposable elements. Two insertions of transposable elements (of 44.5 kb and 9.5 kb, respectively) flanked the four SSP‐encoding genes in 3D7, but not in 3D1, where a different transposable element insertion was present upstream of the QTL (Fig. 3a). The insertion upstream of the SSP genes in 3D7 consisted of an island of 10 different transposable elements, located 1.3 kb upstream of the start codon of Avr3D1. The closest transposable element to Avr3D1 is a DNA transposable element from the Crypton superfamily, which is relatively rare in Z. tritici. Upstream of the Crypton element, three different long terminal repeats (LTRs) from the superfamily Gypsy, the most frequent retrotransposons in Z. tritici (Grandaubert et al., 2015), were inserted. A Gypsy LTR was also inserted only in 3D7 1 kb downstream of the effector cluster. No transposable element insertions occurred in the QTL region of the reference strain IPO323 or the Swiss strains 1E4 and 1A5 (Fig. S8). As in 3D1, transposable element insertions upstream of the QTL were identified in 1E4 and 1A5 (Fig. S8). Although all the insertions were upstream of the gene Zt09_7_00580, they were located at different positions and classified as different superfamilies (Copia in 3D1 and Mutator in 1E4 and 1A5). We extended the analysis of chromosomal rearrangements to the 132 global strains. Remarkably, we observed that 18% of these strains contained at least one transposable element within 6.5 kb upstream of the cluster. Furthermore, seven different insertions were identified between Avr3D1 and QTL7_5. The inserted transposable elements belonged to different superfamilies and were located at various positions (Fig. 3b), suggesting that several different insertion events occurred independently. Thus, the effector cluster resides in a highly dynamic region of the genome, in accordance with what has been previously described for other pathogenic fungi in which effectors reside in fast‐evolving regions of their two‐speed genome (Raffaele & Kamoun, 2012).

Avr3D1 is highly polymorphic in four global Z. tritici field populations

Escape from recognition is often mediated by modifications in avirulence gene sequences. Therefore, we explored sequence polymorphisms of the avirulence gene Avr3D1. In strain 3D1, the avirulent allele of Avr3D1 (Avr3D1 ) encodes a protein of 92 amino acids with a predicted signal peptide of 21 amino acids and a high number of cysteines (eight residues, 11.3%). Avr3D1 has three exons, of which only exon 1 and exon 2 contain coding DNA. The sequence polymorphism of Avr3D1 was analysed in the same four global Z. tritici populations used for transposable element presence/absence analyses. Among these 132 strains, 31 different alleles were identified, encoding 30 different protein variants, all of which were population‐specific (Fig. 4a). Strikingly, the 500 bp upstream of the start codon and the 500 bp following the stop codon showed lower diversity (πUp flanking = 0.0179; πDown flanking = 0.0023) than the CDS (πCDS = 0.067). In addition, nucleotide diversity was much lower in the first intron (πintron1 = 0.0003) and the signal peptide sequence (πsp = 0.0112) compared to the sequence encoding the mature protein (πmature protein = 0.068, Fig. 4). This pattern is consistent with accelerated diversification of the CDS, as confirmed by the high ratio between nonsynonymous and synonymous mutations (dN : dS) in the populations (Notes S1; Fig. S9). According to the codon‐based maximum likelihood approach, 58 of 96 codon sites were estimated to be under purifying selection, three were neutral and 35 were under diversifying selection (Fig. 4b; Notes S1; Table S3), suggesting that strong diversifying selection has led to high sequence polymorphism of Avr3D1. We hypothesise that numerous adaptive mutations have occurred in this avirulence gene, most probably to counteract recognition by a resistance protein.

Figure 4

Avr3D1 is highly polymorphic and exhibits the signature of diversifying selection. (a) Phylogenetic tree of the protein sequence of Avr3D1 generated from 132 Zymoseptoria tritici strains from four populations and the reference strain IPO323. (b) Upper panel: representation of the Avr3D1 gene including two introns (red arrows), the coding DNA sequence (CDS) of the mature protein (blue arrow) and the signal peptide (yellow arrow). Conservation and gap fractions of each nucleotide in 132 global Z. tritici strains is shown for each nucleotide as green and blue vertical bars, respectively. Nucleotide diversity (π) of Avr3D1 is represented with a black line. Lower panel: protein sequence encoded by the avirulent allele Avr3D1 . Conservation is shown as vertical bars in green, orange and red, representing purifying, neutral or diversifying selection, respectively, as determined by analysis of dN : dS ratios. Residues under significant (P < 0.01) diversifying selection are labelled with red asterisks. Gap fractions are shown as blue vertical bars. The consensus sequence and sequence diversity are depicted as sequence logos.

Despite the high protein diversity, the amino acid substitutions did not affect the signal peptide and did not occur in any of the eight cysteine residues, indicating that the overall backbone structure of Avr3D1 is conserved. Remarkably, in the orthologues in Z. pseudotritici and Z. ardabiliae (with 60.2% and 53.5% protein identity, respectively) all the cysteine residues were also conserved (Figs 4, S10, S11). This conservation of the overall protein structure may indicate a general role for Avr3D1 in host colonisation that was preserved after speciation.

Substitutions in Avr3D1 lead to evasion of recognition

The Avr3D1 variants in the avirulent strain 3D1 and in the virulent strain 3D7 share 86% sequence identity as a consequence of 12 amino acid substitutions and one gap in 3D7 (Fig. 1d). To determine the impact of these differences on recognition, we ectopically expressed the 3D1 (Avr3D1 ) and the 3D7 (Avr3D1 ) alleles of Avr3D1 under the control of the promoter from 3D1 in the knockout background and tested the ability to complement the phenotype of 3D1∆avr3D1. Avr3D1 fully complemented the virulence phenotype of 3D1∆avr3D1 with respect to both lesion development and pycnidia production. However, Avr3D1 did not alter the phenotype of 3D1∆avr3D1, indicating that Avr3D13D1 but not Avr3D13D7 triggers an immune response in cv Runal (Fig. 1c). Moreover, expression of Avr3D1 under the control of the promoter from 3D1 in the 3D7∆avr3D1 background led to a significant reduction in disease (avirulence), whereas expression of Avr3D1 did not alter the phenotype in the same genetic background (Fig. 1c). Therefore, Avr3D13D1, but not Avr3D13D7, is recognised in both genetic backgrounds, demonstrating that substitutions in Avr3D1 led to evasion of recognition in the virulent strain 3D7.

Discussion

In numerous plant pathosystems, a key determinant of host specificity is the resistance protein‐mediated recognition of avirulence factors, which are often SSPs. Although 20 race‐specific large‐effect resistance genes against Z. tritici have been mapped in the wheat genome, their fungal interactors remain unknown with the exception of the resistance gene Stb6. Here, we report the discovery of a new Z. tritici gene, Avr3D1, encoding a cysteine‐rich small secreted protein that triggers quantitative resistance in wheat cultivars harbouring the Stb7 locus.

Avr3D1 recognition induces quantitative resistance

Avr3D1 is a candidate effector that is expressed during the latent phase but downregulated upon the onset of the necrotrophic phase, suggesting a function in host colonisation during the latent phase or in the transition to the necrotrophic phase. The recognition of the avirulent allele leads to a dramatic reduction in the amount of infection and pycnidia formation. This demonstrates that Avr3D1 is an avirulence factor that is likely to be specifically recognised by an Stb gene. The fact that only certain wheat cultivars recognise Avr3D1 suggests that recognition follows the gene‐for‐gene model. In contrast to what has been shown for most other avirulence factors, Avr3D1 recognition does not lead to full resistance, but instead to quantitative resistance in which the pathogen is impaired in its ability to infect, but eventually completes its life cycle. The mechanisms through which Z. tritici eventually circumvents the resistance response remain unknown. We hypothesise that the magnitude of the defence response is not strong enough to prevent the progression of the infection and/or that the downregulation of Avr3D1 during the necrotrophic phase substantially decreases the defence response. The induction of this type of partial resistance might be a shared characteristic among nonobligate pathogens that grow in the apoplast and normally do not trigger a hypersensitive response (Stotz et al., 2014). The strength of this partial or incomplete resistance response may still be sufficient to limit propagation of the pathogen under field conditions, in which case the underlying resistance gene could be a valuable source of resistance for breeding programmes. Pyramiding of Stb resistance genes is an objective in several breeding programmes because this approach is thought to be an effective and durable strategy to control STB in the field (Chartrain et al., 2004; Kettles & Kanyuka, 2016). In fact, TE‐9111 and Kavkaz‐K4500 L6.A.4, two of the cultivars that specifically recognised Avr3D1, contain at least three Stb genes and are major sources of resistance to Z. tritici (Chartrain et al., 2004). Our work might contribute to the identification of the corresponding Stb gene in the future.

In this work, we show that asexual reproduction can occur even upon induction of effector‐triggered defence. In the case of AvrStb6 recognition, Stb6 strongly hinders the progression of infection, abolishing the induction of necrosis (Kema et al., 2000; Ware, 2006; Zhong et al., 2017). By contrast, recognition of Avr3D1 triggers a weaker form of resistance that prolongs the asymptomatic phase, while allowing necrotic lesions to develop and pycnidia to form. Although we cannot rule out the possibility that other protein variants of Avr3D1 might trigger a stronger resistance that completely blocks the progression of infection, these findings highlight the continuum between qualitative and quantitative resistance in gene‐for‐gene interactions. Although we identified an avirulence gene that has a large effect on some wheat cultivars, additional factors must contribute to the differences in virulence between the two strains, because the density of pycnidia formed by the Avr3D1 knockout in the avirulent strain was still lower than the pycnidia density produced by the virulent strain. The provided data further highlight the quantitative nature of wheat–Z. tritici interactions.

Chromosome rearrangements contribute to diversification of the effector gene cluster

Avr3D1 and the three other genes in the effector gene cluster are located on the right arm of chromosome 7, which is distinctive because of its low overall expression levels (Rudd et al., 2015) and its enrichment in heterochromatin (Schotanus et al., 2015). In fact, it was postulated that this region originated from a fusion between an accessory chromosome and a core chromosome (Schotanus et al., 2015). Numerous independent insertions of transposable elements surrounding Avr3D1 were identified in 132 global strains of Z. tritici. Transposable elements are frequently described as an evolutionary force shaping adjacent regions by contributing to diversification through nonhomologous recombination or RIPs (Faino et al., 2016; Wicker et al., 2016). Given that four putative effector genes are clustered in this region, transposable elements could play a similar role in facilitating rapid evolution of these effectors, but may also enable concerted expression of effector genes during infection by chromatin remodelling (Soyer et al., 2014; Schotanus et al., 2015). In the case of this effector gene cluster, transposable elements might have contributed to the high diversity of the avirulence gene Avr3D1 and the presence/absence polymorphisms shown for the other effector genes. Sequence diversification is particularly relevant for pathogen effectors, as they are key players in the coevolution with their hosts. Indeed, sequence modifications of Avr3D1 in the virulent strain allowed an escape from recognition by the corresponding resistance protein.

Avr3D1 sequence variation to evade recognition

A common evolutionary strategy for evading recognition is the loss of an entire avirulence gene (Schürch et al., 2004; Mackey & McFall, 2006; de Jonge et al., 2012; Hartmann et al., 2017). However, loss of Avr3D1 was not observed in any of the 132 global strains, despite its location in a highly plastic genomic region, as shown by presence/absence polymorphisms for neighbouring genes and transposable elements. Other deleterious mutations such as frameshifts, premature stop codons and nonfunctional splice sites were not found. Despite the high overall diversity, all the cysteine residues and the signal peptide, two core features of effector proteins, were completely conserved. The absence of any high‐impact mutations suggests that loss of Avr3D1 may impose a significant fitness cost. We therefore hypothesise that Avr3D1 plays a crucial role in the life history of Z. tritici, although we could not demonstrate a contribution of Avr3D1 to lesion or pycnidia formation in susceptible varieties during the seedling stage under glasshouse conditions. It is possible that the role of Avr3D1 is more pronounced under field conditions or at different developmental stages, for example in adult plants. An additional hypothesis to explain the apparent dispensability of Avr3D1 is that functional redundancy masks phenotypic effects in the knockout mutants (Marshall et al., 2011; Win et al., 2012; Mirzadi Gohari et al., 2015; Rudd et al., 2015).

Conclusion

We identified a new major avirulence factor of Z. tritici (Avr3D1) that we postulate to be recognised by Stb7 or Stb12, although further analysis will be required to validate this hypothesis. Unlike what has been described for most avirulence factors, recognition of Avr3D1 does not prevent lesion formation or pathogen reproduction, demonstrating that race‐specific resistance is not always qualitative. Our comprehensive comparative genomic analysis suggests that effectors in Z. tritici are located in dynamic genomic compartments favouring rapid evolution, which may facilitate adaptation to the evolving wheat host.

Author contributions

L.M. and A.S‐V. conceived and designed experiments. L.M., F.E.H., C.P. and P.C.B. performed the experiments. L.M., D.C. and A.S‐V. analysed the data. L.M. and A.S‐V. wrote the manuscript. P.C.B., D.C., F.E.H. and B.A.M. reviewed and edited the manuscript.

Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.

Fig. S1 Manual annotation of putative effector genes in the QTL for virulence.

Fig. S2 The gene QTL7_5 does not contribute to virulence.

Fig. S3 Z09_7_00581 encodes the avirulence factor Avr3D1 and sequence modifications lead to evasion of recognition.

Fig. S4 Avr3D1 does not explain all the differences in virulence between 3D1 and 3D7.

Fig. S5 In vitro growth of the mutant lines was unaltered under several stress conditions.

Fig. S6 Avr3D1 expression peaks at the end of the latent phase.

Fig. S7 Specific recognition of Avr3D1 by certain wheat varieties leads to a reduction in pycnidia formation.

Fig. S8 Synteny plot of the QTL between five Zymoseptoria tritici strains.

Fig. S9 Frequency distribution of dN : dS ratios for all pairwise Avr3D1 haplotype comparisons.

Fig. S10 Orthologous sequences of Avr3D1 identified in five strains of Zymoseptoria pseudotritici.

Fig. S11 Orthologous sequences of Avr3D1 identified in four strains of Zymoseptoria ardabiliae.

Table S1 Primers used in this study

Table S3 Model test and parameter estimates of diversifying selection with Paml based on the total Avr3D1 data set

Notes S1 Population genetic analysis.

Click here for additional data file.

Table S2 Effector gene cluster annotation. gff file of the manually reannotated effector genes identified in the QTL

Click here for additional data file.