Small RNA and Degradome Sequencing Reveal Roles of miRNAs in the Petal Color Fading of Malus Crabapple

The Malus crabapple is an important woody ornamental plant. The fading of petals during its development significantly affects their ornamental value. Petal color is related to anthocyanin content and miRNAs play an important role in the post-transcriptional regulation of anthocyanin synthesis. However, the mechanisms underlying miRNA regulation of petal fading have rarely been studied. Transcriptome and small RNA sequencing of petals from the blooming phases of Malus. ‘Indian Summer’ varieties S1 (small bud), S2 (initial-flowering), and S3 (late-flowering) allowed us to identify 230 known miRNAs and 17 novel miRNAs, including 52 differentially expressed miRNAs which targeted 494 genes and formed 823 miRNA–target pairs. Based on the target gene annotation results, miRNA–target pairs were screened that may be involved in the fading process of Malus crabapple petals through three different pathways: anthocyanin synthesis, transport, and degradation, involving mcr-miR858-MYB1\MYB5 and mcr-miR396-McCHI inhibiting anthocyanin synthesis; mcr-miR167, mcr-miR390, mcr-miR535, and mcr-miR858 inhibiting anthocyanin transport from the cytoplasm to the vacuole by targeting ABC transporter genes (ABCB, ABCC, ABCD, and ABCG); and mcr-miR398 targeting the superoxide dismutase genes (CZSOD2 and CCS) to accelerate anthocyanin degradation. These findings offer a novel approach to understanding the mechanism of petal fading and serve as a reference for other plants with floral fading.


Introduction
Ornamental crabapple is a broad term used for ornamental plants of Malus Mill. of the Rosaceae family. Because of its colorful flowers, it is considered an important woody flowering plant in spring. On one hand, the plant's vibrant flower colors attract insects for pollination, while also increasing the ornamental value of plants. The inhibition of anthocyanin synthesis, transfer processes, and degradation lead to the fading of petals during flower development, thus reducing the ornamental value of plants. This phenomenon has been observed in Malus crabapple [1], Rosa rugosa [2], Paeonia lactiflora [3], Lycoris longituba [4], Silene littorea [5], chrysanthemum [6], and other plants, but is most evident during the flower opening of ornamental crabapple. Transcriptome and metabolomic analyses have been used to study color fading in plant petals. However, the mechanisms by which miRNAs regulate plant petal fading have rarely been reported.
Anthocyanins are flavonoid pigments that are widely distributed in flowers, fruits, leaves, and other plant tissues. They are synthesized via the phenylalanine pathway by enzymes located on the cytoplasmic face of the endoplasmic reticulum. These enzymes are classified into two groups based on the order in which they are synthesized. Early biosynthetic genes (EBGs), which include phenylalanine ammonia-lyase (PAL), 4-coumaryl: CoA ligase (4CL), chalcone synthase (CHS), chalcone isomerase (CHI), and flavanone 3-hydroxylase (F3H). CHS when expressed as an antisense transgene in Petunia hybrida Int. J. Mol. Sci. 2023, 24,11384 3 of 20 did not conduct further studies. Although studies on the regulation of anthocyanins by miRNAs have been reported for many plants, there is still a lack of studies on the color of ornamental crabapples. Systematic studies on miRNA identification during the crabapple petal fading process and miRNA-target regulatory networks have not been reported.
In this study, transcriptome, small RNA, and degradome sequencing were conducted on petals at three different developmental stages of Malus 'Indian Summer,' to explore the miRNA related to the fading process, and the miRNA-target regulatory network during the fading process was constructed via association analysis. These results provide a new idea for studying the fading mechanism of ornamental petals and a reference for other plants with flower fading.

Reference Transcriptome Assembly
Using strand-specific RNA libraries rather than NEB allows for more accurate gene quantitation, localization, and annotation information as well as translation transcripts and single-exon expression levels in each isoform. Transcriptome sequencing yielded 69,581,194 raw reads, with 68,906,804 clean reads (10.34 GB) obtained after strict quality control. The data had an overall error rate of 0.03%. That of Q20 and Q30 were 97.92% and 93.97%, respectively, and their GC content was 46.92%. After assembly and comparison with the reference genome, high-quality reference transcriptome data were obtained for miRNA and target gene prediction.

