Leveraging single-cell genomics to expand the fungal tree of life

Environmental DNA surveys reveal that most fungal diversity represents uncultured species. We sequenced the genomes of eight uncultured species across the fungal tree of life using a new single-cell genomics pipeline. We show that, despite a large variation in genome and gene space recovery from each single amplified genome (SAG), ≥90% can be recovered by combining multiple SAGs. SAGs provide robust placement for early-diverging lineages and infer a diploid ancestor of fungi. Early-diverging fungi share metabolic deficiencies and show unique gene expansions correlated with parasitism and unculturability. Single-cell genomics holds great promise in exploring fungal diversity, life cycles and metabolic potential.

G enomics research enables a well-resolved phylogenetic backbone for the fungal tree of life and describes how the fungal nutritional toolkit has evolved over a billion years 1 . A complete grasp of these evolutionary patterns requires a thorough understanding of the enzymes and metabolites utilized by all fungi for accessing diverse sources of nutrition. However, much of what we know about the evolution of the fungal metabolic repertoire is biased towards fungi that are model systems, enzyme factories or important human or plant pathogens 2 . This emphasis results in an underrepresentation of uncultured lineages in the fungal genomic tree and the inaccessibility of these clades to comparative genomic analysis.
Environmental DNA amplicon surveys provide ample evidence for a high diversity of uncultured fungi that are unable to complete their life cycle in axenic conditions. Most studies recover species that are not represented in genomic databases [3][4][5] . This demonstrates a collectively limited knowledge of the true diversity of fungi, a problem exacerbated when considering early-diverging lineages, which are primarily microscopic and diverse in biotrophic groups 1,6 . Among these uncultured fungal lineages are entire hyperdiverse phyla, such as the Cryptomycota 7 . Even among the culturable earlydiverging fungi (EDF), such as the Chytridiomycota, environmental DNA studies show almost no overlap with known and sampled species 8,9 . Because genome sequencing has thus far been limited to cultured fungi, the phylogeny of EDF remains poorly resolved. A better understanding of their genomes will provide clues regarding evolutionary history, basic biology and metabolism (for example, ref. 10 ).
Cultivation-independent methods for sequencing environmental microbial taxa, driven largely by shotgun metagenomics, have been applied for over a decade. Owing in part to the complexity of metagenomic data, the application of single-cell genomics has sharply increased over the past 5 years [11][12][13] . These methods rely on isolation and lysis of individual cells with subsequent wholegenome amplification and sequencing 14 . Most of the current environmental work has focused primarily on bacterial and archaeal systems 15 (see ref. 13 for a review); however, an increasing number of eukaryotic genomes have been reported from single-cell sequencing [16][17][18][19] . Although recent efforts successfully demonstrated both single-nucleus de novo genome sequencing of the endomycorrhizal fungus Rhizophagus irregularis 20 and high-throughput microfluidics single-cell sequencing 21 , these methods are neither generalizable nor easily adoptable in other laboratories. Fungi pose several challenges for scaling up single-cell genomics. These include structural challenges, such as the presence of a cell wall and diverse cellular morphologies, as well as genomic challenges, such as multiple chromosomes and higher ploidy, mitochondrial genomes, wide GC variation and transposable elements. To develop methods for robust capture and de novo assembly of environmental fungal genomes, we designed a project using known target species to explore the technical challenges prior to isolating and sequencing more complex environmental samples.
Here, we applied single-cell genomics methods for the first time to uncultured mycoparasitic EDF from the Cryptomycota, Chytridiomycota and Zoopagomycota 22 : Rozella allomycis, Caulochytrium protostelioides, Dimargaris cristalligena, Piptocephalis cylindrospora, Syncephalis pseudoplumigaleata and Thamnocephalis sphaerospora. We focused primarily on biotrophic mycoparasites because they are an understudied group of fungi that are widely represented among EDF 23 and they sporulate abundantly when cultured with their hosts. We also included the pollen saprotroph Blyttiomyces helicus (Chytridiomycota) and Daphnia parasite Metschnikowia bicuspidata (Ascomycota) because samples of these fungi also contain DNA from numerous non-target species, similar to environmental samples. In this study, we provide genomic insights into the biology and evolutionary histories of these uncultivated species. We illustrated robust placement of novel lineages among EDF, demonstrated higher than haploid ploidy as a common characteristic of these lineages and revealed interesting gene family evolution patterns outside the Dikarya. We also highlighted common metabolic deficiencies among uncultured lineages and tested whether these deficiencies could be overcome through culturing efforts with addition of limiting reagents. Collectively, these approaches will facilitate further study on diverse and uncultured environmental eukaryotes.

