Bioinformatics analysis of the Microsporidia sp. MB genome: a malaria transmission-blocking symbiont of the Anopheles arabiensis mosquito

Abstract

Background

The use of microsporidia as a disease-transmission-blocking tool has garnered significant attention. Microsporidia sp. MB, known for its ability to block malaria development in mosquitoes, is an optimal candidate for supplementing malaria vector control methods. This symbiont, found in Anopheles mosquitoes, can be transmitted both vertically and horizontally with minimal effects on its mosquito host. Its genome, recently sequenced from An. arabiensis, comprises a compact 5.9 Mbp.

Results

Here, we analyze the Microsporidia sp. MB genome, highlighting its major genomic features, gene content, and protein function. The genome contains 2247 genes, predominantly encoding enzymes. Unlike other members of the Enterocytozoonida group, Microsporidia sp. MB has retained most of the genes in the glycolytic pathway. Genes involved in RNA interference (RNAi) were also identified, suggesting a mechanism for host immune suppression. Importantly, meiosis-related genes (MRG) were detected, indicating potential for sexual reproduction in this organism. Comparative analyses revealed similarities with its closest relative, Vittaforma corneae, despite key differences in host interactions.

Conclusion

This study provides an in-depth analysis of the newly sequenced Microsporidia sp. MB genome, uncovering its unique adaptations for intracellular parasitism, including retention of essential metabolic pathways and RNAi machinery. The identification of MRGs suggests the possibility of sexual reproduction, offering insights into the symbiont’s evolutionary strategies. Establishing a reference genome for Microsporidia sp. MB sets the foundation for future studies on its role in malaria transmission dynamics and host-parasite interactions.

Background

Malaria cases have been on a steady increase since 2016 with a notable 249 million cases reported in 2022 [1]. A total of 608 000 malaria deaths were recorded during this period, with 95.4% of the total global deaths from the WHO Africa region alone [1]. The major impediments to global malaria control have included the development of insecticide resistance in the Anopheles vector, Plasmodium drug resistance, poor detection of the minor malaria parasites, and the recent invasion of a new vector, An. stephensi, in the Horn of Africa and urban environments [1,2,3,4,5]. Due to these growing challenges, innovations in vector control have been the focus to mitigate resistance and improve the effectiveness of the current strategies. This has included the use of biocontrol tools to complement the current vector control tools.

Symbionts, such as bacteria and fungi, have been proposed as potential tools for controlling mosquito-borne diseases, including malaria [6,7,8,9,10]. These symbiotic microorganisms can affect the lifespan, reproductive competence, and vector competence of mosquitoes, which in turn can affect the transmission of mosquito-borne infectious diseases [6,7,8, 11,12,13,14,15]. Several approaches are being explored in the field of symbiotic control for mosquito-borne diseases. These include disrupting microbial symbionts required by insect pests, manipulating symbionts to express anti-pathogen molecules within the host, and introducing endogenous microbes that affect the lifespan and vector capacity of new hosts in insect populations [16]. Wolbachia is one such symbiont used to successfully suppress dengue transmission in Aedes mosquitoes [17,18,19]. Another unique example is the microsporidia species, which has been highlighted to have the capacity to function effectively as a biocontrol tool [9, 10, 20].

Microsporidia are a set of spore-forming obligate intracellular organisms affecting a diverse variety of vertebrates and invertebrates [21]. Microsporidia have an intricate infection mechanism which is characterized by the formation of a unique polar tube to invade the host cell [22, 23]. They have multiple modes of transmission including horizontal transmission via ingestion of spores, and vertical transmission from an infected mother to their offspring [21]. Microsporidia have the smallest genome sizes among eukaryotes, with the smallest genome of 2.3 Mb recorded from E. intestinalis [24, 25]. These enigmatic microorganisms have garnered significant attention due to their unique genomic characteristics, evolutionary adaptations, and varying host-microbe interactions [26].

Over half of the identified microsporidia species infect insects [27, 28]. These microsporidia have varying modes of infection and can have significant impacts on both the economy and the regulation of insect populations [29]. Several microsporidia species are known to be vertically transmitted, where the parasite is passed from parent to offspring. Strictly vertically transmitted microsporidia, like Nosema granulosis, Dictyocoela roeselum and Dictyocoela muelleri, have also displayed feminizing attributes, impacting the population dynamics of their hosts [30]. Other microsporidia display multiple modes of transmission, like Edhazardia aedis affects the mosquito Aedes aegypti, with both horizontal and vertical transmission pathways. Interestingly, the mode of transmission of Edhazardia aedis has been shown to dictate the blood-feeding success of their host with the microorganism being less virulent during its vertical transmission stages [31]. Some microsporidia exhibit relatively low levels of pathogenicity, making them less harmful to their hosts [32]. For example, Vavraia culicis, which infects various mosquito species with minimal pathogenic effects [26]. Agriculturally important pathogenic microsporidia, such as Nosema bombycis and Nosema apis (also known as Vairimorpha apis), affect silkworm and honeybee populations, respectively. Conversely, other microsporidia species have been explored as biocontrol tools against agricultural pests such as the use of the pathogenic Paranosema locustae to control grasshopper populations in the field and Nosema pyrausta in the control of the European corn borer [10, 27].

Mosquitoes are important hosts for microsporidian parasites showing high levels of host specificity and varying effects on their respective hosts [33]. These species belong to the genera Amblyospora, Hazardia, Encephalitozoon, Enterocytospora, Nosema, and Microsporidium [33, 34]. Owing to the extensive effect of microsporidia in mosquitoes, these have the potential for use in the control of vector populations and disease transmission [35]. Pathogenic microsporidia infecting mosquitoes include Amblyospora connecticus isolated from the saltmarsh mosquito [36], and Nosema algerae infecting An. stephensi [37], and Parathelohania anophelis affecting several Anopheles mosquitoes including An. quadrimaculatus [38]. Edhazardia aedis has been isolated from the Aedes species [31] while Vavraia culicis has been shown to affect plasmodium development in An. gambiae [39]. These microsporidia demonstrate significant potential in the control of vector-borne diseases. The most recent of which has been the discovery of Microsporidia sp. MB conferring protection against malaria development [40].

Microsporidia sp. MB (often denoted as Microsporidia MB) is a microsporidian symbiont that has been found in Anopheles mosquitoes in various regions of Africa, including An. gambiae s.l. in Kenya, An. coluzzii in Niger and both An. gambiae s.s. and An. coluzzii in Ghana [40,41,42,43,44]. It has been shown to have a strong malaria transmission-blocking phenotype, making it a potential candidate for the development of a symbiont-based malaria transmission-blocking strategy [40, 45]. The symbiont has multiple transmission routes including horizontal transmission between adult mosquitoes, specifically between male and female An. Arabiensis and maternal transmission [40, 42, 43]. Female An. arabiensis that acquire Microsporidia sp. MB horizontally can transmit the symbiont vertically to their offspring [42]. It has also been found that Microsporidia sp. MB can infect Anopheles funestus s.s., another primary malaria vector in Africa [42]. The prevalence of Microsporidia sp. MB in Anopheles mosquitoes varies across different regions and in different seasons, with a prevalence of 6% in Kenya and 1.8% in Ghana [41]. Further investigations are needed to understand the diversity and range of Microsporidia sp. MB across sub-Saharan Africa. Overall, Microsporidia sp. MB shows promise as a potential tool for controlling malaria transmission, but more research is needed to fully understand its effectiveness and potential applications [45].