Annotation of Small RNA and Prediction of Novel miRNA
In total, 5.398 G of data was obtained from nine small RNA libraries of M. 'Indian Summer' flower at three different development stages (S1-1, S1-2, and S1-3; S2-1, S2-2, and S2-3; S3-1, S3-2, and S3-3). An average of 0.60 G of raw reads was obtained from each small RNA library. By removing low-quality reads, reads with an N (N indicates that base information cannot be determined) ratio greater than 10%, reads contaminated by a 5 adapter, reads without a 3 adapter sequence and inserted fragments, and poly A/T/G/C reads, high-quality clean reads were obtained for subsequent analysis (Table 1).  55.26-59.57% could be mapped to the positive-sense strand of the reference sequence, and the rest could be mapped to the antisense strand (Table 3). Furthermore, all alignments of sRNAs on the reference sequences were classified and annotated ( Figure 1). Repeated sequences accounted for the highest proportion (approximately 31.3%, on average). These repeated sequences may bind to different Argonaute proteins to participate in important biological processes, such as DNA methylation and transposon regulation. In contrast, the number of sRNAs containing miRNAs accounted for only 2.13% (S1), 3.81% (S2), and 4.96% (S3), suggesting that miRNAs are involved in regulating flower development and the fading process.   The statistics of the length and quantity of sRNA obtained demonstrated that the characteristic distribution of sRNA in nine samples at the three stages was similar, and the quantity of sRNA at different lengths was significantly different, among which 24 nt sRNAs were the most abundant, accounting for 62.15%, 55.60%, and 46.71% in the three different developmental stages, respectively, followed by 23 nt and 21 nt (Figure 2). The statistics of the length and quantity of sRNA obtained demonstrated that the characteristic distribution of sRNA in nine samples at the three stages was similar, and the quantity of sRNA at different lengths was significantly different, among which 24 nt sRNAs were the most abundant, accounting for 62.15%, 55.60%, and 46.71% in the three different developmental stages, respectively, followed by 23 nt and 21 nt (Figure 2). The statistics of the length and quantity of sRNA obtained demonstrated that the characteristic distribution of sRNA in nine samples at the three stages was similar, and the quantity of sRNA at different lengths was significantly different, among which 24 nt sRNAs were the most abundant, accounting for 62.15%, 55.60%, and 46.71% in the three different developmental stages, respectively, followed by 23 nt and 21 nt (Figure 2). Figure 2. Distribution of small RNA length in S1, S2, and S3.

Identification of Known miRNA and Novel miRNA
A total of 3614 unique sRNAs were obtained using the miRDeep-P2-v1.1.4 and mirDeep2.0.1.2 software to analyze 2,946,473 sRNA sequences containing miRNAs by removing inaccurate Dicer shear sites (miRNA/miRNA* overhang of 1< or >3 nt, as well as miRNA expression of <10). According to the RNAfold prediction, 194 pre-miRNAs containing hairpin structures were obtained, corresponding to 247 miRNAs ( Table 4, Table  S2). Finally, the miRbase22.0 database was searched (E-value < 10), and 230 known miR-NAs were retrieved under the conditions of a mismatch of <4 and a gap of <2.

