The first chromosome-scale genome assembly of a microcyclic rust, Puccinia silphii

Abstract

Background

Rust fungi are destructive pathogens in crop plants, having led to epidemics and damaging crops all over the world. The rust Puccinia silphii (Basidiomycota) infects Silphium integrifolium, which is a member of the most speciose plant family, Asteraceae.

Results

This study analyzes the first chromosome-scale genome sequence of a rust that infects any dicot, or any species outside of the Poaceae (grasses). We found it to be the smallest genome among the available Pucciniales genome assemblies. Our final assembly was 41.7 Mb in size and consisted of 19 pseudomolecules. The genome had a BUSCO completeness score of 92.1%, with nearly all BUSCO gene losses shared with other Puccinia genomes, and gene losses concentrated in the sulfate assimilation categories. Other gene losses were unique to P. silphii, whose genome contained by far the fewest protein coding genes. A total of 10,399 protein coding genes were predicted, compared to 14,257 in the next smallest genome. Gene losses in P. silphii were concentrated in categories related to meiosis.

Conclusion

This newly assembled genome provides insights into the size limitations of Puccinia genus genomes, as well as important protein gene families that are evolving within the genus including sulphate assimilation genes and DNA replication and its role in the evolution of microcyclic rusts.

Background

Predicted continuing increases in global human population size along with climate-related challenges to agriculture and food distribution are expected to lead to substantial and increasing food shortages and global food insecurity [1, 2], due in part to the consequences of climate change, such as the increase in severe crop disease epidemics, drought and other stresses [3]. This will lead to costly crop yield losses across the planet, with even more losses expected if no crop adaptation or improvements are implemented [4, 5]. Perennial crops are an essential contributor to food production globally and offer advantages over annual species [6, 7]. They can have longer growing seasons, reduced erosion risk, and sequester more carbon in their roots compared to annuals [5]. Perennial grain, oilseed, and legume crops have also been shown to stabilize soil, prevent nutrient losses, reduce surface runoff, and retain top soil, and are thus being developed as alternatives to existing, widely-planted annual species [8,9,10]. Despite these many advantages, perennial crops face increased pressure from diseases because annual crop rotations are impossible, and pathogens can accumulate within tissue and plant debris [11]. As new perennial crops are developed, the pathogens that affect yield and mortality must be identified, characterized, and counteracted.

The rust fungus Puccinia silphii (Pucciniales, Basidiomycota) is a troublesome pathogen of Silphium integrifolium Michx. (Fig. 1), which is a native of the North American prairie, commonly known as silflower, silphium, or rosinweed. Silphium is a close relative of the common sunflower (Helianthus annuus), and is currently being domesticated as a drought tolerant perennial oilseed [12, 13]. Puccinia silphii is a microcyclic rust fungus, which is only known to produce the hyaline basidiospores and over-wintering teliospores and lacks haploid spermatia and spermatagonia, dikaryotic aeciospores, and long-distance dispersing urediniospores of macrocyclic Puccinia spp [14]. This species is only known to infect the genus Silphium and is characterized as producing orangish-brown to black pustules of teliospores on the leaves and stem [13, 15]. In general, rust-causing pathogens like Puccinia spp. are linked to decreased biomass production and yield, but can also cause a reduction in photosynthesis and other physiological changes such as reduced water use efficiency [13, 16]. Rust fungi normally possess larger genome sizes compared to other plant pathogenic fungi [17, 18], with an estimated average size within the genus being 305 Mb. They also contain large proportions of transposable elements and repetitive sequence [19].

Fig. 1

a) Color scale of Puccinia silphii infection on different leaves; b) Infection on the whole Silphium integrifolium plant, showing spots on the stems as well as the leaves; c) A detail of a Silphium leaf infected with P. silphii

Puccinia silphii infestations have been an issue for breeding plots during the Silphium domestication effort [13]. Initial plantings of breeding populations in Kansas, USA between 2003 and 2013 had relatively low incidence of rust fungi, but by 2014 and 2015 onwards there were many plants that failed to flower or produced small and late flowers due to infections from Puccinia and others [12]. In the Kansas Silphium plots, necrosis of the entire stalk was commonly seen when plant leaves and stems had severe disease (> 20% of the total area infected), which highlighted the need to better understand the biology of this disease [12].

