Role of MicroRNA-Like RNAs in the Regulation of Spore Morphological Differences in the Entomopathogenic Fungus Metarhizium acridum

Entomopathogenic fungi are some of the most promising agents for the biocontrol of insect pests. Some species have been used for pest control as biopesticide products. As of 2017, over 200 products based on entomopathogenic fungi and nematicidal fungi were registered for use against various pests (Kumar et al. 2019). As a microbial pesticide, Metarhizium acridum is widely used for locust and grasshopper control in Asia, Africa, and Australia, and its use is based on aerial conidia (Hunter et al. 2001; Lomer et al. 2001; Peng et al. 2008). Conidiation includes a period of vegetative growth (Bosch and Yantorno 1999) and is involved in conidial germination, mycelium formation, thick-walled foot cell formation at the tip of the aerial mycelium, and the production of multinuclear conidiophores.

In contrast with aerial conidia (CO), blastospores (BS) are produced by hyphal constriction, separation at the septa, and yeast-like budding. In insect mycology, BS are commonly referred to as any hyphal bodies produced in insect blood by hyphal budding, so BS are also considered hyphal bodies (Fargues et al. 2002). CO and BS are two different sporulation patterns that occur in the life cycle of M. acridum, and they are significantly different in cell morphology, structure, and activity. Conidiation is of crucial importance because conidia are infectious propagules and active components in mycoinsecticides. However, the high production cost and poor conidia production efficacy have retarded their application as mycoinsecticides (Lacey et al. 2001; Hajek et al. 2007). Thus, the BS were studied because of the lower cost of using industrial fermenting tanks (Adámek 1963). However, their thinner cell walls and lower tolerance to environmental stress and storage stability limit their application in the field (Pereira and Roberts 1990). Although BS have several drawbacks, they have other advantages, such as a fast germination rate and higher spore activity and virulence, so BS are usually used for seeds in solid-state fermentation (Hegedus et al. 1992; Faria and Wraight 2007). Most previous studies on the difference between CO and BS have focused on visual characteristics, virulence, storage conditions, and field application (Hegedus et al. 1992; Stephan et al. 1997; Leland et al. 2005; Wassermann et al. 2016). Nevertheless, little is known about the difference between the morphologies of CO and BS in entomopathogenic fungi on the transcriptome levels.

MicroRNAs (miRNAs) are a class of small non-coding RNA molecules that are approximately 22 nucleotides (nt) in length, which can partially or entirely bind to the 5’-UTR, 3’-UTR, or coding region of target genes to down- or up-regulate the expression of target genes (Grey et al. 2010; Fang and Rajewsky 2011; Reczko et al. 2012; Helwak et al. 2013; Hussain et al. 2013; Asgari 2014). These miRNAs are pervasive in animals and plants and act as posttranscriptional regulators that specifically guide target gene recognition to regulate gene transcriptional start or repression (Carthew and Sontheimer 2009). In animals, miRNAs have been shown to play important roles in cell development, proliferation, and differentiation. Targeted silencing of miRNA-132-3p expression is an advantage for bone marrow-derived mesenchymal stem cells (BMSC) osteogenic differentiation and osteogenesis (Hu et al. 2020). The overexpression of miRNA-324-5p exerts cell growth and migration-promoting effects through activating Wnt signaling pathway and epithelial to mesenchymal transition (EMT) by negatively regulatory suppressor of fused gene (SUFU) in gastric carcinomas (Peng et al. 2020). Chi-miR-199a-5p inhibited TGF-β2 expression at both mRNA and protein translation levels in fibroblasts (Han et al. 2020). In plants, they also play various roles in plant development, stress response, and antibacterial resistance. MiRNA157 regulated floral organ growth and ovule production in Gossypium hirsutum by negatively regulatory promoter-binding protein-like (SPL) gene (Liu et al. 2017), and miRNA156 could increase tolerance to heat stress by downregulating promoter-binding protein-like (SPL) gene in Arabidopsis (Stief et al. 2014). In fungi, miRNA-like RNAs (milR-NAs), with similar characteristics to miRNAs in animals and plants, were discovered in Neurospora crassa. They are produced by at least four diverse pathways that use a distinct combination of factors, including QDE-2, Dicers, the exonuclease QIP, and an RNase III domain-containing protein, MRPL3 (Lee et al. 2010). Subsequently, milRNAs in various species of fungi were discovered, such as Cordyceps militaris, Metarhizium anisopliae, Magnaporthe oryzae, Aspergillus fumigatus, Aspergillus flavus, and Penicillium marneffei (Zhou et al. 2012; Lau et al. 2013; Özkan et al. 2017; Li et al. 2020).

