This study utilized two strains, XY03792 and XY03829, of Rhizopus arrhizus. The XY03792, sourced from soy sauce in Malaysia, is a fermentative strain. The XY03829 is a wild strain obtained from animal dung in Pakistan. Both strains demonstrated the capability to ferment, resulting in the production of various compounds such as maltose, glucose, ethanol, lactic acid, fumaric acid, malic acid, glycerol, among others (Yao et al., 2018; Liu et al., 2022). These strains were preserved at Shandong Normal University under −20°C with 15% glycerine.

To prevent bacterial contamination, sporangiospores underwent a meticulous two-step surface sterilization process following the method outlined by Becard and Fortin (1988). In the initial step, sporangiospores were immersed in a 0.05% Tween 20 solution for two minutes, followed by a 10 min soak in a 2% chloramine T solution. Subsequently, they were thoroughly rinsed three times with sterile distilled water. This soak and wash procedure was repeated once more, after which the sporangiospores were preserved in a sterile solution containing 200 mg/L streptomycin and 100 mg/L gentamicin at 4°C. Moving to the second step, the stored sporangiospores underwent another round of soaking in a 2% chloramine T solution and were washed with sterile distilled water immediately before inoculation. The surface-sterilized sporangiospores, treated through this two-step process, were then cultured on potato dextrose agar (PDA: 200 g/L potato, 20 g/L glucose, 20 g/L agar, and 1,000 mL distilled water) at 30°C.