In particular, genetic resources and knowledge of its evolution compared to other rust fungi could inform management of both agricultural and wild plant populations. To this end, we used Pacific Biosciences SMRT HiFi libraries of Puccinia silphii extracted from heavily infected Silphium integrifolium tissue to assemble and annotate a chromosome-scale genome of this rust pathogen. We asked whether the genome of this microcyclic species was substantially different in size, and whether any changes in genome size and content can explain its unusual life history and other aspects of its biology.

Methods

Tissue cultivation and DNA extraction

Inoculum in the form of Puccinia silphii telia-infected leaves was collected on November 23, 2020 from a research plot located in Gypsum County, KS, USA. Fresh leaf samples were dried thoroughly overnight in a low humidity cool drying room at The Land Institute in Salina, KS after collection. Leaves were then soaked at 20˚C in a lighted growth chamber in plastic bags containing tap water for 6 h before suspending the leaves on fiberglass mesh screens that were adhered to the inside of a lid of a plastic box containing rust susceptible clones of S. integrifolium. Household humidifiers maintained 100% relative humidity conditions for 16 h, after which the inoculated S. integrifolium plants were moved to greenhouse benches for 14 days before infections were visible. A large (1 cm) single telium was extracted from the leaf of one plant; subsequent inoculations were conducted in the same manner serially, using telia from the previous inoculation, all derived from the initial single telium (12/30/2020–2 plants, 1/12/2021–6 plants, 2/2/2021–14 plants, 2/23/2021–35 plants, 3/3/2021–32 plants with no infection, 3/31/2021–22 plants only 1 telium, 7/28/2021–12 plants with no infection, 8/18/2021–12 plants with no infection. Nutrient stress beginning 3/31/2021 may have limited the ability of the plants to become infected; subsequent inoculations in July and August may have been unsuccessful due to the age of the inoculum). Telia that developed after inoculations with the single telium were used in the next of 6 rounds of inoculations 1/12/2021. Following each inoculation, telia were carefully extracted from the leaves using a scalpel and immediately frozen using dry ice and stored in a -80˚C freezer.

Library Preparation and DNA sequencing

Rust samples were sent to the Arizona Genomics Institute at the University of Arizona for DNA extraction using a CTAB based protocol [67], with chloroform extraction prior to DNA precipitation followed by RNAse treatment, with Qubit 4 for quantification using the dsDNA HS Assay Kit, Nanodrop for quality check, then Pulse Field electrophoresis gel for DNA size profiling. The excised teliospore sample combined to 0.3 g of tissue and yielded 83.7 ng/mL for a total of 8.4 µg of DNA. Library construction of the gDNA was performed by University of Minnesota Genomics Center using the PacBio HiFi SPK 3.0 kit and sequenced across 2 Sequel II 8 M SMRT Cells.

Genome assembly, isolation of Puccinia and genome quality assessments

The raw PacBio HiFi library was assembled de novo using Hifiasm version 0.16.1-r375 [20]. The resulting assembly was checked for taxonomic partitioning using Blobtools v1.1.1 [21] and the nt/nr database (accessed on May 20th 2023). Contigs that were categorized as Basidiomycota were considered to be nuclear contigs of the Puccinia silphii genome and retained as the final assembly. Quality checks on the Puccinia assembly were performed with QUAST v. 5.2.0 [22] and BUSCO v. 5.5.0 [23] with the basidiomycota_odb10 database. QUAST assesses assembly statistics such as N50, L50 and assembly length which are measures of the completeness of the contigs in the assembly. The Benchmarking Universal Single-Copy Orthologs (BUSCO) tool contains a list of genes that should be universal among a particular clade. Genes that are missing from this list can usually be attributed to incomplete assemblies– the closer the score to 100%, the more complete the assembly. Finally, we also used TIDK [24] to check for telomeric repeats on the ends of each chromosome.

Repetitive content and gene prediction

