milR20 negatively regulates the development of fruit bodies in Pleurotus cornucopiae

The mechanism underlying the development of fruit bodies in edible mushroom is a widely studied topic. In this study, the role of milRNAs in the development of fruit bodies of Pleurotus cornucopiae was studied by comparative analyses of the mRNAs and milRNAs at different stages of development. The genes that play a crucial role in the expression and function of milRNAs were identified and subsequently expressed and silenced at different stages of development. The total number of differentially expressed genes (DEGs) and differentially expressed milRNAs (DEMs) at different stages of development was determined to be 7,934 and 20, respectively. Comparison of the DEGs and DEMs across the different development stages revealed that DEMs and its target DEGs involved in the mitogen-activated protein kinase (MAPK) signaling pathway, protein processing in endoplasmic reticulum, endocytosis, aminoacyl-tRNA biosynthesis, RNA transport, and other metabolism pathways, which may play important roles in the development of the fruit bodies of P. cornucopiae. The function of milR20, which targeted pheromone A receptor g8971 and was involved in the MAPK signaling pathway, was further verified by overexpression and silencing in P. cornucopiae. The results demonstrated that the overexpression of milR20 reduced the growth rate of mycelia and prolonged the development of the fruit bodies, while milR20 silencing had an opposite effect. These findings indicated that milR20 plays a negative role in the development of P. cornucopiae. This study provides novel insights into the molecular mechanism underlying the development of fruit bodies in P. cornucopiae.

The mechanism underlying the development of fruit bodies in edible mushroom is a widely studied topic. In this study, the role of milRNAs in the development of fruit bodies of Pleurotus cornucopiae was studied by comparative analyses of the mRNAs and milRNAs at different stages of development. The genes that play a crucial role in the expression and function of milRNAs were identified and subsequently expressed and silenced at different stages of development. The total number of differentially expressed genes (DEGs) and differentially expressed milRNAs (DEMs) at different stages of development was determined to be 7,934 and 20, respectively. Comparison of the DEGs and DEMs across the different development stages revealed that DEMs and its target DEGs involved in the mitogen-activated protein kinase (MAPK) signaling pathway, protein processing in endoplasmic reticulum, endocytosis, aminoacyl-tRNA biosynthesis, RNA transport, and other metabolism pathways, which may play important roles in the development of the fruit bodies of P. cornucopiae. The function of milR20, which targeted pheromone A receptor g8971 and was involved in the MAPK signaling pathway, was further verified by overexpression and silencing in P. cornucopiae. The results demonstrated that the overexpression of milR20 reduced the growth rate of mycelia and prolonged the development of the fruit bodies, while milR20 silencing had an opposite effect. These findings indicated that milR20 plays a negative role in the development of P. cornucopiae. This study provides novel insights into the molecular mechanism underlying the development of fruit bodies in P. cornucopiae. KEYWORDS milR20, fruit body development, Pleurotus cornucopiae, comparative transcriptome, MAPK signaling pathway

