Comparative Transcriptome Analysis Uncovers the Regulatory Roles of MicroRNAs Involved in Petal Color Change of Pink-Flowered Strawberry

The pink-flowered strawberry is popular in China due to its high ornamental value. In the present study, sRNAome, transcriptome, and degradome sequencing were performed to understand the functions of microRNAs (miRNAs) and their target genes during flower development in pink-flowered strawberry. Nine small RNA libraries and a mixed degradome library from flower petals at different developmental stages were constructed and sequenced. A total of 739 known miRNAs and 964 novel miRNAs were identified via small RNA sequencing, and 639 miRNAs were identified to cleave 2,816 target genes based on the degradome data. Additionally, 317 differentially expressed miRNAs among the various stages of flower development were identified, which regulated 2,134 differentially expressed target genes. These target genes were significantly enriched in the transcriptional regulation, phenylpropanoid biosynthesis, and plant hormone signal transduction pathways. Furthermore, integrated microRNAomic and transcriptomic analyses suggested that 98 miRNAs targeted several transcription factors, including MYBs (26), bHLHs (12), NACs (14), and SPLs (19), related to anthocyanin accumulation. In addition, 27 differentially expressed miRNAs might affect anthocyanin biosynthesis by regulating 23 targets involved in the hormone signal transduction pathway. The quantitative real-time PCR (qRT-PCR) analysis confirmed the expression changes of 21 miRNA-target pairs. Furthermore, the transient expression of candidate miRNAs was performed in the pink-flowered strawberry cultivar “Fenyun” at the bud stage. Introduction of FamiR156a, FamiR396e, and FamiR858_R-2 in the “Fenyun” increased flower color intensity, while transient expression of FamiR828a decreased flower color intensity. Overall, the present study uncovers the regulatory functions of microRNAs, including anthocyanin biosynthesis, hormone signaling, and regulation factors during flower development and coloration in pink-flowered strawberry. This work expands the knowledge of miRNAs affecting coloration in strawberry and provides rich resources for future functional studies.


INTRODUCTION
Strawberry (Fragaria × ananassa), belonging to the Rosaceae family, is one of the most significant horticultural crops grown globally with high economic value. It is rich in anthocyanins which have strong antioxidant and anti-cancer effects (Garcia-Alonso et al., 2005;Speer et al., 2020). The pink-flowered strawberry, an intergeneric hybrid of Fragaria and Potentilla, was sold as ornamental plants and loved by consumers worldwide due to its colorful petals and fruit production (Xue et al., 2015). Moreover, the plants of pink-flowered strawberry are usually five times costlier than those of white-flowered strawberry because red represents happiness, enthusiasm, and auspiciousness in Chinese traditional culture. Understanding mechanism behind the red petal formation will assist in breeding novel pinkflowered strawberry cultivars.
The petals and fruits of the pink-flowered strawberry are red, but their main anthocyanins are different. The cyanidin-3-glucoside and pelargonidin-3-glucoside are the major anthocyanins in the petals and fruits, respectively (Kosar et al., 2004;da Silva et al., 2007;Xue et al., 2016). Extensive studies have been carried out on the anthocyanin metabolic pathway in strawberry fruits and several genes within the pathway, including FaCHS (Hoffmann et al., 2006), FaCHI (Deng and Davis, 2001), FaPAL (Pombo et al., 2011), FaF3H, FaANR, FaANS, FaFLS, FaLAR (Almeida et al., 2007Jiang et al., 2013), FaUFGT (Salvatierra et al., 2010), FaMYB1 (Aharoni et al., 2001), and FaMYB10 (Lin- Wang et al., 2014) have been identified successively. Transcription factors (TFs) are crucial regulators of anthocyanin accumulation in plants. Among them, MYB/ bHLH complexes play a key role. Other TFs, such as SPL, NAC, ARF, and WRKY have also been reported to play a regulatory role in anthocyanin biosynthesis (Gou et al., 2011;Qi et al., 2011;Maier et al., 2013;Zhou et al., 2015). Several biochemical and molecular biology studies have explored the mechanism of color formation at the transcription level in pink-flowered strawberry (Xue et al., 2019;Hong et al., 2021). Liu et al. (2021) identified three R2R3-FaMYB genes, FaMYB28, FaMYB54, and FaMYB576 based on the genome-wide analysis of the R2R3-MYB gene family in F. × ananassa, which regulated anthocyanin biosynthesis in the petal of pink-flowered strawberry. However, the post-transcriptional regulatory role of miRNAs in flower development and color change remains unclear in pink-flowered strawberry.
Plant miRNAs are small non-coding RNAs that are 20-24 nucleotides long and transcribed by non-coding genes. The miRNAs play a pivotal regulatory role in numerous biological processes, including stress response, hormone signal transduction, and metabolite biosynthesis, by negatively regulating key genes at the post-transcriptional stage (Singh and Sharma, 2017;Cai et al., 2019;Das et al., 2021). They are complementary to target mRNAs and can promote the degradation or translation repression of target mRNAs. In the past few decades, many miRNAs have been reported to regulate anthocyanin biosynthesis in plant, such as miR156 (Gou et al., 2011;Liu et al., 2017), miR7125 (Hu et al., 2021), miR828 (Yang et al., 2013;Zhang et al., 2020), miR477 (Dong et al., 2021), and miR858a . However, studies on the role of miRNAs in regulating strawberry fruit color formation are rare. In order to clarify the posttranscriptional regulation mechanism of flower color formation, flower buds of pink-flowered strawberry cultivar "Sijihong" at three developmental stages were selected to perform a combined analysis of microRNAome, transcriptome, and degradation. These results would provide additional potential targets for the bioengineering and future breeding of new pink-flowered strawberry germplasm.