Repetitive regions and transposable elements belonging to various classes, including LTRs, non-LTRs and DNA transposons, were annotated using the extensive de novo transposable element annotator (EDTA v1.9.3; [25]) with the parameters “--sensitive 1 --anno 1 --evaluate 0 --threads 10 --species others --overwrite 0”. The genome assembly was then masked with RepeatMasker (v4.1.1; [26]), using the soft-masking option (-xsmall) and the gene families library generated by EDTA. Proteins from the genus Puccinia were downloaded as a database to be used as extrinsic evidence in the annotation. Gene prediction was then carried out with both the BRAKER pipeline v. 2.1.6 [27] in the EP-mode and with Funannotate v. 1.8.15 [28]. Functional annotation was performed with EggNOG-Mapper v. 2.1.9 [29], using eggNOG DB version 5.0.2 [30] database of orthologous groups to predict protein functions and gene evolutionary histories. The same annotation pipeline was used for the other five Puccinia species (P. coronata, P. graminis, P. polysora, P. striiformis and P. triticina) utilized on the comparative genomics analyses (described below).

Comparative genomics among Puccinia spp

For the comparative genomics analyses we used other five different species within the genus: Puccinia coronata (ASM287312v1; [31]), P. graminis (GCA_000149925.1; [32]), P. polysora (GCA_025617555.2; [33), P. striiformis (GCA_021901695.1; [34]) and P. triticina (GCA_019358815.1; [35]), which were the other publicly available chromosome-scale assemblies in Puccinia genus. All of them had genomes available on NCBI, which were downloaded as fasta files and annotation files. We performed a comparison among all of them and P. silphii using their genome sizes, transposable elements percentage and gene content number. Because the annotation methods for most of them were different from one another, and thus couldn’t be compared without annotation bias, we re-annotated these five species using the same pipeline we had utilized for P. silphii, generating standardized values for the gene content metrics. We also generated a comparative genomics analysis based on shared orthogroups among the six species, using Orthovenn3 [36] and ProteinOrtho v. 6.3.0 [37]. Orthovenn3 takes the annotated protein sequences of the taxa to be compared and finds orthologous clusters among them. It then takes one ortholog per species and produces a peptide sequence alignment of each ortholog. It outputs a concatenated alignment of all of the orthologs found in the genome. We used this concatenated peptide alignment to infer a phylogeny in RAxML v. 8.2.12 [38] using “-f a -m PROTGAMMAWAG -p 12345 -x 12345 -# 100”, which uses a WAG + Γ substitution model [39], with branch support inferred with 100 bootstraps.

The level of heterozygosity within each assembled genome was calculated by first filtering each assembly for contigs whose best Blast hit to the NCBI nt.gz Database (https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/; Accessed March 15, 2024) belonged to the phylum Basidiomycota (Supplemental Code). The genomic reads were then aligned to the resulting filtered assemblies using bwa mem version 0.7.17-r1188 [40], and variants were called using samtools/bcftools version 1.21 [41] (Supplemental Code). Variant VCF files were then filtered for high quality variants (Q > 99), and polyallelic variants and indels were removed. The resulting number of variants was then divided by the respective filtered genome size to obtain the heterozygosity estimates.

Functional enrichment analyses

We were interested in patterns of gene loss or gain within the Puccinia genus, regarding categorized functions. These analyses were performed in two scales: In the narrow sense, among the BUSCO gene set, and broadly among all annotated protein-coding genes. In the narrow scale analysis, we aimed to test the hypothesis that consistently low (< 95%) BUSCO completeness scores among available Puccinia genomes represents a biological pattern of gene loss related to pathogen evolution rather than technical artifacts of incomplete genome assembly. We therefore identified the BUSCO genes missing in all six, five of six, four of six, and so on, down to the genes missing in one of six Puccinia genomes, and compared the number of observed cases with X missing BUSCO genes to expectations from 10,000 permutations– in which the number of missing BUSCO genes in each taxon was held constant but the identity of these genes was randomly shuffled (see Supplementary R code).

We then examined whether the genes missing in all six genomes were enriched for certain biological functions (i.e. gene ontology [GO] terms) using a Fisher’s Exact test implemented in GOATOOLS v. 1.2.3 [42]. The GO terms for Basidiomycota BUSCOs were obtained using EggNOG-mapper v 2.1.9 [29] using the Eggnog DB version 5.0.2 with the BUSCO consensus ancestral peptide sequences as input for the search algorithm DIAMOND v. 2.0.15 [43]. The background gene set for this test was all 1764 Basidiomycota BUSCOs. The required go-basic.obo file (data-version: releases/2022-03-10) was downloaded from http://geneontology.org.

Our broad scale analysis then focused on patterns of genome-wide loss and gain of genes. For the functional enrichment of all protein-coding genes of each species, we again used EggNOG-mapper 43 to predict Gene Ontology terms. OrthoFinder v. 2.5.4 [44] and ProteinOrtho v. 6.3.0 [37] were then used to identify the orthologous groups of genes and proteins that are shared among all studied species or unique to individual species. In order to understand the functional signatures of genes and test for overrepresentation of functional categories in Gene Ontology (GO) terms, we then used results from EggNOG-Mapper and ProteinOrtho to perform a gene enrichment analysis with the topGO package v. 2.50.0 [45] in R [46]. A Fisher’s Exact test was used to determine significant enrichment compared to the null expectation. We tested for enrichment on two sets of genes - those that were present only in Puccinia silphii (and absent in all other five) and those that were absent only in P. silphii (and present in all other five). Enrichment in the unique P. silphii genes was tested using the P. silphii genome to generate the null expectation while enrichment in the missing P. silphii genes was tested using the P. polysora genome to generate the null expectation.

Results

Genome assembly and genome quality

The initial assembly was 2.37 Gb in size and had 74,691 contigs, with a L90 of 60,390 contigs and a L50 of 23,477 contigs (Table S1; Figure S1). After the taxonomic categorization of contigs with Blobtools, and removal of any contigs but those identified as belonging to Basidiomycetes, the final assembly size was 41.7 Mb in a total of 164 contigs, with a L90 of 19 contigs and a L50 of 9 contigs. The majority of the remaining contigs in the assembly were categorized by Blobtools as plant contigs (Viridiplantae; Figure S2), representing contigs assembled from the host plant genome (Silphium integrifolium). Primary assembly and Haplophase 1 were identical in composition after removing off-target contigs. Since Haplophase 1 was the more complete of the two, we are calling it the primary haploid assembly and report relevant statistics here. A second, nearly complete haplophase was also produced by Hifiasm, which was 37.5 Mb in size, a L90 of 17, and a BUSCO completeness score of 88.1%. Details of the Haplophase 2 assembly are reported in Table S1. Sixteen of the chromosome-scale contigs were capped on both ends by telomeres, suggesting that those contigs represent complete chromosomes (Figure S3). Both graphs with the contig length assessment for before and after filtering are shown in Figure S1 and assembly statistics are shown in Table S1. After filtering and manual curation, our final assembly contained 92.1% of the fungal BUSCO complete genes– the highest among other assembled genomes for the genus Puccinia (Figs. 2 and 3A) --as well as a BUSCO score of 90.1% of single-copy genes.

Fig. 2

a) Comparative genomics and phylogenetic tree for Puccinia silphii and five other species in the Puccinia genus. All genomic assemblies within the Puccinia genus were annotated using a standardized pipeline with EDTA, RepeatMasker and BRAKER. While P. silphii contains the least repetitive content (TEs), it also has much less protein-coding sequence than the other Puccinia genomes. The highest values in each column are bolded. The phylogeny was inferred using maximum likelihood under a WAG + Γ substitution model in RAxML and is rooted in the outgroup Melampsora larici-populina