Over the past few decades, researchers have made significant progress in understanding the genomic structures of microsporidia, an organism that is otherwise difficult to tame in laboratory setups owing to its obligatory intracellular nature. Through genome sequencing, scientists have been able to gain valuable insights into the biology, evolution, gene diversity, and pathogenicity of microsporidia [24, 25, 30, 46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71]. These findings have vastly highlighted the highly specialized lifestyle of microsporidia and their dependence on their host for essential biochemical processes. Furthermore, genome size reduction and gene loss in microsporidia species are closely linked to their reliance on host-derived nutrients and biochemical processes, as many genes involved in essential metabolic pathways are no longer needed when these functions are supplemented by the host[24, 46, 48, 49, 51, 57, 58, 72, 73]. This dependence on the host makes microsporidia an interesting model for studying the intricate interactions between microorganisms and their hosts. Challenges in annotating microsporidia genomes arise from their small size, high gene density, absence of introns, and the need for specialized methodologies to accurately predict and validate gene structures, which are influenced by their unique evolutionary adaptations to obligate intracellular parasitism [74].

To understand the mechanism of protection of Microsporidia sp. MB against malaria transmission, a high-quality reference genome is necessary. The draft genome of Microsporidia sp. MB isolated from Anopheles arabiensis in Kenya was recently sequenced using DNA NanoBalls Sequencing (DNBSEQ) short paired-end sequencing technology from dissected mosquito ovaries and de novo assembled [75].

In this study, we describe the structure of the Microsporidia sp. MB genome and further highlight its gene composition and key conserved proteins involved in important metabolic pathways within the host (75–80) using several computational tools. Microsporidia, characterized by their highly reduced genomes and lack of mitochondria, have lost most of the components of the glycolytic pathway but utilize host resources, indicating their reliance on host energy sources [76, 77]. This study revealed the presence of 10 of 11 key glycolytic genes in Microsporidia sp. MB, that have otherwise been lost in other microsporidia genomes. Additionally, we explored the RNA interference (RNAi) machinery within this symbiont, a complex previously implicated in the regulation of gene expression in microsporidia such as Nosema ceranae [78, 79].

Results and discussion

Comparative analysis of Microsporidia sp. MB and Vittaforma corneae ATCC 50505 genome structures

Microsporidia sp. MB and Vittaforma corneae ATCC 50505 are the closest relatives in Clade IV (Enterocytozoonidae) group based on 18S taxonomic classification [80], yet they differ significantly in their biological roles, transmission dynamics, and effects on their hosts. While Microsporidia sp. MB is a symbiotic organism primarily found in mosquitoes with a strong malaria transmission-blocking phenotype, V. corneae is known for causing ocular infections in humans, particularly in immunocompromised individuals. Based on these key differences, similarities and the availability of these sequence genomes online, a comparative overview of the two complex microsporidia genomes was conducted.

The de novo assembly of the Microsporidia sp. MB genome was based on a total of 3,874,899 reads, with an estimated sequencing depth of 100-fold coverage. The quality of the assembled genome was analyzed using the QUality Assessment Tool (QUAST) revealing that Microsporidia sp. MB had an average size of 5.9 Mb spanning across 2335 contigs compared to 3.2 Mb (314 contigs) in V. corneae ATCC 50505 (GCA_000231115.1; accessed December 14, 2023) (Table 1). Despite having relatively different genome sizes, the number of genes predicted in both genomes was relatively similar (Microsporidia sp. MB = 2247; V. corneae = 2239).

Table 1 Comparison of Microsporidia sp. MB Genome Statistics with its closest relative V. corneae ATCC 50505

The statistics above were further visualized on BlobTools as a snail plot (Fig. 1a). The Benchmarking Universal Single-Copy Orthologs (BUSCO) tool suite was used to assess the level of completeness of the assembled genome by checking the presence of core microsporidia genes that were represented in the genome (Fig. 1b).

Fig. 1

Microsporidia sp. MB genome statistics. a The snail plot created using BlobToolKit highlights the genome metrics. A basic legend of the assembly statistics is provided in the top left corner. The outer blue rings highlight the GC content over the assembly, corresponding to each of the contigs in the grey ring. The overall GC content of the assembly is summarized at the bottom left corner. The light grey rings show the total length of the assembly. The red line highlights the longest scaffold, while the dark and light orange rings show the length of the contigs at which 50% and 90% of the total assembly is represented, respectively. b A pie-chart representation of the BUSCO gene completeness analysis shows a total of 81% single-copy core microsporidia genes present in Microsporidia sp. MB

Several reports suggest that the varying sizes of microsporidia genomes are mostly attributed to the proportion of repeats in their assemblies. To identify whether the size of the assembled genome was due to the presence of such repeats, their composition was analyzed using RepeatMasker v4.0 [81]. The Microsporidia sp. MB genome constituted 99.43% non-repeat regions and 0.57% repeats which included: Small RNA repeats (n = 821; 0.01%), low complexity repeats (n = 10,645; 0.18%), and simple repeats (n = 22,870; 0.38%) (Table 2). Short interspersed nuclear elements (SINE), Long interspersed nuclear elements (LINE), long terminal repeat retrotransposons (LTR), and DNA transposons were notably absent in this assembly. Small RNA repeats refer to repetitive sequences associated with non-coding small RNAs, including microRNAs (miRNAs) and small interfering RNAs (siRNAs). A list of the different types of repeats identified per contig and their positions in the sequences are highlighted in Supplementary Data File 1.

Table 2 Distribution of different types of repeats identified in the Microsporidia sp. MB assembly

Gene prediction and gene organization

A total of 2247 genes were predicted using GeneMark-ES intronless eukaryote mode and each protein was checked against the Microsporidia taxon (Taxon ID: 6029) using NCBI BLASTp (E-value = 0.0001). The predicted proteins in Microsporidia sp. MB range in size from 99 to 263 amino acids. Supplementary Data Files 2–4 contain information on the predicted gene sequences including genbank file, predicted coding sequences and predicted proteins respectively. Moreover, Supplementary Data File 5 contains a highlight of BLAST output results from these sequences.

A set of genes annotated using GeneMark-ES and corresponding to BUSCO identification of single-copy and complete genes were used to plot the genome architecture schematic in Fig. 2 below. The representation was based on the top contigs of lengths greater or equal to 10 kb, ranging from 26.338 – 10.051 kb (n = 48).

Fig. 2