Plant Materials and Total RNA Extraction
The pink-flowered strawberry cultivar "Sijihong" and released by Shenyang Agricultural University in China was used in this study. According to the degree of flower opening and the change in petal color, the developmental stages were categorized into the young bud stage (PF_L), the coloration beginning stage (PF_Z), and the big bud stage (PF_D; Figure 1A; Xue et al., 2019). Petals of these three developmental stages with obviously different colors were collected for miRNA and degradome analyses, with three biological replicates per stage. The pink-flowered strawberry cultivar "Fenyun" was used to perform transient expression of candidate miRNAs and determination of anthocyanin contents. All samples were immediately frozen in liquid nitrogen and stored before −80°C for RNA extraction.
Total RNA was extracted using a modified CTAB method and treated with RNase-free DNase I (Takara, Dalian, China) to remove genomic DNA contamination. The RNA purity was measured by agarose gel (2.0%) electrophoresis. The Nanodrop 2000c spectrophotometer and Agilent 2100 were used to quantify the RNA amount and assess the RNA integrity with RIN number > 7.0, respectively.

Small RNA Libraries Construction and Bioinformatic Analysis
Approximately, 1 μg of total RNA was used to prepare the sRNA library using the TruSeq Small RNA Sample Preparation Kit (Illumina, United States). Different ends of RNA molecules were ligated to 3′ and 5′ adapters. The ligated RNAs were reverse-transcribed to single-stranded cDNA. The small RNA libraries were used for 50 bp single-end sequencing on an Illumina Hiseq 2500 (LC Bio, China), following the vendor's recommended protocol. Nine libraries with three biological replicates were constructed from the "PF_L, " "PF_Z, " and "PF_D" samples. The sRNAome sequencing data can be found in NCBI Gene Expression Omnibus (GEO) database (https://www.ncbi. nlm.nih.gov/geo/) under the accession number of GSE193522.
An in-house program was used to remove adapter dimers, junk, less complex, common RNA families, and repeats from the raw reads following the ACGT101-miR pipeline (LC Sciences, Houston, United States). Subsequently, a BLAST search was performed to map the 18-25 nucleotide long unique sequences to species-specific precursors in miRBase 22.0 and identify known miRNAs and novel 3 and 5p miRNAs. The mapped pre-miRNAs Frontiers in Plant Science | www.frontiersin.org were aligned against the reference Fragaria vesca 1 to detect their genomic locations. Then, Bowtie was used to align the unmapped unique sequences, and the RNAfold software was used to predict the hairpin RNA structure in the sequence from the 120 nt flanks. The differentially expressed (DE) miRNAs identified based 1 https://www.rosaceae.org/species/fragaria_vesca/genome_v4.0.a2 on the normalized deep sequencing count among the PF_L, PF_Z, and PF_D samples were further analyzed using a Student's t-test and considered statistically significant at p < 0.01 and p < 0.05.