Fig. 3

a) Genome size for each Puccinia species in this study, with bars colored according to the BUSCO completeness score. b) Genomic content broken down by element type. c) A comparison of the genome size versus the percent taken up by coding DNA sequence. Note the log scale of the x-axis

Repetitive content and gene prediction

A total of 10,605,109 bp (25.42% of the genome) were masked as repeat elements in the genome of Puccinia silphii, with 4.58% of the total length being LTRs from the class mdg-4/Ty3, 0.42% were unknown LTRs, 2.24% were TIRs from the class Mutator, and 12.31% were unidentified repeat regions (Table 1; Fig. 3B). Some 10,399 protein coding genes were then predicted by BRAKER3, while Funannotate annotated 7,357. The Funannotate annotation contained 29% fewer predicted genes. We found that the number of genes predicted by BRAKER compared favorably to the numbers predicted by the original genome publications (P. coronata: our 20192 vs. their 17877 and 17294; P. striiformis: our 16067 vs. their 15303; P. polysora: our 21827 vs. their 20802 and 21519; P. graminis: our 14257 vs. their 17773). The outlier to the agreement between our annotations and the original is P. triticina, with BRAKER predicting 14,718 vs. the original prediction from 27,678 to 26,384, which we attribute to their usage of an older software which may overpredict genes by over-splitting them [51]. Due to the general concordance between BRAKER predicted numbers and those originally annotated in [31,32,33,34,35], we proceeded with analysing the genes predicted from BRAKER. A BUSCO score of 84% was obtained for complete BUSCO genes in the annotation from BRAKER. Of those, 8,675 of the protein coding genes were functionally annotated by EggNOG-Mapper.

