Integrated Analysis of Small RNA, Transcriptome and Degradome Sequencing Provides New Insights into Floral Development and Abscission in Yellow Lupine (Lupinus luteus L.)

The floral development in an important legume crop yellow lupine (Lupinus luteus L., Taper cv.) is often affected by the abscission of flowers leading to significant economic losses. Small non-coding RNAs (sncRNAs), which have a proven effect on almost all developmental processes in other plants, might be of key players in a complex net of molecular interactions regulating flower development and abscission. This study represents the first comprehensive sncRNA identification and analysis of small RNA, transcriptome and degradome sequencing data in lupine flowers to elucidate their role in the regulation of lupine generative development. As shedding in lupine primarily concerns flowers formed at the upper part of the inflorescence, we analyzed samples from extreme parts of raceme separately and conducted an additional analysis of pedicels from abscising and non-abscising flowers where abscission zone forms. A total of 394 known and 28 novel miRNAs and 316 phased siRNAs were identified. In flowers at different stages of development 59 miRNAs displayed differential expression (DE) and 46 DE miRNAs were found while comparing the upper and lower flowers. Identified tasiR-ARFs were DE in developing flowers and were strongly expressed in flower pedicels. The DEmiR-targeted genes were preferentially enriched in the functional categories related to carbohydrate metabolism and plant hormone transduction pathways. This study not only contributes to the current understanding of how lupine flowers develop or undergo abscission but also holds potential for research aimed at crop improvement.


Introduction
Yellow lupine is a crop plant with remarkable economic potential. Because of the symbiotic bond with nitrogen-fixing Rhizobium bacteria it does not need fertilizers, and its protein-rich seeds may be an excellent source of protein for both human consumption and animal feed [1][2][3]. Lupinus luteus flowers are stacked in whorls along the common stem forming a raceme. Pods are formed at the lowest whorls, while the flowers above them fall off [4]. The estimated percentage of dropped flowers is 60% at the 1st Our observations of L. luteus generative development suggest that the fate of flowers (pod set or shedding) is determined on a molecular level during flower development. This study aims to characterize and investigate the role of these important molecules and their target genes during flower development and abscission.
In order to achieve this goal, an integrated analysis of small non-coding RNAs (sncRNA), transcriptome and degradome sequencing data was performed. We identified both known and presumably new miRNAs and siRNAs from flowers at different developmental stages, specifically the lower flowers (usually maintained and developed into pods) and the upper flowers (usually shed before fruit setting). Moreover, in our comparisons of libraries from the upper and lower flowers, differentially expressed miRNAs were found. In order to identify the miRNAs involved exclusively in flower abscission, we compared sRNA libraries from the pedicels of flowers that were maintained on the plant and those that were shed. A transcriptome-and a degradome-wide analysis was carried out to identify the target genes for the conserved or new L. luteus sRNAs. The targeted transcripts were then functionally annotated to outline the putative regulatory network in which these sRNAs might have a role to play. Our results of next-generation sequencing (NGS) analysis indicate that the identified miRNA-targeted modules may be vital in regulating yellow lupine flower development, both generally and depending on the flower location on the inflorescence. Furthermore, these scnRNA also display differential accumulation during flower abscission in this plant.

Sequencing and Annotation of Yellow Lupine sRNAs from Flowers and Flower Pedicels
Flowers collected from the top and bottom parts of the inflorescence were separated into four categories based on the progression of their development, and thus: Stage 1-closed green buds, parts of which were still elongating. Stage 2-closed yellow buds, around the time of anther opening.  Table S1). Flower pedicels from flowers undergoing abscission (FPAB) or maintained on the plant (FPNAB) were also collected, as they had been in our previous study [17]. This division resulted in ten variants of small RNA libraries, which were subjected to single-end deep sequencing performed on the Illumina HiSeq4000 platform (Illumina, Great Abington Cambridge, United Kingdom). After removing low-quality reads, a total of 303,267,263 reads (from 14,186,278 to 15,504,860 reads per library) and 128,060,403 unique reads (from 5,677,701 to 6,990,061 per library) were obtained (Table S2).The length distribution of the small RNAs (15-30 nt) revealed that a length of 24 nt was the most frequent and that of 21 nt was the second most abundant class of the clean and redundant reads (Figure 2), which was compliant with many other RNA-Seq experiments [46][47][48] and correlated with the abundance of siRNAs and miRNAs, respectively.     The unique reads were annotated against Rfam [49,50] and miRBase [51] databases, and from the latter both mature (named in tables as 'miRBase') and precursor sequences (named as 'Hairpin') were taken into account. However, many of them remained unassigned (Table 1). The unique sequences were annotated into different RNA classes against the Rfam database using BLAST [52] such as known miRNAs, rRNA, tRNA, sn/snoRNA and others ( Table 2). A total of 690,436 sRNAs were annotated into all libraries, with the highest number observed in the upper flowers and abscising pedicles. Between these libraries, the most abundant classes were rRNAs and tRNAs, with average values of 26,390 and 3,726 sequences, respectively, followed by snoRNAs and different subtypes of snRNAs with average values ranging from 863 to 876 sequences ( Table 2).  893  496  561  540  591  1125  683  827  748  Intro  1480  1238  1094  1201  1403  1340  1166  1339  1407  1546  Retro  829  800  681  707  792  742  803  766  751  852  U1  415  100  64  103  83  66  62  84  124  124  U2  620  323  263  261  286  294  275  308  312  346  U3  433  244  150  172  169  163  189  184  245  215  U4  248  61  51  63  54  54  58  67  91  82  U5  69  10  12  15  16  16  13  9  13  21  U6  349  81  52  81  64  58  102  76  90  108  Total  45,858

Identification of Known, Conserved miRNAs
After analyzing the results of the alignment against miRBase [51], 394 unique miRNAs containing 366 conserved miRNAs were identified (Table S2). The number of identified miRNAs in each library is shown in Table S1. The identified miRNAs belonged to more than 67 families (Table S2), while most of them belonged to the MIR156, MIR159, and MIR166 families, with more than 35 members in each (Figure 3a). Each discovered miRNA received an identification number in the following format: Ll-miR(number). In case of miRNAs displaying identity to sequences from miRBase, annotation Ll-miR(number)/miRBase annotation is used, for example, Ll-miR224/miR393. Fabaceae species (lower panel). Known miRNA families of L. luteus identified from small RNAseq are listed in the top row. The colors represent relative miRNA families classified into different groups with similar conservation. Blue, yellow, magenta, green and orange represent relative miRNA families with homologs across more than 20, 10-19, 5-9, 2-4 species and in 1 species, respectively.

Evolutionary Conservation of microRNAs Identified in Lupinus luteus
Since this study is the first wide-scale analysis of yellow lupine miRNAs, we decided to explore the evolutionary characteristics of these sequences when compared to the data of almost all [52] Eudicotyledons species present in miRBase [51]. The same analysis was performed exclusively against nine Fabaceae species. As shown in Figure 3b, the 67 known miRNA families exhibited different Fabaceae species (lower panel). Known miRNA families of L. luteus identified from small RNA-seq are listed in the top row. The colors represent relative miRNA families classified into different groups with similar conservation. Blue, yellow, magenta, green and orange represent relative miRNA families with homologs across more than 20, 10-19, 5-9, 2-4 species and in 1 species, respectively.