Transcriptome Data Analysis
Our transcriptome data (NCBI GEO database under access no. GSE125777) was reanalyzed by mapping to the F. vesca Frontiers in Plant Science | www.frontiersin.org genome using the HISAT2 software. Then, StringTie was used to assemble the mapped reads with default parameters. The DE mRNAs with a fold change >2 or <0.5 were selected, and the nested linear models were compared using a parametric F-test (p < 0.05) by R package. The heatmap of DE genes was constructed using the TBtools software .
Degradome Libraries Construction, de novo Assembly, and Bioinformatic Analysis The poly(A) RNA was obtained from the plant total RNA (20 μg) after two rounds of purification using poly-T oligoattached magnetic beads. Then, a 3′-adapter random primer was used to reverse the first-strand cDNA. The final cDNA library obtained had an average insert size of 200-400 bp. The degradome library was constructed by mixing samples from three stages and sequenced on the Illumina Hiseq 2500 platform to obtain 50 bp single-end reads. The original data obtained by degradome sequencing were decoupled and filtered to get the clean data via Illumina software. Then, the clean data were compared with the transcriptome data using CleaveLand software, and the matched sequence was recorded to generate the degradation density file. For a given read, only perfectly matched alignments were retained for degradation analysis. Target Finder was used to predict the genes targeted by the most abundant miRNAs (Bai et al., 2020). Further, the miRNA-target regulatory network was analyzed following GO annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using R script.

Data Validation and Gene Expression Analysis by qRT-PCR
To validate the high-throughput sequencing results, qRT-PCR was performed to examine the DE miRNAs and their targets. The primers of each miRNA and target gene were designed using the miRNA Designed V1.01 and Primer 3.0 software, respectively and listed in Supplementary Table S1. Meanwhile, miRNA first Strand cDNA Synthesis Kit (by stem-loop) and HiScript first Strand cDNA Synthesis Kit (+gDNA wiper; Nazyme, China) were used to transcribe miRNAs and their targets, respectively. Then, the qRT-PCR for miRNAs and their targets were performed using a miRNA Universal SYBR qPCR Master Mix and ChamQ Universal SYBR qPCR Master Mix (Nazyme, China), respectively, on a QuantStudio™ Real-time PCR system (Life, United States) with a final volume of 20 μl. Three biological replicates were used per sample, and the U6 and FaDBP (Schaart et al., 2002) genes were used as the reference for miRNAs and target genes, respectively. The petals at young bud stage were used as the control and the relative expression levels were calculated following the 2 -ΔΔCT method. Data were analyzed using the multiple range test of Duncan in SPSS v.17.0 and represented as mean ± SE.

Construction of Transient Expression Vectors of Candidate miRNAs
The precursor sequences of FamiR156a, FamiR396e, FamiR828a, and FamiR858_R-2 were amplified and cloned into the pRI-101 AN (TaKaRa, Dalian, China) binary vector digested by BamH I/Sal I using the One Step Cloning Kit (Vazyme, China; primers listed in Supplementary Table S1). All constructed vectors were confirmed by sequencing at Sangon Biotech Co., Ltd. (Shanghai, China). Then, those vectors were transformed into Agrobacterium tumefaciens (strain GV3101) using the heat shock transformation method.

Transient Expression and Determination of Anthocyanin Contents in Petals of Pink-Flowered Strawberry
The pink-flowered strawberry cultivar "Fenyun" at the bud stage was infiltrated by A. tumefaciens for transient expression.
The Agrobacteria suspension carrying target constructs was diluted to OD 600 = 0.5 with suspension buffer (10 mM MES, 10 mM MgCl 2 , and 100 μM acetosyringone) for petal injection. Agrobacteria containing pRI-101 AN-FamiR156a, pRI-101 AN-FamiR396e, pRI-101 AN-FamiR858_R-2, and pRI-101 AN-FamiR828a constructs were individually introduced into five "Fenyun" plants, and the empty pRI-101 AN vector as well. All plants were covered with clear plastic lids to maintain the humidity, and kept in dark for 24 h. Plants were then transferred to greenhouse with a 16 h light/8 h dark photoperiod at 25°C for recovery. Petals were collected within 3-7 days after agroinfiltration. The expression of miRNAs and target genes were analyzed by qRT-PCR. The anthocyanin content was measured by the pH differential method (Chen et al., 2015). The sample extracts were diluted 20 times (V: V) in pH 1.0 and pH 4.5 buffer, respectively. After 15 min equilibrium, each solution was measured at 530 and 700 nm by UV-Vis spectrophotometry (UV2450, ShimadzuLtd., Japan). The anthocyanin content was expressed by milligram of cyanidin 3-glucoside equivalents (CGE) mg/g fresh weight (FW).