Table 1 EDTA annotation of transposable elements in detail for Puccinia silphii, showing the proportion of each TE class found in the genome

Comparative genomics and functional enrichment

In our comparative genomics analysis we compared Puccinia silphii with other five species of Puccinia: P. coronata, P. graminis, P. polysora, P. striiformis and P. triticina. In the phylogeny inferred from the multiple sequence alignment, we found P. silphii to fall outside of a clade containing P. coronata, P. graminis, P. striiformis, and P. triticina, with P. polysora appearing earlier diverged relative to the rest of the taxa. These relationships were observed in [33].

We observed major differences among these species, especially comparing P. silphii’s genome size and gene content to its congeners. P. silphii has the smallest genome size in the comparison: 41.7 Mb in size, compared to 850 Mb of P. polysora, 124 Mb of P. triticina, 99.2 Mb of P. coronata, 80.4 Mb of P. striiformis and 88.6 Mb of P. graminis (Fig. 2). Despite having the smallest genome, which might indicate an incomplete assembly, it had the highest BUSCO score (92.1%) for complete BUSCO genes. The others scored 91.0%, 91.0%, 88.5%, 90.9% and 88.4%, respectively. Our results show that P. silphii has had a reduction in its gene content number, having only 10,399 predicted genes in its genome compared to the other species: 21,827 genes for P. polysora, 14,718 genes for P. triticina, 20,192 genes for P. coronata, 16,067 genes for P. striiformis and 14,257 genes for P. graminis (Fig. 2).

We found a core group of 3937 orthogroups that were shared among all Puccinia analysed (Fig. 4). The next largest category in the Upset Plot are genes that are present in all but P. silphii, totaling 797 genes. P. silphii shares the largest number of unique orthogroups with P. coronata (143), followed by P. graminis (119), then P. striiformis (112), then P. polysora (72), then P. triticina (29).

Fig. 4

Upset plot with the comparison of shared orthogroups among all Puccinia species used in this analysis, orthogroups missing from just P. silphii, and pairwise comparisons between P. silphii and each of the other species

We found that the high number of BUSCO genes missing in the published Puccinia genomes was not due to artifacts of the genome assembly. None of our ten thousand permutations had as many BUSCO genes missing in zero, three, four, five and six species, while all 10,000 permutations had more BUSCO genes missing in one or two species. These results confirm that the majority of missing BUSCO genes are likely real (i.e. informed by phylogeny) rather than a random technical artifact (p < 0.0001) such as assembly incompleteness or errors (Fig. 5). A total of 43 of the 1764 Basidiomycota BUSCOs were missing in all six of the Puccinia genomes analyzed in this study (Table S2). Gene ontology enrichment analysis of these 43 genes revealed a single, highly enriched term: “sulfate assimilation” GO:0000103 (p = 6.8e-09). Specifically, there were five genes annotated with this GO term in the background Basidiomycota BUSCO set, and all five were among the 43 missing genes (Fig. 5). These were: MET3 (BUSCO ID: 14386at5204; ATP sulfurylase catalyzing the first step of intracellular sulfate activation), MET5 (1398at5204; Sulfite reductase beta subunit), MET8 (74445at5204; Bifunctional dehydrogenase and ferrochelatase), MET10 (8706at5204; Subunit alpha of assimilatory sulfite reductase), and MET16 (61744at5204; 3’-phosphoadenylsulfate reductase; PAPS reductase), where “MET” stands for “methionine requiring” (Table S5).

