Genome analysis and avirulence gene cloning using a high-density RADseq linkage map of the flax rust fungus, Melampsora lini

Rust fungi are an important group of plant pathogens that cause devastating diseases in agricultural, silvicultural and natural ecosystems. Throughout history rust epidemics have caused significant losses in cereal production and today they are considered a major threat to global food security due to the recent evolution and spread of the highly-virulent Ug99 strain of the wheat stem rust pathogen (Puccinia graminis f. sp. tritici) and persistent worldwide losses due to wheat stripe rust disease (caused by Puccinia striiformis f. sp. tritici) [1, 2]. In addition, many other cultivated plants, including soybean, coffee, poplar and pine, suffer rust diseases that drastically reduce productivity [3–6] and the recent incursion of myrtle rust (caused by Puccinia psidii) into Australia threatens native ecosystems and forestry production [7]. Consequently, there is great interest in understanding how these fungal pathogens infect their hosts and how plants can be better protected from rust disease.

The infection of flax (Linum usitatissimum) by the flax rust fungus (Melampsora lini) has been an important model system for understanding rust biology since the pioneering studies of Harold Flor in the late 1930’s and early 1940’s. Flor determined that the success or failure of pathogen attack is dictated by a ‘gene-for-gene’ interaction between resistance genes in the plant and avirulence genes in the pathogen [8]. When a pathogen expressing an avirulence gene infects a plant expressing a corresponding resistance gene, the plant is able to recognise the presence of the pathogen and trigger a multifaceted defence response, which includes localised cell death at the point of attempted infection [9]. Most disease resistance genes encode cytoplasmic receptor proteins with a nucleotide binding (NB) site and leucine rich repeat (LRR) domain [10, 11], whereas avirulence genes encode small, secreted proteins that are thought to be part of a diverse array of secreted effectors with virulence functions that assist in the establishment and maintenance of pathogen infection [12–15]. Five avirulence genes (AvrL567, AvrM, AvrP, AvrP123 and AvrP4) have been identified at four different loci in M. lini [16–18]. They encode secreted proteins produced in a specialised fungal structure called the haustorium, which is a major site for signal exchange and nutrient acquisition in rust pathogens and other biotrophic fungi [19, 20].

Genome sequencing of plant pathogenic fungi has become an important approach to identify effector genes that might be important during infection. For species with relatively small genomes (~40 Mb or less) and a low percentage of repetitive DNA, high quality reference genome assemblies have been produced with large scaffolds that often represent complete chromosomes from telomere to telomere [21–25]. By contrast, species with larger genomes or a high percentage of repetitive DNA have proved more difficult to assemble (e.g. [26]). Rust fungi have repetitive genomes [27–29] that are amongst the largest known for plant pathogens, with genome sizes estimated by flow cytometry that range from 77 Mb in Puccinia triticina to 2.5 Gb in Uromyces bidentis [30–32]. In addition, they produce dikaryotic spores that are highly heterozygous [29, 33], which makes genome assembly particularly challenging. The wheat stem rust (89 Mb) and poplar rust (Melampsora larici-populina; 101 Mb) pathogen genome assemblies were produced by Sanger sequencing, yet these still contain 392 and 462 scaffolds, respectively [27]. Rust fungal genome assemblies produced using next-generation sequencing (NGS) technologies are considerably more fragmented [29, 34, 35]. For example, the M. lini genome assembly was produced by Illumina sequencing and contains 21,130 scaffolds, which together span 189 Mb of the predicted ~220 Mb genome. The longest scaffold in the assembly is only 239.7 kb and 50 % of the assembly comprises scaffolds of 31.5 kb or less [28].

One strategy to improve the utility of fragmented genome assemblies is to anchor the scaffolds onto a genetic linkage map. In recent years, NGS technologies have been developed that allow simultaneous identification and scoring of single nucleotide polymorphism (SNP) markers, the two most commonly-used methods being the genotyping-by-sequencing (GBS) technique described by Elshire et al. [36] and restriction-site associated DNA sequencing (RADseq; [37]). GBS and RADseq are reduced representation sequencing technologies that rely on restriction digest of genomic DNA to reduce sequence complexity, allowing the high depth of coverage required for accurate SNP identification to be achieved economically. Using these technologies to sequence mapping families it is possible to produce high-density genetic maps that contain thousands of sequence-characterised markers, which can be used for quantitative trait loci (QTL) analysis, map-based cloning, anchoring genome assemblies and genome-wide analysis of recombination rate [38–44]. In fungi, QTLs linked to temperature ensitivity and melanisation have been identified in Zymoseptoria tritici using RADseq markers [45, 46] and a GBS-based linkage map was used to anchor the Pyrenophora teres f. sp. teres genome assembly [47]. However, these powerful technologies have yet to be applied widely in fungal genetics research.

