Mycoavidus sp. Strain B2-EB: Comparative Genomics Reveals Minimal Genomic Features Required by a Cultivable Burkholderiaceae-Related Endofungal Bacterium

The intracellular features of BRE dwelling in Mortierella parvispora E1425 were characterized by transmission electron microscopy (TEM) (Fig. 1). Bacterium-like cells were observed inside the host-fungal hyphae; they were rod shaped and 0.3 to 0.5 by 0.8 to 1.4 μm in size, with a smooth double-layered cell envelope and a nucleoid surrounded by numerous ribosomes in the cell center. As in previous reports (15, 22), the endobacterial cells aggregated in close proximity to lipid bodies (Fig. 1A and B). In addition, our microscopic observations showed that endobacterial cells were also located near other organelles such as vacuoles and mitochondria, even directly in contact with them (Fig. 1C to E). In contrast to other Mycoavidus-related endofungal bacteria (4, 15), the endobacteria in M. parvispora were not enveloped in fungal vacuoles according to our observations. Together, these morphological features extend the understanding of intimate interactions between the BRE and their host fungus at the subcellular level.

To date, only one endofungal bacterium, M. cysteinexigens B1-EB^T, hosted by fungi in the genus Mortierella has been successfully isolated in pure cultures (20). The BRE present in M. parvispora E1425 was identified as the closest relative of strain B1-EB^T, based on phylogenetic analyses using an almost-full-length 16S rRNA gene (4). In order to isolate the BRE, aliquots from filtered fungal homogenates of M. parvispora E1425 were incubated on B-CYEα agar plates for 30 days at 23°C. Because the bacterial colonies appearing on the plate were homogeneous, a single colony from the plate was subcultured repeatedly on new B-CYEα agar plates to assess its purity. A pure culture of a Gram-negative rod-shaped bacterium was obtained and designated strain B2-EB (Fig. 2A). Although this strain shared 99.0% nucleotide identity in the 16S rRNA gene sequence with strain B1-EB^T, its colonies on B-CYEα agar plates showed a color (translucent to white) different from that of strain B1-EB^T (white to cream) (Fig. 2B). This result implied that a novel cultivable BRE strain had been successfully isolated in pure cultures.

Genomic DNA of strain B2-EB was separately sequenced on Illumina HiSeq2500 and PacBio RS II platforms, and 1.48 Gb of paired-end reads and 0.93 Gb of long reads were obtained, respectively. Hybrid assembly and read correction using SPAdes (version 3.10.1) and Pilon (version 1.22) resulted in a 1,876,900-bp circular chromosome with 1,117-fold genome coverage and 48.9 mol% GC content, which contained 2 rRNA operons (rrn), 42 tRNA loci for 20 amino acids, and 1,627 protein-coding sequences (CDSs) with an average size of 1,027 bp (Table 1). Of these CDSs, 1,319 were identified as single-copy CDSs accounting for 81% of total CDSs. In addition, 1,157 and 1,069 CDSs were assigned to COG and KEGG orthologs, respectively (Table 1); 215 CDSs matched hypothetical proteins of unknown function, and 59 of them were less than 100 amino acid residues in size (data not shown). Compared with other BRE genomes, the B2-EB genome featured the lowest GC content, the highest proportion of single-copy genes, and the smallest genome size among the known cultured BRE.

A reconstructed maximum-likelihood (ML) tree based on concatenated amino acid sequences of 52 single-copy housekeeping genes conserved in the Burkholderiaceae genomes showed that strain B2-EB was phylogenetically closest to the cluster of the species M. cysteinexigens, forming an independent clade neighboring the “Ca. Glomeribacter gigasporarum” clade (Fig. 3). This result was consistent with previous analyses based on the 16S rRNA gene sequence (4, 22). In addition, the whole-genome average nucleotide identity (ANI) between strains B2-EB and B1-EB^T (88.6%) was lower than the ANI cutoff value of same species (94 to 95%) (33), suggesting that strain B2-EB can represent an undescribed species in the genus Mycoavidus. Moreover, our phylogenetic analysis using 49 representative complete genomes in the family Burkholderiaceae showed that Mycetohabitans endosymbionts hosted by Rhizopus microsporus were closer to bacteria with large and multipartite genomes in the genera Caballeronia, Paraburkholderia, and Burkholderia than to Mycoavidus and “Ca. Glomeribacter gigasporarum” endosymbionts in the Mortierella spp. and AMF, respectively (Fig. 3). These results supported that Mycetohabitans endosymbionts have an independent origin, sharing their most recent common ancestor with the free-living bacteria in these genera, whereas Mycoavidus and “Ca. Glomeribacter gigasporarum” endosymbionts seemed to originate from another common ancestor (16). From the cultivability, translucent colony color, and ANI analysis results, strain B2-EB could represent an undescribed species in the genus Mycoavidus.