Introduction
MicroRNAs (miRNAs) are small non-coding RNA molecules that are 18-24 nt long, play important regulatory roles in gene regulation, and influence various biological processes in plants and animals. Numerous miRNAs of plants and animals have been identified to date. The first miRNA-like fungal RNAs (milRNAs) were discovered in Neurospora crassa (Lee et al., 2010) and the milRNAs have been subsequently identified in other filamentous fungi and basidiomycetes (Kang et al., 2013;Wang et al., 2021). The characteristics of fungal milRNAs are similar to those of plant and animal miRNAs. For instance, the miRNA precursors of plants and animals have a typical hairpin structure that is similar to those of fungi. However, the biosynthesis mechanism of fungal miRNAs are different from the animal and plant miRNAs. The milRNAs of N. crassa are produced by four different mechanisms that include a distinct combination of factors, including Dicers, Argonaute protein QDE-2, the exonuclease QIP, and the RNAse III domain-containing MRPL3 protein. While the miRNA in animals and plants were produced by Dicer-like enzymes or Drosha proteins in miRNA maturation (Jones-Rhoades et al., 2006;Lee et al., 2010).
The majority of recent studies on miRNAs are primarily focusing on the miRNAs of animals or plants, and there is a scarcity of research on fungal miRNAs. The miRNA-mediated post-transcriptional regulation of genes is a novel gene regulation strategy that is used to regulate the expression of protein-coding genes by targeting mRNAs via cleavage or translational repression. Our current understanding of target recognition by miRNAs suggests that the mRNA sequence is complementary to bases 2-8 of miRNAs (referred to as the seed sequence) in the majority of miRNA-mediated silencing complexes (miRISCs). The seed region has the highest complementarity to the 3′-untranslated region (UTR) of the mRNA of the target gene, and previous studies have demonstrated that various miRNAs with numerous functions regulate multiple target genes via different mechanisms (Kiriga et al., 2020). The recent advancements in sequencing technologies and bioinformatics tools have facilitated the identification of milRNAs in various fungi Mu et al., 2015;Li et al., 2016). However, there is a scarcity of information regarding the functions and target recognition mechanism of fungal milRNAs.
Understanding the regulatory mechanisms of fungal milRNAs may aid in breeding novel varieties of edible mushroom. Pleurotus cornucopiae is one of the most extensively cultivated mushrooms in China, and has a high nutritional and medicinal value. The mechanism underlying the development of fruiting bodies is a complex process that is regulated by both genetics and environment, and has been a topic of immense interest in recent years. Numerous genes that play an important role in the development of mushrooms have been identified and characterized. For instance, it has been demonstrated that the genes that encode the FvHmg1 and LFC1 transcription factors negatively regulate the fruit body development of Flammulina velutipes Meng et al., 2021). Additionally, genes encoding SsNox2 NADPH oxidases contribute to the generation of reactive oxygen species (ROS), which are essential for the sclerotia development of Sclerotinia sclerotiorum (Kim et al., 2011). Glutathione peroxidase (GPX), which aids in maintaining ROS homeostasis, has a complex influence on the growth of the filamentous fungi Hypsizygus marmoreus . Additionally, the genes encoding protein kinases in the mitogen-activated protein kinase (MAPK) signaling pathway play an important role in cellular regulation in fungi by regulating phosphorylation and dephosphorylation. A previous study on Metarhizium robertsii revealed that the MAPK signaling cascade plays a regulatory role in the formation of conidia and stress tolerance (Chen et al., 2016). The SakA response factor of Aspergillus nidulans can transmit osmotic and oxidative stress signals to the MAPK signaling pathway and regulate the growth, development, and stress response of A. nidulans (Lara-Rojas et al., 2011). Another study reported that the adenosine cyclase of the cyclic adenosine monophosphate (cAMP) signal transduction pathway aids in the transformation of yeast morphology to mycelial morphology, and plays a crucial role in mycelial growth (Rocha et al., 2001).
Recent studies have demonstrated that small RNAs play vital roles in fungal development. For instance, it has been demonstrated that milR4 and milR16 mediate the development of fruiting bodies in Cordyceps militaris. The disruption of milR4 results in the non-formation of fruiting bodies while the disruption of milR16 results in the formation of abnormal fruiting bodies with pale yellowcolored primordia (Shao et al., 2019). In contrast, the overexpression of Po-MilR-1 in P. ostreatus results in slow mycelial growth and formation of abnormal pilei with irregular edges (Xu et al., 2021). However, only one milRNA of P. ostreatus that plays a vital role in mycelial growth has been identified and purified to date. Therefore, the potential roles of milRNAs in the development of fruiting bodies of P. cornucopiae is poorly understood to date owing to the scarcity of information, and further studies are necessary in this regard.
In this study, the genes encoding the Dicer, argonaute (AGO), and RNA-dependent RNA polymerase (RDRP) proteins of P. cornucopiae were identified and their expression profiles were determined at different developmental stages. The milRNAs related to the development of P. cornucopiae were determined by small RNA sequencing and in silico analyses. The potential targets of the milRNAs in the genome of P. cornucopiae were additionally detected, and the expression and functions of these target genes were determined by transcriptome sequencing and bioinformatics analyses. The theoretically predicted milRNAs and their gene targets were experimentally validated by quantitative real-time PCR (qRT-PCR) and dual-luciferase activity assay. Finally, one milRNA of P. cornucopiae was identified, and its functions in the development of fruiting bodies of P. cornucopiae were determined by overexpression and silencing. The results are anticipated to provide a foundation for research on milRNA function and the application of milRNAs in the development of edible mushrooms.

Strains and media
The CCMSSC 00406 strain of P. cornucopiae was obtained from the China Center for Mushroom Spawn Standards and Control. The fungal mycelia were inoculated on potato dextrose agar (PDA) at 28°C for 6 days. The cottonseed hull culture medium was used to the mushroom production experiment according to our previous study (Qiu et al., 2018). The samples of mycelia, primordia, and caps of fruiting bodies, denoted as M, P, and C, were collected and stored at −80°C after freezing with liquid nitrogen.

Identification of RDRP, Dicer, and AGO proteins
The amino acid sequences of the proteins which were functionally annotated as RDRP, Dicer, and AGO proteins were derived from the genome of P. cornucopiae. The sequences were subjected to domain analyses using the conserved domain database of NCBI for determining protein function. The sequences of the RDRP, Dicer, and AGO proteins of other fungi were retrieved from GenBank, and aligned to the corresponding proteins of P. cornucopiae with CLUSTALW (Thompson et al., 1994;Hu et al., 2013;Wang et al., 2021). A phylogenetic tree was constructed using the maximum likelihood method based on the Tamura-Nei model and 1,000 bootstrap replicates with the MEGA software, version 5.0 (Kumar et al., 2016) for analyzing the relationships between the RDRP, Dicer, and AGO proteins of P. cornucopiae and those of other fungi in literature.

