A Near-Complete Haplotype-Phased Genome of the Dikaryotic Wheat Stripe Rust Fungus Puccinia striiformis f. sp. tritici Reveals High Interhaplotype Diversity

Abstract

A long-standing biological question is how evolution has shaped the genomic architecture of dikaryotic fungi. To answer this, high-quality genomic resources that enable haplotype comparisons are essential. Short-read genome assemblies for dikaryotic fungi are highly fragmented and lack haplotype-specific information due to the high heterozygosity and repeat content of these genomes. Here, we present a diploid-aware assembly of the wheat stripe rust fungus Puccinia striiformis f. sp. tritici based on long reads using the FALCON-Unzip assembler. Transcriptome sequencing data sets were used to infer high-quality gene models and identify virulence genes involved in plant infection referred to as effectors. This represents the most complete Puccinia striiformis f. sp. tritici genome assembly to date (83 Mb, 156 contigs, N₅₀ of 1.5 Mb) and provides phased haplotype information for over 92% of the genome. Comparisons of the phase blocks revealed high interhaplotype diversity of over 6%. More than 25% of all genes lack a clear allelic counterpart. When we investigated genome features that potentially promote the rapid evolution of virulence, we found that candidate effector genes are spatially associated with conserved genes commonly found in basidiomycetes. Yet, candidate effectors that lack an allelic counterpart are more distant from conserved genes than allelic candidate effectors and are less likely to be evolutionarily conserved within the P. striiformis species complex and Pucciniales In summary, this haplotype-phased assembly enabled us to discover novel genome features of a dikaryotic plant-pathogenic fungus previously hidden in collapsed and fragmented genome assemblies.IMPORTANCE Current representations of eukaryotic microbial genomes are haploid, hiding the genomic diversity intrinsic to diploid and polyploid life forms. This hidden diversity contributes to the organism's evolutionary potential and ability to adapt to stress conditions. Yet, it is challenging to provide haplotype-specific information at a whole-genome level. Here, we take advantage of long-read DNA sequencing technology and a tailored-assembly algorithm to disentangle the two haploid genomes of a dikaryotic pathogenic wheat rust fungus. The two genomes display high levels of nucleotide and structural variations, which lead to allelic variation and the presence of genes lacking allelic counterparts. Nonallelic candidate effector genes, which likely encode important pathogenicity factors, display distinct genome localization patterns and are less likely to be evolutionary conserved than those which are present as allelic pairs. This genomic diversity may promote rapid host adaptation and/or be related to the age of the sequenced isolate since last meiosis.

Keywords: Dikaryon; basidiomycetes; genomics; plant pathogens.

FIG 1

The Pst-104E genome assembly is highly contiguous and complete. (A) Comparison of the Pst-104E primary and haplotig assemblies with the two most complete publicly available P. striiformis f. sp. tritici genome assemblies, Pst-78 and Pst-130. The histograms and the left y axis show log₁₀ counts of contigs within each size bin. The dots and the right y axis show the cumulative sizes of small to large sorted contig lengths. Each dot represents a single contig of the given size, shown on the x axis. Each plot also shows the number of contigs or scaffolds, total assembly size, N₅₀ of the assembly, and NG₅₀ assuming a genome size of 85 Mb. NG₅₀ is the N₅₀ of an assembly considering the estimated genome size instead of the actual assembly size. This enables comparisons between different-sized assemblies. (B) Genome completeness was assessed using benchmarking universal single-copy orthologs (BUSCOs) for Basidiomycota (odb9) as proxy. The graph shows BUSCO results for Pst-104E primary (p), haplotig (h), and nonredundantly combined (ph) assemblies, in comparison to all publicly available P. striiformis f. sp. tritici genome assemblies with gene models, including Pst-78, Pst-130, Pst-21, Pst-43, Pst-0821, and Pst-887. The analysis was performed on the protein level, using publicly available gene models. An asterisk indicates the actual number of identified BUSCOs for the complete Pst-104E ph assembly before filtering gene models for similarity with genes related to transposable elements.