Fig. 5

Colored ridges represent the expected number of BUSCO genes missing if they were missing at random due to incomplete genomic assemblies. Black diamonds represent the observed number of genes that were missing among the six Puccinia genomes the number of times listed to the left of the plot. For example, the teal curve (top row) represents the expectation that 0–1 genes should be missing in all six Puccinia if they were missing at random, whereas the observed number that were missing in all six genomes was 43

The GO term abundances of annotated genes were compared among Puccinia silphii and the other five Puccinia species in order to understand what functions of genes were acquired or lost in P. silphii. Of the set of genes that were unique to P. silphii, we found 12 GO terms in the ‘biological process’ category that were significantly enriched (Table S3). Among those, the most significant terms were related to macromolecule metabolic process, cellular nitrogen compound metabolic process and cellular biosynthetic process, among others. Of the set of genes that were missing in P. silphii, we found 27 statistically significant GO terms (Table S4). Among those, the most significantly enriched ones were related to meiotic DNA replication initiation, regulation of initiation of premeiotic DNA replication and positive regulation of initiation of premeiotic DNA replication.

Discussion

Genome assembly

The genome assembly of Puccinia silphii is the first chromosome-scale genome of a microcyclic rust, and the smallest genome known in the genus Puccinia [47]. Microcyclic rusts have only the teliospore and basidiospore stages, while macrocyclic rusts have all five spore stages (spermatia, aeciospore, urediniospore, teliospore and basidiospore) and have been better studied genomically [48]. Due to the complexity and the high number of transposable elements and repeats in rust genomes [48], only a few chromosome-scale complete genomic sequences are currently available for the group [48,49,50]. Moreover, despite the complete assembly of many fungal genomes, such sequencing efforts have missed key features of the diversity of life-history, physiological, and phylogenetic diversity of fungi, such as rusts. Here we present the first chromosome-scale assembly of a microcyclic rust fungus– P. silphii (the Silphium rust), which was achieved solely based on PacBio HiFi data without the aid of additional library types for scaffolding, such as Hi-C or Omni-C chemistries. The high quality of this 41.7 Mb genome assembly is evidenced by a low L90 of 19, ~ 92% of BUSCO completeness (the highest currently available within the genus; Fig. 54), and the presence of telomeres at the ends of most of the 19 chromosome-length contigs. Three of the contigs (contigs 7, 14, and 18) only contain a telomere at one end (Figure S3), meaning that we can’t rule out the possibility that two of these contigs may join together, resulting in 18 total chromosomes instead of 19. Other species in this genus that have chromosome-scale assemblies contain 18 chromosomes so the parsimonious expectations for P. silphii would be that it contains 18 as well. Additional data from scaffolding chemistries could resolve this ambiguity in the future. The annotation of this genome adds to a trend within the Puccinia published genomes that smaller genome size correlates with less coding DNA sequence (Fig. 3C). The genomic annotation contained 84% of BUSCO genes, which when compared to the BUSCO of 92% for the assembly leads to a reasonable conclusion that we are underestimating the true number of genes in the genome by around 8%. This is a known issue [51] with genomic annotations and improvements in the annotation software and/or manual curation of gene models [68] will only reduce the disparity between these two numbers in the future.