Deep sequencing of mRNAs and small RNAs
Comparative mRNA and milRNA analyses were performed using the CCMSSC 00406 strain of P. cornucopiae. Samples of the different developmental stages, including M, P, and C, of CCMSSC 00406 were collected in three biological replicates and subjected to mRNA and small RNA sequencing. The total RNA was extracted from all the samples using TRIzol reagent (Invitrogen, Carlsbad, United States), according to the manufacturer's instructions, and subsequently treated with RNase-free DNase I (TaKaRa, Shiga, Japan) for removing the genomic DNA. The concentration of the RNA was evaluated using a NanoDrop 2000 spectrophotometer (ThermoFisher, Waltham, United States), and the integrity of the RNA was detected using an Agilent 2100 Bioanalyzer (Agilent, Palo Alto, United States). The cDNA libraries were constructed according to the protocol for library construction and sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego, United States). The small RNAs were isolated from the total RNA by polyacrylamide gel electrophoresis (PAGE) with a 6% Tris, boric acid, and EDTA (TBE)-urea denaturing gel, and ligated to specific 5′ and 3′ adaptors. The cDNA libraries were sequenced on an Illumina HiSeq2500 platform following reverse transcription and appropriate amplification and purification.

Bioinformatics analyses of mRNAs and small RNAs
The clean data (clean reads) were obtained from the raw RNA-seq data by removing the reads containing adapters, poly-N, and low quality reads. The HISAT2 software was used for mapping the reads to the reference genome. 1 The genes were annotated by BLAST search against using the NCBI non-redundant protein sequence (NR), Gene Ontology (GO; Tatusov et al., 2000), Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa et al., 2004), KOG (Koonin et al., 2004), and protein family (Pfam) databases. The expression levels of the genes were estimated by fragments per kilobase of transcript per million fragments mapped (FPKM). The genes that were differentially expressed among the different developmental stages were determined using the DESeq2 tool, with a false discovery rate (FDR) of <0.05. The differentially expressed genes (DEGs) were subjected to GO and KEGG pathway enrichment analyses, and the 20 most enriched pathways with the lowest Q values were selected (Mao et al., 2005).
1 https://www.ncbi.nlm.nih.gov/nuccore/WQMT00000000.2/ The high-quality small RNA sequence reads (clean reads) were filtered from the total reads by removing the low-quality reads, reads containing adaptor sequences, and sequences smaller than 18 nt or longer than 30 nt. The clean reads were subsequently aligned to the genome of P. cornucopiae using the Bowtie software (Langmead et al., 2009). The different non-coding RNAs, including rRNAs, snRNAs, tRNAs, and snoRNAs, were identified using the Bowtie software, and subsequently removed. The remaining unannotated small RNAs were analyzed for detecting the known miRNAs from miRBase, and the miRDeep2 tool was used for predicting the novel milRNAs (Friedlander et al., 2012). The Randfold software was used for predicting the secondary structures of the novel milRNAs. The expression of the milRNAs in different developmental stages was calculated using the transcripts per million (TPM) normalization method (Fahlgren et al., 2007). The milRNAs that were differentially expressed in the M, P, and C at the different developmental stages were identified using the DESeq tool (Zeng et al., 2018). The target genes of the milRNAs were predicted based on the milRNA sequence information using the TargetFinder, 2 miRanda (Enright et al., 2003), and RNAhybrid (Rehmsmeier et al., 2004) webtools, as previously described.

Analysis of the expression of milRNAs and mRNAs related to the development of fruiting bodies
The expression levels of the milRNAs were quantified by stemloop real-time PCR, using 5S rRNA as the internal control for each sample . The stem-loop primers in the reverse transcription kit and the upstream primers used for qRT-PCR were designed using miRNA design software. 3 The first-strand cDNA was synthesized using a miRNA first Strand cDNA Synthesis Kit (Vazyme, Nanjing, China), according to the manufacturer's instructions. The miRNA Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) and an ABI 7500 real-time PCR amplifier (Applied Biosystems, Foster City, CA, United States) were used for qRT-PCR, as described in our previous study. The expression of the target genes of the milRNAs was detected using glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as the internal control for each sample, as previously described (Hou et al., 2021). The relative expression levels of the milRNAs and their target genes in the different stages were quantified using the comparative threshold cycle (CT) 2 −△△CT method. The primers used for qRT-PCR amplification of the milRNAs are enlisted in Supplementary Table S1.

Dual-luciferase activity assay
A ~ 200 bp sequence near the binding site of the wild-type (WT) g8971 and mutant g8971 genes were synthesized and inserted into the pmirGLO vector. Briefly, HEK293T cells were co-transfected with 0.2 μg of the luciferase reporter vector (g8971-WT or g8971-MT) and Frontiers in Microbiology 04 frontiersin.org 10 ng of milR-20-mimic or mimic NC together with the renilla luciferase construct using lipofectamine TM 2000 (Invitrogen), according to the manufacturer's instructions (Grentzmann et al., 1998). The HEK293T cells were collected 48 h post-transfection, and the activities of luciferase and renilla luciferase were measured using a Dual-Luciferase ® Reporter Assay System (Promega, Wisconsin, United States) according to the instructions provided (Cai et al., 2017). Five biological repetitions of the experiment were averaged and analyzed using Student's t test.