Although the genomic features of Mycoavidus sp. strain B2-EB, including GC content, number of rrn and number of tRNA loci, were close to those of M. cysteinexigens B1-EB^T, the other features of genome size, mobile genetic elements, and proportion of single-copy genes were strikingly different between these two cultured endofungal bacteria (Table 1). Compared to the two complete genomes of M. cysteinexigens B1-EB^T and AG77, strain B2-EB possessed an approximately ∼1.0-Mb-smaller genome with a paucity of transposable elements and a large proportion of single-copy genes (Fig. 4A; see Fig. S1 in the supplemental material). However, phylogenetic orthology analysis using OrthoFinder showed that the B2-EB genome shared 91.3% of the orthologous groups (OGs) with the two M. cysteinexigens genomes, which included 1,395 core OGs plus 13 and 33 OGs shared with the B1-EB^T and AG77 genomes, respectively (Fig. 4B). These results indicated that although the genome size of strain B2-EB was reduced, its genetic elements were highly homologous with those of M. cysteinexigens.

Further assessment of genomic variation on a population scale showed that the ratio of nonsynonymous to synonymous sites (d_N/d_S) for individual genes between the B2-EB and M. cysteinexigens genomes displayed a leptokurtic distribution clustered around the mean value regardless of the d_N/d_S estimation method (Fig. 4C). This pattern was similar to those found in nonessential heritable endosymbionts such as Wolbachia, Hamiltonella/Regiella, and “Ca. Glomeribacter,” as well as those of some free-living bacterial lineages such as Burkholderia, Prochlorococcus, Bradyrhizobium, and Bifidobacterium (10). Additionally, the mean d_N/d_S ratio between B2-EB and M. cysteinexigens estimated by three methods was in the range 0.118 to 0.152, suggesting strong purifying selection and large population size like those detected in some microbial endosymbionts (10, 34). This tendency was even stronger than in the case of the “Ca. Glomeribacter gigasporarum” endosymbionts (d_N/d_S, 0.149 to 0.198) estimated in the previous study (10). Therefore, Mycoavidus sp. B2-EB represents a model of nondegenerative genome evolution similar to that for “Ca. Glomeribacter gigasporarum” endosymbionts. This also created a debate on whether the B2-EB genome is reduced or the B1-EB^T and AG77 genomes are extended (35).

A syntenic comparison with M. cysteinexigens B1-EB^T showed that 936 CDSs related mainly to a region for prophages and transposons and accounting for 40.4% of total CDSs in the B1-EB genome were missing in the B2-EB genome (Fig. 5A). Even so, some gene loci within the translatable regions of the B1-EB^T and AG77 genomes were still found in the B2-EB genome (Fig. 5A and B).These residual genes indicate that genome reduction was likely to have resulted in the relatively small genome of strain B2-EB. Therefore, we proposed that the genome reduction occurring in strain B2-EB may be mediated by the repeated insertion-deletion of prophages and transposons, as previously reported (31, 32, 36). The missing CDSs affiliated with the orthologous genes conserved in the genomes of Burkholderiaceae-related bacteria seemed to be reduced from the B2-EB genome over evolution and were involved mainly in the secretion system, amino acid and glycolytic metabolism, and transcription (see Fig. S2 in the supplemental material). Notably, the other specific CDSs in the B1-EB^T genome might be acquired through horizonal gene transfers, such as genes encoding insecticidal toxins close to orthologs possessed by the genus Xenorhabdus (Gammaproteobacteria) living symbiotically with soil entomopathogenic nematodes (37). However, the process of this gene gain event is still unclear.