Small RNA Sequencing and Transcriptome Profiling in Petals of Pink-Flowered Strawberry
The pink-flowered strawberry showed a remarkable change in flower color during the petal development, with the deepest flower color at the big bud stage. Small RNA sequencing and transcriptome profiling were performed for different flower developmental stages to explore the expression profiles of miRNA and coding genes, and the molecular changes. After removing low quality sequences and contaminated adaptors, 15,145,282, 9,291,511, and 6,985,639 valid reads were obtained from PF_L, PF_Z, and PF_D, respectively (Supplementary Table S2). The sRNAs showed similar total read length in the nine libraries. The majority of sRNAs were 21 and 24 nt long, with 24 nt ones showing the highest abundance ( Figure 1B).
Meanwhile, RNA-seq generated 69.57 Gb of raw reads from nine libraries. After removing low-quality sequences and contaminated adaptors, 47,701,449, 50,216,274, and 54,682,967 valid reads were obtained from PF_L, PF_Z, and PF_D, respectively (Supplementary Table S3). After quality control, the valid reads were mapped to the F. vesca reference genome. The mapping ratios of each sample were between 75.31 and 77.41%, with an average value of 76.42%, suggesting genomic difference between pink-flowered strawberry and the F. vesca reference genome.

Known and Novel miRNAs Identified During Petal Coloring of Pink-Flowered Strawberry
A total of 739 known miRNAs and 964 novel miRNAs were identified from the small RNA libraries of pink-flowered strawberry, respectively. Of these mature miRNAs, 213 miRNAs were expressed only in PF_L, 76 in PF_Z, and 61 in PF_D. Meanwhile, 665 miRNAs (348 known and 317 novels) were expressed in all three samples ( Figure 1C). These known miRNAs in pink-flowered strawberry could be classified into 85 miRNA families, among which the miR482 family had the largest members of 43, followed by miR171-1 and miR156 families of 27 and 25 members, respectively, while 31 families only identified one member, such as miR481 and miR1862 ( Figure 1D).

Target Predictions of Both Known and Novel miRNAs by Degradome Sequencing
The targets of miRNAs identified in pink-flowered strawberry were predicted via degradome sequencing. A total of 23,726,719 raw reads were obtained (Supplementary Table S4). After removing the reads without the adaptors, 23,362,636 reads (98.47%) were successfully mapped to the F. vesca genome, and a total of 14,406 targets matched the annotated transcripts in the degradation library via Target Finder. CleaveLand identified 2,816 target genes for 639 miRNAs (Supplementary Table S5). Target analysis based on the degradome data showed that more than 80% of miRNAs regulated more than one gene, with the number of target genes varying from 1 to 278. Here, zma-MIR396f-p5_2ss19CT20TC cleaved the largest number of transcripts (278), while 114 miRNAs were cleaved a single target. Meanwhile, a gene could be the target of multiple miRNAs. For example, ath-mir858b-p5 and mdm-miR858 were predicted to regulate MYB82-like genes. The functions of miRNAs were predicted by analyzing the GO and KEGG enrichment of target genes. These analyses showed that miRNAs targeted a wide range of conserved regions in TFs including MYBs, SPLs, floral homeotic protein, and growthregulating factor.

