Whole-Transcriptome Sequencing Reveals That mRNA and ncRNA Levels Correlate with Pleurotus cornucopiae Color Formation

: Pleurotus cornucopiae cap color is an important commercial trait. The roles of non-coding RNA molecules (ncRNAs) in fungal fruiting body color formation are unclear. Whole-transcriptome analyses were performed, identifying messenger RNA (mRNA) and ncRNA (including long stranded non-coding RNA (lncRNA), micro RNA-like (milRNA)


Introduction
Pleurotus cornucopiae has significant nutritional and medicinal value [1].Cap color is an important commercial trait of P. cornucopiae and exhibits substantial differences among varieties of this fungus.Consumers typically prefer mushrooms with dark caps; therefore, it is of great significance to explore the genetic regulation mechanisms underlying P. cornucopiae color formation.
In recent years, there has been significant interest in studying the color of edible mushroom caps, diverse color is determined by various pigment costituents.Cordyceps xanthophylls in Cordyceps militaris were identified as potentially useful pigments and carotenoid supplements for the food industry [2].Moreover, understanding of rare virulent strains has been deepened through the identification of the molecular structure of the red pigment, syringone, in Lactarius lilacinus [3].While lack of tyrosine in the culture medium of Auricularia auricula strains can lead to diminished melanin secretion, among other characteristics [4].Further, genes related to cap color have been discovered, cloned, and functionally identified.In the Agaricus bisporus genome, 42 genes were screened for their involvement in melanin synthesis [5], while the investigation of genetic regulation of P. cornucopiae color revealed that cap color is a quantitative trait, controlled by multiple genes, where dark gray is incompletely dominant over white [6].In our laboratory, we localized major quantitative trait loci controlling P. cornucopiae cap color by BSA-seq and XP-GWAS, and found that the tyrosinase gene (PcTYR) can regulate cap color; however, the regulatory networks regulating cap color formation were unclear [7].
Non-coding RNA (ncRNA) molecules do not encode proteins and can perform biological functions by regulating gene expression at transcriptional or post-transcriptional levels [8].ncRNAs include microRNAs (miRNAs), long-stranded non-coding RNAs (lncRNAs), and circular RNAs (circRNAs) [9,10], which have important roles in color formation in plants and animals.By constructing competitive endogenous RNA (ceRNA) networks of lncRNAs and circRNAs in mice with different coat colors, complex interactions between ncRNAs and messenger RNAs (mRNAs) related to skin and melanocyte development were discovered, which co-regulate pigmentation development [11].Wholetranscriptome sequencing and broadly targeted metabolome analysis of circRNAs, mRNAs, and small RNAs (sRNAs) of Melastoma candidum at three different developmental stages led to the construction of nine circRNA-mediated ceRNA regulatory networks that demonstrated important roles in M. candidum flower development and color formation [12].In fungi, ncRNAs, such as lncRNA and micro-RNA-like (milRNA) molecules are involved in cell division, nutrient perception, sexual reproduction, and other processes [13][14][15].CircRNAs regulate gene expression in diverse biological processes in Ganoderma lucidum [16], while milRNAs can regulate P. cornucopiae fruiting body development [17]; however, there are no reports of investigation of the regulation of fruiting body color by ncRNAs.
Whole-transcriptome sequencing can be used to characterize associations among mRNAs, miRNAs, circRNAs, lncRNAs, and other ncRNAs.Unlike transcriptome analysis, which solely considers overall gene expression and structure, this approach can also be used to analyze the expression of lncRNAs, milRNAs, and circRNAs, which can facilitate a deeper understanding of the roles of ncRNAs in biological processes, particularly at the level of post-transcriptional regulation [18][19][20].
The ceRNA hypothesis proposes that mRNAs and ncRNAs can interact by binding to shared milRNA response elements, allowing lncRNAs/circRNAs to act as sponges, absorbing milRNAs, and resulting in alterations in target genes regulated by milRNAs [21,22].Whole-transcriptome sequencing is widely applied in studies of both animals and plants.Complete analysis of the longest dorsal muscles from two distinct pig breeds using whole-transcriptome sequencing revealed that ncRNAs, functioning as ceRNAs of miR-NAs can cooperatively regulate skeletal muscle proliferation and differentiation, significantly impacting its formation [23].Whole-transcriptome analysis of tomato leaf senescence revealed that ncRNAs and mRNAs interact to jointly affect the senescence process [24].Such coordinated regulation of mRNA and ncRNA has also been discovered in studies on the growth and development of other plants, such as cucumber and soybean [25,26].Hence, whole-transcriptome analysis is a feasible approach to examine the regulatory mechanisms underlying ncRNA and mRNA regulation of P. cornucopiae cap color formation.
In this research, we performed whole-transcriptome sequencing of three P. cornucopiae varieties with increasing color depth (white, grayish-white, and grayish-black).Bioinformatics methods were used to investigate variation in mRNA and ncRNA expression levels in strains of different colors.Important pathways for color formation were identified, and mRNAs and ncRNAs associated with these pathways in P. cornucopiae systematically screened.Further, mRNA and ncRNA co-expression networks with the potential to regulate cap color were identified, and the molecular mechanism underlying coordinated regulation of cap color by mRNAs and ncRNAs was elucidated.This research will enable new directions in the study of P. cornucopiae cap color.

Experimental Materials and Growing Conditions
The materials used in this study were P. cornucopiae CCMSSC00358 (P358), CCMSSC00406 (P406), and CCMSSC04656 (PHY), which were cultured in PDA medium for 7 days and subsequently inoculated in cottonseed husk cultivation bottles.Inoculated cultivation bottles were placed in a germination room for dark culture.Once mycelium covered the bottles, mushroom growth experiments were conducted.Young mushrooms were collected as materials, three biological repeats were made for each strain, and they were placed into liquid nitrogen immediately after collection and stored in a laboratory refrigerator at −80 °C.

Strand-Specific Library Construction and Sequencing
For sequencing mRNAs, lncRNAs, and circRNAs, ribosomal RNAs (rRNAs) were removed using an Epicentre Ribo-Zero TM kit.Next, enriched mRNA and ncRNA samples were fragmented into shorter pieces and reverse-transcribed into cDNA using random primers.Second-strand cDNA was synthesized using buffer, dATP, dUTP, dCTP, dGTP, RNase H, and DNA polymerase I.The purified double-stranded cDNA was end-repaired; poly(A) tails were then added and ligated to sequencing adapters.Finally amplified by PCR to produce a cDNA library, the quality of library was checked by Qsep-400 method and then sequenced (NovaSeq 6000 system; paired-end 150 bp cycles) by Biomarker Technologies Corporation (Beijing, China).

Small RNA Library Construction and Sequencing
An sRNA library was constructed using the small RNA Sample Pre Kit (Vazyme.NR801-01).Briefly, the sequencing adapters were added to the 3′ end and 5′ end of small RNAs using the T4 RNA Ligase 1 and T4 RNA Ligase 2 , respectively, the cDNA was obtained and amplified by PCR, and the products were purified by polyacrylamide gel electrophoresis.All libraries generated were validated using the Qsep-400 method.Then, sequenced via the HiSeq 2500 (IIIumina) platform which sequencing read length was single-end (SE) 50 nt.Library construct and sequencing programs were executed by Biomarker Technologies Corporation (Beijing, China).

Data Filtering, mRNA and ncRNA Identification, and Differential Expression Analysis
Raw reads containing adapters, unknown nucleotides, and low-quality bases were removed, and clean reads were obtained.HISAT2 (v2.0.4) [27] was utilized to align clean data with the specified reference genome, P. cornucopiae (NCBI accession no.WQMT00000000), and the mapped data was obtained.Then the mapped reads were used to reconstructed the transcripts using StringTie software (v1.3.1), and novel transcripts were identified using HISAT2.Expression levels of all mRNAs were measured according to fragments per kilobase of transcript per million mapped reads (FPKM) [28].CPC (CPC2-beta) [29], CNCI (v2) [30], CPAT [31], and Pfam (v1.3) [32], with default parameters, were used to identify lncRNAs.Shared lncRNAs identified by these four methods were used for subsequent analysis.
Reference genome read sequences were compared with mature milRNA sequences using miRBase (v21), a database of known milRNAs, to identify established milRNAs.miRDeep2 software [33] was used to predict novel milRNAs.
CircRNAs were detected using find_circ [34].The software extracted 20 base anchor sequences from both ends of reads that did not align with the genome, and compared them with the genome to search for distinct matching sites.If the positions of the two anchors were aligned in a reverse orientation, anchor reads were extended to detect circRNA splice sites; circRNAs were recognized if the sequences on both sides of this point were GT/AG splicing signals.
mRNAs and ncRNAs expression levels were assessed by FPKM.DESeq2 software was used to identify significantly differentially expressed mRNAs and ncRNAs in the sample groups [35], with threshold values of fold-change (FC) ≥ 1.5 and false discovery rate (FDR) ≤ 0.05.

Gene Ontology (GO) and Kyoto Enrichment of Genes and Genomes (KEGG) Enrichment Analyses
GO and KEGG enrichment analyses were performed on differentially expressed mRNAs (DEmRNAs), target genes of differentially expressed lncRNAs (DElncRNAs), target genes of differentially expressed miRNAs (DEmiRNAs), and differentially expressed circRNA (DEcircRNA) source genes, to assess differentially expressed RNA (DERNA) biological functions and enrichment pathways.GO database (https://www.geneontology.org/(accessed on 5 May 2023)) and the database of KEGG (http://www.genome.jp/kegg/(accessed on 5 May 2023)) were used to perform GO annotations and KEGG pathway analysis.

Network Visualization
The ceRNA hypothesis describes a novel mechanism of interaction among RNAs, the TargetFinder software was used for predicting the binding sites of the target RNAs and milRNAs, according to the sequences of milRNA response elements [36,37].Based on the results, the regulation networks of circRNA-milRNA-mRNA was constructed to illuminate the functions and relationships among DERNAs, and the circRNA-miRNA-mRNA networks was depicted using the Cytoscape software.
Genes with similar expression patterns often share similar functions.To identify correlations between different DERNAs, we utilized the pearsonr function of the scipy module in Python.DERNAs with absolute correlation value >0.9 and FDR < 0.005 were filtered out.Then, mRNA-lncRNA-circRNA-milRNA co-expression networks were constructed using Cytoscape, to elucidate the interaction relationships among them.

qRT-PCR
The expression levels of selected DEmRNAs, DElncRNAs, DEmilRNAs, and DE-circRNAs were verified using qRT-PCR, based on identified DERNA co-expression networks.Total RNA samples were used as a template for qRT-PCR assays.Reverse transcription was performed using the Vazyme Reverse Transcription Kit.A mixture of oligo (dT) primers and random 6-mer primers were used to reverse transcribe mRNAs and lncRNAs, and primers for mRNAs and lncRNAs were designed using an online website (http://www.biorun.com/tools/(accessed on 20 July 2023)).CircRNAs were reverse transcribed using random 6-mer primers and validated by designing upstream and downstream primers for qRT-PCR using the back-to-back principle at both ends of the circularization site.For milRNAs, stem-loop reverse primers and qRT-PCR primers were designed using Vazyme miRNA primer design software, and reverse transcription and qRT-PCR performed using a Vazyme kit specific for milRNAs.β-actin was used as an internal reference gene for mRNAs, lncRNAs, and circRNAs, and 5S rRNA was used as the internal reference gene for milRNAs.Primer sequences are listed in Supplementary Table S1.qRT-PCR was conducted using the following conditions: 95 °C for 3 min; followed by 40 cycles of 95 °C for 3 s, 60 °C for 32 s, and 72 °C extension for 30 s. Relative expression levels were then analyzed using the 2 −∆∆CT method.

Statistical Analysis
Statistical analysis was performed using Microsoft Excel and IBM SPSS Statistics 24.One-way ANOVA was used to assess the significance of differences among groups and the threshold for statistical significance was p < 0.05.Graphs were constructed and statistically analyzed using GraphPad Prism 9.

Comparison of Cap Color in Different P. cornucopiae Strains
Significant cap colors differences were observed among the lines, P358, PHY, and P406, in both young and mature stage fruiting bodies; lines P358, PHY, and P406 had white, grayish-white, and grayish-black caps, respectively (Figure 1A).Measurement of cap color brightness in the three strains showed a significant gradual decrease from line P358 (Lightness value [L] = 89.85), to line PHY (L = 62.5) and P406 (L = 45.80)(Figure 1B), indicating a gradual deepening of cap color.Young-stage fruiting bodies, at the earliest time of color difference appearance, were selected for used in subsequent analysis.

Global Characteristics Detected by Whole-transcriptome Sequencing
Whole-transcriptome sequencing of young fruiting body cap samples from the three P. cornucopiae strains, P358, PHY, and P406, was performed, with three biological replicates per strain.A total of 167.68 Gb clean RNA sequencing (RNA-seq) data were generated, after removing adapters and low-quality reads.Clean reads that aligned with the ribosomal RNA (rRNA) and unmapped reads were removed, and the remaining reads aligned to the P. cornucopiae reference genome, of which 71.65% to 76.44% were mapped to the reference genome and used to identify mRNAs, lncRNAs, and circRNAs (Supplementary Table S2).For sRNA sequencing, a total of 116.93 million clean reads were obtained, and more than 10.34 million clean reads were achieved for each sample.Highquality clean reads were used to identify milRNAs (Supplementary Table S3).
A total of 13,704 mRNAs (including 2089 newly identified genes), 1745 lncRNAs, 1102 circRNAs, and 56 predicted milRNAs were detected in the three different strains (Figure 2A).Among them, 11,725 mRNAs, 1066 lncRNAs, 314 circRNAs, and 45 milRNAs were identified in strain P358; 12,472 mRNAs, 1419 lncRNAs, 53 milRNAs, and 401 circRNAs in strain PHY; and 12,530 mRNAs, 1448 lncRNAs, 52 milRNAs, and 541 circRNAs in strain P406 (Figure 2B).The numbers of mRNAs and ncRNAs increased gradually with cap color darkness, indicating that the differentially expressed genes among the three strains were closely related to cap color formation.Circos plots showed the expression levels of mRNAs and ncRNAs, and their locations in chromosomes, with mRNAs exhibiting a more uniform distribution, while ncRNAs were expressed in specific locations in different strains, revealing significant RNA distribution differences in the genomes of the different strains (Figure 2C).

Identification and Functional Analysis of DEmRNAs
Pearson analysis of mRNA expression levels was performed for each of the three strains, and the results showed that the Pearson correlation coefficients among biological replicates were >0.88, indicating sufficient reproducibility for subsequent analyses (Figure 3A).According to the screening criteria, a total of 3604 DEmRNAs were identified in the P358 vs. PHY, P406 vs. PHY, and P358 vs. P406 comparison groups.A total of 1618 mRNAs were differentially expressed in P358 vs. PHY, of which 1130 were up-regulated and 488 were down-regulated in PHY.A total of 2580 mRNAs were differentially expressed in P358 vs. P406, of which 1567 and 1013 were up-and down-regulated in P406, respectively.A total of 1396 DEmRNAs were detected between P406 vs. PHY, of which 721 were upregulated and 675 were down-regulated in PHY (Figure 3B,C).The number of up-regulated genes was greater than that of down-regulated genes in strains with darker color compared to that with white caps, indicating that the up-regulated genes have active roles in the color formation process.A clustering heatmap of DEmRNAs is presented in Figure 3D, illustrating that the white strain, P358, was separated from the dark strains, P406 and PHY, and that some DEmRNAs were highly expressed in the dark strains, P406 and PHY, while they were present at low levels in the white strain, P358, indicating that these genes are important in dark color formation.
GO and KEGG enrichment analyses were performed to explore the potential functions of DEmRNAs.The top 15 most abundant GO terms in each category are shown in Figure 3E, comprising metabolic processes, cellular processes, and signaling in the Biological Process (BP), category, cell and membrane components in the Cellular Component (CC) category, and catalytic activity and binding in the Molecular Function (MF) category.KEGG metabolic pathway analysis showed that DEmRNAs were mainly enriched in steroid biosynthesis, biosynthesis of antibiotics, nitrogen metabolism, starch and sucrose metabolism, sphingolipid metabolism, and tyrosine metabolism, among other pathways (Figure 3F).Ten DEmRNAs, Pleurotus_ostreatus_newGene_65, g1302, g1734, g4221, g4279, g5806, g6479, g7919, g8533, and g9059, were enriched in tyrosine metabolism, which is a pathway important for melanin production.These DEmRNAs may affect melanin production and play vital roles in P. cornucopiae cap color formation.

Identification and Functional Analysis of DElncRNAs
A total of 520 DElncRNAs were identified in the P358 vs. PHY, P406 vs. PHY, and P358 vs. P406 comparison groups.The number of DElncRNAs between P358 and PHY was 223, of which 170 were up-regulated and 53 were down-regulated in PHY.The number of DElncRNAs between P358 and P406 was 404, of which 285 were up-regulated and 119 were down-regulated in P406.The number of DElncRNAs between P406 and PHY was 169, of which 53 were up-regulated and 116 were down-regulated in PHY (Figure 4A,B and Supplementary Table S4).A clustering heatmap of DElncRNAs showed that the three replicate samples of each strain were clustered together and the expression pattern of DElncRNAs in PHY was closer to that of P358 (Figure 4C), suggesting that the DElncRNAs may have been correlated with the formation of light color.
lncRNAs can cis-regulate the expression of neighboring upstream and downstream genes.GO enrichment analysis of DElncRNA target genes showed that highly enriched GO terms included metabolic processes, cellular processes, signaling, and biological regulation in the BP category; cells and membranes in the CC category; and binding in the MF category (Figure 4D).KEGG enrichment analysis showed that DElncRNA target genes were enriched in the glycosphingolipid biosynthesis-globo series, sulfur metabolism, arachidonic acid metabolism, other glycan degradation, thiamine metabolism, synthesis and degradation of ketone bodies, lysine biosynthesis, and lipoic acid metabolism, among other pathways (Figure 4E).These results suggest that DElncRNAs regulating these metabolic pathways may have important roles in cap color formation.

Identification and Functional Analysis of DEmilRNAs
A total of 29 DEmilRNAs were identified among the three P. cornucopiae strains, P358, PHY, and P406.Among them, the number of DEmilRNAs between P358 and PHY was 14 (10 up-and 4 down-regulated in PHY); that between P358 and P406 was 24 (14 up-and 10 down-regulated in P406), and that between P406 and PHY was 9 (6 up-and 3 down-regulated in PHY) (Figure 5A,B and Supplementary Table S5).A clustering heatmap of DEmil-RNAs showed that the three replicate samples from each strain clustered together, and that the white strain, P358, was separated from the dark strains, P406 and PHY (Figure 5C), suggesting that DEmilRNAs may be positively associated with dark cap color formation in P. cornucopiae.To further study the biological functions of DEmilRNAs, TargetFinder software was employed to predict potential target genes.A total of 125 putative target genes, derived from 29 DEmilRNAs, were identified and used for GO and KEGG pathway analyses.GO analysis showed that DEmilRNA target genes were enriched in metabolic process, cellular process, biological regulation and signaling in the BP category; cell and membrane in the CC category; and catalytic activity and binding in the MF category (Figure 5D).KEGG enrichment analysis showed that the DEmilRNA target genes were enriched in RNA polymerase, the phosphatidylinositol signaling system, inositol phosphate metabolism, phagosome, arachidonic acid metabolism, amino acid metabolism, RNA degradation, and endocytosis, among other pathways (Figure 5E).Together, the results of GO and KEGG analyses indicated that DEmilRNAs regulate genes that participate in these signaling systems and may have important roles in cap color formation.

Identification and Functional Analysis of DEcircRNAs
A total of 20 DEcircRNAs were identified among the three P. cornucopiae strains, P358, PHY, and P406.The number of DEcircRNAs between P358 and PHY was 6 (4 up-and 2 down-regulated in PHY), that between P358 and P406 was 16 (9 up-and 7 down-regulated in P406), and that between P406 and PHY was 13 (9 up-and 4 down-regulated in PHY) (Figure 6A,B and Supplementary Table S6).Clustering heatmap analysis of DEcircRNAs showed that the three replicates of each strain were clustered together, while strain P406 was separated from P358 and PHY (Figure 6C), suggesting that DEcircRNAs may be closely correlated with light cap color formation.
GO and KEGG analysis of DEcircRNA targets were conducted to further understand the potential functions of DEcircRNAs.GO enrichment analysis showed that DEcircRNA target genes were enriched in cellular processes, meiotic cell cycle, protein ubiquitination, translation, and protein phosphorylation in the BP category; small ribosomal subunits, cellular components, and cytoplasm in the CC category; and protein tags, structural constituents of ribosomes, protein serine/threonine kinase activity, RNA binding, and ATP binding in the MF category (Figure 6D).KEGG metabolic pathway analysis showed that DEcircRNA target genes were enriched in the ribosomal pathway (Figure 6E).

A ceRNA-miRNA-Gene Regulatory Network
To explore the global regulatory network of ncRNAs and mRNAs involved in cap color formation, a ceRNA network connecting mRNAs and ncRNAs was constructed based on miRNA-mRNA and miRNA-circRNA targeting relationships (Figure 7A); a network including four circRNAs, four milRNAs, and ten mRNAs was generated.milRNAs link mRNAs and circRNAs by repressing mRNA expression and being regulated by circR-NAs [38].Four milRNAs (unconservative_scaffold_84_24490, unconservative_scaf-fold_1_1231, unconservative_scaffold_197_34481, and unconservative_scaffold_39_15572) were the main nodes of the network, acting as response elements which interact with ten mRNAs and connect with four circRNAs (scaffold_42:19565|20053, scaf-fold_3:977102|977467, scaffold_51:7348|8628 and scaffold_11:25806|26198).These mil-RNAs, along with circRNAs, regulate the expression of single or multiple mRNAs.GO enrichment analysis of the mRNAs in the network showed that they were mainly enriched in protein phosphorylation, intracellular protein transport, and transmembrane transport in the BP category; membrane components and intracellular components in the CC category; and ATP binding, protein serine/threonine kinase activity, and Ran GTPase binding in the MF category (Figure 7B).For example, g10683, g9233, and g4098 were enriched in the 'integral components of the membrane' category; g8994, g9127, and g1810 were enriched as 'integral components of membrane and transmembrane transport'; and g3400, g6870, g8994, g9127, and g1810 were enriched for 'protein serine/threonine kinase activity', 'ATP binding', 'protein phosphorylation', and 'Ran GTPase activity', which have important roles in color formation.These analyses suggest that the identified circRNA-milRNA-mRNA regulatory network may be related to color regulation in P. cornucopiae.More studies are needed to elucidate the functions of the ceRNAs, milRNAs, and their target genes in cap color formation.

Construction of Co-Expression Networks
Next, we constructed co-expression networks using all DEmRNA, DElncRNA, DEmiRNA, and DEcircRNA data.The results showed that 215 DEmRNAs, 110 DElncRNAs, 9 DEmiRNAs, and 5 DEcircRNAs were included in mRNA-lncRNA-circRNA-miRNA co-expression networks.The co-expression regulatory networks were classified into three categories.In the first category, 5 miRNAs were in the center of the network, and were co-expressed with 143 mRNAs and 27 lncRNAs, indicating that they potentially play a role in connecting the mRNA and lncRNAs (Figure 8A).In the second category, 4 miRNAs, 16 lncRNAs, 2 circRNAs, and 56 mRNAs formed a complex co-expression regulatory network; lncRNAs and circRNAs in this network may impact the expression and activity of miRNAs through interactions, which in turn co-regulate mRNA expression (Figure 8B).In the third category, 16 mRNAs were co-expressed with 3 circR-NAs and 67 lncRNAs, without any milRNAs (Figure 8C).KEGG pathway enrichment analysis of all DEmRNAs in the co-expression network revealed that g10696, g3311, g7069, and g9452 were enriched for the mRNA surveillance pathway; g4531 was enriched for sesquiterpenoid and triterpenoid biosynthesis; g5381, g6946, and g7651 were enriched for cysteine and methionine metabolism; and g757 was enriched for phenylalanine, tyrosine, and tryptophan biosynthesis (Figure 8D).The phenylalanine, tyrosine, and tryptophan biosynthesis pathways are closely related to color formation.These results suggest that the mRNAs and ncRNAs in the co-expression network may interact with one another to form a complex regulatory network involved in the control of P. cornucopiae cap color formation.

Verification of RNA-Seq Results
To confirm the quality of our RNA-seq data and DEmRNA and DEncRNA expression patterns in P. cornucopiae strains with different cap colors, we selected 6 DEmRNAs (Figure 9A), 5 DElncRNAs (Figure 9B), 4 DEmiRNAs (Figure 9C), and 4 DEcircRNAs (Figure 9D) that may play important roles in cap color for analysis by qRT-PCR (Supplementary Table S7).The results showed that expression levels of these genes were largely consistent with those determined by RNA-seq, suggesting that our RNA-seq data are reliable.Further, DEmRNA and DEncRNA expression levels were related to P. cornucopiae cap color depth, suggesting these RNA molecules may be important in P. cornucopiae cap color formation.

Discussion
P. cornucopiae are cultivated and consumed worldwide because of their delicious taste, nutritious value, and high medicinal value, and consumers prefer them to be darker in color.In our research, we discovered significant differences in color of the P358, PHY, and P406 strains through analysis of mushroom caps and color transparency measurement.Pigments are the foundation of color formation.Zhang et al. [39] extracted and identified pigments from three different colored strains of P. ostreatus and discovered that they contain a combination of true melanin and phaeomelanin, with variation in the relative proportions and content of these two pigments resulting in different colors of P. ostreatus.The differences observed in the colors of P. cornucopiae P358, PHY, and P406 may also be attributable to variation in the distribution or amount of melanin.
Some key genes and enzymes in the melanin synthesis pathway have been identified in P. ostreatus, while their associations with ncRNAs remain unclear.Whole-transcriptome sequencing has been widely used in previous studies of color formation mechanisms in plants and animals [40,41].Genes associated with melanin formation have been demonstrated to be highly expressed in young mushrooms [7]; therefore, we used young mushrooms at the earliest stage when color variation appears as the material in this study to assess how mRNAs and ncRNAs act in concert to regulate P. cornucopiae color at the molecular level.
Melanin is produced by a complex biochemical process that begins with tyrosine and its metabolite, dopa, production of which is catalyzed by specific melanogenic enzymes (TYR, TYRP-1, and TYRP-2).Expression of the genes encoding these enzymes is driven by MITF transcription factors, whose activity is regulated by various signaling pathways, such as protein kinase C (PKC), cyclic AMP, and mitogen-activated protein kinase (MAPK)/ERK kinase [42].The receptor tyrosine protein kinase pathway is major cell signaling pathway involving the activation of MAPK, PKC, and phosphatidylinositol 3-kinase, to trigger downstream biological effects [43].Hence, tyrosine metabolism and protein phosphorylation are closely related to color regulation.Zhang et al. [7] found that overexpression and interference of PcTYR changed the color of the young fruiting bodies in P. cornucopiae.Based on whole-transcriptome sequencing data, we identified 3604 DEmRNAs, 520 DElncRNAs, 29 DEmilRNAs, and 20 DEcircRNAs in young mushrooms from three P. cornucopiae strains.GO and KEGG enrichment analysis showed that the DEmRNAs and target genes of DEncRNAs were enriched in tyrosine metabolism, protein phosphorylation, and other metabolic pathways related to color formation, indicating that these DEmRNAs and DEncRNAs may contribute to the color differences among the P358, P406, and PHY strains.Previous study found that there are two pathways for melain synthesis in fungi [5] and the gene involved in the melanin synthesis pathway were identified in P. ostreatus [44].In this research, we also found that four key genes in the melanin synthesis pathway, namely, g11035 (CM), g11056 (4CL5), g7673 (AT_2), and g9683 (PAL1), were differentially expressed in three differently colored P. cornucopiae strains, suggesting important roles for these genes in color formation.
We also created networks of mRNAs and ncRNAs, to explore the functions of multiple ncRNAs in P. cornucopiae coloration regulation.The ceRNA network contained 18 nodes, and mRNAs in the network were enriched for protein phosphorylation pathways, which transfer phosphoryl groups to specific amino acid residues (serine, threonine, and tyrosine) in substrate proteins.Protein phosphorylation has crucial roles in metabolism of numerous amino acids, and is the most important mechanism for regulating and controlling protein activity [45].The four circRNAs in the network resemble sponges, with adsorption effects for the four target miRNAs, which in turn inhibit the binding of miRNAs to mRNAs; hence, mRNAs, miRNAs, and circRNAs regulate one another to influence enzyme activity and protein degradation, by regulating levels of protein phosphorylation, leading to changes in color.
Additionally, we constructed co-expression networks of mRNAs and ncRNAs.The DEmRNAs in the co-expression networks were enriched in the biosynthetic pathways of phenylalanine, tyrosine, and tryptophan.In plants, the phenylalanine biosynthetic pathway is associated with regulation of plant phytochromes [46], and the phenylalanine can be oxidized to tyrosine by phenylalanine hydroxylase.Further, melatonin, a derivative of tryptophan, can directly inhibit the production of melanin in melanocytes [47].Thus the biosynthesis of tryptophan is involved in color formation.In this study, the expression level of the genes enriched in the phenylalanine, tyrosine, and tryptophan biosynthetic pathways showed an gradually increasing or decreasing trend in samples P358, PHY and P406, which was consistent with or opposite to the cap color depth of these three strains.The expression of milRNAs, lncRNAs, and circRNAs in the regulatory networks was correlated with that of mRNAs.But the ways in which these ncRNAs and mRNAs interact with still required further investigation.
Finally, the expression level of the 6 DEmRNAs, 5 DElncRNAs, 4 milRNAs, and 4 circRNAs in the co-expression network, which may be associated with color formation, were verified by qRT-PCR and their contribution to color formation were analyzed.The results showed that the expression levels of these genes were largely consistent with those determined by RNA-seq.The expression levels of the DElncRNAs, DEmRNAs and circR-NAs in the network, except for scaffold_0:1316855|1317048, showed an gradually increasing or decreasing trend in samples P358, PHY and P406, which was consistent with or opposite to the cap color depth of these three strains.These results indicated that the DEmRNAs, DElncRNAs and circRNAs in the network may play an positive or negative role in the color formation.However, the DEmiRNAs in the network, except for uncon-servative_scaffold_37_15178, were higher expressed in the dark strains than the white strain, but they did not have any correlation with the depth of color, possibly due to the fact that milRNAs have multiple target genes, and have multiple functions more than just in color formation.These results lay the foundation for subsequent studies on the mechanisms of cap color formation in P. cornucopiae.

Conclusions
This research found that DEmRNAs and DEncRNAs existed in the P. cornucopiae strains with different cap color, and the regulatory network of DEmRNAs and DEncRNAs was involved in protein phosphorylation, phenylalanine, tyrosine and tryptophan biosynthesis process that play important roles in color formation, revealing that mRNAs and ncRNAs correlate with Pleurotus Cornucopiae color formation.This research offers a new perspective for comprehending the color regulatory network of this organism; however, additional investigations are required to verify the functions of specific genes and ncRNAs and their interactions in the regulatory networks.This study provides a theoretical framework for future investigations into the molecular mechanisms of color formation in P. cornucopiae.

Figure 1 .
Figure 1.(A) Phenotypic color characteristics of P. cornucopiae strains P358, PHY, and P406.(B) Cap color brightness (L values) of P. cornucopiae strains P358, PHY, and P406.Note: The data are expressed as mean values ± standard deviations, and different lowercase letters indicate significant differences (p < 0.05) as determined by ANOVA.

Figure 2 .
Figure 2. Global characteristics of P. cornucopiae whole-transcriptome sequencing data.(A) Total numbers of mRNAs and ncRNAs.(B) Numbers of mRNAs and ncRNAs in different strains.(C) Global display of mRNAs and ncRNAs in different strains (different colors represent different strains; circles of the same color represent mRNAs, lncRNAs, circRNAs, and milRNAs, from outside to inside).

Figure 3 .
Figure 3. Identification of differentially expressed mRNAs in P. cornucopiae strains P358, PHY, and P406.(A) Sperman correlation analysis showing sample distribution among different groups.(B) Venn diagram showing the DEmRNAs in the three groups.(C) Number of up-regulated and downregulated DEmRNAs.(D) Heatmap analysis of DEmRNAs.(E) GO enrichment analysis of DEmR-NAs.Top 15 GO terms according to annotation analysis.(F) KEGG pathway enrichment analysis of DEGs.Top 20 pathways enriched for DEmRNAs.

Figure 4 .
Figure 4. Identification of differentially expressed lncRNAs in P. cornucopiae strains, P358, PHY, and P406.(A) Venn diagram showing DElncRNAs among the three strains.(B) Numbers of up-and down-regulated DElncRNAs.(C) Heatmap analysis of DElncRNAs.(D) GO enrichment analysis of DElncRNAs target genes.Top 15 pathways enriched for DElncRNA target genes.(E) KEGG pathway enrichment analysis of DElncRNA target genes.Top 20 pathways enriched for DElncRNA target genes.

Figure 5 .
Figure 5. Identification of differentially expressed milRNAs in P. cornucopiae strains P358, PHY, and P406.(A) Venn diagram showing DEmilRNAs in the three groups.(B) Numbers of up-and downregulated DEmilRNAs.(C) Heatmap analysis of DEmilRNAs.(D) GO enrichment analysis of DEmilRNA target genes.Top 15 pathways enriched for DEmilRNA target genes.(E) KEGG pathway enrichment analysis of DEmilRNA target genes.Top 12 pathways enriched for DEmilRNA target genes.

Figure 6 .
Figure 6.Identification of differentially expressed circRNAs in P. cornucopiae strains, P358, PHY, and P406.(A) Venn diagram showing DEcircRNAs in the three groups.(B) Number of up-and downregulated DEcircRNAs.(C) Heatmap analysis of DEcircRNAs.(D) GO enrichment analysis of DE-circRNA target genes.(E) KEGG pathway enrichment analysis of DEcircRNA target genes, the numbers in the right of the columns represent the number of genes enriched.

Figure 7 .
Figure 7. ceRNA regulatory network analysis.(A) ceRNA network of circRNA-milRNA-mRNA. circRNA are represented by green nodes, milRNA are represented by purple nodes and mRNA are represented by blue nodes.(B) GO enrichment analysis of ceRNAs.The numbers in the right of the columns represent the number of genes enriched.

Figure 8 .
Figure 8. Analysis of co-expression regulatory networks.(A-C) DERNA co-expression networks.lncRNA are represented by purple square nodes, circRNA are represented by blue circular nodes, milRNAred are represented by prism nodes and mRNA are represented by green triangles nodes.(D) GO enrichment analysis of co-expressed mRNAs in networks.

Figure 9 .
Figure 9. Gene expression analysis by qRT-PCR.Relative expression levels of DEmRNAs (A), DElncRNAs (B), DEmiRNAs (C), and DEcircRNAs (D).The data are expressed as mean values ± standard deviations, and different lowercase letters indicate significant differences (p < 0.05) as determined by ANOVA.