On the other hand, 1,490 out of 1,627 CDSs in total in the B2-EB genome were well conserved in the B1-EB genome (Fig. 5A), and 60.5% of these commonly conserved CDSs shared more than 90% amino acid identity (AAI) between these two strains (see Fig. S3 in the supplemental material). Orthologous genes sharing low AAI values of <80% accounted for 21.3% of the total common orthologous genes, which encoded mainly proteins involved in RNA processing and modification, transcription, and energy production and conversion. In the case of the obligate endosymbiont Wolbachia associated with arthropod hosts, the specific genes for involvement in host interaction underwent elevated substitution rates that may be due in part to adaptations by Wolbachia to a new host environment (38). It is tempting to speculate that the genes with elevated mutation rates may represent specific gene elements in the two Mycoavidus genomes for adapting to the different intracellular environments of the two Mortierella species. Furthermore, a high degree of genomic rearrangement was identified in the 0.80- to 1.13-Mb downstream region of dnaA, encoding the replication initiator, showing “symmetric inversions” with endpoints that were equally distant from the origin of chromosomal replication (Fig. 5A and B). Genomic rearrangement likely impacts gene expression and results in the loss of gene function when a breakpoint occurs inside a reading frame (39, 40), which not only affects the organismal phenotype but also contributes to phylotype subdivision (41). It has been reported that extensive genome rearrangements probably triggered or enhanced the Wolbachia speciation event, because it could suppress homologous recombination in chromosomal regions that are not colinear (42). Here, we showed that genomic arrangements occurred between the genomes of Mycoavidus cysteinexigens and Mycoavidus sp. B2-EB (the two closest species to one another), suggesting that they may play an important role in the Mycoavidus speciation process.

Additionally, in spite of the high homology with the previously determined M. cysteinexigens genomes, the B2-EB genome also possessed 137 unique CDSs (Fig. 4B), including 89 and 48 CDSs encoding hypothetical and functional proteins, respectively. Most of the CDSs encoding functional proteins were involved in interactions with the host, mainly including 12 CDSs encoding lipopolysaccharide biosynthesis and secretion enzymes, 4 CDSs encoding secondary metabolite biosynthesis enzymes, and 4 CDSs encoding toxins or antitoxins (see Table S1 in the supplemental material). Such genes have been reported for the Burkholderiaceae-related endofungal bacteria in previous studies (22, 30, 43). The other functional CDSs specific to strain B2-EB were related to DNA replication and repair, nutrient metabolism, cell elongation and division, ribosome biogenesis, membrane transport, and signal transduction (Table S1). Notably, two of the unique CDSs (locus tags MPB2EB_0909 and MPB2EB_1044) were likely to be associated with bacterial survival inside the fungal cell, annotated as putative genes encoding β-lactamase and E3 ubiquitin-protein ligase, respectively. β-Lactamase can provide resistance to β-lactam antibiotics produced by fungi, such as penicillins (44), whereas the E3 ligases of some pathogenic bacteria mimic the activity of eukaryotic E3 ubiquitin ligases to benefit themselves (45). However, the roles of the two genes in the endofungal lifestyle should be examined using molecular biological approaches in the future.

The smallest genomes possessed by Buchnera aphidicola, an obligate intracellular endosymbiont in aphids, are shaped by the evolutionary force of genetic drift, in which slightly deleterious mutations fix rapidly, causing the loss of gene functions, even including DNA repair, which further accelerates the accumulation of mutations (8). Notably, the proportion of transposable element accounted for only 2.4% of the total genome length in the B2-EB genome, versus 7.2 to 11.5% in known Mycoavidus genomes. In contrast, the genome reduction of heritable endofungal bacteria is caused by limited recombination in a largely clonal population, which is well supported by genomic information on “Candidatus Glomeribacter” (10). In this study, we estimated the genome completeness across the family Burkholderiaceae using the taxonomic-specific workflow of CheckM (46). Surprisingly, the results indicated that Mycoavidus endofungal bacteria dwelling in Mortierella possessed 88.6% of the checked single-copy genes (Fig. 6A), which was obviously more than those conserved in “Ca. Glomeribacter gigasporarum” genomes (21.8 to 70.2%). Moreover, the full gene set of DNA repair mechanisms, including base excision repair, nucleotide excision repair, mismatch repair, and homologous recombination, were well conserved in the genomes of Mycoavidus strains B1-EB^T and B2-EB and may be responsible for purging slightly deleterious mutations from the genomes. These findings suggested that nondegenerative genome reduction occurred during the evolution of Mycoavidus-related endofungal bacteria adapting to the intracellular environments of fungal hosts.