Overexpression and silencing of milR20
For the overexpression of milR20, a precursor of milR20 was amplified using the WT genomic DNA as the template, and inserted by homologous recombination into a pCAMBIA1300 vector containing the gpd promoter of P. ostreatus. MilR20 was silenced using the short tandem target mimic (STTM) technology; STTM contains two target mimic (TM) sequences and a 48 nt-long specific linker sequence. The TM sequence corresponds to the sequence that is complementary to milR20 and possesses a tri-nucleotide that is inserted between the 10 and 11th bases of milR20. The STTM sequence was subsequently ligated to the pCAMBIA1300 vector.
All the recombinant plasmids were verified by sequencing and transfected into WT cells by Agrobacterium tumefaciens-mediated transformation (ATMT), as previously described. The strains in which milR20 was overexpressed and silenced were detected by PCR for cloning the hpt gene. The expression levels of milR20 following overexpression and silencing were quantified by qRT-PCR. The diameters of the colonies of the WT strain and the strains in which milR20 was overexpressed (OE-milR20) and silenced (STTM-milR20) were measured using the cross method for determining the mycelial growth rate. The WT, OE-milR20, and STTM-milR20 strains were separately inoculated on a culture medium for analyzing the primordial formation time, the developmental cycle of the fruiting bodies, and the spore print.

Statistical analysis
All statistical analyses were performed using the SPSS 26.0 software (SPSS Inc., Chicago, United States). The data are presented as mean ± SEM values. Statistical significance was defined as *(p < 0.05), **(p < 0.01), and ***(p < 0.001). The GraphPad Prism 8.0.1 software (GraphPad Software Inc., San Diego, United States) and Excel 2010 software (Microsoft, Redmond, WA, SA, United States) were used for drawing figure.

mRNA sequencing and analyses
The expression profiles of the genes expressed in the M, P, and C stages across the three different developmental stages were determined using mRNA-seq. Three biological replicates were sequenced for each tissue type and a total of approximately 428 million clean reads were obtained from all the samples after filtering the low-quality reads. The number of reads in the samples ranged from 39 to 65 million. The reads were mapped to the genome of P. cornucopiae; approximately 32-53 million reads (80-82% of the total reads) mapped to the genome of P. cornucopiae (Supplementary Table S2). The Pearson correlation coefficients results indicated that the reproducibility between biological replicates was high enough for subsequent studies (Supplementary Figure S1).
In order to identify the genes that are involved in the development of P. cornucopiae, the DEGs among the different developmental stages were identified using the following criteria: FDR ≤ 0.05 and fold change (FC) ≥ 1.5. The number of DEGs in the M vs. P, M vs. C, and P vs. C comparison groups was determined to be 4,819, 7,017, and 2,917, respectively ( Figure 1A). A total of 7,934 DEGs were identified from all the comparison groups, of which the number of DEGs in the M vs. C comparison group was highest. This indicated that the number of genes differentially expressed between the different intervals was higher than that between adjacent stages, and this finding was consistent with the organizational difference.
The DEGs were subjected to GO and KEGG pathway enrichment analyses. The results of GO enrichment analysis revealed that the DEGs in the M vs. C, P vs. C, and M vs. P comparison were significantly enriched in 46, 54, and 47 GO terms (p < 0.05), respectively. The significantly GO terms in the overlapped groups and the top 10 most significantly enriched GO terms in a single group were present. In biological process, the significantly enriched GO terms shared in M vs. C and P vs. C were gluconeogenesis, glycolytic, nitrogen compound metabolic process, and carbohydrate transport. The significantly enriched GO terms shared in M vs. P and M vs. C were histidine biosynthetic process, translation, and reciprocal meiotic recombination. The significantly enriched GO terms shared in M vs. C, P vs. C and M vs. P were oxidation-reduction process. The significantly enriched GO terms in single group were translational elongation, response to stress, transmembrane transport etc. In molecular function, the significantly enriched GO terms in different group were ATPase activity, oxidoreductase activity, transporter activity, and signal transducer activity etc. In cellular component, the significantly enriched GO terms in different stage were mitochondrion, ribosome, intergral component of membrane etc. These findings indicated that the DEGs involved in energy metabolic process, signal transduce process, and DEGs located to mitochondrion and membrane could play a role in the development of fruiting bodies in P. cornucopiae ( Figure 1B). The total 7,934 DEGs from all three comparison group were used for KEGG analysis and 20 most enriched KEGG pathways were present. The results revealed that the DEGs were enriched in the MAPK signaling pathway, metabolism, cell cycle, and protein processing endoplasmic reticulum terms ( Figure 1C). These findings indicated that the DEGs that were involved in these pathways could play a key role in the development of P. cornucopiae.