We used RADseq to identify and map genetic markers in an F₂ family derived by self-fertilisation of M. lini strain CH5 [48]. The resulting high-density genetic map comprises 27 linkage groups that contain 13,412 RADseq markers and nine of the ten avirulence genes known to segregate in this mapping family. In addition the map includes the I-1 avirulence inhibitor locus, which prevents pathogen detection by flax plants containing the L1, L7, L8, L10 and M1 resistance genes [48, 49]. The RADseq marker sequences were used to anchor 68.9 % of the CH5 M. lini genome assembly onto the genetic map. The genetic map and anchored assembly were then used to examine recombination frequency and distribution across the M. lini genome, show that loss of avirulence in some CH5-derived mutants is associated with large-scale chromosomal duplications and deletions and identify two new avirulence genes (AvrM14 and AvrL2) by map-based cloning.

We have developed a high-density genetic map for the flax rust fungus (M. lini) that contains 13,412 RADseq markers in 27 linkage groups. As 97.7 % of the 8823 recombination events detected were single recombination events, it is likely that the region of the genome covered by the map is close to marker saturation given the amount of recombination information available in the CH5 F₂ family. Despite the high marker density, the map contains 27 linkage groups whereas the M. lini haploid genome contains 18 chromosomes [55]. In part, this may be due to the fact we selected markers for map construction that are codominant, predicted single-copy sequences with no strong segregation distortion. Therefore, linkage groups might be broken in regions where markers could not be developed (due to lack of sequence diversity or excessive diversity between the parental strains) or in regions composed entirely of repetitive DNA, provided that sufficient recombination occurred in the region to produce a genetic gap that exceeded the 20 cM threshold used in map construction. Gaps exceeding 20 cM are also expected in regions of strong segregation distortion such as fungal mating type loci, two of which are expected to be heterozygous in every member of the CH5 F₂ family [56]. Nevertheless, despite M. lini having a repetitive genome [28], at least nine chromosomes are represented by a single linkage group and the high marker density, low level of missing data (0.78 %; Additional file 16) and low estimated marker misplacement rate (1.52 %) makes this genetic map a valuable tool for genomics research.

Genetic maps have been used to anchor a number of fungal genome assemblies, including Fusarium graminearum and Pyrenophora teres f. sp. teres [47, 57, 58]. F. graminearum has a high quality genome assembly [24], 99.8 % of which could be anchored using a genetic map of just 235 markers [57]. By contrast, we were only able to anchor 68.9 % of the M. lini genome assembly onto the genetic map and provide unambiguous orientation to 687 scaffolds. This is due to the fact that the number of genetic markers required to anchor a genome assembly is inversely proportional to the size of sequence scaffolds. Although the RADseq markers are an average of 16.4 kb apart (generating a dense genetic map), 50 % of the assembly is comprised of scaffolds of 31.5 kb or less (a highly fragmented genome assembly), which limits improvement of the anchored assembly.

During the anchoring process we identified 190 scaffolds that contained markers from different linkage groups, or different parts of the same linkage group. The strict alignment parameters used make it unlikely that this is due to incorrect marker placement on the reference genome and none of these markers showed aberrant segregation or high read depth that would indicate they were derived from multi-copy sequences. Instead, the most likely explanation is that the scaffold sequences are chimeric, a common error in assemblies of complex genomes resulting from the inability of sequence assemblers to resolve repetitive DNA stretches from short read sequence data [59, 60]. Scaffolding errors are very difficult to detect by analysis of the assembly alone, prompting some authors to suggest that genetic maps are an essential component of any genome sequencing project [61–63].