In recent years, many studies have reported that milRNAs can affect the diverse physiological process in organisms, such as cell growth, development, virulence, metamorphosis, and metabolism (Mukherjee and Vilcinskas 2014; Ylla et al. 2017; Zhang et al. 2018; Guo et al. 2020). In P. marneffei, milRNAs regulated the growth process of mycelial and yeast phases (Lau et al. 2013). Ssc-milR-240 was shown to potentially regulate sclerotial development by epigenetic regulation of its target histone acetyltransferase in Sclerotinia sclerotiorum (Xia et al. 2020). Overexpression of milRNA-87 exhibited a dramatic increase in the growth, conidiation, and virulence of Fusarium oxysporum f. sp. cubense by silencing target gene (glycosyl hydrolase coding gene, FOIG_15013) expression (Li et al. 2022). In Verticillium dahlia, milRNA-1 regulated the virulence by binding to the 3’-UTR of a hypothetical protein-coding gene (VdHy1) for transcriptional repression (Jin et al. 2019). In Fusarium graminearum, milRNA-2 combined with 3’-UTR of bioH1 involved in biotin biosynthesis to regulate biotin synthesis (Guo et al. 2019). MilR236 was shown to regulate appressorium formation and virulence of M. oryzae by binding to MoHat1, a histone acetyltransferase type B catalytic subunit (Li et al. 2020), and miR4 and miR16, which are involved in mycelium growth and sexual development in C. militaris (Shao et al. 2019). Hence, milRNAs may play an important role in M. acridum. The study of milRNAs and the expression profiles may help us better understand the roles in spore morphological differences in M. acridum.

M. acridum is an unusually effective model for studying spore morphological differentiation. Previous research has mainly focused on the cell structure and application of CO and BS in filamentous fungi. Nevertheless, very little information is known about the roles of milRNAs in spore morphological differences. Thus, in the present study, the cDNA and small RNA libraries from CO and BS of M. acridum were sequenced, and the DEMs and DEGs between CO and BS samples were screened to elucidate the biological function of milRNAs in spore morphological differences in M. acridum. Our study aimed to provide primary data for further research on spore morphological differences in M. acridum.

mRNA expression profiles. To explore the expression patterns and co-expression of differentially expressed genes (DEGs) in different spore types, i.e., conidia (CO) and blastospores (BS), cDNA libraries from BS and CO were sequenced. As shown in Table I, 567,108,002 reads were obtained from six cDNA libraries after filtering, including 287,099,910 and 280,008,092 from BS and CO, respectively. More than 94.06% of raw reads had Q-phred scores at the Q30 level (an error threshold of less than 0.01%), the BS and CO reads had approximately 46.96% GC content, and more than 64.94% of clean reads were mapped to the genome.

Gene expression was normalized to FPKM (Fragments Per Kilobase of Exon Per Million Fragments Mapped) in the present study. A total of 9,692 expressed genes were found in BS and CO in M. acridum, and the gene expression distribution was obtained in different samples (Fig. 1A). A total of 4,646 DEGs were found between BS and CO, containing 2,640 upregulated and 2,006 downregulated genes, with the criteria as follows: p-values ≤ 0.05 and fold change (FC) ≥ 2 (Fig. 1B). A correlation analysis showed that M. acridum exhibited a significant difference between BS and CO, indicating that there was transcriptional strong differentiation between BS and CO (Fig. 1C).

Fig. 1

Fig. 1

