Identification and characterization of miRNAs and PHAS loci related to the early development of the embryo and endosperm in Fragaria × ananassa

The strawberry fleshy fruit is actually enlarged receptacle tissue, and the successful development of the embryo and endosperm is essential for receptacle fruit set. MicroRNAs (miRNAs) and phased small interfering RNAs (phasiRNAs) play indispensable regulatory roles in plant growth and development. However, miRNAs and phasiRNAs participating in the regulation of strawberry embryo and endosperm development have yet to be explored. Here, we performed genome-wide identification of miRNA and phasiRNA-producing loci (PHAS) in strawberry seeds with a focus on those involved in the development of the early embryo and endosperm. We found that embryos and endosperm have different levels of small RNAs. After bioinformatics analysis, the results showed that a total of 404 miRNAs (352 known and 52 novel) and 156 PHAS genes (81 21-nt and 75 24-nt genes) could be found in strawberry seed-related tissues, of which four and nine conserved miRNA families displayed conserved expression in the endosperm and embryo, respectively. Based on refined putative annotation of PHAS loci, some auxin signal-related genes, such as CM3, TAR2, AFB2, ASA1, NAC and TAS3, were found, which demonstrates that IAA biosynthesis is important for endosperm and embryo development during early fruit growth. Additionally, some auxin signal-related conserved (miR390-TAS3) and novel (miR156-ASA1) trigger-PHAS pairs were identified. Taken together, these results expand our understanding of sRNAs in strawberry embryo and endosperm development and provide a genomic resource for early-stage fruit development.