Recombination is an essential part of the meiotic cell division process and at least one recombination event (or crossover) is required per chromosome pair to ensure proper segregation of homologous chromosomes during the first meiotic division. However, as aberrant segregation can lead to severe genetic disorders (such as Down syndrome in humans), in most species recombination is tightly controlled to produce a low number of widely spaced recombination events on each chromosome [64]. In Caenorhabditis elegans, Arabidopsis thaliana and many mammalian species, only 1–2 crossover events occur per chromosome pair [65–68]. By contrast, we saw an average of 114.6 recombination breakpoints per F₂ individual and for the largest linkage group (LG1) up to 22 breakpoints were detected in some F₂ individuals. Based on the former observation, we infer an average of 6.4 recombination events per chromosome pair during meiosis in M. lini, which is in excess of the number of recombination events thought to be required for correct chromosome disjunction alone. High recombination rates may be typical for fungi as an average of 90.5 recombination events per meiosis have been detected across the 16 chromosome pairs of a diploid Saccharomyces cerevisiae strain [69] and many fungi have relatively large genetic maps given the physical size of their genome. For example, a genetic map of 3006 cM has been constructed for the 90 Mb genome of the pine fusiform rust fungus, Cronartium quercuum f. sp. fusiforme [30, 70]. This is the only other rust fungus for which a comprehensive genetic map currently exists and the size of the map suggests a comparable rate of recombination to that seen here in M. lini. Understanding what controls recombination is an active area of research in many species (see [64, 71]), but we know very little about what controls this process in filamentous fungi. The high rate of fungal recombination is of great benefit in experimental studies as it means that high-resolution mapping can be conducted using relatively small population sizes. In a biological sense, the rapid shuffling of alleles and potential for intragenic recombination (including at avirulence loci) might produce new genetic combinations that are better adapted to changing environmental conditions.

Despite an overall high rate of recombination in the M. lini genome, some regions were devoid of recombination and 59 coldspots were identified in which a recombination bin was associated with a physical region of 200 kb or more (Fig. 2c). In many species strong suppression of recombination is observed at centromeres, which are sites of sister chromatid adhesion during cell division and play a major role in ensuring that chromosomes segregate correctly during mitosis and meiosis [72, 73]. Therefore, some of the coldspots observed in the M. lini genome might correspond to centromeric regions. However, as some linkage groups contained more than one coldspot (Fig. 1), other explanations are possible. Low recombination rates have been associated with non-centromeric heterochromatin and other repeat-rich regions [74, 75] and we observed an increase in mean repetitive DNA content as recombination rate decreased (Fig. 2d). Alternatively, lack of recombination could be caused by chromosomal inversions, translocations or regions of low sequence homology that are specific to the parental combination used to construct the CH5 F₂ mapping family.

We used genetic markers to examine five mutants that arose spontaneously during asexual propagation of strain CH5 (see [76]) and, in three cases, found that conversion to virulence was associated with large-scale chromosomal alterations. In MS45 and MS48, over 1 Mb of sequence that contained the AvrM14 gene was deleted from the strain H chromosome corresponding to the end of LG19 (Figs. 3a and 4) and a small duplication of the strain H chromosome corresponding to the end of LG10 was also observed in these mutants. The relationship between these two events is currently unknown, but perhaps deletion of LG19 can only be sustained if there is an associated duplication of a function located on LG10. Alternatively, the tip of the LG10 chromosome could have replaced the region deleted from the LG19 chromosome during a translocation event.

MS51 sustained the most substantial and complex chromosomal rearrangements of the mutants examined. In this mutant, loss of AvrN was associated with deletion of over 4 Mb from the strain C chromosome corresponding to LG27 and half of LG5 (Fig. 3b). The deletion in LG5 began in a large recombination coldspot at 164.69 cM that might represent a centromere, suggesting that the deletion may involve loss of a complete chromosome arm. In addition we observed duplication of the entire strain H chromosome corresponding to LG5/LG27 and duplication of strain C sequences corresponding to the central region of LG2 and the end of LG14 (Fig. 4). Again, the relationship between these changes is unclear. However, the fact that deletion of LG5/27 from the strain C-derived nucleus was accompanied by duplication of LG5/27 from the strain H-derived nucleus suggests that there may be physical association or communication between the two haploid nuclei, which have traditionally been considered separate entities in this dikaryotic life stage. Chromosome length polymorphisms and dispensable chromosomes have been observed in many fungal species [77, 78] and loss of heterozygosity has been associated with changes in pathogenicity in the oomycete pathogen, Phytophthora capsici [79], but to our knowledge this is the first demonstration of gross chromosomal changes affecting virulence in a rust fungus. Large-scale chromosomal changes may impose a fitness cost in the field. However, it will be interesting to determine whether such changes contribute to genome evolution in natural populations of rust fungi, particularly for heteroecious species that require two host species to complete their life cycle in regions where reproduction is necessarily asexual due to absence of the sexual stage host.