Functions of Differentially Expressed Genes and miRNAs During Petal Coloring of Pink-Flowered Strawberry
To better understand the dynamic flower development, the DE genes and miRNAs in the nine libraries of pink-flowered strawberry were analyzed. As shown in Figure 2A, more DE genes were detected in PF_L vs. PF_D with 2,368 upregulated genes and 2,746 downregulated genes. Meanwhile, 1,940 and 3,741 DE genes were detected in PF_L vs. PF_Z and PF_D vs. PF_Z, respectively. To further analyze the biological functions of DE genes in pink-flowered strawberry, KEGG enrichment analysis was performed in PF_L vs. PF_Z, PF_L vs. PF_D, PF_Z vs. PF_D, and PF_L vs. PF_Z vs. PF_D. The DE genes were significantly enriched in plant hormone signal transduction and primary and secondary metabolic pathways, such as biosynthesis of amino acids, glutathione metabolism, terpenoid backbone biosynthesis, and phenylpropanoid biosynthesis (Supplementary Figure S1). Notably, 118 DE genes related to phenylpropanoid biosynthesis were identified in PF_L vs. PF_D, probably involved in petal color change during flower development. Meanwhile, the plant hormone signal transduction was significantly enriched in each stage of flower development.
Meanwhile, 251, 121, and 116 DE miRNAs were identified in PF_L vs. PF_D, PF_D vs. PF_Z, and PF_L vs. PF_Z, respectively ( Figure 2B). According to the changes in the expression levels, these DE miRNAs could be categorized into nine groups (profiles 0-8; Figure 2C), among which profile 1, profile 6, and profile 8 showed significantly difference (p < 0.05). The profile 1 was continuously upregulated including 63 miRNAs, while profiles 6 and 8 were downregulated, including 75 and 65 miRNAs, respectively. A total of 2,134 DE genes targeted by DE miRNAs were identified in pink-flowered strawberry. These target genes were involved in transcriptional regulation, plant hormone signal transduction, phenylpropanoid biosynthesis, spliceosome, and endocytosis (Supplementary Figure S2). Collectively, these results confirm the significant role of miRNAs in petal color formation of pink-flowered strawberry.