Background
Cultivated strawberry (Fragaria × ananassa Duch.) is one of most popular berry crops in the Rosaceae family, and China has been the world's largest strawberry producer since 1994 [1]. In the botanical sense, the strawberry fruit is neither a berry nor a true fruit; it is actually a group of achenes dotting the receptacle's surface. Total or partial removal of the achenes from the receptacle resulted in failure of the fleshy fruit to grow and ripen or an abnormal shape unless exogenous auxin was supplied [2,3]. Achenes are rich sources of auxin, whose concentration peaks during early development and then declines quickly as fruit growth and ripening progress. This is consistent with the levels of the auxin-related amino acid tryptophan (Trp) and the lack of free auxin in the receptacles during all developmental stages [2,[4][5][6]. This auxin level pattern correlates with auxin promoting receptacle enlargement and inhibiting strawberry ripening [7,8]. In contrast to the auxin level, the abundance of abscisic acid (ABA) in achenes remained constant during early development and then increased at the onset of fruit ripening [9]. Moreover, ABA is self-induced in achenes and coordinates with IAA in the regulation of early fruit development and later ripening in strawberry [9]. As a consequence, achenes are essential to strawberry receptacle enlargement and ripening. Each achene consists of four tissues: a hard and relatively thick ovary wall, a thin seed coat, a single-layer endosperm, and a small embryo [6,10]. Among them, the endosperm is a nourishing tissue that supports embryo growth and germination, and failure of endosperm development will cause embryo arrest [11]. In addition, endosperm may play a more prominent role than the embryo in the biosynthesis of auxin and gibberellin for strawberry fruit set [10]. In this regard, the role of endosperm is widely accepted as being altruistic. However, the embryo also regulates the development of the endosperm; for example, scanning electron microscopy (SEM) photos that show the decomposition of starch granules in the endosperm cells facing the scutellum of the embryo [12], and VP1 genes expressed in or GA/ABA signalling released from the embryo modulate endosperm development [13][14][15]. Therefore, successful communication between the endosperm and embryo is necessary for seed development and has a profound impact on fruit set and quality.
At present, many studies have revealed that small RNAs (sRNAs), such as miRNAs and short interfering RNAs (siRNAs), are key genetic and epigenetic regulators in diverse biological processes related to plant growth, development, and stress responses [16]. miRNAs are the most abundant sRNAs and show wide conservation in plants. Mature miRNAs are 20-22 nucleotides in length, are derived from single-stranded RNAs, and regulate their target genes through either transcript cleavage or translation inhibition [17]. Distinct from miRNA biogenesis, siRNAs are generated from long double-stranded RNAs, and distinct siRNA species from both strands depend on the Dicer proteins involved. Accordingly, siRNAs can be further grouped into several subclasses, including heterochromatic small interfering RNAs (hc-siRNAs), phasiRNA and other endogenous siRNAs [18]. Among them, hc-siRNAs are predominantly 24 nt long and are mainly derived from repetitive sequences or transposons. They direct de novo cytosine DNA methylation at homologous loci through an RNA-directed DNA methylation (RdDM) pathway [19]. PhasiRNAs are a special type of siRNA that rely on the function of miRNAs in their biogenesis. After cleavage of the target miRNA, the 5' fragment of the target mRNA is rapidly degraded and the 3' fragment is converted into dsRNA via the activity of RNA-DEPENDENT RNA POLYMERASE6 (RDR6). Then, the dsRNA is iteratively cleaved by a Dicer protein (DCL4 or DCL5) from the strand containing the cleavage site, producing duplexes of phasiRNAs (DCL4 forms 21-nt phasiRNAs, DCL5 forms 24-nt phasiRNAs). The biogenesis of phasiRNA is guided by miRNAs through either one or two miRNA binding sites [17]. PhasiRNAproducing loci are named PHAS loci, and the first PHAS loci, TAS3, was identified in Arabidopsis as a conserved noncoding genes among land plants. The TAS3 gene targeted by miR390 generates several tasiRNAs that target AUXIN RESPONSIVE FACTOR (ARF) genes [20]. The three components, miR390, TAS3 and ARF, are indispensable regulators in auxin signalling, and the miR390-TAS3-ARF pathway plays an important role in leaf morphogenesis and flower development [21,22]. Pha-siRNAs are also generated from protein-coding genes; for instance, the PENTATRICOPEPTIDE REPEAT (PPR), nucleotide-binding leucine-rich repeat receptor (NB-LRR) and v-myb avian myeloblastosis viral oncogene homologue (MYB) families were shown to generate phasiRNAs in Arabidopsis [23], Medicago [24], apple (Malus × domestica 'Golden delicious') [25], litchi (Litchi chinensis 'Huaizhi') [26], woodland strawberry (Fragaria vesca) [27] and Rosa rugosa [28].
To date, research on sRNAs in strawberry has mainly focused on fruit ripening [29,30] and stress [31], but there is little research on seed (embryo and endosperm) development. In this study, we focused on the developing embryo and endosperm at the early development stage (12 days after pollination (DAP)) by genome-wide characterization of miRNAs, their target genes and the PHAS pathway. Our results demonstrate the large miRNA population and some miRNA-PHAS pairs in the embryo and endosperm, which enhances our understanding of plant sRNAs in strawberry seed development.

Library construction and sequencing
To identify known and novel miRNAs in the embryo and endosperm of strawberry, twelve sRNA libraries derived from reciprocally crossed embryo and endosperm tissues of F. × ananassa cultivar 'Benihoppe' and 'Sweet Charlie' at 12 DAP were constructed and sequenced. Initially, a total of 9,935,499 to 54,785,143 raw reads were obtained for each library (Table S1). After filtering the low-quality reads, adaptor and contaminant sequences, reads with lengths < 16 nt and reads containing polyA in the small RNA fragment, 8,100,514 to 44,609,527 clean reads were obtained for each library (Table S1). We analysed the length distribution of the sequences in the 12 libraries and found that the majority of sRNAs ranged from 20 to 25 nt in length. The concentrated length distribution peaked at 24-nt sRNAs, while 23-nt and 22-nt sRNAs were the second and third most abundant populations in each library, respectively. Furthermore, there were more 24-nt sRNAs present in the embryo than in the endosperm, whereas the 21-22-nt sRNA showed a reversed trend in the embryo and endosperm of BS, which may be due to tissue-specific expression of sRNA in strawberry development (Fig. 1). Moreover, all clean reads were aligned against the GenBank database [32], Rfam 11.0 database [33], and cultivated strawberry genome [34]. An average of 37.16%, 0.76%, 0.14%, and 2.58% reads matched other noncoding RNAs, including rRNAs, tRNAs, snRNAs, and snoRNAs, and these sRNA populations contained more reads in the endosperm than in the embryo at the 12 DAP stage. In addition, a small portion of reads was mapped to protein coding sequences and intron sequences, which are likely to be RNA degradation products ( Figure S1).

Identification of known and novel miRNAs in the developing embryo and endosperm
Clean reads (excluding reads mapped to snRNAs, snoR-NAs, tRNAs and coding sequences or intron sequences) from 12 libraries were aligned against the miRbase 22 database [35]. A total of 352 known miRNAs were identified (Table S2); these miRNAs belong to 46 miRNA families. Twenty-two out of the 46 families are highly conserved in plants, and a few less conserved miRNAs, including miR477, miR482/2118, and miR535, were identified. As shown Fig. 2A, the most abundant miRNA families are miR156 (38 members), miR482 (25 members), and miR171 (24 members). In addition, 52 novel mature potential candidate miRNAs from 12 libraries were identified, of which most were 21 nt in length ( Fig. 2B & Table S2). Pri-miRNAs with stem-loop structures are cleaved into miRNA::miRNA* duplexes, and in this study, we also analysed the miRNA* structure (Fig. 2C & Table S2). The abundance of miRNA* was significantly lower than that of mature miRNAs; in contrast, we found that fan-miR408* and fan-miRN1* were more abundant than their references, fan-miR408 and fan-miRN1. Interestingly, miR408* was also found to have a higher abundance than rco-miR408 in caster been [36].

Expression profiling of miRNAs in the developing embryo and endosperm
To explore known miRNA families that are highly expressed in the development of embryos and endosperm in strawberry, the TPM reads of miRNA families were counted by summing. Among them, miR1511 was the most abundant family, followed by miR156 and miR159 ( Figure S2). With these expression data, we can examine how highly expressed miRNA families are involved in the embryo and endosperm. These high expression levels of known highly expressed miRNA families were divided into five groups according to their expression profile (Fig. 3A). Among them, group I and group V were highly expressed in the endosperm of SB and BS, respectively, Fig. 1 Length distribution and abundance of sRNAs in the reciprocally crossed embryo and endosperm of strawberry. BS and SB represent 'Benihoppe' (♀) × 'Sweet Charlie' (♂) and 'Sweet Charlie' (♀) × 'Benihoppe' (♂), respectively. em: embryo; en: endosperm. For each tissue, three biological replicates (em1, em2, and em3; en1, en2, and en3) were performed and we considered that this was related to variety specificity. The expression levels of group III miRNA families, including miR164, miR167, miR168, miR396, miR477, miR1511, miR6275, miR7125, miR11289, miR11292 and miR11312, decreased in the embryo but increased in the endosperm, indicating that their functions are potentially involved in the development of the endosperm. In contrast, group IV had a high expression level in the embryo from the reciprocal cross, which suggested that these miRNA families play an important role in the embryo. Novel miRNA families with more than 50 TPM reads were found to be highly expressed in the development of strawberry seeds. Of them, 26 novel miRNA families were classified, and miRN12 was the most abundant family, followed by miRN10 and miRN8. According to the different expression levels between the embryo and endosperm, these novel families were divided into three groups (Fig. 3B). Ten miRNA families showed upregulated expression in the endosperm from the reciprocal cross relative to the embryo, so we placed these in group I. In group II, miRN23 showed tissue-specific expression in the embryo of BS, and miRN24 showed strain-specific expression in SB. Group III miRNA families have upregulated expression in the embryo from reciprocal crosses, showing their regulatory function in heart embryo development. To verify the results of sRNA sequencing, tailing-reaction quantitative reverse transcription PCR (qRT-PCR) was performed to determine the expression level of 12 differentially expressed miRNAs in the development of strawberry embryos and endosperm. As shown in Fig. 3C, the expression patterns of miRNAs determined by qRT-PCR were broadly consistent with those from sRNA sequencing data.

Prediction and analysis of miRNA target genes in the embryo and endosperm
To analyse the function of miRNAs in strawberry embryo and endosperm development, we predicted targets of the 404 miRNA candidates, including 352 known miRNAs and 52 novel miRNAs; 2984 target genes were detected (Table S3). The majority of conserved miRNA targets are various transcription factors, including SPL, MYB, ARF, NFYA, AP2, TCP2, and HD-ZIP, which play an important role in regulating plant development [16,37]. Other conserved miRNA targets include NAC (miR164), scarecrowlike protein (miR171), TAS3 (miR390), sulfate transporter (miR395), laccase genes (miR397), umecyanin-like (miR398), ubiquitin-conjugating enzyme or inorganic phosphate transporter (miR399), uclacyanin (miR408) and SPX domain-containing membrane protein/protein TORNADO 2 (miR827). miR162 and miR168 were found to target genes encoding DCL1 and AGO1, respectively, which are key components of the general process of miRNA biogenesis and activity. miR391, miR482, miR2109, miR2118 and three novel miRNA families containing miRN19, miRN25 and miRN29 were found to   [16,38] and have been identified to be targeted by the conserved miRNA of miR394, the less conserved miRNAs including miR11288 and miR11315, and two novel miRNA families, miRN1 and miRN20. Two less conserved (miR11285 and miR11301) and two novel (miRN24 and miRN26) miRNA families targeted the PPR protein. Receptor-like serine/threonine protein kinases are a large family of biological enzymes that participate in the phosphorylation of proteins [39]; these kinases were targeted by several miRNA families, including miR390, miR482, miR535, miR7125, miR11284 and miR11308. Moreover, novel miRNAs targeted a number of functional genes, such as NEDD8-specific protease 1 (miRN3), RNA polymerase II degradation factor 1-like (miRN4), sulfate transporter (miRN5), elongation of fatty acids protein (miRN11), ACT domain-containing protein (miRN14), piriformospora indica-insensitive protein (miRN17), protection of telomeres protein (miRN18), and phosphoenolpyruvate carboxylase (miRN31). To further demonstrate that these identified genes were targeted by miRNA, qRT-PCR was performed to analyse the fan-miR159a-FaMYB101, fan-miR160a-FaARF18, fan-miR164a-FaNAC, fan-miR167a-FaARF8, fan-miR167h-FaARF6, and fan-miR172b-FaAP2 pairs. The qRT-PCR results showed that these 6 target gene expression profiles were opposite those of their corresponding miRNAs (Fig. 4), which was consistent with previous studies; for instance, miR167 directly targets ARF6 and ARF8 and decreases their transcript levels, which is essential for the fertility of both ovules and anthers [40]. miR172 inhibits the expression of AP2, which limits cell division and expansion in Arabidopsis [41], and overaccumulation of miR172 leads to the silencing of AP2 in apple [42].

Identification of PHAS loci in the developing embryo and endosperm
To investigate phasiRNA biogenesis and regulation in embryos and endosperm, PHAS loci were identified based on small RNA reads by using sRNAanno [43]. Hundreds of PHAS loci (251 twenty-one-nucleotide PHAS loci and 798 twenty-four-nucleotide PHAS loci) were identified. When using a stringent phase score (phase score > = 10), 81 21-nt PHAS loci could still be identified. We identified 75 24-nt PHAS loci by using a phase score of 5 (phase score > = 5) (Table S4).
To refine the putative annotation of the predicted PHAS loci, we aligned these PHAS loci to the NCBI Nucleotide Collection (nr/nt) database and the Plant-Rep database (Table S4). Among the 21-nt PHAS loci, some PHAS loci were protein coding genes, such as PPR, TAR2, AFB2, NAC and kinase, which are conserved in plants [44,45] (Fig. 5A). We found that six CM3 genes and one ASA1 gene in strawberry also generated 21-nt phasiRNAs, which play an important role in the shikimic acid pathway [46]. In addition, 36% of the PHAS loci originated from ncRNA, and two identified 21-nt PHAS loci were TAS3 loci; we did not find TAS1 and TAS2 loci in strawberry. The 81 PHAS loci were mapped to the genome using TBtools [47] software, and among them, the majority of PHAS loci were clustered on chromosomes Fvb4-3 and Fvb6-1 (Fig. 5C). Up to 55 PHAS loci overlapped with annotated protein-coding genes of cultured strawberry, including CM1/3, PPR, peptidyl-prolyl cis-trans isomerase FKBP62, TAR2, DNA repair protein, E3 ubiquitin-protein ligase, AFB2, and other proteins (Table S4). A large number (33/55) of PHAS loci overlapped with genes of unknown function, many of which were single copies in the genome, suggesting cis-rather than trans-activity. The remaining 26 PHAS loci were found in intergenic regions of the genome (Table S4).

The expression patterns of PHAS loci in the developing embryo and endosperm
We also examined the expression of each PHAS locus ( Fig. 6 and Figure S3). Among 81 of the 21-nt PHAS loci, 56 showed differential expression in the embryo and endosperm, where 25 and 31 loci accumulated at higher levels in the embryo and endosperm, respectively (Fig. 6A). In this study, we found four PPR loci in strawberry, and only one PPR locus was specifically expressed in embryos (Fig. 6B). Conversely, other conserved PHAS loci, including two TAR2 and one AFB2, showed higher expression in the endosperm (Fig. 6C). In addition, among six CM3 loci, strawberry also had a higher level in the endosperm (Fig. 6D). However, we observed that most of the 24-nt PHAS loci showed embryo-specific expression, and only nine had higher expression in the endosperm, indicating that 24-nt phasiRNAs play important roles in embryo growth and development ( Figure  S3A). To understand the specific expression patterns of 24-nt PHAS loci, sRNA abundance and phasing score at an endosperm-specific zinc transporter protein locus and embryo-specific MYB locus were examined, as shown in Figure S3B and C, respectively.

Discussion
Bidirectional dialogue between the endosperm and embryo is important for fruit set and development [10,[12][13][14][15]48]. However, studies that focused on the sRNAs involved in endosperm and embryo development have not been found for Fragaria.

Abundance of sRNA in the developing embryo and endosperm
In this study, we constructed 12 sRNA libraries from reciprocally crossed embryo and endosperm tissues from strawberry and sequenced the libraries using highthroughput sequencing technology to gain insight into the identified miRNAs and PHAS loci that might be involved in strawberry embryo and endosperm development. Here, we found that the majority of the sRNAs ranged from 20 ~ 25 nt in length, with the 24-nt class being the most abundant in the embryo and endosperm (Fig. 1). The result was similar to previous reports from castor bean [36] and maize [49,50], while in the green fruit stage and red fruit stage of cultivated strawberry, 21-nt sRNAs represented the most abundant sRNAs [51], suggesting that 24-nt sRNA biogenesis is conserved in embryo and endosperm development.

Features of the miRNAs in the developing embryo and endosperm
miRNAs are the most functionally important and most widely studied class of sRNA in plants. Our current study confirmed that 352 known miRNAs belonging to 46 miRNA families were correlated with the embryo and endosperm development processes in strawberry. We also found 52 novel miRNAs that have not been reported elsewhere. Among the miRNA families, the miR156 family was the largest family, with 28 members, followed by miR482 and miR171, with 25 and 24 members, respectively ( Fig. 2A). These results were similar to those of previous studies, albeit with some differences. For example, Huang et al. [49] found that the miR166 and miR167 families were the largest families, with 15 members in maize endosperm, followed by the miR156 (13 members) and miR171 (11 members) families. The differences in these miRNA family members in different plants revealed that they could function in common or unique regulatory pathways. Compared to the miRNA target genes identified in other plants, many conserved miRNAs, such as miR156, miR160, miR166, miR167, miR395, miR397, and miR408, etc., showed similar targets as their homologues in strawberry [16,37]. These conserved miRNAs and their targets usually play fundamental regulatory role in plant growth and development. For instance, the conserved miR156 family is involved in regulation of flowering and the vegetative-to-reproductive transition, by cleaving SPL transcripts and promoting their translational repression [52]. Moreover, miR156 and miR172 are positively regulated by the transcription factors they target, participates in the stability of juvenile-to-adult phase transition in plants [53,54]. In plants, miR165 and miR166 share the same function of targeting and cleaving HD-ZIP III genes, which act against meristem maintenance [55]. miR397 regulated fruit cell lignification in pear by inhibiting expression of laccase genes which is useful for improving pear fruit quality [56]. However, some conserved miRNAs and most novel miRNAs exhibited species-specific targets in cultivated strawberry. For example, ARF transcription factors, which are essential in the early auxin response, were targeted by the miR160 and miR167 families, which is consistent with our results [16,57]. Our study showed that miR160 family members were all more highly expressed in the embryo than in the endosperm and targeted ARF17 and ARF18. In the miR167 family, many members were highly expressed in the embryo, while miR167j, miR167k and miR167l were highly expressed in the endosperm. Their target genes were ARF6 and ARF8. Moreover, we found that ARFs were also targeted by miR159. The results of our study showed that miR159 family members, including miR159m, miR159n, miR159o, miR159p, miR159q and miR159r, showed endosperm-specific expression ( Figure  S5) and targeted ARF23 (Table S3). Previous research showed that ARF23 is a pseudogene in A. thaliana, and its biological function and ability to bind DNA is not  [58][59][60][61]. Among the various functional mechanisms of pseudogenes, antisense transcription-related regulation and competition with miRNA and RNA-binding proteins appear to be common in humans [62,63]. The genetic functions and evolutionary trends for ARF23 in the strawberry genome are unknown. Additionally, it is unclear how the miR159 family that was predicted to target ARF23 regulates the development of strawberry endosperm and embryos. These are interesting issues that require further investigation. In addition, we classified the miRNA families according to their expression level profiles to explore diverse molecular and physiological pathways of miRNAs in strawberry embryos and endosperm. We found that some of the conserved miRNA families, such as miR164, miR167, miR168 and miR395, showed conserved expression in the developing reciprocal cross endosperm, and other conserved miRNA families, including miR159, miR160, miR162, miR166, miR171, miR319, miR395, miR399, and miR535, showed conserved expression in the developing reciprocal cross embryo (Fig. 3A). This result is similar to the research on embryos and endosperm in maize [64].

PhasiRNAs involved in the developing embryo and endosperm
In addition to miRNAs, phasiRNAs have been shown to be involved in developmental regulation and stress responses in flowering plants. Genome-wide identification of phasiRNAs is an important tool that can be used to investigate gene regulation involving small RNAs in legumes [24], soybean [65], tomato [66], litchi [26], Panax notoginseng [67], wheat [68], and pummelo [69]. In the woodland strawberry F. vesca, a total of 515 PHAS loci, including 192 21-nt PHAS loci and 323 24-nt PHAS loci, were identified [27]. In our study, 81 21-nt PHAS loci and 75 24-nt PHAS loci were identified in the cultivated strawberry F. × ananassa genome (Table  S4). However, only 49.38% (40/81) of the 21-nt PHAS loci and 16% (12/75) of the 24-nt PHAS loci were identified as targets of miRNA triggers. Many studies have revealed that a few phasiRNA pathways are highly conserved alongside conserved miRNAs in land plants [50]. TAS3 is arguably the most well-conserved PHAS locus, as it has been identified across a broad range of species, and the miR390-TAS3-ARF pathway is the most conserved and archetypal [44,70]. Xia et al. [38] identified two TAS3 genes, named TAS3L and TAS3S, in F. vesca. Cultivated strawberry also has two TAS3 genes targeted by miR390, but we did not find other TAS genes in the embryo and endosperm, which suggested that TAS3 is primarily conserved in Fragaria. A number of PPRs were shown to generate phasiRNAs and to be conserved, especially in eudicots. The pathway for the generation of PPR phasiRNAs includes a small number of diverse trigger miRNAs. For instance, miR161, miR400, and miR173 in Arabidopsis [20,71,72], miR1509 in legumes [24], miR7122 in Rosaceae [73], and miR11285 in F. vesca [38]. In this study, our research showed that PPR generated phasiRNAs by triggering miR11285 or miR11301 ( Fig. 7C and Table 1); this result was similar to that in F. vesca. We also identified some conserved PHAS loci in strawberry, i.e., TAR2, NAC, and AFB2  5A). As an auxin biosynthetic gene, TAR2 plays a role in the indole-3-pyruvic acid (IPA)-dependent route for the conversion of Trp to IAA, which is the main mechanism for the de novo synthesis of auxin in plants and is required for proper embryo development, lateral root growth and several developmental processes in Arabidopsis [74][75][76][77][78]. As previously described in Arabidopsis, FaTAR2 works as a Trp aminotransferase in strawberry and is involved in auxin biosynthesis. Transient silencing of this gene was accompanied by a diminished responsiveness to auxin in ripening receptacles [48]. NACs are a family of genes specific to plants that play a role in a diverse set of developmental processes; their functions are partially similar to that of TAR2. For example, NAC1 is also involved in the auxinsignalling pathway to promote lateral root development, and in Arabidopsis, miR164-directed regulation of the NAC-domain gene CUC1 is necessary for normal embryonic development [79,80]. In addition, NAC1, as a transcriptional activator, regulates the expression of two auxin-responsive genes, AIR3 and DBP [79]. AFB2 has positive regulatory roles as an auxin receptor in the regulation of fruit development compared with the AFB gene family [81,82] and has the highest expression level in the early fruit development stage in cucumber and the diploid strawberry F. vesca [10,83]. Moreover, AFB2 plays a major role in the control of jasmonic acid conjugation through the ARF6/ARF8 auxin signalling module during Arabidopsis adventitious rooting initiation [84]. In addition, the PHAS genes ASA1 and CM3 were newly identified in strawberry and found to play an important role in the shikimate pathway [46]. The shikimate pathway serves as an important precursor for auxin biosynthesis, and Trp accumulation occurs via this pathway to sustain auxin production [76]. Overall, we suggest that the biosynthesis of auxin-related phasiRNAs is an indispensable regulator of early endosperm and embryo development and of receptacle fruit set in strawberry.