Identification and analysis of the genes involved in milRNA biogenesis and function
RNA-dependent RNA polymerase, AGO, and Dicer proteins play a key role in milRNA biogenesis and function in eukaryotes. Therefore, the presence of these proteins could indicate that P. cornucopiae contains miRNAs. Therefore, the genes encoding RDRP, AGO, and Dicer proteins were identified in the genome of P. cornucopiae by  (Figure 2A). The results of phylogenetic analysis demonstrated that the sequences of the RDRP, AGO, and Dicer proteins of P. cornucopiae were highly homologous to those of C. sinensis and C. militaris ( Figure 2B). These findings suggested that P. cornucopiae possesses functional milRNA machinery. In order to further explore the possible role of milRNAs in the developmental process of P. cornucopiae, the expression levels of RDRP, AGO, and Dicer at the different developmental stages were analyzed from the transcriptome data. The results demonstrated that the expression levels of AGO were higher than that of the other genes at each of the developmental stages. The expression levels of the genes encoding RDRP, AGO, and Dicer proteins varied across the different developmental stages of P. cornucopiae ( Figure 2C), which suggested that the expression and functions of milRNAs were various during development.

Sequencing and analyses of the milRNAs in the different development stages
In order to identify the milRNAs that are related to the development of fruiting bodies in P. cornucopiae, the small RNAs in the M, P, and C stages were subjected to sequencing, which was performed in triplicate. The samples of M tissues were denoted as M-1, M-2, and M-3, while the samples of P tissues were denoted as P-1, P-2, and P-3, and the samples of C tissues were denoted as C-1, C-2, and C-3. Approximately 20 million raw reads were obtained from each sample. Approximately 10 million clean reads with lengths varying between 18 and 30 nt were obtained after filtering the low-quality reads and trimming the 3′-specific adaptors, and the remaining small RNAs were annotated. The clean reads were aligned to the genome of P. cornucopiae, and the results demonstrated that the 2-6 million reads included various small ncRNAs, including rRNAs, tRNAs, and snoRNAs, which accounted for 21.35-60.19% of the total clean reads obtained from the different developmental stages. The 4-8 million unannotated clean reads were subsequently analyzed for further prediction of milRNAs, which accounted for 39.64-77.94% of the total clean reads obtained from the different developmental stages (Supplementary Table S3).
A total of 32 milRNAs were finally identified from the different developmental stages of P. cornucopiae, and 31, 26, and 30 milRNAs were identified from the M, P, and C stages, respectively (Supplementary Figure S2). The milRNAs were denoted as milR1-milR32 (Supplementary Table S4). Analysis of the length distribution of the milRNAs revealed that the lengths of the mature milRNAs ranged from 18 to 25 nt. The majority of these milRNAs were 18-22 nt long, and accounted for 93.75% of the total mature milRNAs, which was higher than the number of milRNAs with other lengths ( Figure 3A). The results of nucleotide bias analysis revealed that the nucleotides at the 5′-terminus had a strong preference for uracil (U) in the milRNAs that were 18-22 nt long, which was similar to that observed in animals and plants. However, the milRNAs that were 24-25 nt long were enriched in adenine (A) at the 5′-end ( Figure 3B). The abundance of milRNAs was normalized according to the TPM normalization method. Of the 32 milRNAs, 25 were expressed in all the three developmental stages, while the other seven milRNAs were expressed in one or two of the developmental stages (Supplementary Figure S2). The heatmap in Figure 3C demonstrates that the expression patterns of the milRNAs varied across the different development stages and nine of these milRNAs were highly expressed during the entire development of P. cornucopiae. The findings suggested that the milRNAs that were expressed at high levels in all the three developmental stages could play a crucial role in the development of P. cornucopiae.
To explore the milRNAs that are related to the development of fruiting bodies in P. cornucopiae, we analyzed the differentially expressed milRNAs (DEMs) across the three developmental stages. p < 0.05 was regarded as the threshold for determining the significant differences in milRNA expression. The results demonstrated that 13, 13, and 6 milRNAs were significantly different expressed in the M vs. P, M vs. C, and P vs. C comparison groups, respectively. A total 20 DEMs were identified in the three comparison groups, of which three DEMs were common between the M vs. C and P vs. C comparison groups, and could play an important role in the development of fruiting bodies ( Figure 4A). Analysis of the expression levels of the DEMs in the different comparison groups revealed that the number of DEMs downregulated was higher than that of upregulated in the C vs. P comparison group, while the number of upregulated DEMs was approximately equal to that of the downregulated DEMs in the M vs. P and P vs. C comparison group (Supplementary Table S5; Figure 4B). These results indicated that the DEMs that were downregulated along the development of the fruiting body could play a more crucial role in the development of P. cornucopiae.
In order to elucidate the potential functions of the DEMs in the development of P. cornucopiae, the target genes of milRNAs were predicted using the TargetFinder, miRanda, and RNAhybrid software, as previously described. A total 17 milRNAs in the 20 DEMs were predicted to target 44 genes. The findings revealed that some of the DEMs could regulate several target genes and more than one milRNA could regulate the same target gene (Supplementary Table S5), which was consistent with the reports of previous studies on plant and animal miRNAs. These target genes were subjected to functional enrichment analyses using the GO and KEGG databases. The findings demonstrated that 41 of the target genes of DEMs were functionally annotated. The results of GO enrichment analysis revealed that the target genes were primarily enriched in the transporter activity, Frontiers in Microbiology 08 frontiersin.org receptor activity, and catalytic activity terms in the molecular function category; the signaling pathways, cellular processes, and response to stimulus terms in the biological process category; and the organelles and membranes terms in the cellular component category (Supplementary Figure S3). The results of KEGG enrichment analysis demonstrated that the target genes were enriched in different pathways, including the phosphatidylinositol signaling system, endocytosis, MAPK signaling pathway, protein processing in endoplasmic reticulum, and other metabolic pathways ( Figure 4C).