FIG 2

The Pst-104E genome is characterized by high levels of interhaplotype variation. (A) Summary of interhaplotype variation between primary contigs and their respective haplotigs, analyzed using Assemblytics. Each plot indicates the number of bases that are spanned by the specific variation category, which is illustrated by a cartoon. The number labeling each histogram represents the percentage of the total size of primary contigs with haplotigs that are contained within this variation type and size bin. (B and C) Two representative whole-genome alignments of primary contigs 019 and 028 with their respective haplotigs. This illustrates the large-scale variations summarized in panel A.

FIG 3

Allele transposition in the Pst-104E genome. (A to C) Dot plots of whole-genome alignments generated using the mummer toolset, where the x axis represents primary contig and the y axis shows the haplotig sequence. (A) The whole-genome alignments of haplotigs_027_xxx to primary contig 014. (B) The whole-genome alignment of haplotigs_027_xxx to primary contig 027. (C) The whole-genome alignment of haplotigs_014_xxx to primary contig 014. Black lines indicate alignments in the forward direction, and red lines indicate alignments in the reverse direction in the haplotig sequence. The black rectangles highlight an ~40-kb region in haplotig_027_006 that does not align to primary contig 027 yet aligns to a region in primary contig 014, which is not covered by an associated haplotig of 014. (D) Microsynteny analysis of this extended region, with primary contig 014 on top and haplotig_027_006 on the bottom. Gene models identified as alleles are labeled with their locus tag and shaded with a light blue background. Vertical gray shading illustrates the blastn identity between sequences on both contigs, according to the scale shown in the right bottom corner next to the sequence scale bar. Start and stop positions for each contig sequence are given at the start and the end of each contig.

FIG 4

Identification of candidate effectors based on detailed expression analysis of secreted proteins of both Pst-104E assemblies. (A) Clustering of Pst-104E secretome expression profiles for genes located on primary contigs. Blue color intensity indicates the relative expression level based on rlog-transformed read counts in spores, germinated spores, haustoria, and in wheat tissue at 6 and 9 days postinfection. For example, cluster 8 shows the lowest relative expression in spores and the highest in haustoria, compared to the other clusters. (B) Clustering of Pst-104E secretome expression profiles for genes located on haplotigs.

FIG 5

Candidate effector genes are spatially associated with conserved genes and with each other. (A) Nearest-neighbor gene distance density hexplots for three gene categories, including all genes, BUSCOs, and candidate effectors. Each subplot represents a distance density hexplot with the log₁₀ 3′-flanking and 5′-flanking distance to the nearest-neighboring gene plotted along the x axis and y axis, respectively. (B) Violin plots for the log₁₀ distance to the most proximal transposable element for genes in each category without allowing for overlap. (C) Violin plots for the log₁₀ distance to the most proximal gene in the same category for subsamples of each category equivalent to the smallest category size (n = 1,444). (D) Violin plots for the minimum distance (log₁₀) of candidate effectors and BUSCOs to each other or a random subset of genes (n = 1,444). The P values for panels B, C, and D were calculated using the Wilcoxon rank-sum test after correction for multiple testing (Bonferroni; alpha = 0.05) on the linear distance in bases.

FIG 6

The candidate effector allele status influences association with conserved genes. (A) Violin plots for the log₁₀ distance to the most proximal BUSCO for candidate effectors in each category. The Kruskal-Wallis one-way analysis of variance of all three categories showed a significant difference between the three samples (P, ~2.36e⁻⁰⁶). (B) Violin plots for the log₁₀ distance to the most proximal gene for candidate effectors in each category. The Kruskal-Wallis one-way analysis of variance of all three categories showed no significant difference between the three samples (P, ~0.08). The P values in panels A and B were calculated using the Wilcoxon rank-sum test after correction for multiple testing (Bonferroni; alpha = 0.05) on the linear distance in bases. *, Wilcoxon rank-sum test comparisons with interhaploid genome paralogs lacked statistical power due to the small sample size (n = 28).