Identification of Known miRNA and Novel miRNA
A total of 3614 unique sRNAs were obtained using the miRDeep-P2-v1.1.4 and mirD-eep2.0.1.2 software to analyze 2,946,473 sRNA sequences containing miRNAs by removing inaccurate Dicer shear sites (miRNA/miRNA* overhang of 1< or >3 nt, as well as miRNA expression of <10). According to the RNAfold prediction, 194 pre-miRNAs containing hairpin structures were obtained, corresponding to 247 miRNAs ( Table 4, Table S2). Finally, the miRbase22.0 database was searched (E-value < 10), and 230 known miRNAs were retrieved under the conditions of a mismatch of <4 and a gap of <2. Among the 230 known miRNAs, 128 miRNAs with a length of 21 nt were the largest, followed by 24 nt with 43, 20 nt, 22 nt, and 23 nt with 16, 31, and 12 nt, respectively. Among the known miRNAs, the proportion of the four bases was the highest in U (28.27%), followed by G (24.74%), A (23.88%), and C (23.11%), all with a relatively uniform distribution. The results indicated that the first base at 5 had the highest frequency towards U but the highest resistance to G among 230 known miRNAs, which was consistent with the preference of the first base in plant miRNAs ( Figure 3). the known miRNAs, the proportion of the four bases was the highest in U (28.27% lowed by G (24.74%), A (23.88%), and C (23.11%), all with a relatively uniform di tion. The results indicated that the first base at 5′ had the highest frequency toward the highest resistance to G among 230 known miRNAs, which was consistent w preference of the first base in plant miRNAs ( Figure 3). The novel miRNA prediction method was the same as that used for known mi After comparison in miRbase 22.1, 17 novel miRNAs of Malus crabapple were pre including five novel miRNAs and 12 novel miRNAs*. The lengths of these novel m ranged from 20 to 22 nt, among which the number of miRNAs with a length of 21 the largest (nine), and the number of miRNAs with a length of 20 and 22 nt was sev one, respectively. This result was similar to the length distribution of apple miRN the miRBase 22.1 databases. All miRNAs with a length of 21 nt were dominant, w miRNAs with lengths of 23 nt and 24nt were found.In terms of the distribution of A and C bases, G occupies the highest proportion (29.34%), followed by U (24.5 (23.93%), and A (22.22%). In addition, the first base at the 5′ end of novel miRNA strong base preference for U, which was consistent with the preference of the first plant miRNAs (Figure 4). The novel miRNA prediction method was the same as that used for known miRNAs. After comparison in miRbase 22.1, 17 novel miRNAs of Malus crabapple were predicted, including five novel miRNAs and 12 novel miRNAs*. The lengths of these novel miRNAs ranged from 20 to 22 nt, among which the number of miRNAs with a length of 21 nt was the largest (nine), and the number of miRNAs with a length of 20 and 22 nt was seven and one, respectively. This result was similar to the length distribution of apple miRNAs in the miRBase 22.1 databases. All miRNAs with a length of 21 nt were dominant, while no miRNAs with lengths of 23 nt and 24nt were found.In terms of the distribution of A, U, G, and C bases, G occupies the highest proportion (29.34%), followed by U (24.50%), C (23.93%), and A (22.22%). In addition, the first base at the 5 end of novel miRNAs had a strong base preference for U, which was consistent with the preference of the first base in plant miRNAs ( Figure 4).

miRNA Expression Profile and Identification of Differentially Expressed miRNAs
To determine the expression levels of miRNAs in the three different developmental stages, we created a heat map ( Figure 5) after normalizing the miRNA expression levels in each sample. A total of 23,102 and 112 miRNAs were classified as high-, medium-, and low-abundance according to TPM (Table S3). The expression levels of nine miRNAs, including mcr-miR1511c-3p, mcr-miR396b-5p, and mcr-miR396b-3p, were the highest.

miRNA Expression Profile and Identification of Differentially Expressed miRNAs
To determine the expression levels of miRNAs in the three different developmental stages, we created a heat map ( Figure 5) after normalizing the miRNA expression levels in each sample. A total of 23,102 and 112 miRNAs were classified as high-, medium-, and low-abundance according to TPM (Table S3). The expression levels of nine miRNAs, including mcr-miR1511c-3p, mcr-miR396b-5p, and mcr-miR396b-3p, were the highest.
To explore the expression trend of miRNAs in the flower fading process, we found that 62 miRNAs were differentially expressed in the S1, S2, and S3 stages, among which 35 miR-NAs were upregulated and 28 were downregulated (Table S4). According to the results of the differential expression of miRNA ( Figure 6), there were 26 differentially expressed miRNAs in S1 vs. S2, 39 differentially expressed miRNAs in S1 vs. S3, and 20 differentially expressed miRNAs in S2 vs. S3. There were 24 differentially expressed miRNAs shared by S1 vs. S2 and S1 vs. S3, and the numbers of upregulated and downregulated miRNAs were consistent (12 in both cases). The results showed that miRNAs actively respond to the developmental process of flowers, and different miRNA responses may differ at different developmental stages, suggesting that miRNAs play an important role in regulating the fading process of flowers.
To verify the accuracy of the miRNA sequencing results, a total of 11 miRNAs were randomly selected from DEMs for real-time quantitative PCR (RT-qPCR) analysis. By comparing the RT-qPCR results with the TPM results, it was found that the expression trend of these miRNAs between RT-qPCR and TPM was consistent, which proved the reliability of the sequencing results ( Figure 7). To explore the expression trend of miRNAs in the flower fading process, we that 62 miRNAs were differentially expressed in the S1, S2, and S3 stages, among 35 miRNAs were upregulated and 28 were downregulated (Table S4). According results of the differential expression of miRNA ( Figure 6), there were 26 differentia pressed miRNAs in S1 vs. S2, 39 differentially expressed miRNAs in S1 vs. S3, a differentially expressed miRNAs in S2 vs. S3. There were 24 differentially expressed NAs shared by S1 vs. S2 and S1 vs. S3, and the numbers of upregulated and down lated miRNAs were consistent (12 in both cases). The results showed that miRN tively respond to the developmental process of flowers, and different miRNA resp may differ at different developmental stages, suggesting that miRNAs play an imp role in regulating the fading process of flowers. To verify the accuracy of the miRNA sequencing results, a total of 11 miRNAs were randomly selected from DEMs for real-time quantitative PCR (RT-qPCR) analysis. By comparing the RT-qPCR results with the TPM results, it was found that the expression trend of these miRNAs between RT-qPCR and TPM was consistent, which proved the reliability of the sequencing results (Figure 7).

Target Gene Prediction and Functional Enrichment Based on Degradome Sequencing
To reduce the number of false-positive target genes, the degradome sequencing results were used to remove false-positive targets without mRNA cutting sites. A total of 22,232,259 raw reads and 8,699,925 unique raw reads were obtained via the degradome sequencing. Remove short reads (<15 nt) and 3′ adaptors, 8,660,707 unique reads were mapped to 4,997,426 unique reads, which covered 29,591 transcripts (66.23% of 44,677 reference transcripts) ( Table 5). The target genes were predicted using CleaveLand v4.0, and miRNA-target pairs were divided into zero, one, two, three and four types according to DegradomeCategory, and the third and fourth types were eliminated to improve prediction accuracy. Finally, 1641 target genes of 194 miRNAs were obtained, resulting in 3971 putative miRNA-target pairs, of which there were 300, 80, and 3591 target plots of types 0, 1, and 2, respectively. Target plots (T-plots) of these miRNA-target pairs were plotted (https://figshare.com/articles/dataset/The_T-plot_of_all_the_miRNA-tar-get_pairs/23592387, accessed on 28 June 2023). Functional annotation identified 172 transcription factors (TFs), 30 transcription regulators (TRs), 82 kinases, and 211 other genes and predicted putative proteins in 1148 reference genomes among the 1641 target genes.

Target Gene Prediction and Functional Enrichment Based on Degradome Sequencing
To reduce the number of false-positive target genes, the degradome sequencing results were used to remove false-positive targets without mRNA cutting sites. A total of 22,232,259 raw reads and 8,699,925 unique raw reads were obtained via the degradome sequencing. Remove short reads (<15 nt) and 3 adaptors, 8,660,707 unique reads were mapped to 4,997,426 unique reads, which covered 29,591 transcripts (66.23% of 44,677 reference transcripts) ( Table 5). The target genes were predicted using CleaveLand v4.0, and miRNA-target pairs were divided into zero, one, two, three and four types according to DegradomeCategory, and the third and fourth types were eliminated to improve prediction accuracy. Finally, 1641 target genes of 194 miRNAs were obtained, resulting in 3971 putative miRNA-target pairs, of which there were 300, 80, and 3591 target plots of types 0, 1, and 2, respectively. Target plots (T-plots) of these miRNA-target pairs were plotted (https:// figshare.com/articles/dataset/The_T-plot_of_all_the_miRNA-target_pairs/23592387, accessed on 28 June 2023). Functional annotation identified 172 transcription factors (TFs), 30 transcription regulators (TRs), 82 kinases, and 211 other genes and predicted putative proteins in 1148 reference genomes among the 1641 target genes.     The KEGG metabolic pathway analysis of target genes can be classified into metabolism, genetic information processing, environmental information processing, and cellular. There are five major types of processes and organismal systems involved in 287 KEGG The KEGG metabolic pathway analysis of target genes can be classified into metabolism, genetic information processing, environmental information processing, and cellular. There are five major types of processes and organismal systems involved in 287 KEGG metabolic pathways. Among them, spliceosomes (ko03040), protein processing in the endoplasmic reticulum (ko04141), ribosomes (ko03010), glycolysis/gluconeogenesis (ko00010), and plant hormone signal transduction (ko04075) were significantly enriched (Figure 9). In addition, enrichment analysis of DEMs target genes showed that plant hormone signal transduction (ko04075), spliceosome (ko03040), pyruvate metabolism (ko00620), glycolysis/gluconeogenesis (ko00010), and the citrate cycle (TCA cycle) (ko00020) were significantly enriched. Additionally, pathways involved in anthocyanin synthesis and transport were significantly enriched. metabolic pathways. Among them, spliceosomes (ko03040), protein processing in the endoplasmic reticulum (ko04141), ribosomes (ko03010), glycolysis/gluconeogenesis (ko00010), and plant hormone signal transduction (ko04075) were significantly enriched ( Figure 9). In addition, enrichment analysis of DEMs target genes showed that plant hormone signal transduction (ko04075), spliceosome (ko03040), pyruvate metabolism (ko00620), glycolysis/gluconeogenesis (ko00010), and the citrate cycle (TCA cycle) (ko00020) were significantly enriched. Additionally, pathways involved in anthocyanin synthesis and transport were significantly enriched.

miRNA Regulatory Network and Key Components in the Petal Fading Process
The construction of a DEM-target regulatory network plays an important role in studying the molecular mechanisms underlying the post-transcriptional regulation of flower fading in ornamental crabapples. Among the 494 target genes, regulatory network analysis showed that 52 of the 61 DEMs targeted 494 genes, forming 823 miRNA-target pairs. Functional annotation identified 69 transcription factors, 20 kinases, 69 other genes, and 336 predictive proteins. To further study the DEM-target gene network regulating petal fading, we screened 37 genes that were annotated to the MYB, bHLH, WRKY, and other transcription factors, phenylpropanoid biosynthesis (ko00940), ABC transporter (ko02010), and peroxisome (ko04146) pathways, which were targeted by 20 miRNAs from nine families ( Figure 10, Table S5). For example, two MYB transcription factor genes, HF00466 (MYB1) and HF18993 (MYB5), as well as HF06013 (ABC transporterlike) were targeted by mcr-miR858. HF23861 (McCHI) was targeted by mcr-miR396b-5p. The ABC transporter gene HF14440 (ABCC) was targeted by mcr-miR167g-5p and mcr-miR167h-5p. Similarly, HF37268 (ABCB) was targeted by mcr-miR390a-3p, mcr-miR390c-3p, and mcr-miR390d-3p. Additionally, mcr-miR535a-3p targeted the ABC transporter gene HF14280 (ABCD), while mcr-miR535b-3p targeted HF24977 (ABCG). Additionally, mcr-miR398a-3p targeting three genes were annotated to the peroxisome (ko04146) pathway including HF06091 (CZSOD2), HF08261 (CCS) and HF40618 (CZSOD2), all of them belong to the superoxide dismutase family. These miRNA-target pairs might inhibit anthocyanin synthesis and transportation, and accelerate the anthocyanin degradation rate, thus directly or indirectly leading to the fading of petals. synthesis and transportation, and accelerate the anthocyanin degradation rate, thus directly or indirectly leading to the fading of petals.

Discussion
Flower fading throughout development has a significant impact on the ornamental value of this essential woody ornamental plant in the temperate zones of the Northern Hemisphere. Therefore, it is critical to investigate the fading of petals and reveal the underlying molecular regulatory mechanisms. miRNAs are important regulatory factors in plants that regulate gene expression at the post-transcriptional level and play an important role in regulating anthocyanin synthesis in flowers and fruits [27,30,31]. Petal fading is the opposite phenomenon of flower and fruit coloration and is associated with anthocyanin synthesis, transport obstruction, and anthocyanin degradation. However, it is uncertain whether or not this process is regulated by miRNAs in the same manner as that of miRNAs during coloration. As a result, in this study, small RNA sequencing was performed on the petals of ornamental crabapples at different developmental stages to identify miRNAs involved in the fading process, and combined with degradome sequencing, miRNA target genes involved in the regulation network of petal fading were constructed, providing a basis for further research on the fading of ornamental crabapple petals.

Analysis of miRNA Sequence Characteristics during Petal Fading
With the development of high-throughput sequencing technology, small RNA sequencing has become an important tool for miRNA discovery, and is widely used to identify miRNAs in many plants. Yu

Discussion
Flower fading throughout development has a significant impact on the ornamental value of this essential woody ornamental plant in the temperate zones of the Northern Hemisphere. Therefore, it is critical to investigate the fading of petals and reveal the underlying molecular regulatory mechanisms. miRNAs are important regulatory factors in plants that regulate gene expression at the post-transcriptional level and play an important role in regulating anthocyanin synthesis in flowers and fruits [27,30,31]. Petal fading is the opposite phenomenon of flower and fruit coloration and is associated with anthocyanin synthesis, transport obstruction, and anthocyanin degradation. However, it is uncertain whether or not this process is regulated by miRNAs in the same manner as that of miRNAs during coloration. As a result, in this study, small RNA sequencing was performed on the petals of ornamental crabapples at different developmental stages to identify miRNAs involved in the fading process, and combined with degradome sequencing, miRNA target genes involved in the regulation network of petal fading were constructed, providing a basis for further research on the fading of ornamental crabapple petals.

Analysis of miRNA Sequence Characteristics during Petal Fading
With the development of high-throughput sequencing technology, small RNA sequencing has become an important tool for miRNA discovery, and is widely used to identify miRNAs in many plants. Yu [34]. Currently, 322 miRNAs from apples are included in the miRBase database. Although considered an important ornamental plant in the Malus spp., ornamental crabapples still lack relevant reports on miRNA identification and functional studies.
In this study, during the development of ornamental crabapple flowers, high-throughput small RNA sequencing was performed on the petals of S1, S2, and S3. In total, 247 miRNAs were identified, including 230 known and 17 novel miRNAs. According to the miRNA length statistics, both the known and novel miRNAs had the largest number of 21 nt, and the length of apple miRNA in the miRBase database was also mainly 21 nt, indicating that the length characteristics of miRNAs of ornamental crabapple were consistent with the overall characteristics of apples. The first base preference analysis showed that both novel miRNAs and known first ribonucleic acid bases have a strong preference for U and some resistance to G, which is consistent with the behavior of most plants. This bias is attributed to the stability of RISC, the recognition pattern of the AGO protein for the 5 nucleotide of miRNA, and the specificity of the cleavage site [35,36], proving the accuracy of sequencing and prediction results.

Target Gene Prediction
miRNA-mediated mRNA splicing is the main mechanism by which miRNAs function in plants; therefore, the identification of target genes is an important step in the study of miRNA function. Degradome sequencing technology has been widely used to predict miRNA target genes in plants owing to its high efficiency, speed, accuracy, and high throughput. In this study, 1641 target genes of 193 miRNAs (247 in total) were predicted via degradome sequencing, and 3972 putative miRNA-target pairs were identified. The 494 target genes of the 52 DEMs (61 in total) were predicted to form 823 miRNA-target pairs. In addition, the target genes of the 54 identified expressed miRNAs were not detected via degradome sequencing; therefore, we speculated that these 54 miRNAs might exert their functions through translation inhibition.

miRNA Targeting Anthocyanin Biosynthesis Is Involved in Petal Fading
During flower development, the petal color switches from light to dark or from dark to light, affecting the attractiveness of pollinators [25]. However, the color of the petals changes from dark to light (fading), thus seriously affecting the ornamental value of the plants. Studies have shown that anthocyanin biosynthesis is related to petal color changes. Han et al. used HPLC-DAD to determine the anthocyanin content in the petals of S1, S2, and S3 of Malus hupehensis and found that cyanidin-3-galactoside was the most important pigment, which gradually decreased during floral development. The expression levels of the PAL, CHS, CHI, DFR, and ANS gene related to anthocyanin biosynthesis were consistent with the changes in anthocyanin content. During the development of Malus halliana flowers, the floral fading phenomenon is also attributed to the expression levels of anthocyanin biosynthesis-related genes, and the decrease in their expression may be affected by the methylation of the MYB10 promoter region [1,37]. As important regulatory factors, miRNAs are widely involved in the synthesis of secondary metabolites, and their involvement in the regulation of anthocyanin synthesis has been revealed in recent years. In the results of this study, the upregulated expression of mcr-miR858 targets two MYB transcription factors, MYB1 (HF00466) and MYB5 (HF18993), where MYB1 is a positive regulator of anthocyanin synthesis of bayberry and blueberry [38,39]. MYB5 positively regulates the accumulation of pigments and proanthocyanidins during biosynthesis [40]. Moreover, miR156-SPL9 and miR858-MYBL2 are positive regulatory factors that promote the synthesis and accumulation of anthocyanins in Arabidopsis thaliana [26,41]. Brassica rapa BrmiR828-BrPAP1, BrMYB82, and BrTAS4, apple miR172-MYB10, apple mdm-miR858-MdMYB9, and MdMYBPA1 inhibit anthocyanin biosynthesis [31,42,43]. Therefore, the relationship between miRNAs and anthocyanin synthesis depends on the function of their target genes. Therefore, we speculated that with the opening of flowers, the gradually increased expression of mcr-miR858-3p may inhibit the expression of MYB1 and MYB5, thus inhibiting anthocyanin synthesis and, consequently, the fading of crabapple petals.
Many studies have reported that miRNAs not only participate in anthocyanin synthesis by targeting transcription factors but also regulate anthocyanin synthesis by targeting structural genes involved in anthocyanin biosynthesis. Novel_miR_138 targets PsCHI in Paeonia suffruticosa [44] and miR168 targets CHS in Canna [45]; miR2616 and novel-miR25 target F3GT and F3GT7 in peony [46], respectively. These target genes are the key structural genes involved in anthocyanin synthesis. Inhibition of or interference with their expression can directly affect anthocyanin synthesis. The upregulated expression of McCHI (HF23861), an important structural gene in anthocyanin synthesis targeted by mcr-miR396b-5p, was found in the following quantitative analysis. McCHI expression was the highest in S1 but decreased in S3, which is complementary to the expression trend of mcr-miR396b-5p. mcr-miR396b-5p-McCHI inhibits anthocyanin synthesis and plays a key role in flower petal fading. MdCCR (cinnamoyl-coenzyme A reductase gene) was targeted by miR7125, thus inhibiting lignin synthesis and inducing anthocyanin accumulation in apples [47]. In this study, DEM mcr-miR10996-3p targeted cinnamate dehydrogenase McCAD1 (HF22691), which plays an important role in catalytic lignin synthesis [48]. mcR-mir10996-3p-McCAD1 may affect anthocyanin synthesis by competing with anthocyanins as a substrate, causing the petals to fade.

miRNA-Targeting Transporters Are Involved in Petal Fading
After anthocyanins are synthesized in the cytoplasm, they must be transferred to vacuoles for storage and color formation. The efficiency of anthocyanin transport in plants has an important effect on plant color to a large extent [49]. In the analysis of DEMs, the upregulated expression of miRNAs of the mcr-miR167, mcr-miR390, mcr-miR535, and mcr-miR858 families targeting ABC transporters included ABCC (HF14440), ABCB (HF37268), ABCD (HF14280), ABCG (HF24977), and HF06013. There are four main proteins involved in anthocyanin transport: glutathione transferase (GSTs), multi-drug resistance-related protein (MRP), multidrug and toxic compound excretion (MATE), and BTL-homologue of bilirubin translocation enzyme [50][51][52][53]. Among these, GSTs participate in anthocyanin transport by conjugating to ABC transporters. For example, AtABCC2 and Vitis vinifera VvABCC1 participate in anthocyanin vacuolar transport by conjugating with GSH [54,55]. In this study, the upregulation of miRNAs inhibited the expression of ABC transporters, thereby affecting the formation of conjugates with GSTs, which inhibited anthocyanin transport to the vacuoles and participated in petal fading.

miRNA Is Involved in Anthocyanin Degradation, Leading to Petal Fading
The decrease in anthocyanin content is the main reason for petal fading. This process is related to anthocyanin biosynthesis and transport; anthocyanin degradation is another important factor. Compared to synthesis, there are few studies on anthocyanin degradation. Vaknin et al. were inspired by the results of anthocyanin degradation in wine and fruit juice and found that there was a significant increase in peroxidase activity during anthocyanin degradation in Brunfelsia calycina, which was independent of petal aging [24]. During the development of petals, the downregulated expression of mcr-miR398a-3p resulted in increased expression levels of three superoxide dismutase family genes, HF06091 (CZSOD2), HF08261 (CCS) and HF40618 (CZSOD2), in the peroxisome pathway, which accelerated the degradation of anthocyanins and may be a factor leading to the fading of ornamental crabapple petals. However, further experimental verification is required. Nevertheless, these results provide an important reference for studying the miRNAs involved in anthocyanin degradation.

Plant Materials
The materials used in the experiment were obtained from the national repository of Malus spp. germplasm, which is located in Jiangdu District, Yangzhou City, Jiangsu Province (119 • 55 E, 32 • 42 N), China. Through the phenological observation of 105 varieties in the resource nursery, the cultivar M. 'Indian Summer' whose petals faded obviously during the development stages was selected as the research object. M. 'Indian Summer' was bred by Robert Simpson at the Simpson Nursery in the United States and has been reported in the International Ornamental Crabapple Society Bulletin [56]. The phenotypic characteristics of this variety can be found on the website (https://www.malusregister. org/cn, accessed on 28 June 2023), where it is described as deep pink buds and pink petals. Based on the relevant study [1], petal samples were collected from S1 (small bud, deep), S2 (initial flowering, light), and S3 (late flowering, lightest) in April 2021 ( Figure 11). Three clones of this variety were selected from each period as biological replicates. Nine petal samples were placed in liquid nitrogen and stored at −80 • C. Total RNA from each sample was isolated in accordance with a slightly improved CTAB-LiCl protocol [57] and agarose gel electrophoresis was used to detect the degree of RNA degradation and contamination. In addition, an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) was used for quality control, and samples with an RNA integrity index (RIN) of ≥7 were used for subsequent analysis.

Plant Materials
The materials used in the experiment were obtained from the national repository of Malus spp. germplasm, which is located in Jiangdu District, Yangzhou City, Jiangsu Province (119°55′ E, 32°42′ N), China. Through the phenological observation of 105 varieties in the resource nursery, the cultivar M. 'Indian Summer' whose petals faded obviously during the development stages was selected as the research object. M. 'Indian Summer' was bred by Robert Simpson at the Simpson Nursery in the United States and has been reported in the International Ornamental Crabapple Society Bulletin [56]. The phenotypic characteristics of this variety can be found on the website (https://www.malusregister.org/cn, accessed on 28 June 2023), where it is described as deep pink buds and pink petals. Based on the relevant study [1], petal samples were collected from S1 (small bud, deep), S2 (initial flowering, light), and S3 (late flowering, lightest) in April 2021 ( Figure  11). Three clones of this variety were selected from each period as biological replicates. Nine petal samples were placed in liquid nitrogen and stored at −80 °C. Total RNA from each sample was isolated in accordance with a slightly improved CTAB-LiCl protocol [57] and agarose gel electrophoresis was used to detect the degree of RNA degradation and contamination. In addition, an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) was used for quality control, and samples with an RNA integrity index (RIN) of ≥7 were used for subsequent analysis.

Stand-Specific Sequencing for Reference Transcriptome
Stand-specific RNA libraries were constructed with the total RNA of the nine samples using NEBNext ® Ultra™ Directional RNA Library ep Kit for Illumina ® (NEB, Ipswich, MA, USA), in accordance with the manufacturer's instructions. Subsequently, Agilent 2100 Bioanalyzer System was used to determine the library insert size. After qualified library testing, Illumina Novaseq 6000 (Illumina, San Diego, CA, USA) was used for highthroughput sequencing. Clean reads were obtained by removing low-quality reads, the reads containing poly-N and 5′ adaptors, or the reads without 3′ adaptors [58]. The clean reads were compared with the apple reference genome (https://github.com/moold/Genome-data-of-Hanfuapple, accessed on 28 June 2023) [59] to obtain a standard reference transcriptome.

Stand-Specific Sequencing for Reference Transcriptome
Stand-specific RNA libraries were constructed with the total RNA of the nine samples using NEBNext ® Ultra™ Directional RNA Library ep Kit for Illumina ® (NEB, Ipswich, MA, USA), in accordance with the manufacturer's instructions. Subsequently, Agilent 2100 Bioanalyzer System was used to determine the library insert size. After qualified library testing, Illumina Novaseq 6000 (Illumina, San Diego, CA, USA) was used for high-throughput sequencing. Clean reads were obtained by removing lowquality reads, the reads containing poly-N and 5 adaptors, or the reads without 3 adaptors [58]. The clean reads were compared with the apple reference genome (https: //github.com/moold/Genome-data-of-Hanfuapple, accessed on 28 June 2023) [59] to obtain a standard reference transcriptome.

Small RNA Sequencing and miRNA Identification
The NEBNext Multiplex Small RNA Library Prep Set for Illumina ® (NEB, USA) was used to construct sRNA libraries for the nine samples. LongAmp Taq 2X Master Mix, SR Primer from Illumina, and index (X) primers were used for PCR amplification. The PCR products were purified and recovered on 8% polyacrylamide gel (100 V, 80 min). Library quality was assessed on Agilent Bioanalyzer 2100 System using DNA High Sensitivity Chips and the NovaSeq 6000 platform was used for sequencing. Library construction and high-throughput sequencing were performed by NovoGene Co., Ltd., Beijing, China. Clean reads were obtained after the quality control of raw reads and compared with a standard-specific reference transcriptome to analyze sRNA distribution and classification using the Bowtie2-2.1.0 software [60] at a 0 mismatch base. Small RNAs from mRNA exons, introns, degraded fragments, protein-coding genes, repeat sequences, rRNA, tRNA, snRNAs, and snoRNAs were excluded to improve the accuracy of miRNA prediction.

Identification, Quantification, and Differential Expression Analysis of miRNA
To avoid a high false positive rate in direct miRBase prediction, miRNAs were identified via prediction and then compared ( Figure S1). The first step was the prediction of the pre-miRNAs. miRDeep-P2-v1.1.4 [61] and mirDeep2.0.1.2 [62] were used to predict the secondary structure of the tag sequences containing sRNA, the Dicer cleavage site, and the minimum free energy. In addition, the number and distribution of pre-miRNA tags were determined and non-conforming pre-miRNA tags were removed. Then, according to the annotation standards of plant miRNAs, the pre-miRNAs were manually proofread [63]. Finally, the sRNA tags containing pre-miRNA secondary structures, miRNAs, and miRNAs* were obtained. RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/ RNAfold.cgi, accessed on 28 June 2023) was used to predict the hairpin RNA sequences in the second step. The miRNAs and miRNA* were compared with plant miRNAs in the miRbase22.0 database (http://www.mirbase.org/, accessed on 28 June 2023) [64]. The known miRNAs and pre-miRNAs were classified and named for comparison, and the remaining miRNAs were considered novel miRNAs. Finally, the miDeep2 quantifier.pl was used to obtain the length, base bias on the first position, and base on each position of all the identified miRNAs.
miRNA expression levels were estimated via transcripts per million (TPM). Differentially expressed miRNAs (DEMs) were identified using the DESeq R package (version 1.24.0) with adjusted p-values using the Benjamini and Hochberg methods. A corrected p-value of <0.05 was set as the threshold for screening differentially expressed genes.
Total RNA was extracted from petals in accordance with the method described by Huang et al. (2020) [65]. DEMs were verified using RT-qPCR, and their expression patterns were analyzed. Reverse transcription of miRNAs was performed using Mir-XTM miRNA First-Strand Synthesis and TB Green RT-qPCR (Takara Bio, Inc., Mountain View, CA, USA). Specific primers and the universal primer mRQ3 were used for the RT-qPCR analysis. Specific primers for miRNA RT-qPCR are listed in Table S1, and all primers were synthesized by GenScript Biotechnology Co., Ltd., Nanjing, China. The reaction was performed in accordance with the instructions of PowerUp SYBR Green Master Mix (ABI, Carlsbad, CA, USA). The reference gene, which was selected from small RNA sequencing, was used as a control for normalization with the 2 −∆∆Ct method.

Prediction and Screening of Target Genes via Degradome Sequencing
Poly(A) RNA was purified from total plant RNA (20 µg) using poly T oligo-attached magnetic beads after two rounds of purification. The 5 adapters were ligated to the 5 end of the 3 cleaved mRNA fragments from miRNA-induced cleavage by RNA ligase. Biotinylated random primers and mRNA were mixed for reverse transcription and amplified via PCR to construct a cDNA library with an average insert size of 200-400 bp. Finally, we performed 50 bp single-end sequencing on Illumina HiSeq 2500 (LC Bio, Hangzhou, China) following the vendor's recommended protocol.
Clean reads were used to validate the miRNA-target pairs after quality control and noncoding RNA removal. Targets were predicted using CleaveLand v4.3 program [66] and Oligomap was used to map the degradome data to the reference transcriptome and construct a degradome density file. Standard sequences that were valuable for the degradome density files were compared in the NRPM database (per million reads) to remove redundancy. Target genes paired with miRNA sequences were predicted using GSTAr. Finally, the results of the two software programs were integrated to determine the common mRNA that was the target of the miRNA. The T-plot of the miRNA-target pairs was plotted based on degradome density files. Based on the abundance of reads at the cleavage site (RCSs), the miRNA targets were divided into five categories from high to low reliability: 0, 1, 2, 3, and 4.

Target Gene Function Enrichment
Gene ontology (GO) and KEGG enrichment analysis were used on the target gene candidates of DEMs to investigate the biological significance of miRNA. The GOseq R package was used to realize the Wallenius non-central hypergeometric distribution for the GO enrichment analysis. KOBAS software (version 2.0) was used to test the statistical enrichment of target genes in the KEGG pathway [67]. The p-value (cut-off 0.05) was used to detect the enrichment of GO terms and KEGG pathways.

Conclusions
In this study, miRNAs were found to be involved in regulating the fading of ornamental crabapple flower petals. A total of 247 miRNAs, including 230 known and 17 novel miRNAs, were identified using transcriptome and small RNA sequencing. Combined with degradome sequencing, we identified 194 miRNAs targeting 1641 target genes, forming 3971 putative miRNA-target pairs. The results of differentially expressed miRNA analysis showed that 52 DEMs targeted 494 target genes, among which 20 targeted 38 genes that may play an important role in flower fading. Through the miRNA-target regulatory network, it was found that miRNAs participate in the fading of flower petals by inhibiting anthocyanin synthesis via targeting transcription factors and key structural genes. It is also possible to inhibit anthocyanin transport from the cytoplasm to the vacuole by targeting the ABC transporter gene and accelerate anthocyanin degradome by targeting the peroxidase gene, which participates in the fading of petals. This study provides a novel approach to revealing the fading mechanism of ornamental crabapple petals and provides a reference for other plants exhibiting this fading phenomenon.