Correlation analysis of milRNAs and mRNAs
The Venn diagram depicting the DEGs and the target genes of DEMs indicated that 29 of the 44 target genes of the DEMs were differentially expressed across the different developmental stages ( Figure 5A). Of these, only seven DEGs could be mapped by KEGG pathway analyses, and were found to be enriched in the MAPK signaling pathway, phosphatidylinositol signaling system, protein processing in endoplasmic reticulum, and other metabolic pathways that may play crucial roles in the growth and development of P. cornucopiae (Table 1). For instance, the findings revealed that the g4622 gene was involved in the phosphatidylinositol signaling system, g9630 was enriched in endocytosis, g8971 was enriched in the MAPK signaling pathway, and the remaining genes were involved in other pathways.
MiRNAs are important regulators of gene expression and act via degradation or translational repression of target mRNAs (Zdanowicz et al., 2009). Correlation analysis of expression profiles of milRNAs and their targets showed that only milR20 and milR14 had the relative opposite expression trend to their targets in partial development stage (Figures 5B,C), while other miRNA-target did not ( Figures 5D-F). The results demonstrated that milR14 could play a minor role in the development of fruiting bodies because it was not expressed in the P and C stages. However, the stage-specific expression pattern of milR14 suggest it may have important function during stage M development. However, the expression of milR20 was downregulated from the P to the C stage, while the expression of its target gene, g8971, was upregulated from the P to the C stage. These result indicated that g8971 may be the target of milR20. milR20 was highly expressed in the three development stage with the expression level higher than 40,000, and it was a DEM from the M vs. C and P vs. C comparison group with the expression level decreased from M to C stage. The predict target of milR20, g8971, was DEGs at different developmental stages, and encodes a pheromone receptor that involved in the MAPK signal pathway which plays an important role in the development. These results therefore indicated that milR20 could negatively regulate the development of fruit bodies in P. cornucopiae.
In this study, qRT-PCR analysis was also performed for validating the expression profiles of the miRNA-target modules of interest, namely, milR20 and g8971. The expression pattern of milR20 and its target gene, g8971, obtained by RT-qPCR analysis was similar to that determined by high-throughput sequencing, and the findings revealed that the expression of milR20 tended to decrease from stages P to C (Supplementary Figure S4), while the expression level of g8971 increased from stage P to C. The findings revealed that the milR20 had opposite expression trend to g8971 from stage P to C. So milR20 was selected for further analyses.

milR20 targets g8971 and inhibits its expression
The results of prediction using the RNAhybrid software revealed that milR20 targeted the 848-868 nt region of the g8971 gene ( Figure 6A). Dual-luciferase reporter assays were performed for elucidating the targeting relationship between milR20 and its target gene, g8971. The WT g8971 reporter vector (g8971-WT) and the mutant plasmids (g8971-MT) were constructed, and co-transfection experiments were performed. The results of the dual-luciferase assay demonstrated that milR20 significantly reduced the relative luciferase activity of g8971-WT. However, there was no significant effect on the relative luciferase activity of g8971-MT. These findings therefore indicated that milR20 could negatively regulate the expression of g8791 by directly binding to and targeting g8971 ( Figure 6B).