The second surprising finding was that strain B2-EB represented the most streamlined genomic architecture in the Burkholderiaceae bacteria with sequenced genomes, even beyond the streamlined genomes of Polynucleobacter species (47–49). The B2-EB genome was 1,877 kb with 1,627 CDSs and featured the highest proportion of single-copy genes in the genome (Fig. 6B). Currently, the free-living bacterium with the smallest genome (1,304 kb) is an uncultured but abundant ocean methylotroph of the Betaproteobacteria, HTCC2181, with 1,377 CDSs (50), followed by the hyperthermophilic bacterium Aquifex aeolicus with a slightly larger genome (1,591 kb and 1,613 genes) (51). The smallest photosynthetic cell, represented by the oceanic cyanobacterium Prochlorococcus marinus, has a 1,660-kb genome containing 1,765 genes (52). In the family Burkholderiaceae, the smallest genome is represented by Polynucleobacter necessarius endosymbionts isolated from ciliates living in aquatic ecosystems (53), with 1,560 kb and 1,279 CDSs (47). However, the single-copy genes in this genome accounted for 64.3% of the total genes (Fig. 6B), which was lower than the proportion in the B2-EB genome. Strain B2-EB has a small genome comparable to the minimal genomes of the above-mentioned free-living bacteria and possesses a radically streamlined structure. A reduction in genome size has previously been identified in the M. cysteinexigens genomes (22, 30), and our analysis here showed that the proportions of single-copy genes were 61.9% and 68.5% in the B1-EB^T and AG77 genomes, respectively (Fig. 6B). Further comparisons with their free-living relatives (proportion of single-copy genes, 48.4% to 52.8%) (Fig. 6B) suggested that all three of the Mycoavidus genomes were streamlined, but the most streamlined was the B2-EB genome. Therefore, the genetic elements conserved in the B2-EB genome may represent the minimal genomic features required by a cultivable endofungal bacterium.

Compared with those of close relatives in another fungal hosts and free-living relatives in various habitats, the genomic features of Mycoavidus were characterized by the harboring of more abundant genes responsible for functions including intracellular trafficking, secretion, vesicular transport, RNA processing and modification, and recombination (Fig. 7). Furthermore, the streamlined genome of strain B2-EB harbored almost complete sets responsible for encoding the 30S and 50S ribosomal proteins, tRNA aminoacylation, rRNA and tRNA modifications, DNA replication and repair, transcription, translation, and protein folding (see Tables S2 to S9 in the supplemental material). Although some genes involved in amino acid metabolism (e.g., asnA, asnB, gltB, and hutH) and secretion systems (e.g., yscF, virB5, gspC, gspI, gspL, gspM, vgrG, hcp, impL, impK, and vasG) were absent in the B2-EB genome (see Fig. S4 in the supplemental material), it was similar to the genomes of other BRE that possessed many of these genes (6, 22, 29, 30), suggesting that they have an important role in the endofungal lifestyle. In addition, a large number of genes involved in fatty acid metabolisms were present in the B2-EB genome, indicating that fatty acids are likely to be a major carbon source for Mycoavidus endofungal bacteria (Fig. S4). Together with other essential genes for life, such as those involved in tricarboxylic acid (TCA) cycling, energy generation, and various transport systems, the genetic elements harbored by strain B2-EB showed the minimal genomic features required by both an endofungal lifestyle and artificial culture.

BRE have been detected in several fungal lineages, including Gigaspora, Racocetra, and Cetraspora (Glomeromycotina), Rhizopus (Mucoromycotina), and Mortierella (Mortierellomycotina) (10, 54). The sequenced BRE genomes have revealed common features such as reduced genome size, loss of genes encoding primary metabolism, and specialization in fungal metabolite uptake. In particular, the BRE affiliated with the Mycoavidus-Glomeribacter clade were elucidated as having amino acid auxotrophies that have limited their artificial cultivation in pure cultures. Our previous efforts in reverse genomics have contributed to the isolation of Mycoavidus cysteinexigens from the host fungus Mortierella elongata by supplying exogenous cysteine to substitute for the genetic defect in bacterial cysteine biosynthesis (19, 26).