Functional analysis of the differentially expressed genes (DEGs). To better reveal the function, roles, and biological processes of DEGs in different spore types, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to evaluate the function of 4,646 DEGs using DAVID (the Database for Annotation, Visualization and Integrated Discovery, v. 6.8, https://david.ncifcrf.gov) GO enrichment analyses included categorization into the biological process (BP), cell component (CC) and molecular function (MF). A total of 2,951 DEGs were annotated in 54 GO classes at GO level 2 (Fig. 2A). For BP analyses, most genes were involved in biological (1,919), metabolic (1,750), and single-organism processes (1,494). When GO classification was based on CC, most genes were assigned to the cell (2,094), cell part (2,089), and organelle (1,596). Among the MFs, catalytic activity (1,805) was the most commonly represented, followed by binding (1,571).

Fig. 2

According to KEGG enrichment analyses, a total of 1,076 DEGs were subdivided into 256 KEGG pathways (Table SII). The KEGG database classified these genes into 33 pathways at the second level of classification. Among which the most significantly enriched pathways were carbohydrate metabolism (193), amino acid metabolism (183), translation (171), and global and overview maps (164) (Fig. 2B). The top 20 of 33 significantly enriched pathways showed that the most significantly enriched pathways were biosynthesis of amino acids, carbohydrate metabolism, ribosome, and oxidative phosphorylation (Fig. 2C). These results suggested that the DEGs were mainly involved in metabolism of BS and CO, which might be associated with cell activity and structure.

Analysis of milRNA sequences. To explore the regulation of gene expression by milRNAs in different spore types, i.e., BS and CO, six small RNA libraries were sequenced. A total of 81,627,529 raw reads were generated from six small RNA libraries, and after filtering, 67,072,593 clean reads were obtained from BS and CO samples. Clean reads were aligned against the small RNA database and annotated, and an overview of the small RNA classification annotation results is shown in Table II.

Identification of conserved milRNAs and novel milRNAs in M. acridum. The known milRNAs were identified by alignment against the miRBase v. 21 database. A BLAST search identified 2,350 conserved mil-RNAs in M. acridum (Table SIII). Of these conserved milRNAs, 113 were differentially expressed in BS vs. CO, including 12 upregulated milRNAs and 101 downregulated milRNAs (Fig. 3A and 3B). The length of these milRNAs ranged from 17 nt to 25 nt and were most commonly 22 nt (Fig. 3C). They were classified into 359 conserved milRNA families, and let-7, miR-21, miR-10, miR-30 and miR-26 were the most abundant known milRNA families. The novel milRNAs were predicted by miRDeep2 software. A total of 28 novel milRNAs with lengths between 18 and 25 nt were obtained (Table SIII). Novel 13 mature, novel 4 mature and novel 6 mature milRNAs were mostly enriched novel milRNAs.

Fig. 2

Prediction of differential milRNA target genes and functional annotation. To understand the functions and roles of differential milRNAs target genes, the prediction of target genes was performed using the software miRanda with the following parameters: S ≥ 150, ΔG ≤ –30 kcal/mol, and demand strict 5’ seed pairing. In a comparison between BS and CO, 493 target genes were identified for 54 DEMs (Table SIV), whereas no targets were identified for the other 59 milRNAs. The results indicated that most milRNAs regulated more than one target gene, and different milRNAs could regulate the same target gene. Through GO enrichment analysis, the DEMs target genes were annotated in 43 Gene Ontology (GO) terms. The GO terms with the most significant numbers of enriched target genes were cellular process (261), metabolic process (193), cell (183), cell part (183), catalytic activity (146) and organelle (133) (Fig. 4). According to Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, the 20 most enriched pathways of target genes were mainly RNA transport, purine metabolism, pyrimidine metabolism, glycine, serine and threonine metabolism, and cell cycle (Fig. 5), which were the critical pathways in genetic processing, transcription, posttranscriptional regulation and metabolism.

Fig. 2

Fig. 3

Fig. 3

Fig. 4

Fig. 5

Integration analysis of the milRNAs and mRNAs. To investigate the relationship between milRNAs and target genes, a potential milRNA-mRNA network that might affect spore morphology and structure was established. We speculated that some pathways might be associated with spore activity and characteristics, including transcriptional and posttranscriptional, cellular membrane and wall integrity, cell division, and cellular osmotic pressure. The relationship between milRNAs and target genes is shown in Fig. 6. A total of 18 milRNAs, including two up-regulated milRNAs and 16 down-regulated milRNAs targeted 42 DEGs that included 22 up-regulated genes and 20 downregulated genes. Of the 22 up-regulated target genes, 20 target genes corresponded to downregulated milRNAs, two target genes corresponded to two up-regulated milRNAs. The other 20 down-regulated target genes corresponded to 12 down-regulated milRNAs. Five down-regulated milRNAs corresponded to both up-regulated and down-regulated target genes. Based on these results, we conclude that several miRNA-mRNA pairs indicated negative regulation, and some milRNA-mRNA pairs showed the same trend in expression difference. In transcriptional and post-transcriptional, milRNAs targeted 16 DEGs, of which 10 genes were up-regulated in the BS vs CO, such as the mmu-miR-337-5p and mmu-miR-341-3p targeted MAC_07490 and MAC_06220, respectively. In a cellular membrane and wall integrity, milRNAs targeted 13 DEGs, of which 11 genes were down-regulated in the BS vs. CO, such as bta-miR-328 targeted MAC_09224, Chi-miR-330-5p targeted MAC_02770, and MAC_04818, pal-miR-370-3p targeted MAC_00296. In the cell division pathways, milRNAs targeted 10 DEGs that were significantly up-regulated, while, in cellular osmotic pressure, milRNAs targeted three DEGs that were significantly down-regulated. For all milRNA and mRNA pairs, most genes participated in transcriptional and posttranscriptional, and cell division pathways were significantly up-regulated, while that participated in cellular membrane and wall integrity and cellular osmotic pressure were significantly down-regulated. The results indicated that milRNAs might play an important role in cell growth and cellular morphological changes.

Fig. 6

Validation of milRNA and target gene expression with qRT-PCR. We selected six DEMs (bta-miR-150, bta-miR-328, cgr-miR-132, efu-miR-30b, efu-miR-7a, and sly-miR482c) and six DEGs (MAC_02738, MAC_03363, MAC_03964, MAC_07508, MAC_07780, and MAC_08938) to validate the accuracy of the milRNA and mRNA sequencing data using qRT-PCR. As shown in Fig. 7, all expression patterns showed similar trends. The results indicated that milRNA and mRNA sequencing data were reliable and partially validated the reliability of our findings in this study.

Fig. 7

With people expanding their perception of the environment and health, mycoinsecticides have become increasingly important substitutes for chemical insecticides due to their low toxicity, target specificity, and harmlessness toward non-target organisms. As a microbial pesticide, M. acridum is widely used for locust and grasshopper control in Asia, Africa, and Australia, which is based on aerial conidia (Hunter et al. 2001; Lomer et al. 2001; Peng et al. 2008). CO and BS are two types of spores that occur in different patterns during the life cycle of M. acridum, and they have significant differences in cell morphology, structure, and activity. While most previous studies on the differences in BS and CO have focused on visual characteristics and field application, a few works have analyzed the molecular mechanism underlying the differences between BS and CO in fungi.

In this study, we integrated milRNA and mRNA data to identify and investigate the roles of milRNAs in spore morphological differences in the entomopathogenic fungus M. acridum. 4,646 DEGs were obtained between BS and CO, including 2,640 up-regulated and 2,006 down-regulated genes. GO enrichment analysis showed that the most significantly expressed genes were classified as a cellular process, metabolic process, cell, cell part, and catalytic activity. Studies have shown that BS has a thinner cell wall and lower stability and is maintained for less time than CO (Pereira and Roberts1990). The cell wall is a vital structure for fungal cells that protect against various environmental stresses, such as heat shock, UV-B irradiation, oxidation, and desiccation. Forty cell wall glycoproteins in the filamentous fungus Neurospora crassa were identified by proteomic analyses, and the major cell components included chitin, β-1,3-glucan, mixed β-1,3-/β-1,4-glucans, glycoproteins, and melanin (Patel and Free 2019). Studies have reported that O-mannosyltransferase and β-1,3-glucanosyltransferase mutants reduced fungal tolerance to heat shock and UV-B irradiation, formed thinner cell walls, and significantly reduced total sugar and β-1,3-glucan in entomopathogenic fungi (Luo et al. 2018; Zhao et al. 2019).

In this study, the transcriptome analysis showed that mannosyltransferase (MAC_03068, MAC_08610, MAC_05793, MAC_09434, and MAC_06369) and cell wall glucanosyltransferases Mwg1 and Mwg2 (MAC_02745 and MAC_01181) were significantly down-regulated in BS compared to CO. It might be connected to the thinner cell wall, lower stability, and shorter maintenance of BS in M. acridum. In M. acridum, β-1,3-glucan and chitin are the major polysaccharides. Synthase and hydrolytic enzymes are critical during the synthesis, branching, and cross-linking of polymers (Adams 2004). Interestingly, the expression levels of β-1,3-glucan synthase were not different between BS and CO.

In contrast, the two chitin synthases MaChsII and MaChsVI were significantly differentially expressed, down-regulated and up-regulated in BS compared to CO, respectively. Previous studies identified seven chitin synthases in M. acridum, and the chitin synthase family influences M. acridum growth, stress tolerance, cell wall integrity, and virulence.

Regarding these two chitin synthase genes, ΔMaChsII mutants germinated more rapidly than WT, and MaChsVI participated in the chitin synthesis (Zhang et al. 2019). Chitinases and glucanases have critical roles during cell separation, cell wall modification and cell wall remodeling (Adams 2004). Seventeen and 15 differentially expressed chitinase (three down-regulated and 14 up-regulated), and glucanase (nine down- and six up-regulated) genes were identified, respectively, in this study, which might be associated with morphogenesis. KEGG pathway enrichment analysis showed that the DEGs were enriched in the biosynthesis of amino acids, carbohydrate metabolism, ribosome, and oxidative phosphorylation, resulting in a higher spore activity in BS than in CO. In addition, hydrophobin-like protein ssgA (MAC_07330) was significantly downregulated (log₂fold change = –6.36), presumably due to alterations in spore surface structures that resulted in BS hydrophilicity.

MiRNAs are a class of key regulatory factors that can partially or entirely bind to the 5’-UTR, 3’-UTR, or coding region of target genes to down- or up-regulate the expression of target genes (Grey et al. 2010; Fang and Rajewsky 2011; Reczko et al. 2012; Helwak et al. 2013; Hussain et al. 2013; Asgari, 2014). In this experiment, 2,350 conserved milRNAs were identified in the BS and CO samples of M. acridum. Of these conserved milRNAs, 113 milRNAs were differentially expressed in BS vs. CO, including 12 up-regulated milRNAs and 101 down-regulated milRNAs. However, only 54 DEMs, targeting 493 genes were obtained. The other 59 milRNAs were not targeted, which suggested that most fungal milRNAs were imprecisely complementary to their target genes or that the database had limitations. A similar phenomenon was also found in animals, N. crassa and M. anisopliae (Ambros 2004; Lee et al. 2010; Zhou et al. 2012).

GO classification and KEGG enrichment analyses of the target genes demonstrated that they participated in various essential genetic information processes and metabolic processes, such as translation, carbohydrate metabolism, and nucleotide metabolism, which was consistent with transcriptome analyses, indicating the potential roles of milRNA in cell morphology, structure, and activity.

Previous studies have suggested that a miRNA may not be a one-to-one target gene. A miRNA could target hundreds of genes, and an mRNA could be regulated by multiple miRNAs (O’Day and Lal et al. 2010). A similar phenomenon was found in this study; 34 DEMs were identified to have more than one target, such as btamiR-328 targeting 65 mRNAs and mmu-miR-341-3p targeting 83 mRNAs. Previously, miRNAs acted as negative post-transcriptional regulators that guided target gene recognition to regulate gene transcriptional start or repression (Carthew and Sontheimer 2009). Our study found that 21 DEGs (16 up-regulated and five down-regulated) displayed the opposite trend as their corresponding milRNAs at the expression level. These milRNA-mRNA pairs suggested possible negative regulation. Apart from these, 23 DEGs (22 down-regulated and one up-regulated) showed the same trend, and a similar phenomenon was found in Trichophyton rubrum (Wang et al. 2018). Furthermore, we found 14 milRNAs corresponding to up- and down-regulated target genes, such as down-regulated pal-miR-370-3p, which targeted 18 downregulated mRNAs and 25 upregulated mRNAs.

Previous reports have suggested that miRNAs interact with transcription factors (TFs) to regulate gene expression (Pitto et al. 2008). Our data demonstrated that five DEMs targeted six TFs in the BS vs. CO stages. For example, mmu-miR-341-3p targeted zinc finger transcription factor ace1 and transcription factor TFIIA complex subunit Toa1, while pal-miR-370-3p and novel3_mature targeted HLH transcription factor and C6 transcription factor, respectively. In addition to six TFs, seven DEGs involved in transcription and transcriptional regulation were identified in the BS vs. CO stages, including RNA-directed 5’–3’ RNA polymerase and eukaryotic translation initiation factor 3 subunit EifCa. These results indicated that milRNAs regulated gene expression at the transcription level along with TFs.

Subsequently, the milRNAs that might be involved in the control of target gene expression, related to cell morphology, structure, and activity were analyzed. Among all the DEMs, 18 associated with transcription, cell proliferation, cell wall, member integration, and cell osmotic pressure were selected. As shown in the milRNA-mRNA network, dno-miR-328-3p, novel18_mature, age-miR-127, pal-miR-370-3p and bta-miR-328 regulated eight target genes that participated in cell division and growth. In our study, dno-miR-328-3p, age-miR-127, pal-miR-370-3p and bta-miR-328 were down-regulated, and novel18_mature was up-regulated in the BS stage compared with the CO stage. For example, the dno-miR-328-3p target genes: adenosine deaminase, glycosyl hydrolase and G-patch domain protein were up-regulated. Adenosine deaminase is a well-characterized enzyme involved in the depletion of adenosine, and as an important growth factor, adenosine deaminase could participate in cell differentiation or proliferation (Maier et al. 2005; Sekiya et al. 2013).

Previous studies have suggested that glycosyl hydrolase is involved in fungal morphogenesis and bacterial cell division (Kim et al. 2002; Yakhnina and Bernhardt 2020). The age-miR-127 target gene WD repeat-containing protein pop1 was up-regulated, and WD repeat-containing protein pop1 was a component of ubiquitin-mediated proteolysis that plays a crucial role in the control of the cell cycle (Hershko 1997).

Thus, these related milRNAs might play an essential role in cell differentiation or proliferation, suggesting a probable reason of BS higher activity than CO. In addition to these milRNAs related to cell morphology, the analysis of the correlation between the expression of the milRNAs and the mRNAs showed that 11 mil-RNAs regulated the expression of 13 genes associated with the cell membrane and wall integrity. For example, cgr-miR-598 targeted the GPI-anchored wall transfer protein gene, which was down-regulated in the BS stage compared with the CO stage. GPI-anchored wall transfer protein is involved in cell wall construction and remodeling in Saccharomyces cerevisiae and likely has a role in the integration of chitin polymer within the cell wall matrix (Rodriguez-Peña et al. 2002; Kapteyn et al. 1999). The chi-miR-330-5p targeted phosphate transporter gene (Pho88) was also downregulated in the BS stage, and this phosphate transporter is involved in inorganic phosphate transport and mediates lipid accumulation (James and Nachiappan 2014). These results indicated that milRNAs might be involved in the cell membrane and wall integrity and play a critical role in M. acridum, suggesting a probably regulated function that gave BS a thinner cell wall and higher activity than conidia. Analysis of the results suggested that the milRNAs play critical roles in cell growth, cellular morphological changes, and metabolism in the entomopathogenic fungus M. acridum.

In summary, we first provided insight into the milRNA-mRNA relationship based on a comparison between BS and CO stages and provided useful information on the milRNAs involved in cellular morphological differences of M. acridum. One hundred thirteen mil-RNAs showed altered expression in the BS stage compared to the CO stage. The target genes of milRNAs were classified into a wide range of functional categories, especially those related to cellular processes, metabolic processes, and the cell cycle. These results suggested that milRNAs may play a critical role in cell growth, cellular morphological changes, and metabolism. The transcriptomics data also provided a foundation for further studies aimed at understanding the functions of milRNAs.