Functional analysis of milR20 by STTM-mediated silencing and overexpression
In this study, the copy number of milR20 in the genome of P. cornucopiae was first determined by comparing the milRNA sequence with the genome using the BLASTn program. The results demonstrated that there was only one copy of milR20 in the genome. The precursor sequence of milR20 (pre-milR20) was identified and analyzed using the miRDeep2 software. A 250 bp-long sequence of  pre-milR20 was obtained and represented as hairpin structures, which verified that milR20 was a real milRNA (Supplementary Figure S5). In order to explore whether milR20 has any role in the development of P. cornucopiae, milR20 was separately silenced and overexpressed and the phenotypic effects were analyzed. STTMmediated silencing has been shown to be an effective tool for inhibiting the activity of endogenous mature miRNAs in plants. In this study, we designed milR20 STTM constructs containing two same non-cleavable milRNA binding sites (highlighted in blue in Figure 7A), and linked by a 48-88 nt spacer (colored in yellow). We generated transgenic strains in which milR20 was overexpressed or silenced with STTM. Analysis of the expression levels of milR20 in the M stage of the transgenic strains revealed that the expression of milR20 increased significantly in the OE-milR20 strain compared to that in the WT, while its expression in the STTM-milR20 strain (32-75%) decreased significantly compared to that of the WT ( Figure 7B).
The effect of milR20 on the development of P. cornucopiae was detected by analyzing the mycelial growth of the WT and mutant strains that had been incubated on PDA plates. The rate of mycelial growth in the OE-milR20 strain was significantly lower than that of the WT strain; however, the rate of mycelial growth in the STTM-milR20 strain was significantly higher than that of the WT strain ( Figures 7C,D). These findings indicated that milR20 may play a negative role in mycelial growth.
The effect of milR20 on the development of the fruiting bodies of P. cornucopiae was subsequently detected by cultivating the mutant and WT strains on cultivation substrates. The time required for the formation of primordia was initially analyzed statistically. The primordia appeared on the 28th to the 29th day in all the strains, and there were no significant differences between the mutants and WT with respect to the time required for the formation of primordia. This finding indicated that milR20 had no influence on the formation of primordia (Supplementary Figure S6). Analysis of the morphology of the fruiting bodies revealed that the developmental cycle of the fruiting bodies was prolonged in the OE-milR20 strain that overexpressed milR20, compared to that of the WT. However, the developmental cycle of the fruiting bodies was slightly shortened in the STTM-milR20 strain compared to that of the WT. Analysis of the spore morphology revealed that the spore print density of the OE-milR20 strain did not exhibit significant alterations compared to that of the WT, while the quantity of spores produced by the STTM-milR20 strain increased slightly ( Figure 7C). Detection of the expression levels of the target g8971 gene of milR20 revealed that g8971 was downregulated in the OE-milR20 strain, while there was no significant difference in the expression of g8971 in the STTM-milR20 strain ( Figure 7E). These results indicated that milR20 may negatively regulate the development of the fruiting bodies of P. cornucopiae by inhibiting the expression of g8971.

Discussion
Previous studies have demonstrated that miRNAs have an important role in the development of plant and animals. The development of high-throughput sequencing technologies has enabled the identification of milRNAs from various species of fungi in recent years. However, there is a scarcity of information regarding the functions of milRNAs in fungi. Pleurotus cornucopiae is an important mushroom that has been used for studying the functions of milRNAs in fungal development. AGO, Dicer, and RDRP proteins are key components of miRNA maturation and function in N. crassa and are conserved in C. militaris and other fungal species that have been reported to possess milRNAs (Lee et al., 2010;Yang et al., 2013;Shao et al., 2019;Wang et al., 2021). The present study demonstrated that the AGO, Dicer, and RDRP proteins of these fungal species were closely related to those of P. cornucopiae, which indicated that mechanisms of milRNA biogenesis also exist in P. cornucopiae. The number of genes encoding Dicer, AGO, and RDRP proteins vary across different fungal species, and the expression patterns of these homologs vary across the different developmental stages (Lee et al., 2010;Shao et al., 2019Shao et al., , 2020Wang et al., 2021). This suggests that the homologs of genes encoding Dicers, AGOs, and RDRPs may play different roles during the formation of mature milRNAs from dsRNAs, and the expression and functions of milRNAs also vary during development. In this study, the expression level of genes encoding Dicers, AGOs, and RDRPs were analyzed by transcriptome analysis across the different developmental stages of P. cornucopiae. All the genes were expressed at different developmental stages, and there were Verification of the targeting relationship between milR20 and g8971. (A) The specific binding sites of milR20 in the sequence of g8971 are depicted, and the mutation sites were designed based on this sequence. The sequences in red represent the mutation sites. (B) The activity of luciferase was measured by dual-luciferase reporter assays. The activities of luciferase of g8971-WT decreased markedly in cells transfected with miR-milR20 compared to that of the control. The values are depicted as the mean ± SE; ***p < 0.001. Microbiology  11 frontiersin.org variations in the expression patterns, which suggested that the homologs of genes encoding Dicers, AGOs, and RDRPs may function in a coordinated manner to regulate the expression and function of milRNAs during the development of P. cornucopiae.