In this study, we successfully isolated an endobacterium phylogenetically close to M. cysteinexigens, named Mycoavidus sp. B2-EB, from the host fungus Mortierella parvispora. Moreover, the hybrid assembly of deep-sequencing data sets reconstructed a complete genome sequence of strain B2-EB, revealing a reduced genome size (1.88 Mb), which was less than those of M. cysteinexigens strains (2.64 to 2.79 Mb) and comparable to those of “Ca. Glomeribacter gigasporarum” genomes (1.34 to 2.36 Mb). Further assessment of genomic variation on a population scale indicated that Mycoavidus sp. B2-EB possessed a nondegenerative genome comparable to that found in “Ca. Glomeribacter gigasporarum” endosymbionts. In addition, comparative genomic analysis with other BRE and free-living relatives affiliated with the family Burkholderiaceae showed that the B2-EB genome had a relatively high genome completeness and high proportion of single-copy genes but only a few transposable genetic elements. These results suggested that the B2-EB genome reached a state of relative evolutionary stability. The genetic elements conserved in the stabilized genome appear to have resulted from the adaption of endobacteria to the host metabolism. In our previous study (30) and this study, we found no genes involved in invasion or lysis of fungal cells in the Mycoavidus genome according to comparative genomic analysis between Mycoavidus and Mycetohabitans, indicating that strain B2-EB has already lost a way to get in and out of fungal cells. To the best of our knowledge, the B2-EB genome represents the smallest genome size among cultivable endofungal bacteria. Our findings are expected to fill a knowledge gap in the understanding of genome evolution of Mycoavidus-related endofungal bacteria, which results in a minimal genome of a cultivable endosymbiont.

The information gleaned from this genome revealed the minimal genome features required by both an endofungal lifestyle and artificial cultivation, which furthers our understanding of genome reduction in fungal endosymbionts and extends the culture resources for the biotechnological development of engineering synthetic microbiomes.

Transmission electron microscopy (TEM) with sample preparation using rapid freezing and a freeze-substitution method was performed to observe the endohyphal features in the fungus Mortierella parvispora E1425 (JCM39028) at Tokai Electron Microscopy Inc., Nagoya, Japan. In brief, 7-day-old mycelia grown on sterilized cellophane sheets placed onto glycerol agar (containing, per liter of distilled water, 10 g glycerol, 0.5 g Casamino Acids, 2 g NaNO₃, 1 g KH₂PO₄, 0.5 g KCl, 0.5 g MgSO₄·7H₂O, 0.01 g FeSO₄, and 15 g agar) at 23°C were sandwiched with copper disks and frozen in liquid propane at −175°C, followed by substitution with 2% osmium tetroxide in acetone and 2% distilled water at −80°C for 48 h, −20°C for 4 h, and 4°C for 2 h, dehydrating through anhydrous acetone twice for 30 min each and 3 changes of 100% ethanol for 30 min each at room temperature (RT), infiltrating with propylene oxide (PO) twice for 30 min each and a PO-resin mixture (7:3) for 1 h, and polymerizing in fresh 100% resin at 60°C for 48 h. The polymerized sample was ultrathin sectioned at 90 nm using an Ultracut UCT ultramicrotome (Leica, Vienna, Austria) and stained with 2% uranyl acetate at RT for 15 min and with lead stain solution (Sigma-Aldrich, Tokyo, Japan) at RT for 3 min. The grids were observed using a JEM-1400Plus transmission electron microscope with an EM-14830RUBY2 charge-coupled device (CCD) camera (JEOL, Tokyo, Japan) at an acceleration voltage of 100 kV.