Complete and phased genomic assemblies of rust fungi have also helped decipher many of the reproductive mechanisms these organisms are involved in [52], such as “somatic hybridization”, in which two different strains will exchange complete haploid nuclear genomes [53]. Somatic hybridization has been documented outside of sexual recombination as a powerful generator of genetic variation [54]. In Puccinia coronata, somatic hybridization has been observed across global scales (USA, Australia and Taiwan; 55] and has been shown to help the rust adapt to host pathogen defenses [54]. Also, haplotype-aware genomic assemblies have helped to uncover high heterozygosity within nuclei of P. coronata f. sp. avenae as well as P. striiformis f. sp. tritici [54]. In general, rust fungi that undergo somatic hybridization should contain substantial genetic variation and be able to adapt quickly to sudden environmental changes [52]. Our assembly of P. silphii reveals that somatic hybridization is unlikely given the low heterozygosity contained in it relative to the heterozygosity levels of other Puccinia assemblies (Figure S4).

This newly assembled genome is an important resource for research into Pucciniales rusts, with unique features such as the superlatively small genome size, a small overall gene content, and depletion of genes associated with sulfate assimilation, and will also be essential for future genomic and pathogenicity studies. These resources will be particularly impactful as the first representative chromosome-scale genome of a rust outside of the already very-well studied grass-infecting lineages. Moreover, this is the first chromosome-scale genome assembly of any rust affecting the Asteraceae, the Asterids, and even the Dicots, which will provide insights into the cellular mechanisms of infection by the rusts outside of those infecting grasses.

Genome size and comparative genomics