Frontiers in
The results of mRNA analysis demonstrated that the DEGs among the different development stages were mainly enriched in the signaling and growth terms in the biological process category of GO, and in the MAPK signaling pathway of KEGG. These findings indicated that the genes that were involved in the MAPK signaling pathway could be involved in the development of P. cornucopiae, and this finding was consistent with the results of previous studies on plants (Yi et al., 2016;Xiao et al., 2017;Chen et al., 2021). MAPK cascades are known to transmit extracellular signals to intracellular targets and play a crucial role in regulating several fundamental processes, including proliferation, differentiation, and cellular response to diverse extrinsic stresses (Guo et al., 2020).
In this study, the results of milRNA analysis demonstrated that a total of 32 milRNAs were identified in the M, P, and C stages, and the majority of these milRNAs were expressed at all the stages. However, the expression levels of the milRNAs varied across the three developmental stages. A total of 20 milRNAs were differently expressed in the three M vs. P, M vs. C, and P vs. C comparison groups, of which the number of downregulated DEMs was higher than the number of upregulated DEMs. These results indicated that the downregulated milRNAs could play a vital role in the development of P. cornucopiae. MiRNAs are important regulators of gene expression and function via degradation or translational repression of the target mRNAs (Moran et al., 2017). Integrated analysis of the milRNAs and mRNAs revealed that the target DEGs of the DEMs mapped to the MAPK signaling pathway, and this finding was consistent with the results of mRNA analysis. These results indicated that the milRNAs which regulated the MAPK signaling pathway could play a significant Effect of milR20 on the development of Pleurotus cornucopiae. (A) Schematic representation of the structure of the STTM plasmids that were used for silencing milR20. A target mimic with an unmodified central sequence (highlighted in blue) that was complementary to the central portion of milR20 and had a trinucleotide bulge was inserted at the cleavage site located at 10-11 nt of the miRNA. (B) Analysis of the relative expression levels of the WT and mutants by qRT-PCR. (C) Phenotype of the WT and mutant strains at different developmental stages. (D) Statistics of the colony diameter of the WT and milR20-recombinant strains under normal temperature, determined using SPSS. (E) Detection of the expression levels of the target g8971 gene of milR20 in the WT and mutant strains by qRT-PCR; * and ** indicate significant difference at p < 0.05 and p < 0.01, respectively.
Frontiers in Microbiology 12 frontiersin.org role in the development of P. cornucopiae, and was consistent with the findings of previous studies which reported that miRNAs regulate the activity of the MAPK cascade and influence cellular proliferation in animals (Chen et al., 2017;Xiao et al., 2018;Xu et al., 2018;Safa et al., 2020;. Correlation analysis of the milRNA expression profiles and their target genes revealed that only a small number of the milRNA-mRNA pairs exhibited an opposite expression pattern in the different developmental stages. These results indicated that the expression of the majority of target genes was possibly not negatively regulated by the milRNAs. Previous studies have also demonstrated an incoherent regulation between miRNAs/milRNAs and their target genes (Shao et al., 2019). This could be attributed to the fact that milRNAs primarily mediate gene regulation by repressing mRNA translation in fungi and not via mRNA degradation. Considering the complex regulatory network of gene expression, this could be alternatively explained by the fact that the expression of target genes can also be regulated by other genes, including genes encoding transcription factors, and competing endogenous RNAs that competitively bind to miRNAs (Kartha and Subramanian, 2014).
In this study, the results of the dual-luciferase activity assay and qRT-PCR results revealed that milR20 targeted the g8971 and inhibited the expression of g8971. The expression of milR20 was downregulated in both the M vs. C and P vs. C comparison groups, and the target gene of milR20, g8971, was invovled in the MAPK signaling pathway and could play a vital role in the development of P. cornucopiae. Therefore, the functions of milR20 in the development of P. cornucopiae were subsequently analyzed using overexpressing and silencing by short tandem target mimic (STTM) technology in this study.
The STTM technology mimics the binding of target miRNAs to RNA-induced silencing complex (RISC) to inhibit the functions of target miRNAs (Teotia et al., 2016). This method has been proven to be an effective and stable tool for blocking the activity of endogenous mature miRNAs in plants (Jia et al., 2015;Zhang et al., 2017). In this study, the expression levels of milR20 decreased significantly in the STTM-milR20 strain, indicating that the STTM technology can also be used to effectively silence fungal milRNAs. The mycelial growth rates and the development of fruit body in the OE-milR20 strain were reduced and prolonged, while those of the STTM-milR20 strain were the opposite in comparison to that of the WT strain. The expression level of g8971 in the OE-milR20 were significantly decreased, while no significant difference in STTM-miR20 strains. These could be explained by the expression fold change difference of milR20 in different strains or the regulation of g8971 by other genes. These findings indicated that milR20 could negatively regulate the growth of P. cornucopiae by repressing the expression of g8971 regulating the MAPK signaling pathway.

Data availability statement
The original sequence data of transcriptome can be found at the following link: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA943625 and the milRNA sequence data can be found at the following link: https:// www.ncbi.nlm.nih.gov/bioproject/PRJNA944818. Other data presented in this study are available in Supplementary material.