Conclusions
In this study, a total of 404 miRNAs (352 known and 52 novel) could be found in strawberry embryo and endosperm tissues, of which four and nine conserved miRNA families displayed conserved expression in the endosperm and embryo, respectively. We also identified 81 PHAS loci producing 21-nt phasiRNAs and 75 PHAS loci producing 24-nt phasiRNAs. Among them, some auxin signal-related genes, CM3, TAR2, AFB2, ASA1, NAC and TAS3, were found, which demonstrates that PHAS locus-mediated IAA biosynthesis is important for endosperm and embryo during early fruit growth. Additionally, some auxin signal-related conserved (miR390-TAS3) and novel (miR156-ASA1) trigger-PHAS pairs were identified. The above results improve our knowledge of sRNAs in strawberry embryo and endosperm development and provide candidate PHAS loci for the early developmental stage of fruit in F. × ananassa.

Plant materials and growth conditions
Two cultivated strawberry cultivars, F. × ananassa cultivars 'Benihoppe' and 'Sweet Charlie' , were used as materials for sRNA identification and were grown in the Baima Teaching and Scientific Research Base of Nanjing Agricultural University, Nanjing, China. Developing fruits (receptacles with achenes), BS achenes ('Benihoppe' as the female, pollinated with 'Sweet Charlie') and SB achenes ('Sweet Charlie' as the female, pollinated with 'Benihoppe') from reciprocal crosses between 'Benihoppe' and 'Sweet Charlie' were collected at 12 DAP. The embryo and endosperm of a developing ovary (i.e., achene) were hand-dissected under a stereomicroscope, washed with tissue separation buffer [85], and then frozen immediately in liquid nitrogen. Approximately 2500 embryos and endosperm samples were dissected and pooled into one biological replicate, and three biological replicates were collected.