Identification and Characterization of Genes and miRNAs Related to Petal Color Change of Pink-Flowered Strawberry
The key genes involved in flower color formation of pinkflowered strawberry were summarized from the transcriptome database. The annotations of 56 DE genes were related the anthocyanin biosynthesis, including anthocyanin and flavonol synthesis-related structural genes, anthocyanin modified genes, and anthocyanin transporters (Supplementary Table S6). Majority of DE structural genes showed a gradual increase in expression during flower development ( Figure 3A). Two phenylalanine ammonia-lyases (PALs) phenylpropanoid biosynthesis pathway were DE. Among other DE genes, one chalcone synthase (CHS), two chalcone isomerases (CHI), four flavanone 3-hydroxylases (F3H), three flavonoid 3′-hydroxylases (F3'H), and one flavonoid 3′,5′-hydroxylase (F3'5'H) were DE during the upstream of anthocyanin biosynthesis, and three DE dihydroflavonol 4-reductases (DFR), two DE anthocyanidins synthases (ANS), and four DE UDP-glycosyltransferases (UGTs) were in the downstream of anthocyanin biosynthesis. Three DE genes FvCHS (FvH4_7g01160), FvDFR (FvH4_2g39520), and FvANS (FvH4_5g01170), significantly increased in PF_L vs. PF_D. In addition, two types of anthocyanin transporters were obtained. There were 16 DE genes encoding glutathione S-transferases F (GSTF) and 13 ABCC transporter family members, such as FvGSTF (FvH4_1g27460, FvH4_3g16240) and FvABCC genes (FvH4_3g44230) showing significant difference in PF_L vs. PF_D (log 2 fold change > 3; Figure 3B). Several TFs have been confirmed to play an important role in the petal color change. In this study, 173 DE TFs, such as MYBs (37 genes), bHLHs (36 genes), SPLs (13 genes), and NACs (30 genes), were identified in the PF_L vs. PF_Z vs. PF_D. The expression patterns of the DE TFs in the three developmental stages were shown in Figure 3C. Among the DE TFs, the majority exhibited a lower expression in PF_D than the control PF_L stage.

Screening of Target Genes and miRNAs Related to Hormone Signal Transduction Pathways
According to the KEGG analysis, eight signal transduction pathways, including auxin, cytokinin, gibberellins, and abscisic acid, were annotated in the transcriptome data of pink-flowered strawberry. A total of 177 DE genes were identified, among which 38 were involved in auxin signal transduction, 15 in cytokinin signal transduction, three in gibberellin signal transduction, 36 in abscisic acid signal transduction, 32 in ethylene signal transduction, 27 in brassinosteroid signal transduction, 20 in jasmonate signal transduction, and six in salicylic acid signal transduction (Supplementary Figure S4). In the auxin signaling pathway, six DE genes encoding auxin transporters (Aux1), nine encoding auxin-responsive proteins IAA (AUX1/IAA), 13 encoding auxin response factors (ARF), and seven encoding auxin-responsive GH3 showed distinct expression dynamics across three developmental stages. Multiple FvARF genes, were highly expressed in PF_L, and gradually decreased alone the development of flower; However, two FvAUX1 genes (FvH4_6g11840, FvH4_2g22520) showed high fold change between PF_L and PF_D (>3.5), and their expression levels increased gradually alone flower development. Most DE genes within the cytokinin pathway were expressed higher in PF_L than PF_D, such as one FvAHP gene (FvH4_3g01870) and two FvARR genes (FvH4_7g25970, FvH4_1g01300). Their expression level decreased gradually during petal development.
Three DE genes were involved in gibberellin signaling pathway, in which one transcript encoding FvGID (FvH4_6g04960) was differentially expressed during flower development. A total of 29 transcripts of probable protein phosphatase 2C gene (PP2C) were annotated in abscisic acid signaling pathway, most of which showed increased expression levels; FvH4_7g31810 showed a large fold change of 4.17 between PF_L and PF_D. As shown in Supplementary Figure S4, most of the DE genes of the ethylene signal transduction pathway were expressed higher in PF_L than PF_D. Additionally, seven and 14 transcripts were annotated as serine/threonine-protein kinases (CTR) and ethylene-responsive TFs (ERF), respectively. Three transcripts FvH4_2g26080, FvH4_7g26410, and FvH4_4g03470 encoding FvERF genes showed high accumulation at PF_L. Seven transcripts encoding serine/threonine-protein kinase (BSK), 10 encoding xyloglucan endotransglucosylase/hydrolase (TCH4), and four encoding cyclin-D3 (CYCD3) involving in the brassinosteroid signal transduction were identified. The expression Frontiers in Plant Science | www.frontiersin.org level of the FvTCH4 gene (FvH4_4g22090) decreased first and then increased during development, corroborated with its FPKM value. In addition, the expression level of the JAR1 gene (FvH4_7g16040) within the jasmonic acid pathway increased gradually and reached the maximum at the PF_D stage. We further investigated the target genes of miRNAs involved in plant hormone signal transduction based on the degradome data. Twenty-seven DE miRNAs participated in auxin, GA, ethylene, jasmonate, and abscisic acid signal transduction (Figure 6). Meanwhile, 14 miRNAs targeting eight genes of auxin signal transduction accounted for half of the miRNAs related to hormone signal transduction pathway. Four ARF genes were targeted by seven miR160s, including mdm-miR160a, mdm-miR160a_R + 1, mdm-miR160a_L + 3, mdm-miR160a_ L + 1R-1, mdm-miR160a_1ss21AG, csi-miR160c, and csi-miR160c-5p_R + 1_1ss19CT (TPM > 50). From bud stage to full bloom stage, the expression of mdm-miR160a, mdm-miR160a_ L + 1R-1, mdm-miR160a_R + 1, and csi-miR160c-5p_R + 1_ 1ss19CT gradually decreased; in contrast, the expression level of the candidate target genes gradually increased. Different miR393 members targeted FaAUX1 and FaTIR1 (Supplementary Figure S5). In the abscisic acid signal transduction pathway, six miRNAs targeted six transcripts of FvPP2C gene family. In ethylene signal transduction, three miRNAs (mdm-miR172a_R + 1_1ss1AG, csi-MIR156f-p5_1ss19AT, and zma-MIR396f-p5_2ss2AT18CG) targeted FvERF and FvCTR1. Additionally, mdm-miR172a_R + 1_1ss1AG and csi-MIR156f-p5_1ss19AT showed a progressive increase of expression during flower Frontiers in Plant Science | www.frontiersin.org development, leading to reduction of FvERF expression. On the contrary, zma-MIR396f-p5_2ss2AT18CG showed relatively high expression levels during the initial developmental stages, while its target FvCTR1 showed an opposite expression pattern (Supplementary Figure S5). Our data suggest miRNAs could alter the expression of key genes related to auxin, GA, and ABA pathways during flower development. Furthermore, five identified miRNAs were selected and their expression was measured by the qRT-PCR. Along flower development, a progressive decrease in FamiR160a_L + 1R-1 was observed, while the abundance of its target FaARF (FvH4_6g15960.1) increased significantly. Similar miRNA target expression patterns were also found for FamiR160c-5p_R + 1 and FamiR160a_R + 1 (Figure 7A). Indeed, for these FamiR160 members, a significant decrease (p < 0.05) in abundance was seen in PF_L vs. PF_D, and their target FaARF (FvH4_7g19350.1) showed a progressive increase (p < 0.05). The FamiR172a_R + 1 was found highly accumulated at PF_D, opposite to the target FaERF (FvH4_6g15180.1) gene. The pair of FamiR156f-p5 and its target FaCTR1 (FvH4_3g39380.1) exhibited the same trend ( Figure 7B).

Functions of Candidate miRNAs Involved in Anthocyanin Metabolism
We analyzed transient expression of petal and the functions of candidate miRNAs involved in coloration in pink-flowered strawberry "Fenyun. " The petals injected with pRI-101 AN were used as the control. We observed changes of flower color phenotypes within 3-7 days after injection ( Figure 8A). Compared with the control, the petal color turned darker after injecting pRI-101 AN-FamiR156a, pRI-101 AN-FamiR396e, and pRI-101 AN-FamiR858_R-2, while the petal color turned lighter after injecting pRI-101 AN-FamiR828a. Expression levels of candidate miRNAs and corresponding target genes were determined by qRT-PCR. The result showed that expression levels of overexpressed candidate miRNAs were significantly higher than those of the control, while the expression levels of their corresponding target genes were decreased ( Figure 8B). The anthocyanin contents in agroinfiltrated petals of pRI-101 AN-FamiR156a, pRI-101 AN-FamiR396e, and pRI-101 AN-FamiR858_R-2 were significantly more than those in the control (Figure 8C), indicating that transient expression of these miRNAs activated the anthocyanin accumulation, while the content in pRI-101 AN-FamiR828a was less than that in the control, indicating FamiR828a inhibited the anthocyanin accumulation in "Fenyun. "

DISCUSSION
The present study provides the first comprehensive overview of miRNAs expression dynamics during flower development and coloring in pink-flowered strawberry. High throughput sequencing technology provided us with an efficient method to identify miRNAs. A large number of miRNAs have been identified in many plants, including tomato (Karlova et al., 2013), grapevine (Pantaleo et al., 2010), and sweet potato (He et al., 2019;Tang et al., 2020). Up to date, many Fragaria miRNAs have been identified and deposited to miRBase 22.0 (Ge et al., 2013;Li et al., 2019). Xu et al. (2013) obtained 88 known miRNAs and 1,244 new candidate miRNAs involved in fruit senescence of F. × ananassa. In this study, a total of 739 known miRNAs and 964 novel miRNAs were identified in PF_L, PF_Z and PF_D libraries, expanding the strawberry miRNA database, and facilitate the understanding of the regulatory mechanisms underlying coloration in pinkflowered strawberries. Both structural genes encoding anthocyanin biosynthesisrelated enzymes and TFs controlling the expression of structural genes regulate anthocyanin synthesis. In this study, 56 DE structural genes and 173 TFs, such as MYB (37 genes), bHLH (36 genes), NAC (30 genes), and SPL families (13 genes), were identified during flower development. No miRNAs were found to regulate structural genes related to anthocyanin biosynthesis directly. However, miRNA families of FamiR156s, FamiR159s, FamiR396s, FamiR828s, and FamiR858s potentially targeted FvSPLs, FvMYBs, and FvbHLHs. The MYB genes have been identified as one of the major regulators of anthocyanin biosynthesis pathway. Four subgroups of R2R3-MYBs related to the flavonoid synthesis have been identified in Arabidopsis thaliana, including subgroup 4, subgroup 5, subgroup 6, and subgroup 7 (Dubos et al., 2010). In this study, according to degradome data and expression profiles of DE miRNA-targets, a total of eight FvMYBs were targeted by different miRNAs. Through constructing phylogenetic tree with AtR2R3-MYBs in A. thaliana, FvMYB308 and FvGAMYB were found to belong to subgroup 4, which inhibited anthocyanin synthesis, FvMYB114 in subgroup 6 promoted anthocyanin synthesis, and FvMYB4 and FvMYB5 in subgroup 7 controlled the synthesis of flavonoids in petals of pink-flowered strawberry (Supplementary Figure S6A). Therefore, those corresponding miRNAs might have a regulatory function in petal coloring of pink-flowered strawberry including FamiR828a (target FvMYB114 FvH4_2g31020.1), FamiR159d-R-1 (target FvGAMYB FvH4_7g04470.1), and FamiR858_R-2 (target FvMYB308 FvH4_2g01320.1).
Previous studies indicated that AtMYB114, PyMYB114, and MdMYB114 were active regulators of anthocyanin biosynthesis in Arabidopsis (Gonzalez et al., 2008), pear (Yao et al., 2017), and apple . Other studies suggested that MdMYB12 and MrMYB6 negatively regulate anthocyanin accumulation in apple (Mao et al., 2021) and Chinese bayberry (Shi et al., 2021). Here, our degradome analysis showed that FamiRNA828a and FamiRNA858_R-2 targeted FaMYB114 and FaMYB308, respectively. The expression of FamiR828a was opposite to that of FvMYB114, suggesting a negative correlation between FamiR828a and anthocyanin accumulation in the petal of pink-flowered strawberry. Similarity, FamiR858_R-2 negatively regulated FaMYB308, which might effectively influence the target genes involved in anthocyanin biosynthesis in the petal of pink-flowered strawberry. In addition, a previous study demonstrated that miR828 and miR858 promoted anthocyanin and flavanol accumulation in grape by regulating VvMYB114. And in apple, miR858 and miR828 interacted with multiple MYBs to form a complex regulatory network involved in regulating anthocyanin metabolism (Xia et al., 2012;Guan et al., 2014). Therefore, FamiR828a and FamiR858_R-2 identified in this study might play an important role in anthocyanin biosynthesis in the petal of pink-flowered strawberry.
Previous studies have demonstrated that miR156s positively regulate anthocyanin biosynthesis by targeting SPL TFs, and SPLs negatively regulate anthocyanin accumulation through MYB-bHLH-WD40 complex (Gou et al., 2011). Wang et al. (2020) reported that three microRNA-target modules, including miR156-SPLs, miR160h-ARF18, and miR858-MYB39, may control anthocyanin accumulation in miR156 overexpressed transgenic poplar. In our study, the expression levels of FamiR156s gradually increased during flower development, in contrast to their target genes. These results suggest that FamiR156s-SPLs may be involved in anthocyanin accumulation. Hormones are important factors inducing anthocyanin accumulation. The ABA, ethylene, and jasmonate promote the accumulation of anthocyanin, while auxin inhibits anthocyanin accumulation (Loreti et al., 2008;Gao et al., 2021). Wang et al. (2018) reported that MdARF13 as a repressor directly bound to the promoter of MdDFR to regulate anthocyanin biosynthesis in apple. Based on these reports, we hypothesize that some growth regulators might promote flower coloration, and affect anthocyanin accumulation. In our study, FvARF was downregulated while FvERF was upregulated during flower development, and their corresponding miRNAs showed a reverse trend, suggesting the role of miR160s and miR172s in regulating anthocyanin accumulation through the auxin and ethylene signaling pathways. These data collectively supported that miRNAs could alter the expression of key genes related to hormone biosynthesis and were involved in the regulation of flower development and coloring in pinkflowered strawberry.

CONCLUSION
Our comprehensive analysis of the sRNAome, transcriptome, and degradome indicated the regulatory functions of microRNAs in flower development and color changes of pink-flowered strawberry. We analyzed the changes in miRNAs expression level and predicted the target genes of the candidate miRNAs. Our analysis revealed that 98 miRNAs targeted several TFs, including MYBs, bHLHs, and SPLs, related to anthocyanin accumulation. Additionally, 27 DE miRNAs may affect anthocyanin metabolism by regulating 23 genes in the hormone signal transduction pathway. Given the key role of plant microRNAs in various processes, further genetic and molecular evidence based on loss of function and overexpression analyses is needed to verity miRNA functions and its regulatory role in anthocyanin biosynthesis in petal. These results expand the pool of identified miRNAs in strawberry, which may serve as important resources for future research and also provides a foundation for using miRNAs to change strawberry flower color.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
JY and LX conceived the study and drafted the manuscript. ZL and CZ analyzed the sRNAsome data. JZ and YZ analyzed the transcriptome data. HZ, CT, and ZZ analyzed the degradome data. JL revised the manuscript and provided guidance on the whole study. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the National Natural Science Foundation of China (31701964).
Frontiers in Plant Science | www.frontiersin.org