Endohyphal bacteria living in the fungus M. parvispora E1425 were isolated using the method reported previously (20), with a minor modification. In brief, fungal strain E1425 was cultivated for 7 days at 23°C on half-strength CMMY agar. The cultivated mycelia were homogenized using a sterilized pestle in a 1.5-ml microcentrifuge tube at maximum rotation speed for 2 min and centrifuged at 2,000 × g for 10 min. The supernatant was filtered through 8-μm- and 3-μm-pore-size membrane filters to remove fragmented hyphae and sporangiospores. One-milliliter aliquots of the filtered suspension were spread on a B-CYEα agar plate [containing, per liter of distilled water, 10 g yeast extract, 2 g charcoal powder, 0.4 g l-cysteine hydrochloride, 0.25 g ferric pyrophosphate (soluble), 10 g N-(2-acetamido)-2-aminoethanesulfonic acid (ACES), 1 g potassium 2-oxoglutarate, and 15 g agar, pH 6.9] (Eiken Chemical, Tokyo, Japan) and incubated for 30 days at 23°C. When the bacterial colonies grew on the plate, a single colony from the plate was subcultured repeatedly on a B-CYEα plate in order to assess its purity.

Genomic DNA of endobacterial isolates grown on a cysteine-containing buffered charcoal-yeast extract agar (B-CYEα) plate was extracted using a modified lysozyme buffer method, as in our previous study (30). The obtained DNA (39.2 μg in total) was used to build a pair-end library with an insert size of ∼550 bp and a 20-kb SMRT library and sequenced on the Illumina HiSeq2500 and PacBio RS II platforms, respectively, at GeneBay, Yokohama, Japan. Sequence processing and hybrid assembly were performed in accordance with our previous strategy (55). Briefly, the single-end reads were trimmed out from the HiSeq paired-end reads using sickle v1.33 with default settings. The qualified paired-end reads and PacBio long reads were hybrid-assembled using SPAdes v3.10.1 with the pipeline option of –careful. Genome polishing and finishing were performed using pilon v1.22 and GenoFinisher v2.1 with the BWA-MEM v0.7.12 algorithm.

The genome sequence was basically annotated using the DDBJ Fast Annotation and Submission Tool (DFAST v1.1.0) pipeline (56), followed by determination of the metabolic pathways using the KEGG online services BlastKOALA v2.1 (57) and KEGG mapper v3.2, identifying the COG function category using EggNOG-mapper v2.0 (58), and searching the secondary metabolic gene clusters using antiSMASH v4.0 (59). Prophage sequences and transposons were identified using the PHASTER web server (60) and the previously reported method (30), respectively. The complete genome sequences of the type species for each genus in the family Burkholderiaceae were obtained from the GenBank database to perform the comparative genomic analysis. Orthologous clusters across the selected genomes were grouped using OrthoFinder v2.3.3 (61), and the COG profiles were characterized using the EggNOG-mapper. Whole-genome phylogenetic analysis was performed using ezTree v0.1 (62) with the complete genome sequences of type strains for the species in the family Burkholderiaceae. Genome completeness was estimated using CheckM v1.0.13 with the option of taxonomy_wf to calculate the conservation rate of the 568 single-copy housekeeping genes found in the genomes of free-living bacteria in the family Burkholderiaceae (46). To estimate the dN/dS ratio across the Mycoavidus genomes, the orthologous groups detected between the B2-EB and M. cysteinexigens genomes (B1-EB^T and AG77), three different methods were used, NG86 (63), LWL85m (64), and YN00 (65). The R package orthologr was used to perform the d_N/d_S calculation (66). In addition, these orthologous groups were sorted based on the physical locations and directions in the genomes and used to estimate the genome rearrangement by syntenic comparison.

Circos v0.69-6 (67) was used to generate a schematic spherical illustration of the complete B2-EB genome sequence and visualize the syntenic comparative analysis between the B2-EB and M. cysteinexigens genomes. The R package ggplot2 was used to generate a heat map showing the genome profiles of representative Burkholderiaceae genomes and violin plots of the distribution of the d_N/d_S values and the amino acid identity of individual genes in the Mycoavidus sp. B2-EB genome compared with those in the M. cysteinexigens B1-EB^T and AG77 genomes. The metabolic features of B2-EB were determined using the KEGG mapper with the option of “Reconstruct Pathway” (https://www.genome.jp/kegg/tool/map_pathway.html) and visualized using Microsoft PowerPoint.

The complete genome sequence of Mycoavidus sp. B2-EB has been deposited in DDBJ/ENA/GenBank under the accession number AP021872. The raw sequence reads for the Illumina paired-end and PacBio RS II libraries are available under accession number DRA009105.