Puccinia silphii has a comparatively small genome for the first published microcyclic rust genome, unlike macrocyclic rust genomes (Fig. 3A. In general, rust fungi (order Pucciniales, Basidiomycota) have genomes that are larger in size compared to other pathogenic fungi, with 305.5 Mb as the estimated average genome size for the group [47, 48]. However, there is great variation in genome size among families, genera, species and isolates, that in fact could be caused by different sequencing technologies or pipelines for assembly and not represent true differences in genome content [47]. While Puccinia silphii contained a genome size of ~ 41 Mb, Puccinia polysora is a gigabase-scale fungal pathogen, with its haploid genome size estimated at around 850 Mb [33]. Another rust with a gigabase-sized genome is Austropuccinia psidii, which causes myrtle rust, and has a haploid genome size of 1018 Mb [55]. It has been suggested that the expansion of the genome in rust fungi, when compared to other pathogenic fungal species, could be associated with the biotrophic and parasitic life cycle, and could be a sign of convergent evolution in rust fungi alongside Ascomycetes and Oomycetes, which also have had genome expansions [33]. As an example, the smallest genomes in the phylum Ascomycota are both free-living and parasitic but have lost most of their introns and transposable elements [56, 57], while the largest ones are obligately parasitic and present large genomes with a great content of transposable elements [55, 56]. Obligate biotrophy has been associated with gene loss and genome size expansions in many scenarios and species [19, 50]. Primary and secondary metabolism genes have been lost in some obligate biotrophs such as the oomycete Hyaloperonospora arabidopsis (downy mildew disease-causing oomycete) and some rust species (Basidiomycota) [50, 58]. However, the Parauncinula polyspora genome, a powdery mildew fungus species, contains a very compact genome (< 30 Mb) and has possibly diverged early during evolution [59]. Despite the relatively large genome size of powdery mildews in general, the P. polyspora is an example of a small genome size in a genus of very large genomes, a similar situation to what we are observing here. Therefore, it is plausible that the genome size of P. silphii resembles the ancestral genome size of the clade, and that the larger sizes of the published Puccinia genomes represent expansions possibly due to the parasitic lifestyle.

Other common features that have been documented for rust fungal genomes are having large proportions of transposable elements and repeats (usually more than 30%) and also encoding a large number of proteins (14,000–32,000) [47, 48]. As presented here, the Puccinia silphii genome is comprised of only 25.42% repetitive content and only 10,399 protein coding genes. In general, there is a great range of genome sizes in fungi species [60, 63], with genome reductions having occurred in many different lineages [61]. Yeast-forming lineages are one of the most prominent examples [64, 65] and in these cases, major gene losses could be related to the transition to single cellularity [58]. Changes in genome content and dynamics can be influenced by selective pressures but could also be the result of random genetic drift. Because small effective population sizes (N_e) cause drift to be more powerful than very small selective forces that would act on genome size, build-up of mildly deleterious “junk” DNA could increase genome size during population bottlenecks [62]. In contrast, theory predicts that species with large effective population sizes can evolve smaller genomes due to the more effective removal of deleterious accumulations of extra DNA, such as transposable elements, introns and gene duplications that are due to genetic drift [62]. We do observe significantly higher transposable elements content in the larger Puccinia genomes (Fig. 3B). Changes in effective population size could also be resulted from a newly acquired plant-associated lifestyle, imposing limits to reproductive modes and dispersal, or even could be due to constraints and selective forces on cell size.

Another possible explanation for the trend we see is that if P. silphii represents a reduction in genome size, it could have resulted from the simplified (microcyclic) life cycle in this species. For example, de novo origin of life-stage-specific genes has been observed in brown algae [64] to resolve intra-locus generational conflict (or in this case, conflict between the selective benefit of genes at different stages of the life cycle). It stands to reason that loss of stages of the life cycle could result in the loss of selective pressure to maintain those extra genes. We see tentative evidence in support of this hypothesis in the depleted GO categories related to the regulation of DNA replication (Table S4) which may have been lost after the transition to a life cycle that skips the aecial, uredinial, and pycnial life stages. For example, genes related to making spore-producing structures (spermagonia, aecia and pycnia) will face reduced selective pressure for retention in the genome. Microcyclic rusts may also begin to lose the regulatory framework (i.e. transcription factors, replication machinery) to initiate the development into those life stages [66].

Finally, in addition to the reduced complexity of the life cycle of P. silphii, the species has not been observed infecting a secondary host, unlike most other rusts. If it is the case that this species completes its life cycle on a single host, then the genomic underpinnings of secondary host infection will likely be undergoing losses, leading to a smaller genome and fewer coding genes. However, it is equally plausible that this species follows the two-host paradigm, and the secondary host has simply not yet been observed.

Functional evolution within Puccinia

One striking genus-wide pattern was repeatedly observed in genes that in most fungi are highly conserved. Loss of sulfate assimilation function among BUSCO genes across all six Puccinia genome assemblies suggests this genus may have followed a similar evolutionary trajectory as other obligate biotroph lineages. That is, gene losses in the metabolism of inorganic sulfur, nitrogen, and thiamine [61, 63,64,65] are among the top metabolic pathways missing among unculturable fungi (as all rusts are). Interestingly [32], found Puccinia graminis f. sp. tritici lacks the alpha and beta subunits of sulfite reductase (SiR) and the set of genes missing from P. silphii contains sulfite reductase subunit beta, but not subunit alpha, and is also missing four other genes in the sulfite assimilation pathway. This suggests that genes in the sulfate assimilation pathway are being lost differently among different species in Puccinia. This finding provides an intriguing hypothesis to test in future studies, which could take a whole-genome approach to comprehensively characterize gene loss across Puccinia as it relates to such functions. In contrast to these genus-wide patterns, only P. silphii is missing genes involved in meiotic DNA replication initiation and regulation of initiation of premeiotic DNA replication (Table S3). As this is the first microcyclic Puccinia species with a chromosome-scale genome assembly, these observed gene losses may be related to this life cycle shift. Interestingly, P. silphii appears to have gained some genes and functionality related to nitrogen cycling, metabolism of macromolecules, and cellular biosynthesis. Future research on this taxon should explore the metabolic carbon degradation capabilities of this fungus, which might be relevant for crop pathogenicity and ecosystem function.

Conclusions

Puccinia silphii has a unique genome compared to other known rust genomes, in terms of its small size and unusual genome content and is also important as the first microcyclic rust with a full chromosome-scale genome sequence, as it may yield insights into the genomic mechanisms required for the dramatic transition in life cycle. Gene loss and lower repetitive sequence make it novel and of interest for comparative genetics to other rusts. These unique characteristics also provide potential targets for controlling this pathogen in the emerging perennial crop, silflower, or indeed other dicots. The decreased representation of gene ontologies associated with meiosis poses questions for the species’ ability to undergo meiosis.