Results
Single-cell pipeline successfully captures fungal genomes with high completeness. We developed and applied our single-cell pipeline ( Supplementary Fig. 1) to eight diverse target species. We recovered genomes from individual cells ('1-cell'), as well as pools of multiple cells sequenced as one library ('10-/30-/50-/100-cells'). These individual libraries (Supplementary Table 1) were combined in separate co-assemblies (a technique routinely used in microbial single-cell projects 15,24,25 ) to maximize completeness. Figure 1b summarizes genome completeness (per cent CEGMA (Core Eukaryotic Genes Mapping Approach) 26 ), assembly size and total gene content of co-assemblies and individual assemblies. Generally, CEGMA percentages increased when more cells were incorporated into a library, ranging from 2.8% (1 cell, C. protostelioides) to 96.7% (100 cell, M. bicuspidata). For all species, co-assemblies were most complete (relative to the 1-cell or 50-cell/100-cell libraries), ranging from 73.36% (P. cylindrospora) to 99.34% (D. cristalligena). Assembly sizes of the 1-cell libraries range from 0.5 Mb (C. protostelioides) to 21.1 Mb (B. helicus). The largest single-library assembly was 30.1 Mb (100 cell, B. helicus). Co-assembly sizes ranged from 10.6 Mb (C. protostelioides) to 46.5 Mb (B. helicus). Predicted gene counts ranged from 111 (1 cell, C. protostelioides) to 6,941 (100 cell, D. cristalligena) for single-library assemblies. Co-assembly gene counts ranged from 3,328 (C. protostelioides) to 12,167 (B. helicus). Overall, there is a strong positive correlation between estimated genome completeness and assembly size (Fig. 2a). Within a single taxon, the trend suggests that more cells sorted per library results in more complete libraries. With the exception of R. allomycis, > 75% of reads of a given library were incorporated into an assembly (Fig. 2b).
Single-cell assemblies provide robust fungal phylogeny. To help resolve fungal phylogeny, a maximum likelihood tree was built from whole-genome clustering data from the eight co-assemblies along with major fungal and deep-branching eukaryotic representatives (Supplementary Table 2). The phylogenetic placement of our coassembled genomes agrees with current taxonomic understanding: D. cristalligena groups with Kickxellomycotina; B. helicus and C. protostelioides with Chytridiomycota; R. allomycis with Microsporidia as sister to the rest of the fungi; and P. cylindrospora, S. pseudoplumigaleata and T. sphaerospora form a monophyletic group sister to the Kickxellomycotina within the Zoopagomycota (Fig. 1a).
Because single-cell libraries will only recover a partial genome, we determined the effect of genome completeness and input cell number on generating an accurate phylogeny. We annotated all libraries from D. cristalligena, C. protostelioides, R. allomycis and M. bicuspidata and constructed individual phylogenetic trees using each library (Supplementary Data). These topologies were compared to the co-assembly-based tree in Fig. 1a. Only three libraries (7.7%) were incongruent for that particular taxon, despite libraries having as few as 15 genes belonging to orthologous gene clusters (median: 434): one R. allomycis (1 cell, 1,095 clusters) and two C. protostelioides (100 cell and 10 cell, 1,091 and 1,090 clusters, respectively) (Supplementary Table 3). These results show that less-complete single-celled libraries can be used to generate accurate phylogenies.
Single-cell genomes indicate that higher ploidy is common in EDF. Single-cell genomes provide a unique opportunity to separate polymorphisms among cells from polymorphisms within cells. The ploidy of most EDF species is largely unknown but a few have been shown to be higher than haploid [27][28][29] . Analyses with k-mer graphs and allele frequency spectra can successfully distinguish between haploid and non-haploid organisms [29][30][31] . For C. protostelioides and R. allomycis, our results indicate haploid and triploid patterns, respectively ( Fig. 3a-d). We were able to identify putative singlenucleotide polymorphisms (SNPs) in single-cell libraries of both species, but only for R. allomycis do single-cell SNPs match those identified in isolate sequencing (Fig. 3e).
Single-cell sequencing can also identify heterozygous SNPs in diploid cells 32 . By using C. protostelioides and R. allomycis as haploid and non-haploid models, respectively, we identified five higherthan-haploid species (Fig. 3f). Moreover, because we are able to show that most variants identified in our five non-haploid species are present in multiple libraries (Fig. 3g), our patterns are consistent with heterozygous SNPs expected in all cells of an isolate. These results show that single-cell sequencing can elucidate the ploidy of fungi and suggest that the majority of EDF are non-haploid.
We also identified cell-to-cell polymorphisms among one-cell libraries of D. cristalligena ( Supplementary Fig. 2). We determined that six SNPs for a particular gene were exactly shared between two libraries (AHPZW and AHPZP), whereas one was unique to a third library (AHSAA).
Single-cell genomes highlight common deficiencies in primary metabolism of uncultured fungi. Genome sequencing of uncultured EDF allows us to explore how metabolic content has changed during the course of fungal evolution. Compared to the common ancestor of fungi ( Supplementary Fig. 3), D. cristalligena shows gains in high-level metabolic categories, whereas B. helicus has mainly losses.
Examination of the primary metabolism of genomes can reveal deficiencies responsible for unculturability, ideally resulting in successful culturing through supplementing media with missing metabolites 33 . One challenge with predicting missing function using single-cell genomes is the presence of false negatives from missing data, which we address by searching for commonalities in missing pathways across our target taxa. The conserved pathways found in 75% of 'free-living' fungi, overlaid with those missing from ≥ 5 target fungi, reflect consistent enzymatic losses (Fig. 4a).
The absence of spermidine synthase (enzyme classification (EC) 2.5.1.16) and S-adenosylmethionine decarboxylase (EC 4.1.1.50) suggests that almost all target taxa are unable to make spermidine and/or homospermidine (Fig. 4d). Spermidine is involved in the regulation of processes such as virulence and sporulation 34,35 , and deficiency in polyamines leads to auxotrophy and attenuation of virulence 36,37 .
Another deficiency common to almost all of the target genomes is in the assimilatory sulfate reduction pathway 38 (Fig. 4b) Metabolism of biotin is another common deficiency among target genomes (Fig. 4e) (Fig. 4c).
Entrance into the tricarboxylic acid cycle is accomplished through citrate synthase (EC 2.3.3.1), which is found in all freeliving and target fungi. This enzyme facilitates the conversion of acetyl-CoA to citrate. ATP citrate synthase (EC 2.3.3.8) facilitates the same reaction but also generates ATP in the process. It is absent from all target fungi except for the chytrids C. protostelioides and B. helicus (Fig. 4f).
Based on our primary metabolism results, we attempted preliminary media-supplementing axenic culturing experiments for three mycoparasitic taxa: D. cristalligena, S. pseudoplumigaleata and P. cylindrospora. We tested the efficacy of five supplements to produce axenic growth of these fungi and obtained mixed results. S. pseudoplumigaleata either did not grow or was contaminated by the host fungus. P. cylindrospora grew axenically but weakly on all media treatments, including the control plates, but was not able to complete its life cycle on any of the media formulations. D. cristalligena responded best to our experiments; this fungus grew faster and more abundantly on the media that included all five supplements when compared to all other treatments. However, even on the fully supplemented axenic media, D. cristalligena was not able to complete its life cycle ( Supplementary Fig. 4). Further experiments will be necessary to fully characterize axenic growth of these and related mycoparasitic EDF (see the Supplementary Information

Proteinase and CAZymes reveal ecology of uncultured fungi.
To explore parasitism strategies, we focused broadly on abundance patterns of enzymatic protein domains. Comparing carbohydrate-active enzymes (CAZymes; http://www.cazy.org) 40 and proteases (MEROPS; https://www.ebi.ac.uk/merops/) 41 , we found that facultative mycoparasites have a CAZyme-to-peptidase ratio most similar to plant pathogens, whereas that for obligate mycoparasites is most similar to animal pathogens (Fig. 5). B. helicus, the only saprotroph in this data set, has the most CAZymes among our single-cell genomes, with expansions of families that target pollen polysaccharides (Supplementary Table 4). To complement these broad observations, we focused narrowly on families often associated with mycoparasitism: subtilases, metallopeptidases and chitinases.

Subtilases in Zoopagomycota.
Subtilases, one of the largest clans of serine endopeptidases, are found in fungal entomopathogens, mycoparasites and plant pathogens 42,43 . Recent work has characterized novel categories of serine proteases in fungi 44,45 in the Dikarya. We explored subtilase abundance in EDF mycoparasites and found that D. cristalligena subtilase sequences form a group distinct from others, whereas the Zoopagomycotina (P. cylindrospora, T. sphaerospora and S. pseudoplumigaleata) sequences split roughly equally between known proteinase K-like proteins and the group formed by D. cristalligena (Fig. 6a). Figure 6b shows the domain architecture of orthologous subtilase proteins predicted among the Zoopagomycota and illustrates expansions specific to D. cristalligena, T. sphaerospora and S. pseudoplumigaleata.

Metallopeptidases in Zoopagomycota.
Class M36 metallopeptidases, also known as fungalysins, hydrolyze laminins, elastin, collagen and keratin 46 . The amphibian pathogen Batrachochytrium dendrobatidis has an expansion of this metallopeptidase 47 , which is presumably used in the degradation of keratin-rich amphibian skin. Orthologous metallopeptidase proteins were identified in the three species of Zoopagomycotina and D. cristalligena. A phylogenetic tree ( Supplementary Fig. 5) highlights their relationship to other fungalysin proteins and illustrates a unique expansion among the Zoopagomycotina.

Chitinases in EDF.
Chitin is a definitive component of the fungal cell wall [48][49][50] . It is broken down by chitinases of the glycoside hydrolase family 18 (GH18) and GH19 families and a recently described family (AA11) 51 . The two glycoside hydrolase families do not share similarities in protein sequence, structure or mechanism of action 52 . Furthermore, GH18 chitinases are widely distributed, whereas GH19 chitinases are described mainly from plants and act as defences to insect or fungal invaders.
There are a few known fungal representatives of GH19 chitinases, exclusively from the Cryptomycota and the Microsporidia. We also found GH19 chitinases for the first time in the Chytridiomycota and the Zoopagomycota. A small expansion was found in Rhizoclosmatium globosum, a saprotrophic chytridiomycete. Putative GH19 chitinases were also identified in D. cristalligena, Linderina pennispora and Coemansia reversa (Kickxellomycotina), and in P. cylindrospora (Zoopagomycotina). Only one protein from D. cristalligena had an additional CBM19 domain ( Supplementary  Fig. 6). D. cristalligena also has a unique expansion of the noncatalytic CBM18 family. Finally, the mycoparasites exclusively have AA11 genes, with clear expansions in D. cristalligena and T. sphaerospora (Supplementary Table 4).

Secondary metabolite expansion in D. cristalligena.
Secondary metabolites are non-essential compounds produced by fungi for various purposes, including antagonism of other microorganisms, pathogenesis and iron chelation. Few secondary metabolites have been identified from EDF 53 . Our results support this observation (Supplementary Table 5), with most of the species containing one or fewer non-ribosomal peptide synthetases (NRPSs) or NRPSlike proteins and two or fewer polyketide synthases or polyketide synthase-like proteins. However, D. cristalligena possesses 27 NRPS genes divided between two lineage-specific expansions in two specific clades (Supplementary Fig. 7). Several of these homologous NRPS proteins also show modular synteny, suggesting multiple  Fig. 7), these secondary metabolites are most closely related to epipolythiodioxopiperazine toxins that disrupt cell membranes. Trichoderma mycoparasites are known to produce several epipolythiodioxopiperazine toxins, including gliotoxin 54 . Future work is necessary to determine whether these epipolythiodioxopiperazine-like expansions are related to mycoparasitism in D. cristalligena. Most of the NRPS genes reside alone on relatively short contigs (Supplementary Table 6), precluding identification of a traditional fungal secondary metabolite gene cluster.

Hydrophobins in C. protostelioides.
Hydrophobins are small cysteine-rich proteins involved in the development of aerial hyphae in certain filamentous fungi 55 . These proteins are currently only described in the Dikarya, with no evidence of their presence in EDF 56 . Caulochytrium is the only zoosporic genus known to produce aerial stalks and sporangia reminiscent of filamentous fungi 57 . We found 14 putative hydrophobins in the single-cell C. protos-telioides proteome, all of which contained the highly conserved 8-cysteine marker region. Phylogeny and hydropathy profiles both suggest that the C. protostelioides proteins are more closely related to the group 1 proteins found in the Dikarya (Supplementary Fig. 8). Furthermore, we found one putative group 2 hydrophobin in the Mortierella elongata proteome. These findings represent the first examples of hydrophobins outside the Dikarya.

Discussion
In this study, we used single-cell genomics to create near-complete assemblies of uncultured fungi. This approach allowed us to capture an estimated 73-99% of the genome in multiple-cell co-assemblies ranging from 478 to 8,398 scaffolds. Our analysis shows that genome completeness from a single cell ranges from 6% to 88%, and that combining multiple cells can considerably increase assembly completeness from 6% to 80% (worst case) and from 88% to 96.7% (best case). Co-assemblies of genome data from different single cells further increase genome recovery while losing single-cell resolution for other analyses (for example, heterozygosity).   The co-assemblies allowed annotation of 3,328-12,167 proteins in each species, and common orthologous proteins were used to create a phylogenetic tree, placing these uncultured lineages among their cultured counterparts with robust bootstrap support. We placed Zoopagomycotina as a sister branch to Kickxellomycotina, although with minimal (69% bootstrap) support. A characteristic uniting these subphyla is that certain taxa in both groups produce merosporangia, that is, linear sporangia with few spores 58 . Although Dimargaris has been difficult to confidently place in ribosomal RNA gene phylogenies owing to a rapid rate of sequence evolution 59 , here, D. cristalligena is placed with strong (100%) support as a sister to the Kickxellales. An rRNA-based phylogeny placed B. helicus with the order Spizellomycetales and Rhizophlyctidiales but without statistical support 60 . In our maximum likelihood tree, B helicus groups with the Spizellomycetales representative Spizellomyces punctatus. Caulochytrium was placed in the Spizellomycetales by Barr 61 based on ultrastructural data. However, Barr's concept of the order has been radically reshaped by molecular phylogenetics 60 . Until this study, no molecular phylogenetic analyses have included Caulochytrium; thus, there was no clear null hypothesis for this chytrid. Although C. protostelioides grouped with the Chytridiales representative R. globosum (97% support), the precise relationship of C. protostelioides to other chytrids will require additional genome sequencing.
Single-cell sequencing is uniquely capable of addressing questions in fungal genetics regarding the organization of genetic variation. Analysis of single cells allows testing of whether detected SNPs are specific to an individual or due to variation among multiple genotypes in a variable population. The species targeted by singlecell sequencing were found to correspond to two different groups: non-haploid with high levels of heterozygosity and haploid with negligible heterozygosity (Fig. 3f). Importantly, taxa that we suspected of being diploid had > 10,000 SNPs shared across libraries, whereas taxa that were haploid typically had < 2,000 SNPs shared across libraries. One caveat is that fungi could be genetically diploid but show low heterozygosity, as is the case for Saccharomyces spp. 62 . Surprisingly, we found that R. allomycis is triploid, which is rare in fungi (for example, in Epichloë 63 ) and usually reflects a sexually sterile condition. This observation of greater than diploidy has precedent in the microsporidia (for example, the tetraploid Nosema 29 ), and the chytrid B. dendrobatidis has diploid to tetraploid nuclei 64,65 . Four of the six basal fungi are not haploid, including three of four Zoopagomycota. The preponderance of heterozygous species at the base of the fungal tree implies a diploid (or higher ploidy) ancestor of the fungi. Beyond estimating ploidy, the single-cell methods presented here could be applied broadly to other fungal taxa where culturing has been unsuccessful (for example, certain rust fungi and ectomycorrhizae), to facilitate genetic mapping, test for genetic segregation and establish mating systems using cohorts of meiotically produced spores.
The uncultured fungi in this study have diverse phylogenetic backgrounds and nutritional strategies. Most are parasites, presumably dependent on hosts for certain nutrients, which may explain their unculturability. By mapping the proteome to primary metabolism pathways, we observed major deficiencies in the ability to synthesize the full set of amino acids and polyamines, which probably result in auxotrophies. Overall, there is a mosaic of losses of essential metabolism genes across the uncultured taxa, consistent with patterns seen in the parasitic Cryptomycota 66 . The deficiencies in the biotin and thiamine pathways in mycoparasites are intriguing given that Syncephalis species have been successfully grown using beef liver 67,68 , which contains measurable quantities of both biotin and thiamine 69 . However, Syncephalis do not sporulate nor complete their life cycle under these conditions 68 . Most chytrid pollen saprotrophs are culturable, yet the unculturability of B. helicus may be explained, as it shares deficiencies in spermidine and sulfur metabolism. M. bicuspidata from Daphnia is similar to strain NRRL YB-4993 from brine shrimp (although they are not conspecific 70 ) and lacks 15 enzymes that are broadly linked to aspects of urea, sulfate and thiamine metabolism. Given that the majority are also missing from various other target species, this may help to explain its unculturability. R. allomycis is unlikely to be axenically cultured because it is missing a large number of genes for critical pathways, losses that potentially cannot be corrected for as seen in microsporidia 71 .
Mycoparasitic fungi face unique challenges because they must parasitize a fungal host using fungal enzymes without disrupting their own cells. However, the specificity of this antagonism is not fully characterized 72 . Genomic analysis of the facultative mycoparasite Trichoderma spp. revealed several key features, including gene expansions and numerous antifungal-producing secondary metabolite gene clusters 73 . Nevertheless, many fungal lineages contain mycoparasites and there is no generalized analysis of the traits that are selected for in this ecological transition. We did not find any orthologous genes that were a signal of mycoparasitism across all taxa. In obligate mycoparasites, such as those studied here, we found a CAZyme-to-protease ratio more similar to that of animal pathogens rather than plant pathogens. An expansion of subtilases was observed among mycoparasitic Zoopagomycota, forming a group that was distinct from known subtilase families and among subtilases from other fungi. The GH19 chitinases are exclusively restricted to EDF lineages, which suggests that this gene family is ancestral but has been lost.
Having genomes from these EDF enabled us to uncover traits that were previously only described in Dikarya. We observed a vast expansion of secondary metabolism NRPS genes unique to D. cristalligena among other EDF. Similarly, we provide evidence for multiple hydrophobins in C. protostelioides, the only member of the Chytridiomycota known to form aerial hyphal-like sporeproducing structures.
Here, we show that single-cell methods can dramatically expand our understanding of fungal biology. Given the large numbers of as-yet-unsequenced EDF genomes, many of which have never been cultured, these methods help to clarify the expectations for future assembly and functional annotation of genomes from those lineages. Furthermore, given the large numbers of uncultured Dikarya, these methods can be effectively applied across the phylogeny. As the scale of single-cell genomics increases 13 , we will be able to generate a full picture of the fungal tree of life that precisely reveals the complete extent of fungal diversity, regardless of the culturability of the individual taxa.

Methods
Strains and sample preparation. A dual culture of parasite R. allomycis CSF55 with its host fungus Allomyces sp. was established and used previously to sequence the genome of the parasite 74 . For this study, spore suspensions of R. allomycis were obtained under optimal conditions by washing the plates with a dilute Tween solution. An estimated 10 6 -10 7 spores of the parasite with up to 5% of host spores were obtained. The sample was preserved in 10% sterile glycerol solution, shipped on dry ice and stored at − 80 °C. A dual culture of C. protostelioides ATCC 52028 with its host Sordaria was used to isolate parasitic zoospores at 2.5 × 10 6 per ml. The zoospore suspension was preserved in 10% dimethylsulfoxide with 10% FBS, shipped on dry ice and stored at − 80 °C.
B. helicus was grown to a high density of cells through enrichment methods using spruce pollen in bog water. The sample was obtained from Perch Pond Fen near Old Town, Penobscot County, Maine, USA, in June 2014. This enrichment culture was filtered through a 40-μ m mesh (removing pollen and sporangia) and concentrated by centrifugation to about 5 × 10 4 zoospores per ml. The sample was preserved in 10% glycerol, shipped on ice and stored at − 80 °C.
M. bicuspidata was isolated from an infected population of the water flea Daphnia dentifera grown under laboratory conditions. D. dentifera samples were rinsed repeatedly with deionized water. Then, insect pins were used to puncture the carapace and a micropipette was used to collect haemolymph, which contained a mixture of M. bicuspidata yeast cells and ascospores. Cells were preserved in 10% glycerol at a concentration of 10 5 spores per ml and stored at − 80 °C.
T. sphaerospora RSA 1356 was grown on V8 juice agar in dual culture with C. recurvatus and harvested from petri plates. The sample was stored in 50% glycerol at − 80 °C. P. cylindrospora RSA 2659 was cultivated on potato dextrose agar with the host Umbelopsis isabellina. The culture was grown on many petri dishes and the spores of both the fungus and the host were removed from the culture by washing the plates with 0.2% Triton X-100. An estimated 2.5 × 10 7 spores per ml of parasite with host were obtained and preserved in 10% glycerol at − 80 °C.
All mycoparasites described above are considered obligate mycoparasites. For benchmarking purposes, we used non-single-cell approaches to sequence material gathered from genomic DNA extracted from enrichment cultures of C. protostelioides. We also took advantage of genomic resources from M. bicuspidata NRRL YB-4993, a parasite of brine shrimp 70 related to but not considered conspecific with the Daphnia parasitic M. bicuspidata targeted in this study, and the published genome of R. allomycis 74 . These three genomic DNA isolates were used to compare with their respective single-cell genomes.
Single-cell genomics pipeline. The pipeline schema is shown in Supplementary  Fig. 1. After sample collection, as described in the preceding section, individual cells were isolated using a one-step or two-step FACS (BD Influx Cell sorter) procedure. Target population enrichment in the original sample, determined by microscopy, dictated the use of one-step or two-step FACS. Cells were sorted into 384-well plates at counts of 1, 10, 30, 50 or 100 cells per well. Sorting accuracy was verified microscopically (Zeiss Axio Observer D1) and according to morphology.
As a quality-control step, amplified cells were 18S rDNA screened using the following primers: M13CRYPTO2-2F (5′ -GTTTTCCCAGTCACGACCACAGG GAGGTAGTGACAG-3′ ), M13AU4v2 (5′ -CAGGAAAAGCTATGACGCCTC ACTAAGCCATTC-3′ ), M13DPD360FE (5′ -GTTTTCCCAGTCACGACCGGA GARGGMGCMTGAGA-3′ ), M13DPD1492RE (5′ -CAGGAAACAGCTATGACA CCTTGTTACGRCTT-3′ ), M13SR1RFor (5′ -GTTTTCCCAGTCACGACTA CCTGGTTGATYCTGCCAGT-3′ ) and M13NS4Rev (5′ -CAGGAAACAGCT ATGACCTTCCGTCAATTCCTTTAAG-3′ ). Amplified 18S rDNA fragments were treated with ExoSap-IT (Thermo Fisher Scientific) and Sanger sequenced. Sequences were queried using BLAST against the National Center for Biotechnology Information non-redundant (NCBI nr) and Assembling the Fungal Tree Of Life (AFTOL) 76 databases. Verified target wells were prepared using the TruSeq library method. Between three and seven libraries were pooled together and sequenced on one lane of Illumina HiSeq 1T v4-v6 2 × 150 bp. All pooled libraries belonged to the same species. Prior to establishing this protocol, we tested for crosstalk between TruSeq libraries on the MiSeq platform during pre-assembly library screening and found no such occurrence for the fungal single-cell pipeline. However, owing to previous observations of library crosstalk in the Joint Genome Institute (JGI) bacterial single-cell pipeline, the fungal single-cell pipeline decided on the cautionary same-species pooling until further examination was completed. Fungal single-cell pipeline TruSeq libraries were prepared on eight tube strips with individual caps using automated multichannel pipettes, which excludes plate-based robotic potential cross-contamination between samples.
In the fungal single-cell pipeline, there are several steps of rigorous filtering of potential contaminating sequences, at both the read level (pre-assembly) and the contig level (post-assembly). All Illumina reads were run through a filtering pipeline prior to assembly. BBDuk (https://sourceforge.net/projects/bbmap/) was used with the settings filterk = 27 and trimk = 27 to remove Illumina adapters, known Illumina artefacts, phiX and quality trim both ends to Q12. Resulting reads containing more than one 'N' , with quality scores (before trimming) averaging < 8 over the read or length under 40 bp after trimming were discarded. Remaining reads were mapped to a masked version of human hg19 with BBMap, discarding all hits over 93% identity. Reads mapping to mouse were also similarly discarded. Reads were also normalized to make the coverage more uniform, including removing reads that have very high coverage (> 100). Individual libraries from each species were assembled with SPAdes 77 (v2.6 or v3.0), using the 'careful' and 'singlecell' options, and k-mer sizes of 21, 33 or 55. After assembly, scaffolds of < 2 kb in length were removed as, in our experience, they are phylogenetically ambiguous. Assembled contigs were also compared using BLAST to a set of contaminant databases: bacterial, mouse, human, feline and canine. Finally, we performed tetramer principle component analysis and removed any outlier contigs. This extra cautious read usage for assembly probably leads to reduced completeness and removes symbiotic occurrences but guarantees one species genome instead.
Supplementary Table 1  Annotation completeness. CEGMA 26 as a general measure of completeness has been deprecated in favour of BUSCO 79 . However, because the underlying data for BUSCO relies heavily on Dikarya fungi, we have observed ( Supplementary  Fig. 9) that BUSCO dramatically underestimates the coverage of EDF lineages, particularly within the Chytridiomycota, Blastocladiomycota and Mucoromycota. As such, we continue to use CEGMA metrics for this study as it focuses primarily on EDF.
Heterozygous polymorphism discovery. Paired-end reads were aligned to assembled draft genomes using BWA 80 (v0.7.12-r1044) using default parameters. Variants were identified using FreeBayes 81 (v1.0.2-58-g054b257) with parameters '-pooled-continuous -min-coverage 5' . As non-haploid organisms have multiple genome copies per cell, increased counts of higher-frequency SNPs (allele frequency > 0.25) in their single cells are indicative of heterozygosity. Only SNPs with an allele frequency > 25% were kept for further analysis. Comparisons of ploidy levels between species were performed on four single-cell libraries with the highest completeness and quality. K-mer graphs were plotted using the kmercountexact.sh script of the BBTools package (http://jgi.doe.gov/data-andtools/bbtools/).
For phylogenetic analysis of identified NRPS genes, adenylation domains were mined from proteomes using HMMER 86 (v3.1b2) with a hidden Markov model based on the fungal and bacterial domains identified by Bushley and Turgeon 53 . MUSCLE 87 (v3.8.31) was used to align extracted domains with those used for the hidden Markov model creation and gaps were removed manually. RAxML 84 (v8.2.8) was used for phylogenetic reconstruction using the gamma model of rate heterogeneity and the RTREV substitution matrix with 100 bootstrap replicates.
Phylogenetic analysis. All-v-all blastp 88 (v2.2.26) was run using a cut-off of 1 × 10 −5 on the predicted proteomes of the representative set of fungal and outgroup species provided in Supplementary Table 2. Clusters were predicted using MCL 89 (v1.008) with an inflation value of 2. A python script identified clusters containing at most one gene copy per genome, allowing for up to 8 missing taxa per cluster, which resulted in 805 total clusters selected. Each cluster was aligned using MAFFT 90 (v7.047) on each cluster, with a multiple sequence alignment algorithm detected automatically (using the -auto flag). Alignments were trimmed with Gblocks 91 (v0.91b) and concatenated, followed by rapid bootstrap maximum likelihood tree building with RAxML 84 , 27,751 distinct alignment patterns, 100 bootstrap replicates and using the GAMMA + WAGF protein model. For all 29,255 positions in the alignment, the target genomes had a median of 14% missing, compared to the median of 2.5% missing of others.
Orthogroup reconstruction. Using MCL clustering 89 , orthologues were collected using a minimum of three taxa per orthogroup, ignoring single-cell C. protostelioides and R. allomycis in favour of the respective co-assemblies. The R package APE 92 was used for ancestral-state reconstruction. Protein functions predicted using PRIAM 93 were mapped to each internal node of the phylogenetic tree presented in Fig. 1. Gains and losses in single-cell lineages were considered as relative to the ancestral fungal node (Supplementary Fig. 3): D. cristalligena shows noticeable gains in all high-level metabolic categories, whereas B. helicus has mainly losses.
Metabolic reconstruction. Enzyme classifications were predicted for all 60 proteomes in Supplementary Table 2 using PRIAM 93 . A subset of 24 fungal, nonsingle-cell proteomes were defined as 'free living' based on culturability. The presence/absence counts were converted to a binary presence/absence matrix for each enzyme classification found in each species and filtered according to 'found in free-living' threshold of present in 18 out of 24 (75%) free-living species and absent in at least 5 single-cell species. 'Single cell' includes C. protostelioides and R. allomycis culture isolate genomes in place of those respective single-cell genomes. The resulting pattern of conserved losses among 'core' metabolic pathways was summarized and displayed using IPath 94 .
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The co-assembled genomes and annotations of the target species are available through MycoCosm (https://genome.jgi.doe.gov/fungi) and Genbank using the following MycoCosm URLs and NCBI accessions, respectively: R.

Statistical parameters
When statistical analyses are reported, confirm that the following items are present in the relevant location (e.g. figure legend, table legend, main text, or Methods section).

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement An indication of whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistics including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers upon request. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability All studies must disclose on these points even when the disclosure is negative.

Sample size
The fungal and eukaryotic taxa used for comparative analyses were selected to best represent the breadth of fungal diversity. In general, sample sizes represent the extent of available data.
Data exclusions Exclusions were limited to institutional standard quality control filtering of sequencing data.

Replication
For single-cell sequencing, multiple genome libraries were constructed for each species considered. Multiple individual libraries were constructed from singular individual cells. Additionally, multiple individual libraries were also constructed from pools of individual cells. Cells from a given species were isolated from the same culture.
Randomization Group allocation was not a part of this study, therefore there was no randomization.

Blinding
Group allocation was not a part of this study, therefore there was no blinding.
Reporting for specific materials, systems and methods