A major aim of our research was to identify genes controlling avirulence in M. lini and to this end we identified co-segregating markers for six avirulence genes and the I-1 inhibitor gene, none of which had been cloned. In addition, we successfully identified the AvrM14 and AvrL2 genes by map-based cloning. Notably, although recognised by allelic host resistance genes that encode highly homologous proteins [80–82], AvrL2 is unrelated to AvrL567 and similarly AvrM14 is unrelated to AvrM. AvrM14 is a dual specificity avirulence protein that is recognised by both the flax M1 and M4 resistance proteins. Significantly, AvrM14 has homology to proteins of the nudix hydrolase superfamily and is the first rust avirulence protein for which a prediction of biochemical function can be made from the protein sequence alone. Nudix hydrolases are commonly found in both prokaryotes and eukaryotes and are responsible for hydrolysing a wide range of organic pyrophosphates [83]. Several nudix hydrolases have been implicated in plant-pathogen interactions. The Phytophthora sojae Avr3b avirulence protein has proven nudix hydrolase activity [84] and effectors from Colletotrichum truncatum (CtNUDIX; [85]) and Ralstonia solanacearum (Hpx26; [86]) are predicted nudix hydrolase proteins. An Arabidopsis thaliana nudix hydrolase has also been identified that negatively regulates plant defence signalling (NUDT7; [87]). Therefore, secreted pathogen nudix hydrolases might play an important role in manipulating defence activation and cell death responses by the plant, thereby facilitating infection.

AvrL2 is located in the region of the M. lini genome with the lowest recombination rate, which may be a centromere or other heterochromatic or repeat rich region. A number of avirulence genes have been identified in repetitive regions with unique characteristics, such as Leptosphaeria maculans effectors in AT-rich isochores [22] and Phytophthora infestans effectors in gene-sparse genomic regions [88]. The region surrounding AvrL2 shares neither of these characteristics (Additional file 17). However, we found at least two independent mutants in which conversion to virulence was associated with deletion of the AvrL2-A gene and another individual that contained duplication of the AvrL2 locus and surrounding region. This suggests that the AvrL2 locus is prone to rearrangement, reminiscent of L. maculans AvrLm1 and several Magnaporthe oryzae avirulence genes that are located in unstable genomic regions [89–91]. By contrast, virulence alleles at the AvrL567, AvrM, AvrP123 and AvrP4 loci all result from sequence diversification leading to the production of avirulence protein variants that are not recognised by flax plants containing the corresponding resistance genes [92]. Of these cloned avirulence genes, only AvrL2 is located in a recombination coldspot, therefore the genomic context that influences recombination rate may also shape evolution at different M. lini avirulence loci.

The AvrL2 family proteins have no homology to proteins of known function and this is a common feature of pathogen avirulence proteins. However, they do contain several small predicted intrinsically disordered regions, which are known to provide structural flexibility in many proteins involved in cell signalling pathways by allowing them to adopt different conformations in response to post-translational modification and nucleic acid or protein binding [93, 94]. The predicted 114 amino acid, mature AvrL2-A protein contains three small intrinsically disordered regions that together comprise 36 % of the protein sequence (Fig. 6d). Disordered regions have been identified in a number of pathogen effectors. Structural studies have shown that the oomycete Avr3a and Avh5 effectors contain N-terminal disordered regions [95, 96] and the bean rust RTP1p effector contains a predicted N-terminal disordered domain [97]. Moreover, intrinsically disordered regions are enriched in predicted bacterial effectors, leading Marin et al. [98] to propose these may be commonly occurring regions in pathogen effectors that provide the protein flexibility required to enter plant cells, mimic host functions or evade detection by the plant immune system.

The generation of a high-density genetic map has greatly improved the utility of the fragmented M. lini genome assembly and we have used these genetic resources together to identify two new avirulence genes by map-based cloning. We have also identified RADseq markers that co-segregate with six other avirulence genes and the I-1 inhibitor, which will greatly facilitate future identification of these important genes. The genetic map and anchored assembly have been used to examine recombination rate and frequency across the M. lini genome and the ability to use marker read depth and allele frequency information to calculate DNA copy number has provided novel insights into chromosome rearrangements that occur during asexual reproduction. Taken together, our data suggest that variability in the M. lini genome can be achieved through a high overall recombination rate, consistent with reports of high genetic diversity in natural, sexually reproducing M. lini populations [99], and an ability of the asexual stage to withstand small- and large-scale chromosomal duplications and deletions, which can both impact upon the virulence phenotype. Several mapping families have been produced for other rust fungi, including the wheat stem rust, poplar rust and pine fusiform rust fungi [70, 100, 101]. Approaches like those described here may facilitate genome analysis and avirulence gene cloning in these rust fungi and other species.