Small RNA library construction and sequencing
After total RNA of embryos and endosperm from each reciprocal cross achene was extracted by a TRIzol reagent kit (Invitrogen, Carlsbad, CA, USA), the small RNAs in a size range of 16-30 nt were enriched by 15% denaturing polyacrylamide gel electrophoresis (PAGE). Then, 3' adapters were added, and the 36-44 nt RNAs were enriched. Next, 5' adapters were ligated to the RNAs. The ligation products were reverse transcribed by PCR amplification, and the 140-160 bp PCR products were enriched to generate a cDNA library and sequenced using Illumina HiSeq ™ 2500.

Small RNA sequencing analysis
Clean tags were obtained by removing the adaptor sequences, filtering out the low-quality tags and eliminating short lengths (< 14 nt). Then, all the clean reads were aligned with small RNA sequences in GenBank databases [32] and the Rfam 11.0 database [33] to identify and filter the ribosomal RNA (rRNA), transfer RNA (tRNA), small cytoplasmic RNA (scRNA), small nuclear RNA (snRNA) and small nucleolar RNA (snoRNA). After removing noncoding RNAs, the clean small RNA sequences were mapped to the F. × ananassa 'Camarosa' genome [34] to obtain unique reads with abundance and location information based on the genome using the SOAP 2.0 program (http:// soap. genom ics. org. cn/). The unique RNA sequences that perfectly matched the genome were subjected to subsequent analysis. Sequence reads that mapped to exons and introns of mRNA were removed to avoid mRNA degradation products, and the tags that mapped to repeat sequences were also removed.

Identification of known and novel miRNAs and their targets
To identify the known miRNAs, the unique RNA sequences that perfectly matched the F. × ananassa 'Camarosa' genome [34] were searched against the mature or precursor miRNAs of cultivated strawberry and mature plant miRNA sequences in the miRBase 22 database [35] by using the sRNAanno database [43]. The small RNA sequences that perfectly matched known miRNAs were considered to be known miRNAs. Potential miRNA precursor sequences were obtained by aligning the reads to the strawberry reference genome, and their secondary structure were predicted by using the RNAfold [86]. The results of secondary structure were analysed to verify characteristic miR-NAs based on the previous criteria [28]. Sequences that did not find matches in the miRBase were identified as novel miRNAs by sRNAanno [43], which uses a set of well-established criteria, as previously reported [87]. Secondary structures of miRNAs are shown by RNAfold [86]. Based on the sequences of the known miRNAs and novel miRNAs, the software PatMatch (Version 1.2) [88] was used to predict target genes. The following criteria were applied: (1) no more than four mismatches between the sRNA and the target; (2) no more than two adjacent mismatches in the miRNA/target duplex; and (3) no mismatches in positions 10-11 of the miRNA/target duplex. (4) The minimum free energy (MFE) of the miRNA/target duplex should be ≥ 74% of the MFE of the miRNA bound to its perfect complement.

Identification of differentially expressed miRNAs in the embryo and endosperm of strawberry
To screen differentially expressed miRNAs during embryo and endosperm development in strawberry, the expression level of miRNAs in each sample was calculated and normalized to transcripts per million (TPM) using miRNA counts. The significant differentially expressed miRNAs were determined by a fold change ≥ 2 and P value < 0.05 for each comparison. The heatmap of miRNA expression was generated using TBtools [47].

qRT-PCR analysis of differentially expressed miRNAs and their target genes
For quantification of mature miRNAs, we used the tailing-reaction qRT-PCR method. Tailing-reaction RT-PCR was performed using a miRNA 1st strand cDNA synthesis kit (Accurate Biotechnology, China) with general RT tailing-reaction primers. TB Green ® Premix Ex Taq ™ (TaKaRa, Japan) and tailing-reaction 3' primer (Accurate Biotechnology, China) were used for qPCR on a QuantStudio 5 system (ABI, USA). For expression analysis of target genes, qRT-PCR was performed using a QuantStudio 5 system (ABI, USA), with a PrimeScript ™ RT reagent Kit (TaKaRa, Japan) for reverse transcription and TB Green ® Premix Ex Taq ™ (TaKaRa, Japan) for qPCR. Three independent biological replicates were used, and 26S rRNA and EF-1α were used as the reference genes for the analysis of the expression of miRNAs and target genes, respectively. The relative expression levels of miRNAs and target genes were calculated by using the 2 −ΔΔCt method. All primers used in this experiment are listed in Table S5.

Identification and analyses of putative PHAS loci and miRNA triggers
To investigate phasiRNA biogenesis and regulation in embryos and endosperm, a genome-wide search for PHAS loci was conducted based on small RNA reads by using the sRNAanno database [43]. Loci with a P value < 0.001 were considered to be reliable PHAS loci. These PHAS loci were refined for putative annotation by aligning their sequences to the NCBI Nucleotide Collection (nr/nt) database (https:// www. ncbi. nlm. nih. gov/) and the PlantRep database [89]. A potential miRNA trigger of PHAS loci was also predicted using the sRNAanno database [43].