Genome structure of Microsporidia sp. MB, showing contigs larger than 10 kb to focus on the significant portions of the assembly. The figure illustrates gene density and organization across these contigs, with the Y-axis representing the names of the contigs. Each gene is color-coded based on its functional classification, grouped into 12 distinct categories as indicated in the legend. The gene annotations are derived from GeneMark-ES and validated using BUSCO, ensuring high confidence in the predicted gene functions. This schematic helps visualize the distribution and functional diversity of genes across the larger contigs in the genome.The plot was created in RStudio [82] using the gggenes package (https://cran.r-project.org/web/packages/gggenes)

Analysis of common orthogroups in Microsporidia sp. MB and the Enterocytozoonida group of Microsporidia

The comparative phylogenomic analysis of Microsporidia sp. MB against the Enterocytozoonida group revealed that this symbiont is closely related to V. corneae ATCC 50505 as has been previously seen using the 18S region [80]. Ninety-two common orthogroups were found in all 8 microsporidia species and were used to infer phylogeny within the group (Fig. 3).

Fig. 3

Phylogenomic analysis of common BUSCO orthogroups found in Enterocytozoonida (left) and their respective genome assembly completeness using BUSCO (right). Microsporidia sp. MB is bolded. Nosema granulosis was used as an outgroup. Alignments of common orthogroups were performed using MUSCLE and trimAl [83, 84]. The maximum likelihood tree (1000 bootstraps) was constructed using IQ-TREE v1.6.12 [85] based on the 92 single-copy orthologous genes (with the best model selection (LG + F + I + G4) from ModelFinder [86]. Bootstrap percentage values are shown on each node and the scale bar at the bottom of the tree indicates total changes per site

In addition to comparing the genome of Microsporidia sp. MB with V. corneae ATCC 50505, recently published genomes from the Enterocytozoonida group were also considered, including Agmasoma sp. [87], Enteropsectra and Pancytospora [88] and, Alternosema astaquatica [89]. The nematode-infecting microsporidia, Enteropsectra possesses a genome size of approximately 6.0 Mb, encoding around 3531 predicted proteins [88]. In contrast, the two Pancytospora species exhibit smaller genomes of 4.2 Mb and 4.6 Mb, with predicted protein counts ranging from 2300 to 2800 [88]. When compared to Microsporidia sp. MB, which has a genome size of 5.9 Mb and 2247 predicted proteins, Enteropsectra displays a larger genome and a substantially higher number of predicted genes, suggesting greater genomic complexity. Conversely, the Pancytospora species show reduced genomic size and fewer predicted proteins, highlighting the genomic diversity within this group.

Sequence-based protein analysis

Protein functional analysis

Protein family classification

InterProScan [90] was used to classify the superfamilies of the putative Microsporidia sp. MB proteins predicted using BUSCO and GeneMark-ES tools. A total of 2,247 proteins were predicted with 2235 of these classified into their respective superfamilies (Supplementary Data File 6). Over 30% (n = 707) of these proteins were associated with enzymes and only 12 proteins (0.53%) were left unclassified. The other protein family groups identified included: cell-cycle (n = 530; 23.6%), regulatory/cellular (n = 364; 16.2%), binding (n = 353; 15.7%), transport (n = 178; 0.9%), and structural proteins (n = 103; 4.6%) (Fig. 4a, Supplementary Data File 2). The large distribution of enzymes in this genome is comparatively similar to a report on the protein family classes of V. corneae ATCC 50505 [91].

Fig. 4

Protein family classification using InterProScan. a Seven main protein superfamily groups were identified, with most of the proteins classified as enzymes. b Zooming into the classification of the primary enzyme groups using InterProScan, hydrolases and transferases are highlighted as the most abundant enzyme classes in Microsporidia sp. MB. These two enzymes function in energy metabolism and transport of nutrients, respectively

Distribution of enzyme classes

Since most of the predicted proteins were classified as enzymes, further analysis of the different groups of enzymes was essential to determine the most common enzyme class represented in the dataset. Seven main enzyme classes were identified, with transferases appearing to be the most common category among the annotated enzymes Fig. 4b. These were divided into: hydrolases (n = 237; 38.47%), transferases (n = 228; 37.01%), ligases (n = 45; 7.31%), isomerases (n = 35; 5.68%), oxidoreductases, (n = 33; 5.36%), translocases (n = 31; 5.03%) and lyases (n = 7; 1.14%).

Each enzyme class was further grouped into their respective subclasses based on their different specific functions as shown in Fig. 5 below. Hydrolases are involved in the carbohydrate and energy metabolism of intracellular stages of the microsporidian Nosema grylli [92]. Furthermore, the localization and functionality of microsporidian iron-sulfur cluster assembly proteins, which are likely to be hydrolases, have been investigated by Goldberg [93]. However, the specific functions of hydrolases in microsporidia, particularly in the context of their pathogenicity and host interactions, remain to be fully elucidated. Here, we identified 38.47% (n = 237) of the classified enzymes in Microsporidia sp. MB as hydrolases with most of these acting on acid anhydrides (Fig. 5).

Fig. 5

Distribution of different enzyme classes and subclasses in Microsporidia sp. MB as predicted through InterProScan. “n = ” denotes the number of proteins in each respective subclass while the x -axis represents the proportion of each subclass in their respective enzyme class in percentage

Transferases play a crucial role in microsporidia intracellular lifestyle, particularly in the acquisition of essential nutrients and energy from their host cells. These enzymes also facilitate the import of nucleotides for nucleic acid biosynthesis. Nine main subclasses of transferases were identified and those involved in transferring phosphorus-containing molecules were the most common (Fig. 5).

The role of ligases in microsporidia is further supported by the identification of several components of the mitochondrial iron-sulfur cluster assembly machinery in these organisms. These components, which participate in the biosynthesis of iron-sulfur proteins, are essential for various cellular processes, including DNA repair [94]. The presence of these components suggests that microsporidia may use ligases to maintain their genome integrity and ensure their survival within the host cell [55, 95]. Microsporidia sp. MB contained 7.31% ligases of the predicted enzymes, which were further classified into 5 subclasses depending on the type of bonds they create with different molecules (Fig. 5 above).

A total of eleven oxidoreductase subclasses were observed. Isomerases, oxidases, and thiol-sulfide reductases are important in redox reactions and useful in the survival of the organism in new host environments. Here, 5 isomerase subclasses were identified, some of which had direct links to oxidoreductases and transferases. Lyase proteins were the least-represented enzyme classes in the genome consisting of carbon–carbon, carbon–oxygen, and phosphorus-oxygen subgroups (Fig. 5 above).

Distribution of key domains and repeat regions

Domains, repeats, and conserved sites are key features found in proteins. Domains in proteins are distinct regions that have specific structural and functional properties and are typically composed of several secondary structure elements, such as alpha helices and beta sheets, and can fold independently. The term 'key domains' refers to conserved regions within proteins that are essential for their functional properties, such as enzyme activity or interaction with host cells. These domains include well-characterized motifs like reverse transcriptase domains, leucine-rich repeats (LRR), and others identified through InterProScan. Repeats, on the other hand, are segments of protein sequences that are repeated multiple times within the same protein or across different proteins. Repeats can be indicative of functional motifs or structural elements within proteins.

The predicted proteins contained a total of 659 domains and 19 repeats. A total of 1456 domains were identified in the putative proteins with a majority being the reverse-transcriptase domains (n = 179) as shown in Fig. 6a. Microsporidia lack nucleotide biosynthetic pathways and rely on "stealing" nucleotides from host cells [96, 97]. Reverse transcriptase domains are crucial components of microsporidia proteins as they allow the proteins to perform reverse transcription, a process in which RNA is converted into DNA [98]. Repeat regions within these proteins were also identified through InterProScan as shown in Fig. 6b below.

Fig. 6

Sequence distribution of the different domains and repeats predicted using InterProScan. a A bar plot of the domains identified from the genome shows a large set of domains involved mainly in enzymatic activities. A zoomed image of the main domains excluding the outlier (where less than 10 sequences were represented in each domain) is shown on the right panel. b The distribution of repeat regions is depicted in this bar plot. Leucine-rich repeat regions (LRR) and WD40 were the most abundant. WD40 proteins named for their characteristic WD40 repeat motifs which consist of 40 amino-acids typically ending with tryptophan (W) and aspartic acid (D), are involved in a wide variety of cellular processes due to their ability to mediate protein–protein interactions

Most notable was the presence of leucine-rich repeat regions (LRR) with a total of 17 proteins containing these repeats (Fig. 6b above). The LRRs found in microsporidia proteins play a crucial role in the recognition and interaction with host cells and may represent pathogenicity factors. These regions are often involved in mediating protein–protein interactions, allowing microsporidia to bind to and invade host cells [98]. Microsporidian LRR-proteins secreted into the host have the potential to interact with host protein function. LRR-proteins in S. lophii were shown to have N-terminal signal peptides indicating potential secretion or targeting to the parasite cell surface for interactions with the host [60]. The presence of similar sequences in V. corneae also demonstrates the significance of LRR-proteins in mediating protein–protein interaction [60].

An investigation of the distribution of LRRs in the proteome using MEME Suite showed that 9 of the assembled contigs contained key LRR motifs (Fig. 7). The secretome, which refers to the collection of proteins secreted by an organism, plays an important role in host-parasite interactions by influencing immune responses and cellular processes in the host. To identify whether the LRR-containing proteins could be part of the Microsporidia sp. MB secretome, an important feature in host-parasite interactions, further computational assessment of the presence of signal peptides and transmembrane localization prediction was done. Three of the identified contigs contained putative signal peptides suggesting they could be part of the secretome, in direct contact with the host cell (Supplementary Table 2).

Fig. 7

a Organization of leucine-rich repeats (LRR) across 9 contigs. MEME motif output showing conservation levels across the motif. Larger symbols indicate more common amino acids. Y-axis (bits) shows conservation levels—higher bits reflect higher conservation. X-axis denotes positions along the motif, with higher Y values indicating consistent amino acids at specific positions across analyzed sequences. Two LRR motifs, Motif 1 (red, 17 amino acids) and Motif 2 (blue, 20 amino acids), were identified using MEME Suite. Moreover, the prediction of signal peptides using TargetP v2.0 and SignalP v6.0 (green and brown respectively) in 2 of these sequences (NODE_405 and NODE_1665). b A 3D model of gene_679 (NODE_405) with motif mapping. The structure was modeled on PRIMO [99]. The identified LRR motifs are numbered based on MEME output

One contig (NODE_405) contained 2 main LRR genes (gene_679 and gene_680) with confirmed signal peptides using both TargetP and SignalP prediction tools while contig NODE_1665 contained one gene (gene_1719). Gene_679 was identified to be homologous to the Anopheles plasmodium-responsive leucine-rich repeat protein APL1B which has been shown to play a role in An. gambiae immune response to invasion by plasmodium [100] and could suggest a probable function of gene_679 in blocking malaria in its host.

The 3D model of gene_679 was designed using template APL1 (PDB ID: 3O6N_A; 1.85 Å) and was characterized with 8 alpha helices, 11 beta sheets, and loops (Fig. 7b, Supplementary Fig. 1). QMEAN score of 0.5, Verify3D overall score of 77.98%, and Z-DOPE score of −3.83. Moreover, gene_679 contained both signal peptides and transmembrane segments. No good-quality template was identified from HHSearch to model gene_680 and gene_1719.

Gene ontology

Gene ontology (GO) analysis is frequently used to annotate genes and gene products and to determine the biological features of large-scale genome or transcriptome data [101]. Here, we illustrate using GO that most of the putative proteins in Microsporidia sp. MB are involved in metabolic processes and transport attributed to the bulk of enzymes found in the predicted dataset. Most of these proteins were found in the nucleus and cytoplasm (Fig. 8).

Fig. 8

Gene Ontology Analysis. A highlight of the protein sequence distribution in each primary GO term (Biological Process, Cellular Component, and Molecular Function) shows that the bulk of the putative proteins are involved in metabolic processes, localized within the intracellular anatomical structure (nucleus and cytoplasm), and involved in enzymatic activities

KEGG pathway analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a collective database integrating genomes, biological pathways, diseases, drugs, and chemical substances [102]. A total of 257 pathways were identified with 741 sequences linked to these pathways (E-value = 0.0001). KEGG further categorized the proteins into five distinct groups with most of the proteins involved in metabolism (n = 83; 32.3%) and organismal systems (n = 83; 32.3%) (Fig. 9a). The other categories included environmental information processing (n = 35; 13.6%), cellular processes (n = 30; 11.7%), and genetic information processing (n = 26; 10.1%).

Fig. 9

KEGG metabolic pathway analysis. a A barplot representation of the main KEGG categories shows most of the pathways identified within the genome take part in the organismal systems and metabolism. b Sequence distribution in each KEGG pathway highlighting the top 25 pathways represented in the dataset. A large portion of the sequences were involved in purine and thiamine metabolism, cellular processes, and genetic information processing. “n” represents the total number of proteins in each category

Sequence distribution in each represented pathway showed that most of the predicted proteins are involved in metabolism (purine, thiamine, nucleotide metabolism and biosynthesis of cofactors), genetic information processing (ribosome biogenesis, protein processing, nucleotide excision repair, RNA degradation and RNA polymerase, among others), and general cellular processes (cell cycle, meiosis, cellular senescence and nucleocytoplasmic transport) (Fig. 9b).

Physicochemical properties

The physicochemical properties of proteins, such as surface hydrophobicity, structural stability, and surface charge, play a critical role in determining their biological functions and interactions. These properties are influenced by the protein structure, including the content of α-helices and β-sheets, which in turn affect the flexibility of the protein and functional characteristics. Here, we characterized the physicochemical properties of Microsporidia sp. MB using the EXPASY-ProtParam tool and observed that 1911 (85%) proteins were hydrophilic (Grand Average of hydropathy (GRAVY) score less than zero) while 336 (15%) proteins were hydrophobic (GRAVY score greater or equal to zero). Hydrophilic proteins are more likely to be connected to interactions with the host and thus makes the hydrophobicity test relevant for predicting potential host-symbiont interacting proteins in the dataset. Many of the proteins were identified to be hydrophilic suggesting the heavy reliance of microsporidia on host interactions.

Transmembrane regions and signal peptides prediction

Secretory proteins commonly contain signal peptides or span across membrane regions to allow direct interaction with the host cell. Therefore, identifying signal peptides and transmembrane regions is important in highlighting probable microsporidia-host interactions. The secretion of parasite proteins from different functional categories suggests their active role in manipulating host cell machinery shedding light on the mechanisms employed by microsporidia to modify host metabolic and regulatory pathways [98]. Here, we highlight 84 proteins from different functional groups containing signal peptides using SignalP v6.0 suggesting that these proteins could be part of the Microsporidia sp. MB secretome. Moreover, an overview of transmembrane regions using DeepTMHMM [103] identified 328 (14.6%) proteins spanning these regions while 1919 (85.4%) of these proteins did not have any transmembrane segments (Supplementary Table 3).

Eighty-four proteins had signal peptides and 23 of these contained transmembrane regions (Supplementary Table 3). Moreover, out of the 84 proteins, 69 were hydrophilic based on GRAVY scores. Thirteen of these were characterized as hydrophilic and contained both signal peptides and transmembrane segments.

Microsporidia sp. MB has retained most of the components involved in the glycolytic pathway

Microsporidia sp. MB contained 10 of the 11 proteins involved in the glycolytic pathway (Fig. 10). Among the Enterocytozoonida species, E. cuniculi has been observed to have retained all the proteins in this pathway while E. bieneusi seems to have lost all of these. The presence of all glycolytic genes in E. cuniculi and their absence in E. bieneusi likely reflect biological differences, such as adaptations to varying levels of host dependency and metabolic requirements, rather than being solely due to differences in genome annotation completeness. While E. cuniculi has been widely studied and annotated, the retention or loss of glycolytic genes is likely driven by the specific biological needs of the organism. The missing protein in Microsporidia sp. MB (3-phosphoglycerate kinase) could also suggest that this gene has not been sequenced yet. Comparably, its closest sister V. corneae and Nosema granulosis had the same set of proteins involved in the glycolytic pathway (Fig. 10).

Fig. 10

Overview of the loss and retention of key genes involved in the glycolytic pathway genes across different microsporidia species. Retained genes are highlighted in light blue in the heatmap, while those that have been lost are highlighted in gray. Microsporidia sp. MB has retained most of the genes involved in glycolysis except for Pyruvate kinase, similar to V. corneae and N. granulosis. Comparatively, E. aedis, E. cuniculi, and E. hepatopenaei appear to have lost all these genes, suggesting a total reliance on the host for energy production. The maximum likelihood phylogenomic tree on the left was generated using IQ-TREE v1.6.12 [85] based on the 196 common BUSCO single-copy orthologous genes found in all species. Percent bootstrap support values are shown on each node and the branch length scale is bare at the bottom of the tree

The presence of Meiosis-related Genes (MRG) in Microsporidia sp. MB indicates possible mode of sexual reproduction

Polyploidy, which involves having multiple sets of chromosomes, is common in microsporidia and shows dynamic changes across different species, suggesting its significant role in their adaptability and evolution [104]. Meiosis is a crucial process in sexual reproduction, leading to the formation of haploid gametes. In microsporidia, meiosis-related genes and their regulation have not been fully studied. However, meiosis plays a pivotal role in the formation of spores, which are essential for their transmission and infection [105]. The identification of meiosis-related genes in microsporidia is essential to understanding the mechanisms underlying their sexual reproduction and transmission.

The process of meiosis occurs in microsporidia during the early phases of sporogony [106]. The merging of haploid nuclei within a diplokaryotic sporont results in the production of a single diploid nucleus. A new diplokaryon with two diploid nuclei is then created by the mitotic division of this single nucleus. These diploid nuclei split, causing the formation of synaptonemal complexes within each nucleus. They then fuse to form a single tetraploid nucleus, where chromosomal mingling occurs. Eight haploid nuclei are finally produced by this tetraploid nucleus by two consecutive reductive divisions, two diploid nuclei, four haploid nuclei, and a final mitotic division [107, 108]. Despite the reduced complexity of microsporidia genomes, investigations have revealed genes that show homology or similarity to known meiotic genes in other eukaryotes [109]. These genes include homologs of meiosis-specific genes found in other organisms, such as Dmc1, Msh4, and Msh5, among others (Table 3) and are frequently involved in basic meiotic processes such as recombination, synapsis, chromosome segregation, and spore production. Their presence and functionality indicate that, despite their reduced genomes, microsporidia have kept critical genetic components required for meiosis, albeit with certain modifications and adaptations peculiar to their parasitic lifestyle and evolutionary history. Identifying and understanding these meiosis-related genes aids in the study of microsporidia reproductive strategies and evolutionary biology.

Table 3 A list of MRGs found in Microsporidia sp. MB

This study highlights the presence of several MRGs in the newly assembled genome. Some of the genes involved in meiosis and identified in this study are listed in Table 3 below. The presence of these genes in Microsporidia sp. MB suggests the putative mode of reproduction of this organism.

Microsporidia sp. MB contains the three main components of the RNAi complex

The RNA-induced silencing complex (RISC) is a conserved molecular machinery responsible for post-transcriptional gene regulation through small RNA molecules, particularly microRNAs (miRNAs) and siRNAs. While the RNAi pathway and RISC are well-studied in many eukaryotes, including animals, plants, and fungi, the presence, and mechanisms of RNAi in microsporidia, have only been recently highlighted [61, 78, 79, 110]. The presence of RNAi machinery in microsporidia varies among species (Table 4). While some microsporidia have retained elements of the RNAi pathway, others have lost or significantly modified these components. This variation is attributed to differences in evolutionary history and the degree of genome reduction among microsporidian lineages. If functional, the RNAi pathway in microsporidian potentially plays various roles, including the regulation of host-parasite interactions, suppression of host immune responses, or even self-regulation of microsporidian gene expression. RNAi could also be involved in the maintenance and evolution of the reduced microsporidian genome [13, 54]. Microsporidia sp. MB was observed to contain homologs of key RNA interference proteins. These included the dicer (gene_297), RNA-dependent RNA polymerase (gene_163), and argonaute (gene_271) genes (Table 4). All three of these genes were also present in V. corneae ATCC 50505 but were noticeably missing in the other microsporidians of the Enterocytozoonida group (Table 4).

Table 4 An overview of the retention of the key RNAi components in Enterocytozoonidae and mosquito microsporidia

Argonaute

The phylogenetic classification of argonaute is shown in Fig. 11a below with Microsporidia sp. MB consistently showing high similarity in protein sequence with V. corneae ATCC 50505. The argonaute protein is characterized by two main conserved domains: the PIWI-like and PAZ argonaute-like domains. These were also identified in the Microsporidia sp. MB homolog (Fig. 11b, Supplementary Fig. 2) and the conservation of the domains across several microsporidia species highlighted through multiple sequence alignment.

Fig. 11

Analysis of argonaute orthologs. a A rooted phylogenetic tree of argonaute orthologs based on 17 microsporidian species. The tree was constructed on protein sequences of the argonaute orthologs based on maximum clade credibility with Mr Bayes. Percent bootstrap values are indicated at the nodes. The scale for the branch length is shown at the bottom. b A schematic representation of the position of conserved domains Piwi-like and PAZ are mapped onto the multiple sequence alignment

Dicer

Phylogenetic analysis of the dicer orthologs found in 15 microsporidia species was done using (Fig. 12a) which showed high similarity between Microsporidia sp. MB and V. corneae ATCC50505. An analysis of the conserved domains of the dicer gene on CDD revealed that this contained two main conserved domains including the Ribonuclease III C terminal domain (RIBOc) and the dsRNA-specific ribonuclease domain (Rnc) (Fig. 12b). Multiple sequence alignment showed the presence of inserts between the RIBOc and Rnc domains was found in Microsporidia sp. MB, E. breve, P. epiphaga, and P. philotis. Supplementary Fig. 3 highlights the conservation of these genes in the multiple sequence alignment.

Fig. 12

Analysis of the dicer ortholog in Microsporidia sp. MB. a A rooted phylogenetic tree of argonaute orthologs based on 20 microsporidian species. The maximum likelihood tree was constructed using IQTREE. The scale for the branch length is shown at the bottom and bootstrap support values shown on each node as a percentage (1000 bootstraps). b A schematic representation of the position of conserved domains in the multiple sequence alignment is shown. An insert was observed between the RIBOc and Rnc domains in 4 of the species including Microsporidia sp. MB, E. breve, P. epiphaga, and P. philotis

RNA-dependent RNA polymerase

RNA-dependent RNA polymerase is involved in post-transcriptional gene silencing, and we observed high similarity between the orthologs of this gene found in V. corneae and Microsporidia sp. MB (Fig. 13a). RNA-dependent RNA polymerase was identified to include the RdRP conserved region (Fig. 13b, Supplementary Fig. 4) across all species.

Fig. 13

Analysis of the RNA-dependent RNA polymerase in Microsporidia sp. MB. a Phylogenetic analysis. Tree construction was done using IQTREE with 1000 bootstraps. Bootstrap support values are shown on each node and the branch length scale is indicated at the bottom of the tree. b The identified conserved domain, RdRP, is mapped onto the schematic representation of the multiple sequence alignment

Conclusion

Microsporidia, a group of intracellular parasites, have been found to secrete a range of proteins into infected host cells, potentially influencing host metabolic processes and molecular programs [98]. The genome of the microsporidian parasite E. cuniculi contains a smaller number of protein families, reflecting the streamlining of its genome due to its adaptation to an intracellular parasitic lifestyle, with most orthologs being the smallest among eukaryotes [111]. These genomes consist of a set of host-exposed proteins that are rapidly evolving and belong to large, species-specific gene families, suggesting a common strategy for microsporidia to interact with their hosts [112].

Microsporidia sp. MB, a symbiont of the Anopheles mosquito, exhibits significant potential to interfere with malaria transmission by effectively inhibiting the development of Plasmodium parasites within its host [45, 113]. This study analyzes the key components of the Microsporidia sp. MB genome and focuses on its gene content and structure. This report demonstrates the close relationship between Microsporidia sp. MB and V. corneae ATCC 50505 genomes among Clade IV of microsporidia consistent with previous reports on their taxonomic classification [45, 80, 113]. Notably, the symbiont has a larger genome size of 5.9 Mb compared to the rest of the species in the same clade, with the highest content of core BUSCO genes in the group. In silico methods have been used to annotate hypothetical proteins in V. corneae, revealing that many are involved in binding, enzymatic activity, and regulatory functions, despite its compact genome [91]. Here, a total of 2247 putative proteins were predicted and annotated. Most of these were further classified as enzymes involved in metabolic processes within the symbiont, notably similar to its sister V. corneae.

Microsporidia have been found to possess a range of metabolic pathways. These include polyamine metabolism, which is essential for their growth and development, and carbohydrate metabolism [114, 115]. Interestingly, some microsporidia have lost the ability to synthesize ATP and rely on host cell metabolism for energy [116]. Despite the absence of a recognizable mitochondrion, elements of mitochondrial metabolism, such as the pyruvate dehydrogenase complex, have been retained in these parasites [117]. Here, we analyzed the glycolytic pathway in Microsporidia sp. MB and compared this to other species in the Enterocytozoonida group. The Enterocytozoonida, apart from V. corneae ATCC 50505, have largely lost most of the proteins involved in this pathway except for hexokinase. This study demonstrates that Microsporidia sp. MB has retained all but one of the proteins involved in the glycolytic pathway complex like its closest relative, V. corneae. In addition to identifying a partial glycolytic pathway, which has retained 10 of the 11 key enzymes, we have also uncovered evidence for the functional presence of the pentose phosphate pathway (PPP). The PPP plays a compensatory role in generating NADPH for biosynthesis and redox balance, crucial for the survival of Microsporidia sp. MB in the absence of mitochondria. The interconnectedness of the glycolytic and PPP pathways allows this symbiont to adapt to its intracellular, host-dependent environment, supporting essential metabolic processes that contribute to its ability to block malaria transmission.

The Microsporidia sp. MB genome also showed evidence of ploidy through the identification of MRGs. Our analysis of the putative proteins shows that Microsporidia sp. MB harbors genes linked with a sexual mode of evolution (MRG). To further assert if this microorganism is a diploid or haploid, future analysis of the presence of regions with loss of heterozygosity (LOH) within this genome is necessary.

The study of RNA interference (RNAi) pathways in microsporidia highlights the evolutionary aspects of polyploidy. For instance, the retention or loss of Dicer and Argonaute orthologs across different microsporidian species suggests that these genes are selectively maintained or discarded during lineage divergence. This evolutionary flexibility might be influenced by polyploidy, affecting the parasite's adaptability and interaction with its host [79]. While some microsporidian species appear to retain elements of the RNAi pathway, the functional significance and mechanisms may vary considerably among different species. The presence of the Dicer, Argonaute and RNA-dependent RNA polymerase genes – the three key components of the RNAi pathway, in Microsporidia sp. MB is highlighted here. Despite having little information on the functioning of this complex in most microsporidia, a study of the bee pathogen N. ceranae has in the past confirmed its regulatory functioning. Future studies of the RNAi machinery in Microsporidia sp. MB promise to provide more insight into its regulatory role and effect on the survival of the symbiont in the host.

This study provides a foundational analysis of the genome of Microsporidia sp. MB, but several limitations should be considered when interpreting the results. First, the absence of specific genes, such as pyruvate kinase (pyk) or 6-phosphogluconolactonase (PGLS), may not necessarily indicate their true absence from the genome. These absences could be due to gaps in the genome assembly, low sequencing depth in certain regions, or sequence divergence that complicates gene identification. Additionally, functional predictions based on homology (e.g., from InterProScan and BLAST) are valuable for identifying potential protein functions but may not fully capture the biological roles of these proteins. Thus, experimental validation is needed to confirm the functionality of these predicted genes.

Furthermore, while our study highlights key metabolic pathways, such as glycolysis and the PPP, the relationship of these pathways in Microsporidia sp. MB remains speculative in some respects. The possibility of alternative metabolic pathways or compensatory mechanisms that were not detected in this genome assembly must be considered. Future studies utilizing more advanced genome sequencing technologies, such as long-read sequencing, and functional assays will be critical to address these gaps and provide a more comprehensive understanding of the biology of Microsporidia sp. MB.

Notwithstanding, annotating and analyzing the first reference genome of Microsporidia sp. MB is paramount in understanding its biology of infection and symbiosis within the Anopheles mosquito. The availability and analysis of these sequences act as a backdrop to important biological studies to understand the mechanism of this symbiont including optimizing the efficiency and stability of parasite transmission, understanding the ecological impact of introducing Microsporidia sp. MB into mosquito populations, and assessing any potential negative effects on non-target organisms. Altogether, this study lays the foundation for upcoming studies focusing on host-symbiont molecular interactions, strain variation analysis, and further implementation of RNASeq and other NGS tools in characterizing this important obligate intracellular organism.

Methods

A list of tools and resources used in this study are highlighted in Supplementary Table 1. Figure 14 further outlines the workflow for genome assembly, annotation, and comparative analysis of Microsporidia sp. MB adapted in this study.

Fig. 14

Workflow for genome assembly, annotation, and comparative analysis of Microsporidia sp. MB. This flowchart outlines the step-by-step process used in the genome annotation pipeline. The workflow begins with sample collection and sequencing of Anopheles arabiensis mosquito ovaries using DNBSEQ technology, followed by genome assembly using SPAdes. The quality of the assembled genome was assessed using QUAST and visualized with BlobToolKit to check for contamination. Genome completeness was evaluated using BUSCO by identifying conserved single-copy orthologs specific to microsporidia. Repeat sequences in the genome were detected using RepeatMasker and RepeatModeler2 to annotate known and novel repeats. Gene prediction was performed using GeneMark-ES, optimized for eukaryotes with intronless genomes. Functional annotation of predicted genes was carried out with InterProScan, assigning proteins to families and domains from databases such as Pfam and CDD. Phylogenomic analysis of the genome was conducted using OrthoFinder for ortholog identification and IQ-TREE for maximum likelihood phylogenetic tree construction. Finally, predicted proteins were mapped to biological pathways using KEGG, and physicochemical properties, including transmembrane regions, were characterized using ProtParam and DeepTMHMM. This comprehensive workflow provides a detailed overview of the computational tools used in the experimental design and their roles in the analysis of the Microsporidia sp. MB genome

Genome annotation – quality/completeness/repeats

Sample collection, preparation, genome sequencing, and genome assembly are described in [75]. Highly infected mosquito ovaries were sequenced using BGI DNB technology, followed by stringent quality assessment, de novo genome assembly and gene prediction. The quality of the genome was assessed using QUAST v5.0.2 [118]. Moreover, genome completeness estimation was confirmed using BUSCO v5.4.3 [119] against the microsporidia_odb10 set of 600 core genes. Repeat regions within the genome were identified and annotated using RepeatModeler2 [120] and RepeatMasker v4.0 [81]. BlobToolKit was used to visualize the genome structure [121].

Gene Prediction and Organization

GeneMark-ES v2.5p [122] was used for gene prediction with the intronless eukaryotic option. A comparison of BUSCO and GeneMark-ES output was used to determine common genes across both tools. Gene organizations across the contigs larger than or equal to 10 kb were visualized in RStudio 2023.09.1 + 494 [82] using the library package gggenes v0.5.1 [123].

Analysis of common orthogroups in Microsporidia sp. MB and the Enterocytozoonida group of Microsporidia

OrthoFinder v2.5.4 [124] and BUSCOPhylo v1.3 [125] was used to perform a BUSCO-based phylogenomic analysis of Microsporidia sp. MB and the clade IV microsporidia genomes. Nosema granulosis was used as an outgroup. Protein alignments were conducted using MUSCLE [83]. The maximum likelihood tree was inferred using IQTREE v1.6.12 [85] with the best-fit model LG + F + I + G4 chosen according to the Bayesian Information Criterion (BIC) using ModelFinder [86]. The tree was visualized using Geneious Prime® 2023.2.1 [126].

Protein characterization

Protein function analysis was done using the command line InterProScan tool [90] [accessed 04 October 2023] which includes annotations from multiple freely available tools including Pfam, CDD, Gene3D, HAMAP, MOBIDB, PANTHER, PIRSF, PRINTS, PROSITE, SFLD, SMART, SUPERFAMILY and NCBIFAM [127,128,129,130,131,132,133,134,135,136,137,138]. Gene Ontology (GO) analysis was performed by Blast2GO. The physicochemical properties of the putative proteins were determined using a batch ProtParam in-house script in Jupyter Lab [139]. Prediction of signal peptides and transmembrane regions was done using TargetP v2.0 [140], SignalP v6.0 [141], and DeepTMHMM v1.0.24 [103] respectively. The annotation of metabolic pathways was conducted using KEGG [102]. Visualization of the biological processes (BP), molecular functions (MF), cellular components (CC), and KEGG pathways was done using RStudio 2023.09.1 + 494 [82].

Homology modeling and motif analysis

Protein homology modeling was done using PRIMO (accessed 23 November 2023) as described in [99]. Briefly, HHPred [142] was used for the automatic identification of templates. This was followed by template selection where templates with higher resolution and optimal coverage were selected. Pairwise alignment of the target sequence to the template was done using T-COFFEE. Four models using the "Very slow" refinement method were designed. The quality of the models was assessed using multiple tools including the Qualitative Model Energy Analysis (QMEAN) Server, Protein Structure Analysis (ProSA-web) and Verify3D webserver [143,144,145]. Models with the best Z-dope score and high quality were selected.

LRR motifs were identified from the sequences using MEME and MAST online tools using default parameters [146, 147]. Briefly, MEME was set to identify motifs between 6 and 20 amino acids in length, with a focus on discovering the top three motifs across the sequences. The zoops model (zero or one occurrence per sequence) was applied for motif distribution. MAST was subsequently used to search for these motifs in the sequence database, using an E-value threshold of 10 to report significant matches.

Characterization of the glycolytic pathway

Microsporidia proteins involved in the glycolytic pathway were identified from past literature [46, 49, 76, 77, 148] and 12 Encephalitozoon cuniculi genes involved in this pathway were retrieved from NCBI. These included: pyruvate kinase (ECU09_0640), enolase (ECU10_1690), fructose 1,6-bisphosphate aldolase (ECU01_0240), phosphoglycerate mutase (ECU10_1060), hexokinase (ECU11_1540), pyruvate DHase E1, a sub (ODPA SCHPO) (ECU09_1040), pyruvate DHase E1, b sub (ODPB RICPR) (ECU04_1160), phosphofructokinase (ECU03_0680), phosphoglucose isomerase/glucose-6-phosphate isomerase (ECU05_0650), 3-phosphoglycerate kinase (ECU05_0320), glyceraldehyde-3-phosphate dehydrogenase (ECU07_0800) and triose phosphate isomerase (ECU11_0230).

A homology search of these sequences on NCBI BLASTp [149] was conducted (E-value = 0.001; bit-score ≥ 200) and the top 10 hits for each protein were retrieved. A reverse BLAST of the top 10 hits was used to confirm homology and all hits which included the query sequence were retained. A custom local database was created using the BLAST + makeblastdb command [150] and the retained sequences which included proteins from the Enterocytozoon consisting of V. corneae ATCC 50505, E. bieneusi, E. hepatopenaei, E. canceri, H. eriocheir and other microsporidia such as N. granulosis, V. ceranae, A. locustae, O. bayeri. Putative proteins from the genome assembly were blasted against this local database (E-value = 0.001; bit-score ≥ 200). Microsporidia sp. MB genes that were identified to be homologous to any sequence within the custom glycolysis database were selected and a further homology search against the remote BLASTp was done to confirm identity.

Multiple sequence alignment of each set of proteins was done using ClustalO. MEGAX was used to identify the best-fit model for phylogenetic analysis for each protein set. Tree calculations were done using the best model using ModelFinder and visualization was done in Geneious Prime® 2023.2.1.

Key conserved domains in the putative genes predicted to be involved in the glycolytic pathway were identified on the Batch Web CD-Search tool (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) (accessed 03 October 2023) (E-value = 0.001).

In addition to the partial glycolytic pathway identified in Microsporidia sp. MB, evidence suggests that the pentose phosphate pathway (PPP) may play a compensatory role in the organism’s metabolism. The PPP is closely linked to glycolysis and provides crucial metabolic intermediates while also generating NADPH, which is vital for biosynthetic reactions and maintaining redox balance, particularly in organisms like Microsporidia sp. MB that lack mitochondria.

Several key enzymes of the PPP were identified in the genome of Microsporidia sp. MB, located on specific scaffolds within the genome assembly (Table 5 below). This table lists the scaffolds (nodes) containing predicted genes involved in the pentose phosphate pathway. The table includes descriptions, sequence lengths, enzyme functions, Gene Ontology (GO) terms, and relevant enzyme activities. Key enzymes such as glucose-6-phosphate dehydrogenase (G6PD), 6-phosphogluconate dehydrogenase (PGD), transketolase (TKT), and ribose-5-phosphate isomerase (RPI), which contribute to both the oxidative and non-oxidative branches of the PPP, are highlighted. These enzymes play essential roles in generating NADPH and intermediates necessary for biosynthesis and redox balance in Microsporidia sp. MB.

Table 5 Pentose phosphate pathway (PPP) genes identified in Microsporidia sp. MB

Notably, glucose-6-phosphate dehydrogenase (G6PD), which catalyzes the first step of the oxidative branch of the PPP, was found on NODE_723 (GO:0004345). This enzyme converts glucose-6-phosphate into 6-phosphoglucono-δ-lactone, producing NADPH in the process, a critical function for maintaining cellular redox balance [151]. Additionally, 6-phosphogluconate dehydrogenase (PGD), responsible for converting 6-phosphogluconate to ribulose-5-phosphate in the oxidative branch, was identified on NODE_1158 (GO:0004616). This enzyme further supports the generation of NADPH for biosynthesis and oxidative stress response.

The non-oxidative branch of the PPP also appears to be functional in Microsporidia sp. MB. Transketolase (TKT), found on NODE_1037 (GO:0004802), facilitates the transfer of two-carbon units between sugar phosphates, regenerating glucose-6-phosphate and providing intermediates for nucleotide biosynthesis. Additionally, ribose-5-phosphate isomerase (RPI), located on NODE_1109 (GO:0004751), catalyzes the interconversion of ribose-5-phosphate and ribulose-5-phosphate, a critical step for ensuring the availability of ribose for nucleotide synthesis [152].

However, 6-phosphogluconolactonase (PGLS), which is responsible for converting 6-phosphoglucono-δ-lactone into 6-phosphogluconate, was not identified in the current genome annotation. This gap suggests either incomplete annotation or the possibility that Microsporidia sp. MB uses alternative pathways to process this intermediate. Despite this, the presence of key PPP enzymes indicates that the pathway is likely functional and plays a compensatory role in Microsporidia sp. MB’s metabolism, supplementing the partial glycolytic pathway.

Together, the glycolytic and pentose phosphate pathways form an interconnected system that enables Microsporidia sp. MB to adapt to its host-dependent intracellular environment. These pathways provide metabolic flexibility, ensuring the production of essential biosynthetic precursors and redox cofactors that are crucial for survival [95].

Inference of the putative mode of reproduction in Microsporidia sp. MB: characterization of Meiotic-Related Genes (MRG)

A set of MRGs were identified from a vast literature search using keywords [MRG + Microsporidia] on Google Scholar (https://scholar.google.com/ accessed 3 October 2023). Nineteen microsporidian MRGs were documented [106] and an NCBI search was used to retrieve these sequences. A custom database was created using the retrieved sequences and a local blast search of these genes against the custom database was conducted as described above (E-value = 0.0001). Multiple sequence alignments and tree calculations were done as described in the section above.

Analysis of the RNA interference mechanism in Microsporidia sp. MB

Microsporidia exhibit oxidative stress response mechanisms, including the expression of antioxidant enzymes, which help neutralize reactive oxygen species generated within the host environment. The Dicer, RNA-dependent RNA polymerase, and Argonaute genes have been well-documented to be involved in this system [79]. Their respective orthologs were retrieved from MicrosporidiaDB [153] and homologs of the RNAi complex were retrieved from NCBI (accessed 4 October 2023). Search and analysis of the homologs from the Microsporidia sp. MB genome was done as described in the sections above. Moreover, these genes are known to contain unique conserved domains including PAZ Argonaute-like, DUF, and PIWI-like domains. Therefore, a CD-Search on the Microsporidia sp. MB genes homologous to Argonaute and Dicer genes was done using default parameters to determine the presence of these key domains (E-value = 0.001).