Evolutionary Conservation of microRNAs Identified in Lupinus luteus
Since this study is the first wide-scale analysis of yellow lupine miRNAs, we decided to explore the evolutionary characteristics of these sequences when compared to the data of almost all [52] Eudicotyledons species present in miRBase [51]. The same analysis was performed exclusively against nine Fabaceae species. As shown in Figure 3b, the 67 known miRNA families exhibited different numbers of homologous sequences in both of the comparisons. Twenty of them were the most conserved ones, i.e., had homologues in over 20 species (Figure 3b, shaded in blue). Our comparison across legumes revealed that 8 miRNA families were highly conserved in this taxon, i.e., had homologues in 5-9 species out of 9 ( Figure  A surprisingly high number (39) of miRNA families identified in yellow lupine flowers were not conserved across Fabaceae, probably due to a still incomplete list of miRNAs in these taxa.

Identification of Novel miRNAs
With the use of the ShortStack software (https://github.com/MikeAxtell/ShortStack/) [53], 28 candidates for novel miRNAs were identified (Table 3). This tool identifies miRNAs based on their mapping against a reference genome. Since there was no genome available for the studied species, we used a transcriptome instead (statistical data on de novo assembly is shown in Table S3). The results obtained were filtered against mature miRNAs from miRBase, and unique sequences received names in the following format: "Ll-miRn(number)", (for example, Ll-miRn1). All of these 28 sequences were 21-24 nt in length, with 68% of them being 21 nt long ( Table 3).
The expression of novel miRNAs was also highly diversified across all the libraries. Ll-miRn26 was present only in the LF1 sample, while Ll-miRn21 was present in all the sRNA libraries and had an expression ranging from 3,982.82 to 11,421.55 RPM (Table S4).

Analysis of the Expression Abundance of Known miRNA Families
Since miRNA expression across all libraries displayed high variation, we put the data into five categories based on the maximum value ( Figure 4). Two miRNAs, namely Ll-miR341/miR319 and Ll-miRn21, showed expression maxima of over 10,000 RPM. The maximum expression of another two, Ll-miR260/miR166 and Ll-miR384/miR396, ranged from 2000 to 10,000 RPM. Thirteen miRNAs showed expression maxima ranging from 100 to 2000 RPM. The most numerous category, with 33 elements was the one for miRNAs with expression maxima ranging from 10 to 100 RPM. Another 24 miRNAs were expressed with the maximal RPM values between 1 and 10. The expression value of the five least abundant miRNAs did not exceed 1 RPM (Figure 4, Table S5).
The expression of novel miRNAs was also highly diversified across all the libraries. Ll-miRn26 was present only in the LF1 sample, while Ll-miRn21 was present in all the sRNA libraries and had an expression ranging from 3,982.82 to 11,421.55 RPM (Table S4).

Analysis of the Expression Abundance of Known miRNA Families
Since miRNA expression across all libraries displayed high variation, we put the data into five categories based on the maximum value ( Figure 4). Two miRNAs, namely Ll-miR341/miR319 and Ll-miRn21, showed expression maxima of over 10,000 RPM. The maximum expression of another two, Ll-miR260/miR166 and Ll-miR384/miR396, ranged from 2000 to 10,000 RPM. Thirteen miRNAs showed expression maxima ranging from 100 to 2000 RPM. The most numerous category, with 33 elements was the one for miRNAs with expression maxima ranging from 10 to 100 RPM. Another 24 miRNAs were expressed with the maximal RPM values between 1 and 10. The expression value of the five least abundant miRNAs did not exceed 1 RPM (Figure 4, Table S5). Numerous reports and studies indicate the importance of phased siRNA not only in stress response mechanisms but also in growth regulation [54]. Therefore, we decided to investigate the role of siRNAs during yellow lupine inflorescence development. To achieve this, ShortStack (https://github.com/MikeAxtell/ShortStack/) [53] was used to identify small RNAs that were being cut in phase from longer precursors. We identified 316 siRNA ranging from 21 to 25 nt in length, of which 71% were 24 nt long (Table S6, Figure S1). The identified siRNAs received names in the following format: "Ll-siR(number)", (for example Ll-siR1) and displayed a highly differential expression pattern (Table S7). Some of the sequences showed organ-specific expression, for example, Ll-siR4, -13, -173 were present only in the pedicels of abscising flowers (FPAB), while Ll-siR308 showed an elevated expression in the pedicels (FPAB and FPNAB). On the other hand, Ll-siR246, -291 and -56 were present almost exclusively in the youngest flowers in the lower part of inflorescence (LF1) ( Table S7). Numerous reports and studies indicate the importance of phased siRNA not only in stress response mechanisms but also in growth regulation [54]. Therefore, we decided to investigate the role of siRNAs during yellow lupine inflorescence development. To achieve this, ShortStack (https://github.com/MikeAxtell/ShortStack/) [53] was used to identify small RNAs that were being cut in phase from longer precursors. We identified 316 siRNA ranging from 21 to 25 nt in length, of which 71% were 24 nt long (Table S6, Figure S1). The identified siRNAs received names in the following format: "Ll-siR(number)", (for example Ll-siR1) and displayed a highly differential expression pattern (Table S7). Some of the sequences showed organ-specific expression, for example, Ll-siR4, -13, -173 were present only in the pedicels of abscising flowers (FPAB), while Ll-siR308 showed an elevated expression in the pedicels (FPAB and FPNAB). On the other hand, Ll-siR246, -291 and -56 were present almost exclusively in the youngest flowers in the lower part of inflorescence (LF1) ( Table S7).

Analysis of the Expression Profile of the Identified sRNAs During Yellow Lupine Flower Development
To gain better insight into the dynamics expression of the identified sRNAs during floral development in yellow lupine, we established a wide scope comparison of the following growth stages of flowers from the upper (UF2 vs UF1, UF3 vs UF2 and UF4 vs UF3) and lower (LF2 vs LF1, LF3 vs LF2 and LF4 vs LF3) parts of the inflorescence (Table 4, Figure 5).
The analyses resulted in the identification of 30 differentially expressed miRNAs (DEmiRs) in the lower and 29 in the upper flowers (Table 4), as well as 14 and 7 DE siRNAs, respectively (Table  S8). Between UF2 and UF1, there was a change in the expression of 8 miRNAs, 2 sequences belonging to MIR359 and MIR166 families each, as well as one representative of each of the MIR159, MIR167, MIR396 families and novel Ll-miRn10. Ten DE miRNAs were identified in a comparison of UF3 vs UF2, of which only one Ll-miR258/miR166 was up-regulated. The remaining miRNAs were downregulated and consisted of 3 sequences belonging to the MIR390 and MIR396 families each, and single miRNAs from the MIR168, MIR408, and MIR396 families. A comparison of the UF4 vs UF3 libraries revealed 12 DEmiRs. The most numerous group were members of the MIR390 family, followed by 2 members of MIR167 and MIR319, and singular representatives of MIR398, MIR164, and MIR858, with one novel Ll-miRn11 (Table 4).
During the development of flowers from the lower part of the inflorescence, the miRNAs accumulation dynamics were different. The highest number of the identified DEmiRs was found comparing the youngest flowers (LF2 vs LF1), while, interestingly, a complete lack of DE miRNAs was found when comparing the oldest flowers: LF4 vs LF3 (Table 4). In our comparison of LF2 vs LF1, among the 18 DEmiRs, the most numerous group were novel miRNAs, followed by members of the MIR396 family. Between the LF2 and LF3 stages we confirmed that there was a change in the accumulation of 11 miRNAs, and this pertained to two members of the MIR166 and MIR399 families each, Ll-miRn1 and Ll-miRn22, which were followed by single representatives of the MIR390, MIR395, MIR858, MIR398, MIR408 families (Table 4).
In order to identify miRNAs the presence of which is either common or unique depending on the developmental stage of the upper and lower flowers in lupine, Venn diagrams were constructed ( Figure 6a and Table S9) using Venny 2.1 (https://bioinfogp.cnb.csic.es/tools/venny/) [55]. The results of these analyses revealed that approximately 70% of the identified miRNAs were common in all developmental stages of both the upper ( Figure 6a) and lower flowers ( Figure 5b). However, miRNAs unique to certain developmental stages were also found ( Figure 6 and Table S9).
In regard to siRNAs during flower development in yellow lupine, almost every differentially expressed siRNA was up-regulated. In the lower part of the inflorescence, similarly to miRNAs, there were no differences between the LF4 and LF3 stages. During the upper flower development, most DEsiRs were identified in a comparison of UF2 vs UF1, and the least (only one) when comparing UF3 vs UF2. One noteworthy observation was the presence of the same siRNAs in the comparisons of UF2 vs UF1 and LF2 vs LF1, namely Ll-siR281, -308. and -249, which suggests that an increase in their accumulation is important during phase 1 to phase 2 transition in the development of yellow lupine flowers, regardless of their position on the inflorescence. The complete dataset can be found in (Table S8).
Regarding the phased siRNAs, only 4 of them displayed differential expression, namely Ll-siR119 at stage 1 and Ll-siR224, -100 and -146 at stage 4. These results might suggest that, firstly, miRNAs display differential expression in each and every stage of flower development, regardless of flower position on the inflorescence, and secondly, that miRNAs seem to be much more impactful in comparison with phased siRNA in regards to yellow lupine flower differentiation.
Regarding the phased siRNAs, only 4 of them displayed differential expression, namely Ll-siR119 at stage 1 and Ll-siR224, -100 and -146 at stage 4. These results might suggest that, firstly, miRNAs display differential expression in each and every stage of flower development, regardless of flower position on the inflorescence, and secondly, that miRNAs seem to be much more impactful in comparison with phased siRNA in regards to yellow lupine flower differentiation. Analyses of the Venn diagrams we created (Figure 6c), displaying the presence profiles for the library miRNAs, revealed that in each comparison between the upper and lower flowers (UF1 vs LF1, etc.) around 80% of the identified sequences were common for both the upper and lower flowers ( Figure 6c). However, in each comparison, we were able to identify miRNAs unique to each stage of the development and each flower position. For example, 20 miRNAs were exclusively present in LF1, while 12 miRNAs were unique to UF1. The detailed information on these comparisons can be found in Table S9.
Based on the data received, we suggest that differences in miRNA expression between lower and upper flowers may be related to the fate of these organs (pod formation/flower abscission). To further confirm this function, we performed an experiment in which flowers were removed from the lower whorls, leaving only flower buds from the last, top whorl ( Figure S2). Removing the lower flowers causes maintenance of flowers on the plant and their development into pods, unlike flowers from this whorl in control plants. Thus, their fate seems to be associated with the location in the inflorescence changes. Then, the expression of selected lupine DEmiRs and their target genes were compared during the development of upper flowers after removal of the lower flowers (UFR) in the development stages of S1-S4, with control upper (UF S1-S4) and lower (LF S1-S4) flowers, respectively ( Figure S3). The obtained results show that the removal of lower flowers caused a change in the levels of chosen sRNAs in upper flowers and it similar in this respect to flowers from the lower part of raceme. This indicates a link between these genes and the fate of the flowers.

Comparison of Differentially Expressed sRNAs Between Flower Pedicels with Active And Inactive Abscission Zones
To identify sRNAs possibly involved in yellow lupine flower abscission, mi-and siRNA expression patterns for flower pedicels with an active abscission zone (AZ) (FPAB) and inactive AZ (FPNAB) were compared. As a result, 34 DE miRNAs (including 5 novel ones) ( Table 6) and 20 DE phased siRNAs (Table S8) were identified. 14 miRNAs and 9 siRNAs were up-regulated, while the rest remained down-regulated in FPNAB. Among the up-regulated miRNAs, the most numerous family was MIR167 (5 members), followed by MIR398 (3 members). Among the down-regulated miRNAs, the most abundant were MIR390, MIR396 and MIR395 families with 3 members each (Table 6, Figure 5b). With regard to siRNAs, the most up-regulated in FPANB were Ll-siR173, -4 and -13, and the most down-regulated was Ll-siR208 (Table S8). An analysis of the Venn diagrams based on the presence of the identified miRNAs revealed that approx. 80% of the miRNAs were present in both abscising and non-abscising flower pedicles (Figure 6d). However, 23 miRNAs remained unique to FPAB and 17 to FPNAB (Figure 6d, Table S9).

Validation of the Identified sRNAs in RNA-seq
Stem-loop RT-qPCR technique [56,57] was employed in order to validate the data generated using deep sequencing technology and to confirm the expression patterns of the identified sRNA. Eight identified sRNAs (six conserved miRNAs, one novel miRNA, and one siRNA) were used for this task (Table S10). The qPCR results were similar to sRNA-seq data (Figure 7). For example, in the RT-qPCR analysis, the Ll-siR254 expression increased as the flower developed, showing a positive correlation with the deep sequencing results. Ll-siR249 was preferentially accumulated in yellow lupine pedicels, both in qPCR and RNA-seq. The results of the expression analysis of these sRNAs supported the validity of our sRNA-Seq. An analysis of the Venn diagrams based on the presence of the identified miRNAs revealed that approx. 80% of the miRNAs were present in both abscising and non-abscising flower pedicles ( Figure  6d). However, 23 miRNAs remained unique to FPAB and 17 to FPNAB (Figure 6d, Table S9).

Validation of the Identified sRNAs in RNA-seq
Stem-loop RT-qPCR technique [56,57] was employed in order to validate the data generated using deep sequencing technology and to confirm the expression patterns of the identified sRNA. Eight identified sRNAs (six conserved miRNAs, one novel miRNA, and one siRNA) were used for this task (Table S10). The qPCR results were similar to sRNA-seq data (Figure 7). For example, in the RT-qPCR analysis, the Ll-siR254 expression increased as the flower developed, showing a positive correlation with the deep sequencing results. Ll-siR249 was preferentially accumulated in yellow lupine pedicels, both in qPCR and RNA-seq. The results of the expression analysis of these sRNAs supported the validity of our sRNA-Seq.

Identification of sRNA Target Genes using Degradome and psRNATarget Analysis
In order to estimate accurately the biological function and impact of certain miRNAs, their target genes needed to be identified. To achieve this, we constructed degradome libraries from pooled samples of stage 3 upper and lower parts of the inflorescence. Through total degradome library sequencing, 19,353,278 raw reads were obtained (Table S11). After quality filtering, the degradome data were aligned to the reference transcriptome with CleaveLand 4 [58] to find sliced miRNA and siRNA targets. After processing and analysis, a total of 14,077 targets were identified, and after filtering with a p-value < 0.05, 538 targets emerged (501 targets for 178 known miRNAs and 37 targets for 13 novel miRNAs) (Table S12). For the phased siRNAs, 3,340 targets were initially identified, and after similar filtering, their number dropped to 89 targets for 46 siRNAs (Table S13). Exemplary target t-plots and sequences of the miRNAs and target mRNAs are shown in Figure 8.

Identification of sRNA Target Genes using Degradome and psRNATarget Analysis
In order to estimate accurately the biological function and impact of certain miRNAs, their target genes needed to be identified. To achieve this, we constructed degradome libraries from pooled samples of stage 3 upper and lower parts of the inflorescence. Through total degradome library sequencing, 19,353,278 raw reads were obtained (Table S11). After quality filtering, the degradome data were aligned to the reference transcriptome with CleaveLand 4 [58] to find sliced miRNA and siRNA targets. After processing and analysis, a total of 14,077 targets were identified, and after filtering with a p-value < 0.05, 538 targets emerged (501 targets for 178 known miRNAs and 37 targets for 13 novel miRNAs) (Table S12). For the phased siRNAs, 3,340 targets were initially identified, and after similar filtering, their number dropped to 89 targets for 46 siRNAs (Table S13). Exemplary target t-plots and sequences of the miRNAs and target mRNAs are shown in Figure 8.  As expected, many of the targets for evolutionarily conserved miRNAs were compliant with literature data. For example, Ll-miR329/miR160-5p targeted ARF16 and ARF18, the Ll-miR415/miR171b targeted SCL6, Ll-miR341/miR319q targeted TCP2, Ll-miR224/miR393a-5p targeted TIR1, etc. (Table S12).
A comparison of the expression of four exemplary miRNAs and their target genes (Figure 8) confirmed the reverse-correlation in the accumulation of miRNAs and an abundance of mRNA target genes, especially in flower pedicels. In the flowers, this correlation was not so obvious, presumably As expected, many of the targets for evolutionarily conserved miRNAs were compliant with literature data. For example, Ll-miR329/miR160-5p targeted ARF16 and ARF18, the Ll-miR415/miR171b targeted SCL6, Ll-miR341/miR319q targeted TCP2, Ll-miR224/miR393a-5p targeted TIR1, etc. (Table S12).
A comparison of the expression of four exemplary miRNAs and their target genes (Figure 8) confirmed the reverse-correlation in the accumulation of miRNAs and an abundance of mRNA target genes, especially in flower pedicels. In the flowers, this correlation was not so obvious, presumably because of the organ's more complex nature (with its various elements, such as the stamen and the pistil), where regulation could be tissue specific.
In the case of some of the identified mi-and siRNAs, we were unable to determine the targets with a degradome analysis, which might have been caused by the lack of a sufficient amount of cleavage products ensuing from using only stage 3 flowers to construct the library. In order to find the putative missing target genes, the psRNATarget tool [59] was employed, which rendered plausible target genes through a comparison of the sRNAs with the reference transcriptome containing data from all of the samples. Using this method, we managed to establish putative target genes for most of the mi-and siRNAs, obtaining 66,102 miRNA and 32,725 siRNA targeted transcripts. A full list of the targets identified using the psRNATarget or degradome analysis for DE miRNAs, siRNAs, and novel miRNAs is contained in Table 3, Table 4, Table 5, Table 6 and Table S13. Targets for all of them are shown in Tables S14, S15, S16, S17 and S18.

Function of the miRNAs Potential Targets
Gene Ontology (GO) analysis was performed in order to investigate the functions of the miRNAs targets identified in yellow lupine flowers. Among the 27,547 targets of known and novel miRNAs identified with psRNATarget 26,230 targets exhibited GO terms (Table S17). 23,092 genes were categorized into 'Cellular component', 23,501 into 'Molecular function', and 22,939 into 'Biological process'. Figure 9 shows target gene percentages for each GO category. The largest number of targets classified as 'Cellular component' was attributable to 'cell', 'cell part' and 'organelle'. The majority of targets of the 'Molecular function' category were classified as 'binding' and 'catalytic activity'. Within the 'Biological process', most of the targets were categorized as 'cellular' and 'metabolic process' (Figure 9, Table S17). Within the 'Flower development' category, the targets of 37 miRNAs fit within GO terms related to phytohormones ( Figure S4a), and the targets of 69 miRNAs were placed into the category of GO terms related to the development of flower parts ( Figure S4b). because of the organ's more complex nature (with its various elements, such as the stamen and the pistil), where regulation could be tissue specific.
In the case of some of the identified mi-and siRNAs, we were unable to determine the targets with a degradome analysis, which might have been caused by the lack of a sufficient amount of cleavage products ensuing from using only stage 3 flowers to construct the library. In order to find the putative missing target genes, the psRNATarget tool [59] was employed, which rendered plausible target genes through a comparison of the sRNAs with the reference transcriptome containing data from all of the samples. Using this method, we managed to establish putative target genes for most of the mi-and siRNAs, obtaining 66,102 miRNA and 32,725 siRNA targeted transcripts. A full list of the targets identified using the psRNATarget or degradome analysis for DE miRNAs, siRNAs, and novel miRNAs is contained in Tables 3, 4, 5, 6 and S13. Targets for all of them are shown in Tables S14, S15, S16, S17 and S18.

Function of the miRNAs Potential Targets
Gene Ontology (GO) analysis was performed in order to investigate the functions of the miRNAs targets identified in yellow lupine flowers. Among the 27,547 targets of known and novel miRNAs identified with psRNATarget 26,230 targets exhibited GO terms (Table S17). 23,092 genes were categorized into 'Cellular component', 23,501 into 'Molecular function', and 22,939 into 'Biological process'. Figure 9 shows target gene percentages for each GO category. The largest number of targets classified as 'Cellular component' was attributable to 'cell', 'cell part' and 'organelle'. The majority of targets of the 'Molecular function' category were classified as 'binding' and 'catalytic activity'. Within the 'Biological process', most of the targets were categorized as 'cellular' and 'metabolic process' (Figure 9, Table S17). Within the 'Flower development' category, the targets of 37 miRNAs fit within GO terms related to phytohormones ( Figure S4a), and the targets of 69 miRNAs were placed into the category of GO terms related to the development of flower parts ( Figure S4b). Our miRNAs targets analysis against the Kyoto Encyclopedia of Genes and Genomes (KEGG) revealed that most of the sequences in the main KEGG categories belonged to Metabolism (15,856), followed by Genetic information processing (5,267), Environmental information processing (1,517), Cellular processes (1,326) and Organismal Systems (800) ( Figure S5). A full list of KEGG pathways and numbers of assigned sequences is shown in (Figures S6 and S7). One of the most represented Our miRNAs targets analysis against the Kyoto Encyclopedia of Genes and Genomes (KEGG) revealed that most of the sequences in the main KEGG categories belonged to Metabolism (15,856), followed by Genetic information processing (5,267), Environmental information processing (1,517), Cellular processes (1,326) and Organismal Systems (800) ( Figure S5). A full list of KEGG pathways and numbers of assigned sequences is shown in (Figures S6 and S7). One of the most represented sub-categories was Signal transduction (1,484), with over 700 putative targets in the Plant Hormone Signal Transduction pathway, where almost every sequence was frequently targeted with multiple miRNAs (Figure 10). The second most notable pathway was mitogen-activated protein kinase (MAPK) signaling, which was associated with different abiotic and biotic stress factors, with 350 putative targets distributed across every described stress response ( Figure S8). A complete dataset on the KEGG analysis can be seen in (Table S16). sub-categories was Signal transduction (1,484), with over 700 putative targets in the Plant Hormone Signal Transduction pathway, where almost every sequence was frequently targeted with multiple miRNAs (Figure 10). The second most notable pathway was mitogen-activated protein kinase (MAPK) signaling, which was associated with different abiotic and biotic stress factors, with 350 putative targets distributed across every described stress response ( Figure S8). A complete dataset on the KEGG analysis can be seen in (Table S16).

Discussion
Yellow lupine has great potential to become one of the leading legumes in Europe in both animal and human nutrition. Reduction the economic drawbacks resulting from excessive flower abscission would be the most convincing argument for lupine cultivation. However, this can only be achieved if we gain a deeper understanding of the plant's biology and insight into the molecular basis for the development and maintenance of lupine flowers. Therefore, we believe that the pathways controlling these processes deserve intensive research focus. Our previous analyses of yellow lupine transcriptomes resulted in the identification of transcripts of many genes involved in flower and pod abscission and suggested sRNA involvement in this process [17]. Notably, our observations of L. lupinus floral development indicate that their fate (abscission or pod formation) is determined prior to AZ activation. Therefore, we decided to perform comparative analyses between sRNAs from flowers developing on the upper and lower parts of the raceme. Identifying the miRNAs and their target genes involved in the above-mentioned processes will further our knowledge of the biology of not only lupines but plants in general since the role played by sRNA in organ abscission is still obscure.
Our sRNA-seq analyses shed more light on the molecular mechanisms that control flower development of L. luteus and confirmed the involvement of known miRNAs, such as miR159, miR167 or miR172, in this process [60], but we have also explored the roles of sRNAs in flower abscission and identified species-specific miRNAs.

Known sRNAs and Their Target Genes Are Involved in Regulating Flower Development in Yellow Lupine
Among the known and conserved miRNAs a number of miRNAs commonly associated with flower morphogenesis and development, belonging to, inter alia, the MIR156/157, MIR159, MIR165/166, MIR167 and MIR172 families [10] were spotted.
Studies have shown that miR156 is necessary for maintaining anther fertility in Arabidopsis, by orchestrating the development of primary tapetum cells and primary sporogenous cells [61]. In A. thaliana, SPL13B expression is strictly limited by miR156 to anther tapetum in young buds, while SPL2 is weakly expressed in parietal and sporogenous cells and the surrounding cell layers in young flowers [61], where it is targeted by miR156 to regulate pollen maturation [62]. MiR159 was shown to target the conserved GAMYB-like genes that are a part of the GA signaling pathway [63,64]. In A. thaliana miR159 regulates the morphogenesis of the stamen, and male fertility [65]. Two transcription factors involved in pistil and stamen development in various plant species, ARF6 and ARF8, contain the target site for miR167 [66][67][68]. For Arabidopsis, it has been proven that both these genes are involved in stamen filament elongation, anther dehiscence, stamen maturation and anthesis [69]. In tomato, a reduction in the accumulation of the miR167-targeted ARF6 and ARF8 leads to the lack of trichomes on the style surface, failed pollen germination and, consequently, sterility [11]. Recent research into multiple plant species has shown that miR172 targets genes belonging to the APETALA2 (AP2, TOE1, TOE2, TOE3) family. MiR172 is part of the photoperiodic flower induction pathway and is associated with the functioning of the ABCDE model of floral development [70]. Overexpression of MIR172 causes formation of a phenotype characterized by the absence of perianth, transformation of sepals into pistils and early flowering [70].
Our study showed the presence of at least one member of all these families in flowers (Figure 3, Table S5), which indicated that in lupine how crucial the families were for generative development in lupine, as well. MIR156 and MIR159 are the most numerous families in L. luteus, which suggests they play fundamental roles in its flower development processes.
The differentially expressed miRNAs identified in yellow lupine flowers were clustered by the dynamics of their expression ( Figure 5). The first cluster comprised miRNAs, the accumulation of which increased as the flowers developed, and contained miRNAs belonging to the MIR166, MIR167, MIR319, MIR390, and MIR395 families. The first of these families include Ll-miR177, which guides the cleavage of RADIALIS, a transcription factor from the MYB family that controls the asymmetric flower shape in Antirrhinum majus [71,72], as well as Ll-miR258 and Ll-miR265, which probably target the Homeobox-leucine zipper protein ATHB-15. In Arabidopsis, both miR165 and miR166 target the same HD-ZIP III genes: ATHB15, ATHB8, REVOLUTA (REV), PHABULOSA (PHB), and PHAVOLUTA (PHV) to regulate gynoecium and microspore development [28,73]. In lupine the MIR167 family members that accumulate in larger quantities during flower development are Ll-miR280, Ll-miR281, and Ll-miR285, which probably target ARF6 and ARF8. Ll-miR445 and Ll-miR130 are members of the MIR319 family, while their putative target genes are TCP4 and MYB33, respectively. In Arabidopsis, the miR319a/TCP4 regulatory module is necessary for petal growth and development. Moreover, the overexpression of MIR319 reduces male fertility, and this defect is hypothesized to be caused by the cross-regulation of MYB33 and MYB65 by miR319 and miR159. As the miR319 target site within the MYB33 and MYB65 transcripts exhibit a lower match with miRNA than the miR159 target site, the latter is more efficient at regulating these genes and miR319 is their secondary regulator [74]. This regulatory network is even more complex. In A. thaliana, cooperation of three miRNAs and their target genes, namely miR159/MYB, miR167/ARF6/ARF8, and miR319/TCP4, is a prerequisite for proper sepal, petal and anther development, and maturation. miR159 and miR319 influence the expression of MIR167 genes, which in turn affect each other. These miRNAs orchestrate plant development by regulating the activity of the phytohormones GA, JA, and auxin [75]. Increased accumulation of miR167 and miR319 in the late stages of yellow lupine flower development could also be associated with regulating the growth and development of petals and anthers. Another miRNA showing a similar expression profile is Ll-miR9/miR390-5p. In lupine, it targets the TAS3 transcript, which in turn is a source of tasiR-ARF, a negative regulator of ARF2, ARF3 and ARF4 activity. This regulatory cascade plays a vivid role in development of many plant species [76]. The expression level of miR390 derived from MIR390b reflects auxin concentration in organs, while the repression of ARF2, ARF3, and ARF4 by tasiR-ARF are important for lateral organ development [18,77], and flower formation [78]. Ll-miR118 and Ll-miR119, which target ATP sulfurylase (ATPS) according to our degradome data, belong to the MIR395 family. In Arabidopsis, miR395 targets two gene families, ATP sulfurylases and sulfate transporter 2:1 (SULTR2:1), which are elements of the sulfate metabolism pathway [79]. ATPS regulates glutathione synthesis and is an essential enzyme in the sulfur-assimilatory pathway [80]. In cotton, the miR395-APS1 module is engaged in drought and salt stress response [81]. Sulfate is the main source of sulfur and is taken up by roots, transported throughout the plant and used for assimilation. Sulfate limitation forces a significant up-regulation of miR395 expression [82]. Presumably, during yellow lupine flower development, the demand for sulfur increases and the plant activates mechanisms for its efficient uptake.
Within the cluster of miRNAs, the expression of which decreased as the flowers developed, there were homologues of miR390-3p, miR858, miR396-3p, miR168, miR408-3p and miR398 ( Figure 5). Ll-miR99, Ll-miR100, and Ll-miR102 are identical to miR390-3p (the so-called passenger strand, former star strand). However, their expression showed an opposite trend to that of miR390-5p. The differential expression and functioning of passenger miRNAs have already been described. The research carried out by Xie and Zhang in 2015 on cotton showed that the formation of some miRNA*s, such as miR172* and miR390*, was associated with the phases of the plant's growth [83]. Therefore, miRNA*s can be specifically expressed in various tissues to maintain the steady state of the organism. Our degradome analysis for yellow lupine showed that Ll-miR9/miR390-5p was able to guide the cleavage of the TAS3 transcript. There is no certainty as to the status of its passenger strand, which suggests its locally limited activity or its involvement in regulation of other targets and further research is required to identify its accumulation and function in the organs concerned. Another miRNA from the cluster is Ll-miR155/miR396-3p (passenger strand), which guides cleavage of JMJ25 demetyhylase mRNA (confirmed in degradomes), involved in preserving the active chromatin state [84]. ECERIFERUM1 (CER1), the target gene in lupine for another two homologues of miR396-3p, Ll-miR199 and Ll-miR200, is a homologue encoding an enzyme involved in alkane biosynthesis, and in cucumber is engaged both in wax synthesis and ensuring pollen viability [85]. This cluster also included a miRNA that negatively regulates elements involved in miRNA and ta-siRNA functioning, namely Ll-miR247/miR168 targeting AGO1 mRNA [86]. Another miRNA clustered here was the highly conserved Ll-miR60/miR408-3p, which guides the processing of the mRNA of the copper-binding Basic Blue protein homologue (plantacyanin, PC). In Arabidopsis, PC plays a role in fertility, exhibiting the highest expression in the inflorescence, especially in the transmitting tract. [87]. Transgenic Arabidopsis plants over-expressing MIR408 displayed altered morphology, including significantly enlarged organs, resulting in enhanced biomass and seed yield. Plant enlargement was shown to be primarily caused by cell expansion rather than cell proliferation, and in transgenic plants it was correlated with stronger accumulation of the myosin-encoding transcript and gibberellic acid [88]. It seems that high expression levels of miRNAs grouped in the cluster are correlated with intensive growth and differentiation of young floral tissues.
Among the miRNAs identified in yellow lupine several that seemed to be crucial in particular stages of the plant's development were spotted ( Figure 4, Table 4, Table S7). For example, the largest quantities of miR159 (Ll-miR452 and Ll-miR454) were accumulated in stages 2 and 3 of the plant's development. According to degradome data they targeted GGP-5 (GAMMA-GLUTAMYL PEPTIDASE 5) of an undefined function in plants, and an evolutionarily conserved target for GAMYB, respectively. As already mentioned, this could be associated with miRNA family cooperating with miR167 and miR319 in regulating L. luteus anther maturation. The accumulation of Ll-miR251/miR5168-3p, Ll-miR92/miR1861b, Ll-miR229/miR369-5p, and Ll-miR311/miR5794 increased in stage 2 upper and lower flowers, while -interestingly -in the later stages these miRNAs were only present in lower flowers. According to degradome analysis, Ll-miR251/miR5168 guides cleavage of the mRNAs of the genes encoding the Homeobox-leucine zipper protein ATHB-14 and the chaperone protein dnaJ 13. The miR5168 sequence displays a great similarity to that of miR166, thanks to which they may perhaps share the same target gene ATHB-14, the putative transcription factor engaged in the adaxial-abaxial polarity determination in the ovule primordium in A. thaliana [89]. As confirmed by yellow lupine degradome sequencing, Ll-miR229/miR396-5p targets GROWTH-REGULATING FACTOR 5 (GRF5) and GRF4 transcripts. In Arabidopsis, GRF5 is expressed in anthers at early stages of flower development and in gynoecia throughout the whole flower development, and transcripts of GRF4 accumulate later in sepals, tapetum, and endocarpic tissues of ovary valves [90]. Transgenic rice with Os-miR396 overexpression and GRF6 knock-down suffers from open husks and sterile seeds [91]. GRF6 cooperates with GRF10 to transactivate the JMJC gene 706 (OsJMJ706) and CRINKLY4 RECEPTOR-LIKE KINASE (OsCR4) responding to GA, which is a prerequisite for the flower to successfully develop into a normal seed [91]. An increased share of miRNAs involved in cell division, namely miR396, miR319, and miR164, in NGS analyses was also observed in early grain development in wheat [92].The presence of these miRNAs in yellow lupine flowers suggests that their regulation of cell proliferation also plays an important role in development of generative organs.

Involvement of New miRNAs in L. luteus Flower Development
Using ShortStack [53] software we predicted 28 candidates for new miRNAs (Table 3). Interestingly, many of these novel miRNAs showed similarity to precursor miRNAs from miRBase, which leads to the conclusion that they might be new members of the already known families, for example MIR167 (Ll-miRn12 and Ll-miRn27), MIR172 (Ll-miRn4), MIR393 (Ll-miRn19) or MIR169 (Ll-miRn3, Ll-miRn11, and Ll-miRn15) ( Table S6).The other 13 had no homologues among known miRNAs and were recognized as lupine-specific miRNAs. Some of the new miRNAs displayed differential expression during L. luteus flower development. Ll-miRn3, which shows similarity to pre-miR169, displayed differential expression in UF1 vs LF1 and LF2 vs LF1 library comparisons, wherein it is the most accumulated in LF1, and in flower pedicels (up-regulated in FPNAB). According to degradome data, this miRNA targets SCARECROW2 (SCR2) homologue, a putative activator of the calcium-dependent activation of RBOHF that enhances reactive oxygen species (ROS) production and may be involved in cold stress response [93]. In rice SCR2 expression is relatively high in flower buds and flowers, and after flowering rises in the leaves and roots [94]. In yellow lupine, this gene may be involved in intense cell divisions during early flower development and is down-regulated in the pedicels with an active AZ to stop its growth. Another frequently encountered novel DEmiR was Ll-miRn22, which shows sequence similarity to pre-miR1507, is up-regulated in LF3 vs LF2 and LF2 vs LF1 library comparisons, and its expression escalates with flower development in the bottom whorl. The MiR1507 family is annotated as legume-specific [95]. Through analyses of our degradome data we have not found its target gene, and the psRNATarget hit was the putative disease resistance RPP13-like protein 1. Unfortunately, this protein has been poorly described, therefore it is difficult to determine its function in yellow lupine flowers. Noteworthily, the target genes of Ll-miRn1 and Ll-miRn30 identified through degradome sequencing are SGS3 and DCL2, respectively, and the miRNAs are up-regulated in LF3 vs LF2 comparisons and down-regulated in UF1 vs LF1 comparisons, respectively. SGS3and DCL2-encoded proteins are involved in sRNA biogenesis [96]. Importantly, novel miRNA identified in soybean Soy_25 displays high sequence similarity to Ll-miRn1 and also targets SGS3, which indicates that this regulatory feedback loop for sRNA biogenesis is common for Fabaceae [97]. These results indicate that L. luteus miRNAs play a regulatory role in siRNA biogenesis in early flower development.

miRNA Accumulation Varies in Lower and Upper Flowers in Different Stages of Development
One of our goals was to identify the sRNAs engaged in yellow lupine flower development, with a particular emphasis on the differences between flowers from lower and upper parts of the inflorescence, in order to gain an insight into how early the flower fate is determined.
In our study, we spotted differences in miRNA accumulation patterns as early as the first stage of flower development.
Flowers collected from the lower whorls displayed higher accumulations of sequences corresponding to miR5490, miR5794, miR1861, miR396-5p, miR395, miR166, and miR159-3p (Table 5). miR1861 and miR396 were recognized as positive cell proliferation and development regulators [98][99][100]. In rice, for example, miR1861 exhibited differential expression during grain filling [101], and its expression was higher in superior grains in comparison to inferior ones [102]. This is consistent with our hypothesis, that a higher occurrence of miR1861 and miR396 in lower flowers may be an indication of the plant investing more supplies in this part of the inflorescence.
From the second stage until the end of their development, upper flowers accumulated more miRNAs corresponding to miR319, miR394, miR160, and miR393 ( Figure 4, Table 5). MiR393 regulates the accumulation of transcripts encoding auxin receptors belonging to the TAAR family. Changes in receptor abundance affect the sensitivity of the given tissue to auxin and this is how this molecule influences plant development [102]. In A. thaliana, miR160 directly controls three ARF genes, namely: ARF10, ARF16 and ARF17 [103]. In tomato, sly-miR160 is abundant in ovaries, and changes in its expression affect plant fertility [12]. Down-regulation of sly-miR160 caused improper ovary patterning and thinning of the placenta already prior to anthesis [12]. In view of these facts, higher expression of miR160 in lupine upper flowers in their development means that a slightly different organization of the gynoecia may be one of the crucial determinants of flower fate. Additionally, the elevated expression levels of miR160 and miR393 in upper flowers of lupine suggest a reduction in the abundance of the transcripts of their target genes encoding auxin receptors and auxin response factors. This, in turn, may have led to a reduction in auxin sensitivity. Decreasing the number of transcription factors belonging to the TCP family (targeted by miR319), probably caused different cell proliferation profiles in flowers collected from the upper whorls.
Additional expression studies of selected miRNA (Ll-miR281/miR167, Ll-miR224/miR393, Ll-miR333/miR160, Ll-miR329/miR160) carried out in the upper flowers of yellow lupin developing after removal of the lower ones (UFR) (Figure S2), and consequently with a changed, when compared to the original, fate, provide additional confirmation of the results obtained from RNA-seq analysis ( Figure S3).

sRNAs Are Involved in Flower Abscission in L. luteus
Little is known about sRNA engagement in flower abscission. Research on the involvement of miRNAs in this process has been already carried out in cotton [104], tomato [12,105], and sugarcane [106].
For a genome-wide investigation of miRNAs involved in the formation of the abscission layer in cotton, two sRNA libraries were constructed using the abscission zones (AZ) of cotton pedicels treated with ethephon or water. Among the 460 identified miRNAs, only gra-MIR530b and seven novels showed differential expression in abscission tissues [104], and these miRNAs have no homologues in our dataset.
Besides ovary patterning in tomato, sly-miR160 regulates other two auxin-mediated developmental processes: floral organ abscission and lateral organ lamina outgrowth [12]. In that study, down-regulation of sly-miR160 and the resulting higher expression of its target genes, transcriptional repressors of auxin response ARF10 and ARF17, also resulted in the narrowing of leaves, sepals and petals and an impeded shedding of the perianth after successful pollination [12]. This was consistent with the higher accumulation of Ll-miR329/miR160-5p, Ll-miR332/miR160-5p, and Ll-miR333/miR160-5p in upper flowers designated to fall off in yellow lupine. As these miRNAs showed no differential expression in flower pedicels, it probably does not play a role in the executory module of abscission itself but is rather a part of a mechanism that determines flower fate.
Another research on tomato using sRNA and degradome sequencing libraries explored the roles of sRNAs in AZ formation in the early and late stages of the process additionally accelerated or not by ethylene or control treatment [107]. The study showed that in tomato pedicels, the accumulation levels of, inter alia, miR156, miR166, miR167, miR169, miR171, and miR172 rose in late stages of abscission, while the abundance of miR160, miR396 and miR477 dropped [107]. Although it is difficult to compare ethylene-treated tomato pedicel results to our data, it is worth noting that in the corresponding FPAB vs. FPNAB comparison in our study, the accumulation of some miRNAs was similar: miR396 level was lower, and the levels of miRNAs annotated as miR167 and miR166 were higher in FPAB (Table 6).
It has been proven for sugarcane that among others both mature (5p) and passenger (3p) miRNAs from MIR167 family were up-regulated in 'leaf abscission sugarcane plants' comparing to 'leaf packaging sugarcane plants' (which corresponds to the FPAB vs. FPNAB comparison in our study) [106]. In our study, both mature and passenger members of the MIR167 family were leaders among DEmiRs, too, ( Table 6) pointing to their crucial role in both vegetative and generative organ abscission. Significantly, this applies to evolutionarily distant taxa: both monocots and dicots.
In our paper, among the up-regulated miRNAs, the most numerous family besides already mentioned MIR167 was MIR398 with 3 members being among top-regulated ones. Among the down-regulated miRNAs, the members of MIR390, MIR396 and MIR395 families were most abundant. It was shown for other plant species, that these miRNAs are engaged in the regulation of auxin signal transduction pathway (miR167 and miR390 [108]), regulation of cell division (miR396 [100]) and stress response (miR395 [81,82]).
It is worth noting, that in comparisons of Lupinus pedicel libraries there are novel miRNAs: three are down-regulated in FPAB and one is up-regulated. Furthermore, Ll-miRn3 is up-regulated in both, young flowers designated to be maintained on the plant (LF1) and pedicels with inactive AZ (Table 6), which may indicate its role in preventing flower abscission. In the future, it is worth examining the role of its target gene, which encodes a protein that does not resemble any known protein.
With regard to siRNAs, the most up-regulated ones in FPNAB were: Ll-siR173, Ll-siR4 and Ll-siR13, and the most down-regulated one were Ll-siR208. Unfortunately, the lack of literature data on their targets makes it impossible for the specifics of their function in the studied process. However, it is worth mentioning, that in pedicels high levels of accumulation are displayed by siR249/tasiR-ARF and siR308/tasiR-ARF, which target transcripts encoding ARF2, ARF3 (confirmed in degradomes). These results strongly suggest the involvement of siRNAs in the functioning of lupine pedicels.

Possible miRNA-dependent Regulatory Pathways That Participate in Development and Abscission of Yellow Lupine Flowers
Recent studies have shown that sRNA activity is associated with the hormonal regulation of plant development through influencing the spatio-temporal localization of the hormone response pathway [109].
The auxin signal transduction pathway mainly consists of three elements. Auxin is perceived by members of the TAAR family. There are AUX/IAA repressor proteins and ARF transcription factors downstream of these receptors [110][111][112]. The expression of TAAR receptors is regulated by miR393 and secondary ta-siRNA derived from their own transcripts [20]. miR167 and miR160 affect the ARF6, ARF8 [67] ARF10, ARF16 and ARF17 [113] transcript accumulation, respectively. It has been proven that the expression of ARF2, together with ARF3 and ARF4, is regulated by the ta-siRNA/miR390 module [114]. In the two-hit model, ta-siRNA-containing the TAS transcript is recognized by two miR390 molecules, one of which guides its cleavage, and the other, in a complex with AGO7, serves as a primer for complementary strand synthesis, with its subsequent processing ultimately resulting in ARF-targeting siRNA biogenesis [115].
In our study, among the differentially expressed sRNAa in flowers and flower pedicels, there were members of the MIR167, MIR160, MIR393 and MIR390 families, as well as phased siRNAs targeting ARF2, ARF3, and ARF4. This fact suggests a vivid role of auxin-related sRNAs in flower development and abscission in L. luteus and confirms our previously published results of transcriptome-wide analyses, where we observed differences in expression levels of genes encoding several elements of the auxin signal transduction pathway [17]. The relatively high number of members of the MIR167 family showing differential expression in the studied variants indicates that miR167 is one of the key regulators of flower development and abscission in yellow lupine.
Lupinus LlARF2, LlARF3, and LlARF4 transcripts are possibly down-regulated in the processing that is guided by Ll-siR249 and Ll-siR308 (Table 4), which are identical to tasiR-ARFs in many plant species according to the tasiRNAdb database [116]. These tasiR-ARFs probably originate from TAS3 transcript (TRINITY_DN55534_c4_g1) containing two binding sites for miR390 ( Figure S9a). Ll-miR9/miR390and surprisingly also Ll-siR240, guide the cleavage of another TAS3 mRNA (TRINITY_DN54998_c6_g5_i2) ( Figure S10) which contains only one target site for miR390 ( Figure S9b). This is the first report on TAS3 processing regulated by siRNA. The target site for Ll-siR240 is shifted by 10 nucleotides relative to the target site for Ll-miR9/miR390 ( Figure S10). The expression of Ll-siR249, Ll-siR308, and Ll-miR9 showed a similar profile, as it rose during flower development and was the highest in the pedicels (Figure 7). Ll-siR240 accumulated proportionally to TAS3 with only one target site for miR390, which means that it was least expressed in the pedicels, while in flowers its expression increased with time (Table S18). The identified target transcripts belonging to the ARF2, ARF3, and ARF4 gene families showed differential expression but with no clear trend (Table S18). This may indicate that these siRNAs act locally, repressing only a pool of transcripts expressed in a given tissue, while in other flower parts activity of these genes is regulated in other ways. The presence of all the elements of the miR390/TAS3/tasiR-ARF module among the DE sRNAs in yellow lupine suggests that alterations in its functioning have a great impact on L. luteus flower development. The additional element in the form of siRNA that processes TAS3 mRNA seems to be a new species-specific adjuster of this regulation module.
We have also performed GO enrichment analysis of the target genes for sRNAs identified in flowers of yellow lupine (Figure 9, Figure S4a,b, Table S10). What is most interesting is that quite a considerable number of target genes fell within the 'response to stimulus' and 'signaling' categories, which means that miRNAs modulated the way the plant adapted to environmental stimuli ( Figure 9). An in-depth analysis of GO terms concerning plant hormones ( Figure S4a) showed that most of the miRNAs identified in yellow lupine modulated more than one hormone signaling pathway. For example, Ll-miR181 belonging to the MIR166 family modulated processes associated with four hormones, namely auxin, gibberellin, jasmonic acid, and salicylic acid, by targeting not only transcription factor AS1, a central cell division regulator [117], but also Cullin-3A, an element of the ubiquitination complex [118]. Another two members of this family, Ll-miR173 and Ll-miR177, targeted the same gene, 26S PROTEASOME NON-ATPASE REGULATORY SUBUNIT 8 HOMOLOG A (RPN12A), involved in the ATP-dependent degradation of ubiquitinated proteins during auxin and cytokinin response [119]. Our GO analysis for yellow lupine flowers additionally showed that miRNAs were responsible for guiding the processing of genes simultaneously involved in multiple processes associated with flower development ( Figure S4b). For example, in many plants AP2 is involved in the specification of floral organ identity [120], as well as ovule [121] and seed development [122,123], and in our study, it was targeted by ten lupine miRNAs. On the other hand, seven of these miRNAs additionally targeted a homologue of negative flower development regulator, LIKE HETEROCHROMATIN PROTEIN 1 (LHP1) [124]. This highly degenerated and ambiguous model of gene regulation by lupine miRNAs shows that in this plant the adjustment of key biological processes related to fertility is a complex network of interconnected factors.
We have also conducted KEGG functional analysis of the putative targets identified for miRNAs in lupine which indicated their engagement in regulating a number of metabolic pathways-especially 'carbohydrate metabolism' and 'nucleotide metabolism' ( Figure S5). 'Carbohydrate metabolism' was also one of the most enriched KEGG pathways in our previous L. luteus transcriptome analysis [17], and its activation may be an indication of cell walls being rebuilt or changes in nutrient supply. The next most numerous group of miRNA targets was categorized into the 'Genetic information processing' KEGG pathways, namely, 'spliceosome', 'RNA transport', and 'ubiquitin proteolysis'. This suggests that in yellow lupine flowers most miRNAs regulate processes related to post-transcriptional events and protein degradation. Three KEGG categories within the 'Environmental information processing' category is extremely important in terms of plant development, and they are 'Signal transduction pathways' comprising the MAPK cascade, 'phosphatidylinositol' and 'plant hormone' signaling pathways ( Figure 10, Figure S6, S7, S8). The MAPK pathway is involved in regulating several processes, such as biotic and abiotic stress response (reviewed in [125,126]), and associated with the functioning of hormones such as ethylene [127] and abscisic acid, engaged in organ abscission and other processes (reviewed in [128,129]). The MAPK cascade is also an element of the positive feedback loop amplifying the abscission signal [130]. Auxin seems to be major target of sRNAs in yellow lupine. However, KEGG enrichment analyses of the identified target genes for lupine miRNAs indicated that the signal transduction pathways of gibberellin, cytokinin, the already mentioned ethylene, and ABA were potentially modulated by miRNAs in L. luteus, as well, but in less extent ( Figure 10).
Interestingly, like in the case of GO analysis, KEGG analysis for the MIR166 family showed that it was involved in the auxin, cytokinin, and brassinosteroid signal transduction pathways ( Figure 10). These data show again how the fine-tuning of expression of phytohormone-related genes by sRNAs is important for growth and development regulation.

Plant Material
Commercially available seeds of yellow lupine cv. Taper were obtained from the Breeding Station Wiatrowo (Poznań Plant Breeders LTD. Tulce, Poland). Seeds were treated with 3,5ml/kg Vitavax 200FS solution (Chemtura AgroSolutions, Middlebury, United States) to prevent fungal infections and inoculated with cultures of Bradyrhizobium lupine contained in Nitragina (BIOFOOD s.c., Walcz, Poland) according to seed producer's recommendations [131]. All the research material used for RNA isolation was collected from field-grown plants cultivated in the Nicolaus Copernicus University's experimental field (in the area of the Centre for Astronomy, Piwnice near Torun, Poland, 53 • 05'42.0"N 18 • 33'24.6"E) according to producer's agricultural recommendations [131] until the time of flowering. Flowers were collected 50 to 54 days after germination from the top and bottom parts of the inflorescence and were separated into four categories based on the progression of their development. Flower pedicels from flowers undergoing abscission or maintained on the plant were also collected, as in our previous study [17].
Plants with the same number of flower whorls were selected for the flower removal experiment and control. All plants were grown as described above up to the flowering stage. When plants reached the stage at which the top-most flowers were in the developmental stage S1, other flowers were removed from the inflorescence (UFR samples). The samples were collected for the gene expression analysis in the stages S1-S4. As a control, flowers from stages S1-S4 from upper (UF) and lower (LF) part of the inflorescence were collected.

RNA isolation, Library Construction, and RNA Sequencing
Total RNA from at least 5 plants (25 flowers or pedicels) for each variant was performed using the miRNeasy Mini Kit (Qiagen, Venlo, the Netherlands) including on-column DNA digestion with the RNase-Free DNase Set (Qiagen, Venlo, the Netherlands). The total RNA quality and quantity were evaluated with agarose gel electrophoresis and Nanodrop ND-1000 spectrophotometer (Thermo Scientific Waltham, MA, USA). Both the RNA Integrity Number (RIN), and RNA concentration were measured with the 2100 Bioanalyzer (Agilent Santa Clara, CA, USA) using the Small RNA Kit (Agilent Santa Clara, CA, USA). All the samples had adequate concentrations of RNA and RIN ranging from 8.9 to 10.0 and were sent for library construction to Genomed S.A (Warszawa, Poland) and sequencing BGI (Shenzhen, China).
Small RNA libraries were prepared from the total RNA using the NEBNext Multiplex Small RNA Library Prep kit for Illumina (New England Biolabs, Ipswich, MA, USA) and subsequently sequenced on the HiSeq4000 platform (Illumina, San Diego, CA, USA) in the 50 single-end mode. All libraries were constructed in two biological replications resulting in a total number of 20 sRNA libraries.
The total RNA extracted from pooled material derived from three biological replicates was used to prepare ten transcript libraries using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) and sequenced on the HiSeq4000 platform in the 100 paired-end mode.
For degradome sequencing, total RNA from three biological samples of UF3 and LF3 was pooled to maximize the amount of required material. The protocol for degradome library preparation comprised the following steps: (i) mRNA isolation, where poly (A)-containing mRNA molecules are purified from total RNA using poly(dT)oligo-attached magnetic beads, (ii) synthesis and subjugation of cDNA to ligate 5 adaptors, and purification of the resulting products with TAE-agarose gel electrophoresis, (iii) PCR amplification to enrich the final products, (iv) library-quality validation on the Agilent Technologies 2100 Bio-analyzer and using the ABI StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA), and (v) sequencing of the prepared library on the HiSeq4000 platform in the 50 single-end mode.

De Novo Transcriptome Assembly and Gene Expression Analysis
The transcriptome was assembled de novo from RNA-Seq data using Trinity v 2.4.0 (https: //github.com/trinityrnaseq/trinityrnaseq/releases) with default settings as in our previous study [17]. The expression level was estimated at both the unigene and isoform levels and described by FPKM (Fragments Per Kilobase Of Exon Per Million Fragments Mapped): the number of reads per unigene normalized to the library size and transcript length using RSEM [132] as previously described [17].
Short reads were compared against noncoding RNAs from Rfam [49,50] and both mature miRNAs and their precursors from miRBase [51]. The comparison was performed with Bowtie [133] allowing for no mismatches.
To identify phylogenetically conserved mature miRNAs with sequences and lengths identical to known plant miRNAs we searched miRBase for similarity at the mature miRNA level.
To predict potential novel miRNAs we applied ShortStack [53] with default settings. This tool identifies miRNAs based on their mapping against a reference genome. Since no genome was available for the studied species, we used de novo approach for transcriptome assembly instead. The latter method allowed for identification of miRNAs that showed no similarity to miRNAs annotated in miRBase and these miRNAs were assigned as new.
ShortStack [53] was used to identify small RNAs that were being cut in phase from longer precursors (phased siRNAs) with transcriptomes used as references. The top 200 candidates were selected from each sample, based on the phased score value provided by ShortStack. Finally, lists of such sRNAs from all samples were merged into a single dataset of non-redundant phased siRNAs (Table S6). The expression values were calculated as in the case of miRNAs.

Small RNA Expression Quantifications and Analysis of Differentially Expressed si-and miRNAs
MiRNA counts within each sample were first normalized to RPM values (reads per million values) and then a differential expression analysis was performed with the DESeq2 R package [134]. The files containing raw read counts for miRNAs/siRNAs from treatment and control replicates were used as input, and only candidates with an adjusted p-value (q-value) below 0.05 were retained for further analysis.

Identification of sRNA Targets
For target prediction using degradome analyses after sequencing, the reads were filtered using fastq_quality_filter from the Fastx-Toolkit package (http://hannonlab.cshl.edu/fastx_toolkit/) with at least 95% of nucleotides in each read demonstrating quality >= 20 (Phred Quality Score) with -p 95 and -q 20. The filtered Degradome-seq data, sequences of mature miRNA/siRNA and the assembled transcriptomes were processed with the CleaveLand4 package [58] to determine the cleavage sites for sRNA using default program settings. The final results were filtered based on the p-value < 0.05.
To predict targets for known or novel miRNAs, and phased siRNAs, we used also the psRNATarget tool [59] querying the assembled transcriptomes with the default Schema V2 (2017 release) and an expectation score of up to 4.

Evolutionary Conservation of miRNAs
L. luteus miRNAs were assigned to miRNA families based on miRBase classification, and the same was done for the sequences of all Eudicotyledons species present in miRBase, with the exclusion of Gossypium arboretum (which has only one sequence deposited in the database that cannot be classified as belonging to any known miRNA family). miRNAs from 52 species were compared against L. luteus miRNAs in order to count the numbers of miRNA family members shared amongst the species. The same analysis was performed with data narrowing to nine Fabaceae species.

Expression analysis with RT-qPCR
MiRNAs and siRNAs expression analysis was performed using the Stem Loop RT-qPCR technique, according to [56] with some modifications. An RT primer specific for each sRNA was used for the reverse transcription using total RNA of each sample and the SuperScript III Reverse Transcriptase (Thermo Fisher Scientific, Waltham, MA, USA) in a 20 µL reaction volume. To increase the accuracy and efficiency of the reaction, the pulse RT approach [57] was applied to the reverse transcription which consisted of two steps: 30 min of pre-incubation at 16 • C, followed by 60 cycles at 30 • C for 30 s, 42 • C for 30s and 50 • C for 1 s. qPCR was subsequently performed using specific primers designed according to [57], modified so that the UPL9 hydrolysis probe (Roche, Basel, Switzerland) was used for maximization of accuracy and background reduction. This reaction was performed using the SensiFAST Probe No-ROX kit (Bioline meridian bioscience Cincinnati, OH, USA) and the LightCycler 480 (Roche, Basel, Switzerland). Each 20 µL reaction contained: 1 µL cDNA template (transcribed from 100 ng of total RNA for less expressed miRNAs and 25 ng of total RNA for more expressed miRNAs), 1 µL of 10 µM qPCR specific forward primer, 1 µL of 10 µM Universal-qPCR primer, 10 µL of 2× SensiFAST Probe No-ROX Mix, 0.2 µL of 10 µM UPL9 probe and 6.8 µL ddH 2 O. qPCR was executed by pre-incubation at 95 • C for 10 min, followed by 45 cycles of 95 • C for 10 s, 59 • C for 30 s, and 72 • C for 1 s. Target gene expression was performed as in [17]. Each experiment consisted of three biological and technical replicates. The relative expression levels were calculated using the 2 −∆∆Ct method, and the data were normalized to the CT values for the LlActin reference gene (according to [17]). All primer sequences are given in Table S10.

Gene Ontology (GO) and KEGG Analysis of Target Genes
In order to estimate the potential roles of L. luteus sRNAs in biological processes, GO annotations of their target genes were downloaded from the Gene Ontology using NCBI or UniProt identifiers The Bioconductor GOseq package [135] was used for GO enrichment analysis. KEGG annotation and enrichment analysis were performed to determine the metabolic pathways. The GO terms and KEGG pathways were considered to be significantly enriched with the corrected p-value of 0.05, which was calculated using a hypergeometric test [136].

Data submission to Sequence Read Archive NCBI
The RNA-Seq and small RNA-Seq data have been uploaded to the SRA database and are available under BioProject ID PRJNA419564 and Submission ID SUB3230840.

Conclusions
In this paper, we present the first case of identification and integrated analysis of small ncRNA, transcriptome, and degradome sequencing data, which allowed us to identify known and novel miRNAs, siRNAs and their target genes probably involved in regulating yellow lupine flower development and abscission. These miRNAs and siRNAs, by controlling the expression of their target genes, may have an impact on the development and fate of flowers growing in particular parts of the inflorescence (Figure 11). There appear to be microRNAs controlling auxin signal transduction elements and proliferation regulators in n the central node of the regulatory network controlling flower development or abscission. In addition to the purely cognitive aspects of describing the evolutionary conservation and the species specificity of important mechanisms regulating plant development, this work may contribute to the optimization of field crops and to monitoring the impact of various factors on flowering in yellow lupine. The use of the NGS technique allows for a detailed analysis of the regulatory networks which include sRNAs and their target genes. However, the results of sRNA-seq also contain a large number of uncharacterized sRNAs, the function of which may also have significance for the studied processes. More experimental and bioinformatic research is needed to fully describe the complex mechanisms of plant development regulation by low-molecular-weight regulatory RNAs. Figure 11. MiRNAs and siRNAs participating in yellow lupine flower development and abscission. Scheme based on the results of the current analysis. Arrows pointing upwards and downwards represent sRNAs that are up or downregulated in the transition between two developmental stages, respectively. The plus sign marks significantly expressed sRNAs. Colored circles represent targets found in the degradome, colored squares represent targets found using psRNATarget, as listed below. Multiple markers indicate that the sequence has multiple targets. Abbreviated gene names were acquired from UniProt database, where full target names can be found. Pictures from left to right are as follows: 4-whorl inflorescence of yellow lupine, flower cross-sections and isolated pistils for each stage of development, cross-sections of abscissing and non-abscissing flower pedicels stained with phloroglucinol-HCL solution.  Scheme based on the results of the current analysis. Arrows pointing upwards and downwards represent sRNAs that are up or downregulated in the transition between two developmental stages, respectively. The plus sign marks significantly expressed sRNAs. Colored circles represent targets found in the degradome, colored squares represent targets found using psRNATarget, as listed below. Multiple markers indicate that the sequence has multiple targets. Abbreviated gene names were acquired from UniProt database, where full target names can be found. Pictures from left to right are as follows: 4-whorl inflorescence of yellow lupine, flower cross-sections and isolated pistils for each stage of development, cross-sections of abscissing and non-abscissing flower pedicels stained with phloroglucinol-HCL solution.