To avoid bacterial contamination during incubation, the mycelia cultivated on PDA for 5 days underwent surface sterilization using 30% hydrogen peroxide, following the procedure outlined by Izumi et al. (2006). Metagenomic DNAs of Rhizopus arrhizus and its EFBs were extracted using the GOMag Rapid Plant DNA Kit (GO-GPLF-400, GeneOn BioTech, China). According to the manufacturer’s instructions, approximately 30 mg of thalli were successively lysed, adsorbed, washed, and eluted for metagenome extraction. An empty centrifuge tube served as a negative control. The 16S rDNA was amplified using primers 27F (5′-AGA GTT TGA TCC TGG CTC AG-3′) and 1541R (5′-AAG GAG GTG ATC CAG CC-3′). The PCR mixture (25.0 μL) included 1.0 μL of template DNA (10.0 ng/μL), 1.0 μL of the two primers each (10.0 μM), 12.5 μL of 2 × Taq PCR Master Mix (Biomed Diagnostics Pte Ltd., Singapore), and 9.5 μL of sterile deionized water (Caporaso et al., 2012). PCR amplification involved an initial step at 94°C for 5 min, followed by 35 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 1 min, with a final extension at 72°C for 10 min. Sanger sequencing was carried out using the same primers (27F and 1541R) as used in PCR. Phylogenetic reconstruction employed the maximum likelihood (ML) and Bayesian inference (BI) methods through RAxML version 8.1.5 and MrBayes 3.2.7a, respectively (Ronquist et al., 2012; Stamatakis, 2014). Bootstrap supports (BS) for branches were obtained through 1,000 replicates (Estrada-De Los Santos et al., 2018). The resulting tree was edited online using the interactive Tree of Life platform (iTOL, https://itol.embl.de/itol.cgi; Letunic and Bork, 2019).

To isolate endosymbiotic bacteria, the host fungi were cultivated on PDA at 28°C for 3 days. A pellet of thalli, sterilized glass beads, and 1 mL lysogeny broth (LB: 10 g/L tryptone, 5 g/L yeast extract, 10 g/L NaCl) were added to a 2 mL centrifuge tube. The mixture was homogenized using a high-throughput tissue grinder (SCIENTZ-48) at 45 Hz for 30 s. The homogenized tissue fluid was then filtered through a 5 μM membrane and spread on lysogeny agar (LA: 10 g/L tryptone, 5 g/L yeast extract, 10 g/L NaCl, agar 30 g/L). Finally, it was incubated at 30°C for 7 days, and a single colony was transferred to another LA plate.

A single colony of EFBs grown on LA plates was inoculated into 15 mL of LB medium and shaken at 180 rpm at 37°C for 18 h. Bacterial cells were collected by centrifugation, and genomic DNAs were extracted using the Wizard Genomic DNA Purification Kit (A1120) following the manufacturer’s instructions. DNA integrity was verified on an agarose gel. Whole-genome resequencing was performed on the BENAGEN platform using Nanopore and Illumina NovaSeq PE150 platforms.

Raw data were assessed using FastQC 0.11.8 (Andrews, 2010) and Trimmomatic 0.39 (Bolger et al., 2014) for the filtration of low-quality reads, resulting in clean reads. The clean reads were assembled using MaSuRCA 3.4.3b (Zimin et al., 2013) and SPAdes 3.14 (Bankevich et al., 2012). Gene-coding models were predicted with Prokka (Seemann, 2014). For gene functional annotation, the predicted models were compared to various databases, including UniProt,¹ NR,² Pfam,³ KEGG,⁴ GO,⁵ CAZy,⁶ COG,⁷ CARD,⁸ and VFDB.⁹ Prophages, insertional sequences, and gene islands were predicted using PHASTER,¹⁰ ISFinder,¹¹ and IslandViewer¹², respectively.

Transposase and integrase sequences were sourced from the NCBI protein database, followed by clustering and classification using transposon and integron annotation databases, respectively. Initially, the Diamond blast+ software (version 0.9.31; Buchfink et al., 2015) was employed to compare genome and protein sequences against the Uniprot database. The outcomes of this comparison were integrated with the pre-constructed transposon and integron annotation databases to facilitate the prediction of transposons and integrons. Genomic collinearity analysis was conducted using the MAUVE (Darling et al., 2004). Subsequently, based on the results obtained from the Comprehensive Antibiotic Resistance Database (CARD) and the Virulence Factor Database (VFDB), a Venn diagram illustrating differences in gene numbers was generated using the Venny website.¹³

Two free-living Pandoraea sputorum strains, NCTC13161 (BioProject ID: PRJEB6403) and DSM21091 (BioProject ID: PRJNA262705), were downloaded from the NCBI database as references for P. sputorum EFB03792. The genome of EFB Mycetohabitans endofungorum HKI456 (BioProject ID: PRJNA370785) was used as the reference for M. endofungorum EFB03829.

The strains, preserved at −20°C with 15% glycerine, were inoculated on PDA plates supplemented with 100 μg/mL ampicillin, 50 μg/mL kanamycin, 10 μg/mL tetracycline, and 40 μg/mL ciprofloxacin. The plates were then incubated at 28°C for 36 h. Subcultures were performed under the same conditions for 30 generations. In each generation, the strain’s morphology was documented through photography and verified using live/dead staining. The cured fungal strain and free-living EFB isolated from the corresponding wild-type strain were simultaneously inoculated at the same position on lysogeny agar (LA) plates to observe whether the cured fungi resumed sporulation.

The 16S rDNA was successfully amplified from the metagenome of the two strains of the fungus Rhizopus arrhizus, resulting in a target fragment with a length of approximately 1.5 kb. These sequences were deposited in GenBank under the accession numbers OL413494 and OL413496. Identical 16S rDNA sequences were also annotated from the metagenome of the corresponding fungal strains. The maximum likelihood phylogenetic tree of EFBs based on 16S rDNA sequences is presented in Figure 1. In this phylogram, the two EFBs individually residing in the R. arrhizus strains XY03792 and XY03829 were grouped into the clades Pandoraea sputorum and Mycetohabitans endofungorum, respectively. Specifically, P. sputorum exhibited a close relationship with P. apista and P. norimbergensis, while M. endofungorum was closely related to M. rhizoxinica.

Following live/dead staining, EFBs with green fluorescence were observed within the hyphae, columellae, and sporangiospores of both fungal strains (Figures 2A–I). Fluorescence in situ hybridization (FISH) revealed red fluorescence in the hyphae of Rhizopus arrhizus XY03792 (Figures 2K–N), while green fluorescence was observed in the hyphae of R. arrhizus XY03829 (data not shown), confirming the presence of specific EFBs. Transmission electron microscopy (TEM) images showed clear transverse sections of bacterial rods (Figure 2J), indicating the localization of EFBs within the cytoplasm of fungal mycelia rather than in vacuoles. The bacteria within the cytosol were distinguishable from fungal organelles due to their visible cell walls, nucleoids, reserve materials, and morphology. The cross-section size of these cells ranged from 0.6 to 0.8 μm in diameter.

The genome of Pandorea sputorum EFB03792 comprised approximately 1.00 Gb of clean data from the Nanopore sequencing platform, encompassing 22,138 reads with an average sequence read length of 45,173 bp, achieving full coverage (100%) and an average sequencing depth of around 171×. The Illumina-filtered clean data amounted to approximately 1.19 Gb, encompassing 7,962,964 reads with an average sequence read length of 150 bp, achieving full coverage (100%) and an average sequencing depth of approximately 203×. The assembled genome contained one circular chromosome spanning 5,845,363 bp with a GC content of 62.63%. A total of 5,215 gene models were encoded, including 5,099 coding sequences (CDS), 75 transfer RNA (tRNA), 12 ribosomal RNA (rRNA), and two transfer-messenger RNA (tmRNA) (Figure 3A and Table 1).

The genome of Mycetohabitans endofungorum EFB03829 consisted of approximately 1.00 Gb of clean data from the Nanopore sequencing platform, encompassing 30,159 reads with an average sequence read length of 33,158 bp, achieving full coverage (100%) and an average sequencing depth of around 273×. The Illumina-filtered clean data amounted to 1.14 Gb, encompassing 7,598,064 reads with an average sequence read length of 150 bp, achieving full coverage (100%) and an average sequencing depth of approximately 303×. The assembled genome comprised one circular chromosome spanning 2,660,040 bp with a GC content of 61.28%. A total of 3,359 gene models were encoded, including 3,268 CDS, 48 tRNA, nine rRNA, and one tmRNA (Figure 3B and Table 1). Additionally, two plasmids were assembled, with Plasmid 1 measuring 800,149 bp long and exhibiting a GC content of 59.61%, and Plasmid 2 measuring 181,953 bp long with a GC content of 57.63%.

Among the 5,215 gene models of the strain EFB03792 of Pandorea sputorum, 3,365, 4,612, 5,076, 2,239, 2,721, 3,275, 84, 39, and 455 genes were annotated with the UniProt, Pfam, NR, COG, KEGG, GO, CAZy, CARD, and VFDB databases, respectively (Table 1 and Supplementary File S1). In the case of Mycetohabitans endofungorum EFB03829, out of 3,359 gene models, 1,659, 2,687, 3,159, 1,267, 1,761, 1,852, 65, 35, and 299 genes were annotated with the UniProt, Pfam, NR, COG, KEGG, GO, CAZy, CARD, and VFDB databases, respectively (Table 1 and Supplementary File S2).

The COG annotation results (Supplementary Figure S1) indicated that M. endofungorum EFB03829 had fewer genes in all groups compared to P. sputorum EFB03792, except for mobile genes, which were more in EFB03792 than in EFB03829 (35 vs. 12, Supplementary Figure S1).

In terms of CAZy annotation, P. sputorum EFB03792 possessed 12 auxiliary activity genes (AAs), five carbohydrate-binding module genes (CBMs), 13 carbohydrate esterase genes (CEs), 19 glycoside hydrolase genes (GHs), 34 glycosyl transferase genes (GTs), and one polysaccharide lyase gene (PL; Supplementary File S3). M. endofungorum EFB03829 completely lost CBM and PL genes but had four AAs, five CEs, 18 GHs, and 38 GTs (Supplementary File S3).

CARD annotation revealed that P. sputorum EFB03792 possessed a unique drug resistance gene, specifically the FAD-containing monooxygenase EthA, which confers resistance to ethionamide through antibiotic target alteration. There were no differences in drug resistance genes between M. endofungorum EFB03829 and HKI45.

VFDB annotation results showed that P. sputorum EFB03792 and M. endofungorum EFB03829 had 656 and 309 virulence factor-related genes, respectively (Supplementary Files S4, S5). Among these genes, 28 were specific in P. sputorum EFB03792 and 25 were specific in M. endofungorum EFB03829 (Figures 4C,D and Supplementary Table S1).

Additionally, two prophages were predicted in P. sputorum EFB03792, and one was predicted in M. endofungorum EFB03829 (Table 1). Nine and 25 genomic islands were annotated in P. sputorum EFB03792 and M. endofungorum EFB03829, respectively (Table 1).

For Pandorea sputorum, compared with the free-living strains NCTC13161 and DSM21091, the endosymbiotic strain EFB03792 exhibited an inversion in the structure for more than half of its genes (Figure 5A) and contained a higher number of mobile elements (111 vs. 60; Table 1).

Concerning Mycetohabitans endofungorum, in comparison with the obligate endosymbiotic strain HKI456, the facultative endosymbiotic strain EFB03829 displayed a significant number of genes involved in inversion and/or translocation (Figure 5B) and contained a much larger number of mobile elements (1,517 vs. 392; Table 2).

In this study, CAZy annotation was performed on five strains, namely Mycetohabitans endofungorum EFB03829, HKI456, Pandoraea sputorum EFB03792, DSM21091, and NCTC13161. The newly assembled genomes of M. endofungorum EFB03829 and P. sputorum EFB03792 were predicted with three and six CE10 family genes, respectively, while the other three genomes lacked (Figure 6). The CE10 family genes encoded some enzymes that catalyzed the hydrolysis of carboxylic ester bonds, such as acetyl-hydrolase, monoterpene epsilon-lactone hydrolase, acetyl esterase/lipase, and carboxylesterase.

Fewer CAZy genes were identified in M. endofungorum (65–67) than in P. sputorum (84). Mycetohabitans endofungorum exhibited a complete loss of CBM (CBM32 and CBM73) and PL5 family genes compared to P. sputorum (Figure 6 and Supplementary File S3). The CBM32, CBM73, and PL5 families encoded beta-galactosidase, chitin binding, and alginate lyase, respectively. All strains possessed a rich abundance of CEs, GHs, and GTs (e.g., CE1, GH23, GT4, GT83, and GT9 family) genes (Figure 6). GH108 family genes were enriched in M. endofungorum (seven in EFB03829 and three in HKI456), but none in P. sputorum. These GH108 family genes were presumed to have a putative peptidoglycan binding domain and a predicted peptidoglycan domain.

In this study, 24 and 74 genes related to the type III secretion system (T3SS) were predicted from Mycetohabitans endofungorum EFB03829 and Pandorea sputorum EFB03792, respectively (Supplementary File S6), and M. endofungorum EFB03829 and P. sputorum EFB03792 have completely T3SS. The T3SS of M. endofungorum EFB03829 have five ATPase complexes, two basal bodies, four cytoplasmic rings, three export apparatuses, eight regulators, and two invasion protein genes. The T3SS of P. sputorum EFB03792 exhibited five ATPase complexes, eight basal bodies, six cytoplasmic rings, six export apparatuses, 45 regulators, and five invasion protein genes.

Rhizopus arrhizus underwent continuous sub-culturing on a PDA plate containing antibiotics. Throughout the subculturing process, the sporangiophores of the strain XY03829 exhibited increased bending, and the production of sporangiospores gradually decreased. Starting from the 22nd generation, no sporangiospores were formed, and the strain could not recover to produce any sporangiospores (Figure 7). In contrast, the mycelial morphology of R. arrhizus XY03792 remained unchanged during the subculturing process.

With the progress in exploring EFBs, an increasing number of Burkholderia-related endobacteria (BRE) and Mycoplasma-related endobacteria (MRE) have been identified in fungi, especially within the phylum Mucoromycota (Bianciotto et al., 2003; Partida-Martinez et al., 2007b; Sato et al., 2010; Okrasińska et al., 2021; Uehling et al., 2023). However, Pandoraea sputorum, a member of the family Burkholderiaceae, has never been previously detected within fungal hosts.

Pandoraea sputorum represents an emerging human pathogen known for inducing a pro-inflammatory response that can lead to lung dysfunction in individuals with cystic fibrosis (Xiao et al., 2019). This pathogenic microorganism has been exclusively isolated from respiratory tract sources (Pimentel and Macleod, 2008; Martínez-Lamas et al., 2011; Fernández-Olmos et al., 2012; Pugès et al., 2015; Kwizera et al., 2017) and blood samples (Xiao et al., 2019). Alongside this species, the pathogenic genus Pandoraea includes ten other species (Xiao et al., 2019). While these species have been identified in various specimens, such as sputum, blood, urine, lung tissue, and wounds, they have not been observed within fungi. Therefore, this study marks the initial proposal of the pathogenic bacterium Pandoraea as an EFB genus, particularly thriving within a potential pathogenic fungus of Rhizopus arrhizus. This underscores the heightened relevance of this genus in the field of medicine.

It has been reported that the EFB Mycetohabitans rhizoxinica plays a role in enhancing sporulation in the host fungus Rhizopus microsporus (Partida-Martinez et al., 2007c). In this study, we detected EFB M. endofungorum in the host fungus R. arrhizus, a species closely related to R. microsporus (Liu et al., 2008). The cured R. arrhizus exhibited impaired growth and an inability to produce sporangiospores, underscoring the essential role of EFB M. endofungorum in the growth and reproduction of R. arrhizus. Attempts to restore sporulation through co-culturing on LB plates were unsuccessful, likely attributed to the limited infectivity of EFB03829 on the host. Laser-mediated microinjection (Partida-Martinez et al., 2007c) emerges as a potential superior method for the reintroduction of M. endofungorum into its host R. arrhizus.

Based on the complete assembly genome sequences available in the NCBI database¹⁴ for the Burkholderia genus, it is observed that their genome sizes span a range from 5.23 Mb to 10.63 Mb. Notably, the genome of Mycetohabitans endofungorum EFB03829 is markedly reduced, measuring only 3.64 Mb (composed of a 2,660,040 bp of chromosome, 800,149 bp of plasmid 1, and 181,953 bp of plasmid 2). This represents a significant reduction when compared to its Burkholderia spp. counterparts. While the genome of Pandoraea sputorum EFB03792 (5.85 Mb) does not exhibit notable streamlining when compared to other strains of the same species of P. sputorum, which range from 5.74 Mb to 6.45 Mb (see footnote 14). The prevailing consensus is that obligate endosymbionts undergo genome reduction as an adaptation to their reliance on host-derived resources (Uehling et al., 2023). In contrast, facultative endosymbiotic bacteria are generally not subject to genome reduction (Baltrus et al., 2017). Both M. endofungorum and P. sputorum are facultative endosymbiotic bacteria, demonstrating the ability to thrive not only within fungal mycelia but also on artificial media. Surprisingly, M. endofungorum EFB03829 exhibited a remarkable genome reduction.

According to Salvioli et al. (2017), the genome of Mycetohabitans rhizoxinica includes toxin-antitoxin modules (TAs), which involve in modulating growth under stress conditions and promoting survival in host cells. M. endofungorum also harbors TAs in its genome, potentially influencing the regulation of fungal endobacteria life (Lackner et al., 2011b; Salvioli et al., 2017). The genomes of both M. endofungorum and Pandorea sputorum strains encompass genes associated with virulence factors. And M. endofungorum has demonstrated the ability to produce the toxin rhizonin, exhibiting significant nonspecific hepatotoxicity (Partida-Martinez et al., 2007a). Previous research has indicated that the two R. arrhizus strains under investigation can undergo fermentation to generate glucose and lactic acid for food fermentation purposes (Liu et al., 2022). Consequently, the potential EFBs may pose a risk of toxin production in fermented foods.

The analyses of collinearity indicated a more pronounced change in the gene structure of M. endofungorum EFB03829 compared to P. sputorum EFB03792 (Figure 5). Furthermore, M. endofungorum EFB03829 possesses a significantly higher number of mobile elements (1,517) than M. endofungorum HKI456 (392) and P. sputorum (60–111; Table 1). The repeated insertion and loss of mobile elements, including prophages, can contribute to genome reduction and alterations in gene structure (Vale et al., 2022). Thus, the substantial presence of mobile genetic elements in M. endofungorum EFB03829 is implicated in its genome reduction and structural changes.

CAZymes, or Carbohydrate-Active Enzymes, play a key role in metabolism, involved in the synthesis, modification, and degradation of carbohydrates, including polysaccharides, glycoproteins, and glycolipids (Drula et al., 2022). Our results suggested that the newly sequenced genomes were annotated several CE10 family genes, suggesting their involvement in the metabolism of various compounds in the host R. arrhizus, including drugs, pesticides, and lipids.

The type III secretion system (T3SS) plays a vital role in maintaining symbiosis and is highly conserved in the genomes of Gram-negative pathogenic or symbiotic bacteria (Deng et al., 2017), such as endosymbiont Candidatus Glomeribacter gigasporarum associated with the arbuscular mycorrhizal fungus (AMF) Gigaspora margaritain (Ghignone et al., 2012) and Burkholderia rhizoxinica in the zygomycetous fungus Rhizopus microsporus (Lackner et al., 2011a). In our study, we predicted 74 and 24 genes related to T3SS in P. sputorum EFB03792 and M. endofungorum EFB03829, respectively, indicating their role as symbiotic bacteria with Rhizopus arrhizus. The identification of these specific genes in our study contributes to a better understanding of the mechanisms underlying the actions of symbiotic bacteria in fungi during invasion and colonization stages.