Integrated analysis of small RNA, transcriptome and degradome sequencing provides new insights into floral development and abscission in yellow lupine (Lupines luteus L.)


 Background

Yellow lupine (Lupinus luteus L., Taper c.) is an important legume crop. However, its flower development and pod formation are often affected by excessive abscission. Organ detachment occurs within the abscission zone (AZ) and in L. luteus primarily affects flowers formed at the top of the inflorescence. The top flowers’ fate appears determined before anthesis. The organ development and abscission mechanisms utilize a complex molecular network, not yet not fully understood, especially as to the role of miRNAs and siRNAs. We aimed at identifying differentially expressed (DE) small ncRNAs in lupine by comparing small RNA-seq libraries generated from developing upper and lower raceme flowers, and flower pedicels with active and inactive AZs. Their target genes were also identified using transcriptome and degradome sequencing.

Results

Within all the libraries, 394 known and 28 novel miRNAs and 316 phased siRNAs were identified. In flowers at different stages of development, 30 miRNAs displayed DE in the upper and 29 in the lower parts of the raceme. In comparisons between upper and lower raceme flowers, a total of 46 DE miRNAs were identified. miR393 and miR160 were related to the upper and miR396 to the lower flowers and pedicels of non-abscising flowers. In flower pedicels we identified 34 DE miRNAs, with miR167 being the most abundant of all. Most siRNAs seem to play a marginal role in the processes studied herein, with the exception of tasiR-ARFs, which were DE in the developing flowers. The target genes of these miRNAs were predominantly categorized into the following Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways: ‘Metabolism’ (15,856), ‘Genetic information processing’ (5,267), and ‘Environmental information processing’ (1,517). Over 700 putative targets were categorized as belonging to the ‘Plant Hormone Signal Transduction pathway’. 26,230 target genes exhibited Gene Ontology (GO) terms: 23,092 genes were categorized into the ‘Cellular component’, 23,501 into the ‘Molecular function’, and 22,939 into the ’Biological process’.

Conclusion

Our results indicate that miRNAs and siRNAs in yellow lupine may influence the development of flowers and, consequently, their fate by repressing their target genes, particularly those associated with the homeostasis of phytohormones, mainly auxins.

flowers and, consequently, their fate by repressing their target genes, particularly those associated with the homeostasis of phytohormones, mainly auxins.

Background
Yellow lupine is a crop plant with remarkable economic potential. Thanks to its ability to fix nitrogen, it does not need fertilizers, and its protein-rich seeds may be an excellent source of protein in animal feed [1][2][3]. The crop yield is determined by the number of formed pods, and in yellow lupine this value is significantly limited by its high flower abscission level [1]. 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 1 st (and lowest) whorl, 90% at the 2 nd whorl, and ~100% at the whorls above them. Thus, the problem of flower abscission generates large economic losses in agriculture [1].
Plant organ abscission is an element of the developmental strategy related to reproduction, defense mechanisms or disposal of unused organs [5,6]. In most species, the key players involved in activation of the abscission zone (AZ) are plant hormones, among which auxin (IAA) and ethylene (ET) play the key roles [7,8]. An increased concentration of auxin in the distal region of pedicel (between the flower and the AZ), as compared to its concentration in the proximal region (between the AZ and the peduncle), is a prerequisite for organ abscission [5,6,9]. Our previous transcriptome-wide study [10] proved that the abscission of yellow lupine flowers and pods is associated, inter alia, with intensive changing of auxin catabolism and signaling. Genes encoding auxin response factors ARF4 and ARF2 were more expressed in generative organs that were maintained on the plant, in contrast to the mRNA encoding auxin receptor TIR1 (Transport Inhibitor Response 1), which is accumulated in larger quantities in shed organs [10].
Since (i) some micro RNAs (miRNAs) and small interfering RNAs (siRNAs) restrict the activity of certain ARFs [11,12] and members of the TAAR (TIR1/AFB AUXIN RECEPTOR) family encoding auxin receptors [13], and since (ii) we proved that the precursor of miR169 is accumulated in increased quantities in yellow lupine's generative organs undergoing abscission [10], we predict that sRNAs play significant roles in orchestrating organ abscission in L. luteus. miRNA is a cluster of 21-22-nt-long regulatory RNAs formed as a result of the activity of MIR genes in certain tissues and at certain developmental stages [14][15][16], and also in response to environmental stimuli [17][18][19]. MIR genes encode two consecutively formed precursor RNAs, first pri-miRNAs and then pre-miRNAs, which are subsequently processed by DCL1 (Dicer-like enzyme) into mature miRNAs [20,21]. MIR genes are often divided into small families encoding nearly or completely identical mature miRNAs [22]. miRNA sequences of 19-21 nucleotides are long enough to enable binding particular mRNAs by complementary base pairing, and allow either for cutting within a recognized sequence or for translational repression (reviewed in [23]). Plant miRNAs are involved in, for instance, regulating leaf morphogenesis, establishment of flower identity, and stress response [17][18][19][24][25][26].
Some of them also form a negative feedback loop by influencing their own biogenesis, as well as the biogenesis of some 21-nt-long siRNAs called trans-acting siRNAs (ta-siRNAs). ta-siRNAs are processed from non-coding TAS mRNAs, which contain a sequence complementary to specific miRNAs [27,28].
There is also a large group of plant sRNAs that are referred to as phased siRNA, which are formed from long, perfectly double-stranded transcripts of various origins, mainly processed by DCL4 [29,30].
Numerous studies on sRNA in legumes (e.g. Glycine max, G. soja [ focused on stress response and nodulation (also reviewed in [37]), or have been entirely carried out in silico [38]. In lupines, the knowledge on the roles of regulatory RNAs function, is still negligible [39].
Importantly, the involvement of regulatory sRNAs in mechanisms behind the maintenance/abscission of generative organs has never been explored before. This study aims to investigate the role of low molecular regulatory RNAs (miRNAs and phased siRNAs) and their target genes in Lupinus luteus flower development and abscission.
To achieve this goal, ten variants of sRNA libraries were prepared and NGS sequencing 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 to 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-wide analysis was carried out to identify the target genes for the recognized sRNAs, and a degradome-wide analysis to confirm the processing of the target sequences for a number of conserved or new Lupinus luteus sRNAs. Our results of an NGS analysis of small RNA, transcriptome and degradome libraries indicate that miRNAs play a vital role in regulating yellow lupine flower development, both generally and depending on the flower's location on the inflorescence. Furthermore, these molecules also display differential accumulation during flower abscission.

Plant Material and Total RNA Extraction
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 live cultures of Bradyrhizobium lupine contained in Nitragina (BIOFOOD s.c., Walcz, Poland) according to seed producer's recommendations [187].
All the research material (flowers and flower pedicels) used for RNA isolation was collected from fieldgrown 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 [187]  were also collected, as in our previous study [10]. The full list of sample names with descriptions is available in the Additional File 1: Table S1.
Total RNA extraction was performed using the miRNeasy Mini Kit (Qiagen, Venlo, the Netherlands) with a modified protocol including an increased volume of the QIAzol lysis reagent in the initial phase of extraction (due to the exceptionally high level of nucleic acids in some of the samples), and oncolumn DNA digestion was carried out with the RNase-Free DNase Set (Qiagen, Venlo, the Netherlands) to avoid any DNA contamination in the samples. Total RNA quality and quantity was

RNA Library Construction and NGS Sequencing
Small RNA libraries were prepared from the total RNA using the NEBNext Multiplex Small RNA Library Prep kit for Illumina (New England Biolabs Ipswich, Massachusetts, U.S.A) and subsequently sequenced on the HiSeq4000 platform (Illumina, San Diego, CA, U.S.A) in the 50SE (single-end) mode.
All libraries were constructed in two biological replicates.
Transcript libraries were prepared from the total RNA using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs Ipswich, Massachusetts, U.S.A) and sequenced on the HiSeq4000 platform in the 100PE (paired-end) mode.

Degradome Library Construction and Target Identification
For degradome sequencing, total RNA from UF3-LF3 samples was pooled. 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 (T) 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, California U.S.A.), and (v) sequencing of the prepared library on the HiSeq4000 platform in the SE50 mode.
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 demonstratinga 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 [56] to determine the cleavage sites for sRNA using default program settings.
The final results were filtered based on the P-value <0.05.
The construction and sequencing of all the above-mentioned libraries were performed by Genomed S.A (Warszawa, Poland) and BGI (Shenzhen, China).
Then, redundant and counting read occurrences (i.e. raw expression values) were identified with the fastx_collapser from the same package.
The nucleotide length distribution of small RNAs from all the ten libraries ( Fig. 1) was plotted using the ggplot2 R package [189]. The miRNA accumulation in all the ten libraries expressed in reads per milion (RPM) (Fig. 4) was plotted using the R packages: ggplot2 [189] and gridExtra [190].

Identification of Known and Novel miRNAs
A search for miRNAs was done in the two following ways. Mature miRNAs with sequences and lengths identical to known plant miRNAs from miRBase were identified using a similarity search at the mature miRNA level. In an alternative approach, we applied ShortStack [52] 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 a de novo approach for transcriptome assembly instead.
The latter method allowed for identification of miRNAs that were not similar to anything within miRBase, and these miRNAs were assigned as new.

Identification of ta-siRNAs (phased siRNAs)
ShortStack [52] was used to identify small RNAs that were being cut in phase from longer precursors.
Similarly to miRNAs, the transcriptomes were used as references, since no genome sequence was available. In order to obtain the highest quality set of these sRNAs, only 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 nonredundant phased siRNAs (Additional File 1: Table S8). The expression values were calculated in the similar manner as in the case of miRNAs.

Analysis of Differentially Expressed si-and miRNAs
A differential expression analysis was performed with the DESeq2 R package [191].
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 kept. The heatmaps displaying the expression profiles and hierarchical clustering of miRNAs, the fold change of differentially expressed miRNAs and p-values (Fig. 4) were plotted using the following R packages: ComplexHeatmap [192], circlize [193] and dendextend [194]. For flowers (Fig. 4a), the expression matrix was z-scaled prior to clustering using a scale() R function.
Due to the ineffectiveness of the Venn diagram approach in the visualization for more than 5 sets, common miRNAs in each intersection of 8 datasets were visualized using the UpSetR R package [55].
The frequency cutoff was set at zero (Additional File 3: Figure S2).

Evolutionary Conservation of miRNAs
As a result of adopting the similarity search approach in our miRNA search strategy, all the identified miRNAs were evolutionarily conserved, i.e. they could be found in at least one other plant species. In order to get more insight into the evolutionary relationship of miRNAs between species (at the miRNA family level), an additional analysis was performed. In detail, 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. Same analysis was performed again for 9 Fabaceae species.

Validation of sRNA-seq Expressed miRNAs
Validation of differentially expressed (DE) miRNAs and siRNAs was performed using the Stem Loop RT-qPCR technique, according to [195]. To increase the accuracy and efficiency of the reaction, the pulse RT approach was applied to reverse the transcription step based on [196] using total RNA and the SuperScript III Reverse Transcriptase (Thermo Fisher Scientific Waltham, Massachusetts U.S.A) in a 20 µl reaction volume. The reverse transcription consisted of two steps: 30 min of pre-incubation at 16°C, followed by 60 cycles at 30°C for 30s, 42°C for 30s and 50°C for 1s. The subsequently performed qPCR was performed using specific primers designed according to [196], modified so that the UPL9 hydrolysis probe (Roche, Basel, Switzerland) could be used for maximum accuracy and background reduction. This reaction was performed using the SensiFAST Probe No-ROX kit (Bioline meridian bioscience Cincinnati, Ohio, U.S.A) 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 2x 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. The signal from the probe was detected during the amplification step. The threshold was set automatically and the threshold cycle (CT) was recorded. Each experiment consisted of three biological and technical replicates. The relative expression levels of the miRNAs were calculated using the method, and the data were normalized to the CT values for the β actin gene (expression level of this reference gene was measured according to [10]). All primer sequences are given in the Additional File 2: Table S13.

RNA-Seq: de novo Transcriptome Assembly and Target Gene Expression Analysis
The transcriptome was assembled de novo from RNA-Seq data using Trinity v 2.4.0 with default settings as in our previous study [10]. The expression level was estimated at both the unigene and isoform levels and described by FPKM: the number of reads per unigene normalized to the library size and transcript length (Fragments Per Kilobase Of Exon Per Million Fragments Mapped) using RSEM [197] as previously described [10].

Finding Targets for sRNAs
To find targets for known and novel miRNAs, we used the psRNATarget tool [57] with the default Schema V2 (2017 release) and an expectation score of up to 4. Here, the miRNA targets were searched for among the assembled transcriptomes. The same analysis was performed for phased siRNAs.

Gene Ontology (GO) and KEGG Analysis of Target Genes
In order to estimate the potential roles of the identified sRNAs and transcripts in biological processes, GO annotations of genes targeted by the identified miRNAs were downloaded from the Gene Ontology10, NCBI11, and UniProt12. The Bioconductor GOseq package [198] was used for a GO enrichment analysis. A KEGG annotation and enrichment analysis was 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 [199].

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.

Sequencing and Annotation of Yellow Lupine sRNAs
Ten variants of small RNA libraries from developing flowers and flower pedicels were constructed and sequenced as described in the material and methods section and Additional File 1: Table S1. Singleend deep sequencing was performed on the Illumina HiSeq4000 platform. 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 (Additional File 1: 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 class amount of the clean and redundant reads (Additional File 1: Table S2, Fig. 1), which was compliant with many other RNA-Seq experiments (e.g. [40][41][42]) and correlated with the abundance of siRNAs and miRNAs, respectively. The unique reads were annotated against Rfam [43,44] and miRBase [45][46][47][48][49][50] 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, Additional File 1: Table S2). The unique sequences were annotated into different RNA classes against the Rfam database using BLAST [51].
The small RNA sequences have been clustered into several RNA classes, such as known miRNAs, rRNA, tRNA, sn/snoRNA and others (  Table S3).

Identification of Known, Conserved miRNAs
After analyzing the results of the alignment against miRBase [45][46][47][48][49][50], 394 unique miRNAs containing 366 conserved miRNAs were identified in the ten samples (Additional File 1: Table S4). The number of the identified miRNAs in each library is shown in Table S. The largest number of known miRNAs (181) was observed in the Lower Flowers stage 3 (LF3) and the smallest (148) in the Upper Flowers stage 1 (UF1) (Additional File 1: Table S1). The identified miRNAs belonged to more than 67 families (Additional File 1: Table S4), while most of them belonged to the MIR156, MIR159 and MIR166 families, with more than 35 members in each (Fig. 2a).

Evolutionary Conservation of microRNAs Identified in Lupinus luteus
To explore the evolutionary roles of the known miRNAs identified in yellow lupine, these RNAs were compared with the sequences of almost all (i.e. 52) Eudicotyledons species present in miRBase [45][46][47][48][49][50]. Gossypium arboretum was excluded from this comparison as there is only one (and, additionally, unclassified) sequence deposited in the database for it. The same analysis was performed exclusively against nine Fabaceae species. As shown in Figure 2b

Identification of Novel miRNAs
With the use of the ShortStack software [52], 28 candidates for novel miRNAs were identified (Table   3). This tool identifies miRNAs based on their mapping against a reference genome. As no genome was available for the studied species, we used a transcriptome instead (statistical data on de novo assembly is in Additional File 1: Table S5). 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 21nt long (Table 3). Interestingly, many of these unique miRNAs showed similarity to precursor miRNAs from miRBase, which led to the conclusion that they might be new members of the already known families: MIR167 (Ll-miRn12 and -n27) and MIR169 (Ll-miRn3, -n11 and -n15) ( Table 3) and hinted that they might be exclusive to yellow lupine. The length of their precursor sequences ranged from 71 to 246 nt and their free energy (MFE) ranged from -32,3 to -144,9 kcal/mol for the secondary hairpin structure. The expression of these 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 (Additional File 1: Table S6).

Analysis of the Expression Abundance of Known miRNA Families in Lupinus luteus
Since miRNA expression across all libraries displayed high variation, we put the data into five categories based on the maximum value (Fig. 3). Two miRNAs, namelyLl-miR341 and Ll-miRn21, showed expression maxima of over 10,000 RPM. Maximal expression of another two, Ll-miR260 and Ll-miR384, ranged from 2,000 to 10,000 RPM. Thirteen miRNAs showed expression maxima ranging from 100 to 2,000 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 (Fig. 3, Additional File 1: Table S7).

Identification of Phased siRNA in Yellow Lupine
Numerous reports and studies point to the important role of phased siRNA not only in stress response mechanisms, but also in growth regulation (reviewed in e.g. [53]). Therefore, we decided to investigate the role of siRNAs during yellow lupine inflorescence development. To achieve this, ShortStack [52] 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 24nt long (Additional File 1: Table S8, Additional File 1: Figure S1). The identified siRNAs displayed a highly differential expression pattern. The highest level of expression was observed in Ll-siR308, -274, -249, -206, -191, -161 and -108 (Additional File 1: Table S9), and ranged from 39 to 579 RPM in all the samples. 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) (Additional File 1: Table S9). On the other hand, Ll-siR246, -291 and -56 were present almost exclusively in the youngest flowers in the lower part of inflorescence (LF1). The least abundant siRNAs were Ll-siR50, -175, -250 and-263, which were present only in one library (Additional File 1: Table S9).

Development
In order to investigate the dynamics of the identified sRNAs' expression 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, Additional File 1: Table S10 (Table 4), as well as 14 and 7 siRNAs, respectively (Additional File 1: Table S11).
Between UF2 vs UF1, there was a change in the expression of 8 miRNAs, of which 4 were upregulated. Among these DE sequences, there were 2 sequences belonging to MIR359 and MIR166 families each, as well as one representative of each of the MIR159, MIR167, MIR396 and novel Ll-miRn10 families. Ten DE miRNAs were identified in a comparison of UF3 vs UF2, of which only one Ll-miR258 (miR166 homologue) was up-regulated. The remaining miRNAs were down-regulated 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 DE miRNAs, of which 7 were up-regulated. The most numerous group was that of 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 was different. The highest number of the identified DE miRNAs were 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 DE miRNAs (6 of which were up-regulated) the most numerous group were novel miRNAs, with 5 sequences, namely Ll-miRn22, -27, -10, -3 and -31, followed by members of the MIR396 family (4 sequences, all down-regulated). Subsequently, among the DE miRNas shown in this comparison we found members of the MIR156 and MIR168 families (2 each), followed by single sequences from the MIR166, MIR159, MIR390, MIR858 and MIR393 families. Between the LF2 and LF3 stages we confirmed that there was a change in the accumulation of 11 miRNAs (6 of which were up-regulated), 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).
As regards siRNAs during flower development in yellow lupine, almost every differentially expressed siRNAs were up-regulated. In the lower part of the inflorescence, similarly to miRNAs, there were no differences between the LF4 and LF3 stages. During upper flower development, the most DE siRNAs 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's flowers, regardless of their position on the inflorescence. The complete dataset can be found in (Additional File 1: Table S11).
In order to identify miRNAs the presence of which is either common and unique depending on the developmental stage of the upper and lower flowers in lupine, Venn diagrams were constructed (Fig.  5 a and a, Additional File 1: Table S12) using Venny 2.1 [54]. The results of these analyses revealed that approximately 70% of the identified miRNAs were common in all developmental stages of both the upper (Fig. 5a) and lower flowers. (Fig. 5b). However, miRNAs unique to certain developmental stages were also found ( Fig. 5 and Additional File 1: Table S12). The findings concerning the abundance of miRNAs in various combinations of 8 datasets are displayed in Additional File 3: Figure   S2 as an UpSet R plot [55]. The total numbers of miRNAs detected in each experimental variant are shown as horizontal bars in the bottom left corner. The histogram demonstrates the number of miRNAs corresponding to a given combination, indicated by the dots connected by lines below each bar. The most numerous set with 155 elements was the one with an intersection of all 8 variants. As many as 17 miRNAs were assigned exclusively to stage 3 flowers collected from the lowest whorl (LF3), while for the other 7 variants this number ranged between 6 and 2.

Comparison of Differentially Expressed sRNAs Between Developing Flowers From the Lower and Upper Whorls of the Raceme
In order to determine the differences in sRNA expression in developing yellow lupine flowers, comparative analyses of both the upper and lower flowers was performed for each developmental stage of the inflorescence (LF1 vs UF1, LF2 vs UF2, LF3 vs UF3 and LF4 vs UF4) ( Table 5, Fig. 4, Additional File 1: Table S10). In general, 46 DE miRNAs were identified, of which twelve DE miRNAs were from stage 1, eleven from stage 2, eleven from stage 3 and 2 from stage 4 of the flower development process. Interestingly, all UF1 miRNAs were down-regulated, contrary to stage 2 miRNAs all of which were up-regulated. As forUF3 miRNAs, 5 of them were up-regulated, with 6 downregulated. In UF4, one miRNA (Ll-miR445, homologous to miR319) was up-regulated, while Ll-miR100 homologous to miR390 was down-regulated (Table 5). In the first stage of development, the most numerous group of DE miRNAs was that of the novel sequences (Ll-miRn3, -25, -29 and -30), followed by sequences annotated as miR396 (3 miRNAs). The remaining miRNAs were single sequences annotated as miR167, miR395, miR390 and miR171. In the second stage of flower development,miRNAs belonging to MIR319 family were identified as the largest group (5 sequences), followed by two DE miRNAs annotated as miR160 and miR396, respectively, and singular ones annotated as miR394 and miR393. The third stage turned out to be the most diverse, with 2 representatives of the miR160 family, followed by single sequences annotated as miR394, -393, -858, -159, -396, -5490, -1861 and miR5168. As regards phased siRNAs, only 4 of them displayed differential expression, namely Ll-siR119 in stage 1 and Ll-siR224, -100 and -146 in 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 important than phased siRNAs are with regard to yellow lupineflower differentiation.
Analyses of the Venn diagrams we created (Fig. 5c), 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 (Fig. 5c).
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 Additional File 1: Table S12. It is worth noting, however, that Ll-miR224/miR393 were only present in the upper flowers in every stage of development, while in the lower flowers no common miRNAs were found. Despite this, in the lower flowers the presence of members of the MIR396 family (two in LF1, four in LF3 and two in LF2), 3 members of MIR171 in LF1, 3 members of MIR395 in LF2, and two members of the MIR159 family in LF3 and LF4 each, was observed. Interestingly, the presence of miR319 (two members for UF1, UF2 and UF4 each) could be found in the upper flowers apart from Ll-miR224. Additional analyses of miRNAs common for specific developmental stages regardless of the flower's position on the inflorescence (UF1LF2, UF2LF2, UF3LF3, UF4LF4) (Additional File 1: Table S12, Additional File 4: Figure S3a) revealed that for stage 1 there were no unique and common miRNAs (that would not be present in the other stages). For stages 2 and 3, 6 miRNAs each were identified, with 3 unique miRNAs for stage 4.

Inactive Abscission Zones
In order to identify sRNAs possibly involved in yellow lupineflower 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 phased siRNAs (Additional File 1: Table S11) 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 downregulated miRNAs, the most abundant were MIR390, MIR396 and MIR395 with 3 members each ( Table   6, Fig. 4b). With regard to siRNAs, the most up-regulated in FPANB were Ll-siR173, -4 and -13, and the most down-regulated was Ll-siR208 (Additional File 1: Table S11).

Validation of the Identified sRNAs in RNA-seq
In order to validate the data generated using deep sequencing technology and to confirm the expression patterns of the identified sRNA, we employed stem-loop (SL) qRT-PCR of the 8 identified sRNAs (six conserved miRNAs, one novel miRNA and two siRNAs) (Additional File 2: Table S13). The qPCR results were similar to sRNA-seq data (Fig. 6). For example, in the qRT-PCR 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 distinct and strong expressions of miR155 in LF1 and miR281 in FPAB were also similar in both the experiments (Fig. 6). 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 first. To this effect, we constructed degradome libraries from pooled samples of stage 3 upper and lower parts of the inflorescence (Additional File 2: Table S14). Through total degradome library sequencing, 19,353,278 raw reads were obtained (Additional File 2: Table   S14a). After quality filtering, the degradome data were aligned to the reference transcriptome with CleaveLand 4 [56] to find sliced miRNA and siRNA targets. After processing and an 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) (Additional File 2: Table S15). For the phased siRNAs, 3,340 targets were initially identified, and after filtering their number dropped to 89 targets for 46 siRNAs (Additional File 2: Table S16). The results obtained through CleaveLand4 were divided into categories from 0 to 4 in the following manner: 4 -only one read at this position, 3 ->1 read, but below or equal to the average* depth of coverage on the transcript, 2 ->1 read, above the average* depth, but not the maximum on the transcript, 1 ->1 read, equal to the maximum on the transcript when there was >1 position at the maximum value, 0 ->1 read, equal to the maximum on the transcript when there was just 1 position at the maximum value [56]. Through our analyses, for the categories ordered from 0 to 4, we obtained (after p-value filtering) 353, 54, 65, 17, 47 miRNA targets and 33, 18, 11, 19, 8 siRNA targets, respectively (Additional File 2: Table S14b). Exemplary target t-plots and sequences of the miRNAs and target mRNAs are shown in Figure 7.
As expected, many of the targets for evolutionarily conserved miRNAs were compliant with literature data. For example, Ll-miR329/miR160-5p (sRNA-seq ID/miRBase annotation) targeted ARF16 and ARF18, the Ll-miR415/miR171b targeted SCL6, Ll-miR341/miR319q targeted TCP2, Ll-miR224/miR393a-5p targeted TIR1, Ll-miR401/miR156k targeted SPL2 and SPL12, and Ll-miR229/miR396a-5p targeted GRF3, GRF4, and GRF5 (Additional File 2: Table S15). A comparison of the expression of 4 exemplary miRNAs and their target genes (Fig. 7) confirmed the reversecorrelation 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 plausible missing target genes, the psRNATarget tool (2017 update) [57] was employed, which rendered putative 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 (with expectation from 0 to 4) 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 Additional File 2: Table S16. Targets for all of them are shown in Additional File 2: Tables S15, 16, 17 and 18.

Function of the Potential miRNAs Targets
Our transcriptome analysis against the Kyoto Encyclopedia of Genes and Genomes (KEGG) revealed the putative function of many targets of the miRNAs identified in yellow lupine. 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) (Additional File 5: Figure S4). A full list of KEGG pathways and numbers of the assigned sequences is shown in (Additional File 6: Figure S5 and Additional File 7: Figure 6). 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 (Fig. 8). The second most notable pathway was mitogenactivated protein kinase (MAPK) signaling, which was associated with different abiotic and biotic stress factors, with 350 putative targets distributed across every described stress response (Additional File 8: Figure S7). A complete dataset on the KEGG analysis can be seen in (Additional File 2: Table S19).
In order to investigate the functions played by the miRNAs identified in yellow lupine flowers, we submitted miRNA targets to a Gene Ontology (GO) category analysis. Among the 27,547 targets of known and novel miRNAs identified with psRNATarget, 26,230 targets exhibited GO terms (Additional File 2: Table S20) Table S20). Within the 'Flower development' category, the targets of 37 miRNAs fit within GO terms related to phytohormones (Fig. 9b), and the targets of 69 miRNAs were categorized to GO terms related to the development of flower parts (Fig. 9c).

Discussion
Yellow lupine has a great potential as a substitute for soybean in animal feed production in Europe.
An ability to limit the economic drawbacks resulting from excessive flower abscission would be the most convincing argument for lupine cultivation. This, however, can only be achieved if we understand the plant's biology, especially the molecular basis for the development and maintenance of lupine flowers, of which little is known. Therefore, we believe that the pathways controlling these processes deserve intensive research focus. The commonly used NGS technique allows for comprehensive analyses of sRNAoms, transcriptoms, and degradomes, as well as identification of sRNAs and their target genes in both model and non-model plants, for example: tomato [58], soybean [59], Arabidopsis [60,61], oilseed rape [62], Brassica juncea [63], lily [64,65], peanut [66], orchardgrass [67] or wheat [68,69].
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 [10]. What is puzzling, our observations of L. lupinus floral development indicate that their fate (abscission or pod formation) is determined prior to AZ activation. Therefore, we performed comparative analyses between sRNAs from flowers developing on 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.

Functions of miRNA/Target Modules During Yellow Lupine Flower Development
Involvement of sRNAs in flower development has undergone intense investigation and already been explored in many plant species, such as wild rice (Oryza rufipogon) [70], hickory (Carya cathayensis Sarg) [71], yellow-horn (Xanthoceras sorbifolia) [72], spring orchid (Cymbidium ensifolium) [73], and roses (Rosa sp.) [74]. Our sRNA-seq analyses have not only shed light on the molecular mechanisms that control flower development in one more plant species and confirmed the involvement of known miRNAs, such as miR159, miR167 or miR172, in this process [75], but we have also explored the roles of sRNAs in flower abscission and identified species-specific miRNAs.
In order to identify miRNAs and their target genes involved in flower development and abscission in L.
luteus we conducted comprehensive analyses of ten sRNA and ten transcriptome libraries, and

Yellow Lupine
We explored the large set of data concerning flower development that we obtained through sRNA-Seq analysis, focusing on three aspects: (i) conserved miRNAs involved in flower morphogenesis and development, (ii) expression dynamics during flower development, (iii) miRNAs specific for particular developmental stages.

Conserved miRNA families associated with flower morphogenesis and development
Among the known and conserved miRNAs we spotted 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 (reviewed in [26]).
Within the degradome data sets, we found target genes of two out of four differentially expressed L.
lupinus miRNAs belonging to the MIR156/157 family (Ll-miR124 and Ll-miR401), encoding squamosa promoter-binding-like proteins 13B (SPL13B) and SPL2, respectively. 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 [76]. In this plant, 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 [76], where it is targeted by miR156 to regulate pollen maturation [77].
The net of factors that regulate pistil development include miR156-targeted members of the SPL gene family, such as SPL3 and the miR156-resistant SPL8 gene, with both of them redundantly controlling gynoecium patterning through interfering with auxin homeostasis and signaling. What is striking, none of miRNAs identified in yellow lupine targets SPL3, which positively regulates the MADS box genes APETALA1 (AP1), FUL and LFY, central regulators of flowering [78]. In Arabidopsis, downregulation of miR156-targeted SPLs resulted in a shortened style and deformed gynoecium [79].
Other studies on mutants have shown, that the aforementioned SPL genes cooperate with ARF3, a master regulator of gynoecium development and morphogenesis [80], which directly regulates transcription of several genes involved, among others processes, in auxin synthesis, transport, signaling and response [81].
miR159 was shown to target the conserved GAMYB-like genes that are a part of the GA signaling pathway [82,83]. In A. thaliana, by affecting the expression of MYB33 (MYB DOMAIN PROTEIN 33) and MYB65, miR159 influences the LFY transcription level and gibberellin-dependent flowering induction.
Two transcription factors involved in pistil and stamen development, ARF6 and ARF8, contain the target site for miR167 [87,88]. The roles of miR167 and ARFs are conserved across different species [89]. For Arabidopsis, it has been proven that both these genes are involved in stamen filament elongation, anther dehiscence, stamen maturation and anthesis [90], and that arf6 arf8 double mutants suffer from poor pollination [91]. By inducing jasmonic acid biosynthesis, ARF6 and ARF8 initiate anther dehiscence and flower bud opening and regulate perianth growth [91][92][93] However, these genes may also act in a jasmonic acid-independent way, by influencing the transcriptional activity of GH3, which results in a shift in the free auxin level [94,95]. 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 [96].
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 ABCE model of floral development [97].
Overexpression of MIR172 causes formation of a phenotype characterized by the absence of perianth, transformation of sepals into pistils and early flowering [97].
Our study showed the presence of at least one member of all these families in flowers (Fig. 2, Additional File 1: Table S7), which indicated that in lupine the families were crucial in generative development, as well. MIR156 and MIR159 are the most numerous families in Lupinus luteus, which suggests they play fundamental roles in its flower development processes.

Clusters of miRNAs depending on expression dynamics
The differentially expressed miRNAs identified in yellow lupine flowers were clustered by the dynamics of their expression (Z-scaled expression, to be precise) (Fig. 4). The first of the four clusters comprised miRNAs the accumulation of which increased as the flowers developed. The second cluster grouped all miRNAs the expression of which differed only in connection with the flower's position on the inflorescence. The third cluster comprised miRNAs the expression of which dropped as the flowers developed. The fourth cluster included miRNAs characterized by a heterogeneous trend in expression (Fig. 4).
The first cluster contained miRNAs belonging to the MIR166, MIR167, MIR319, MIR390 and MIR395 families. The first of these families includes Ll-miR177, which guides the cleavage of RADIALIS, a transcription factor from the MYB family that controls the asymmetric flower shape [98,99], as well as Ll-miR258 and Ll-miR265, which probably target the Homeobox-leucine zipper protein ATHB-15. In Arabidopsis line causes a decrease in the expression of all of the aforementioned HD-ZIP III target genes and results in SAM disruption, which may take the form of an enlarged apical meristem, a shortened carpel and female sterility [104]. miR166 has also been shown to control embryonic SAM development in rice and maize [105,106]. The MIR167 family members that accumulate in larger quantities during flower development are Ll-miR280, Ll-miR281 and Ll-miR285, targeting ARF6 and ARF8. In the thermo-sensitive wheat line, an increased tae-miR167 expression in spikelet tissue generated the phenotype of cold-induced male sterility [107]. Similarly in citrus, nineteen members of the MIR167 family, including the diverse isomiRs, have been identified, whereas most of them displayed a significant up-regulation in the anthers of male sterile mutants [42]. 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 exhibits 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 [108]. 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 [109]. An 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, as anthers in stage 3 of their development are already mature and split. Another miRNA showing a similar expression profile is Ll-miR9/miR390-5p. 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 plant development [110,111]. The expression level of miR390 derived from MIR390b reflects auxin concentration in organs, while the repression of ARF2, ARF3 and ARF4 by tasiR-ARF is important for lateral organ development, including roots and leaves [11,[112][113][114], and flower formation [115]. Down-regulation of ARF2 using RNA interference (dsRNAi) in Arabidopsis leads to delayed flowering, rosette leaf senescence, flower abscission and the opening of siliques [116].
Ll-miR118 and Ll-miR119, which target ATP sulfurylase 1, belong to the MIR395 family. In Arabidopsis, miR395 targets two gene families, ATP sulfurylases (ATPS) and sulfate transporter 2;1 (SULTR2;1), which are elements of the sulfate metabolism pathway [117]. ATPS regulates glutathione (GSH) synthesis and is an essential enzyme in the sulfur-assimilatory pathway [118,119]. In cotton, the miR395-APS1 module is engaged in drought and salt stress response [120]. In B. juncea [63], miR395 presumably down-regulates APS1 to cause an increase in GSH expression and suppress the oxidative stress. 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 [121]. Consequently, during yellow lupine flower development the demand for sulfur probably increases, and the plant activates mechanisms for its efficient uptake.
Within the third 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 (Fig. 4). Ll-miR99, Ll-miR100 and Ll-miR102 are identical to miR390-3p (the so-called passenger strand, former star strand) and, probably similarly to miR390-5p, are able to guide TAS3 cleavage. However, their expression showed an opposite trend to that of miR390-5p. A target gene of Ll-miR115/miR858 is MYB5, engaged in regulating seed coat development and trichome branching [122], although in Arabidopsis miR858 modulates another set of MYB genes, which are flavonoid-specific and control pathogen defense [123]. Another miRNA from the third 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 [124]. ECERIFERUM1 (CER1), the target gene 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 [125]. This cluster also included a miRNA that negatively regulates elements involved in miRNA and ta-siRNA functioning, namely Ll-miR247/miR168 targeting AGO1 mRNA [126]. 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. Transgenic A. thaliana plants overexpressing the gene encoding this protein suffered from severely reduced fertility due to the lack of anther dehiscence and an improper twisted path of wild-type pollen germination [127]. What is puzzling is that in yellow lupine the expression of miR408 was high in stage 1 lower flowers and stage 2 lower and upper flowers, where both the anthers and the transmitting tract were immature, only to decrease in the subsequent stages. It is possible that in lupine PC is suppressed in an alternative way, or miR408 expression actually rises in the anthers and the pistil, although in a very local manner, and this trend is unnoticeable due to the high background from other tissues where the expression profile of this miRNA may be opposite.
Especially so, since it has been shown that miR408 plays a role in other Arabidopsis organs, as well.
Transgenic Arabidopsis plants overexpressing 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 (GA). In addition, this plant line exhibited a higher concentration of copper in chloroplasts and an increased level of plastocyanin (PC) expression, resulting in a significantly higher photosynthetic rate [128].

miRNAs specific for particular developmental stage
Among the miRNAs identified in yellow lupine we spotted several that seemed to be crucial in particular stages of the plant's development ( Table 4, Additional File 1: Table S9). For example, the largest quantities of miR159 (Ll-miR452 and Ll-miR454) were accumulated in stages 2 and 3 of the plant's development. They targeted 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 early on in flower development. The accumulation of homologues of miR5168, miR1861, miR369-5p and 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 [129]. As confirmed by yellow lupine degradome sequencing, Ll-miR229/miR396-5p targets Growth-regulating factor 5 (GRF5) and GRF4 transcripts. In Arabidopsis, miR396 plays an essential role in leaf growth by down-regulating GRF, which in turn interacts with GRF-interacting factor 1 (GIF1) to take part in cell proliferation [130]. In this plant, 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 [131].
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 [133]. The presence of these miRNAs in yellow lupine flowers suggests that they regulate cell proliferation mainly in the pistil and the developing ovules. 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 [69].
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 (upregulated 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 [134]. In rice roots, the SCARECROW genes play a role in cellular asymmetric division [135]. In this plant, SCR2 expression is relatively high in flower buds and flowers, and after flowering rises in the leaves and roots [136]. 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 DEM, Ll-miRn10, shows differential expression in LF2 vs LF1 and UF2 vs UF1 comparisons, being up-regulated in stage 1 flowers. Its target gene is CIPK6 that is engaged in calcium-dependent CBL kinase activation, which in turn regulates potassium channel activity [137].
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 [138]. Through analyses of our degradomic data we did not find its target gene, and the psRNATarget hit was the putative disease resistance RPP13like 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. SGS3-and DCL2-encoded proteins are involved in sRNA biogenesis [139,140]. Importantly, novel miRNA identified in soybean Soy_25 displays a high sequence similarity to Ll-miRn1 and also targets SGS3, which indicates that this regulatory feedback loop for sRNA biogenesis is common for Fabaceae [141]. 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. Lower flowers in this developmental stage exhibited higher expression of miR396-3p, miR393-3p, miR858, miR390-3p, miR167-5p, and miR319. From the second stage until the end of their development, upper flowers accumulated more miR319, miR394, miR160, and miR393 (Fig. 4, Table 5). The elevated expression levels of these miRNAs suggests 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, and through decreasing the number of transcription factors belonging to the TCP family, probably caused different cell proliferation profiles in comparisons of flowers collected from the upper and lower whorls. miR393 participates in regulating several developmental processes in the leaves [142], and contributes to root architecture development [143], pathogen defense [144][145][146] and plant response to nitrogen limitation or changes in salinity [13,143,147,148]. miR393 regulates the accumulation of transcripts encoding auxin receptors belonging to the TAAR family, which, in the presence of auxin, guide proteasomal degradation of AUX/IAA repressors, thus enabling transcriptional regulation by ARF. Changes in receptor abundance affect the sensitivity of the given tissue to auxin and this is how this molecule influences plant development [149].
In A. thaliana, miR160 directly controls three ARF genes, namely: ARF10, ARF16 and ARF17 [150]. It is involved in many processes in plants, e.g. flower identity specification, leaf development or fruit formation [150,151]. In tomato, sly-miR160 is abundant in ovaries, and changes in its expression affect plant fertility [152]. Down-regulation of sly-miR160 by introducing a target mimic caused improper ovary patterning and thinning of the placenta already prior to anthesis. Changes in gynoecium shape were followed by changes in fruit morphology [152]. Moreover, the ectopic expression of SlARF10A caused by a mutation in the target site for miR160 (mSlARF10A) resulted in a similar phenotype, with an additional strictly limited number of seeds that failed to germinate [153].
In view of these facts, the higher expression of miR160 in lupine lower flowers early on in their development means that a slightly different organization of the gynoecia may be one of the crucial determinants of flower fate.
Flowers collected from the lower whorls displayed higher accumulations of miR5490, miR5794, miR1861, miR396-5p, miR395, miR166, and miR159-3p (Table 5). miR1861 and miR396 were recognized as positive cell proliferation and development regulators [154][155][156]. In rice, for example, miR1861 exhibited differential expression during grain filling [157], and its expression was higher in superior grains in comparison to inferior ones [154]. This is consistent with our results, and 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.

sRNAs are Involved in Flower Abscission in L. luteus
Little is known about sRNA engagement in flower abscission. Research into the involvement of miRNAs in this process has been already carried out in cotton [158], tomato [152,159], and sugarcane [160]. Our investigation provided a more comprehensive view of this process and allowed us to explore cues for abscission present long before AZ formation.
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 novel showed differential expression in abscission tissues [158], and these miRNAs have no homologues in our dataset.
In tomato, sly-miR160 regulates three auxin-mediated developmental processes, namely ovary patterning, floral organ abscission and lateral organ lamina outgrowth [152]. 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 [152]. This was consistent with the higher accumulation of miR160 in upper flowers designated to fall off. As this miRNA showed no differential expression in flower pedicels, it probably does not play a role in the executory module of abscission itself, but is rather 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 early and late stages of the process, additionally accelerated or not by ethylene or control treatment [161]. The study showed that in tomato pedicels, when compared to the early stages of abscission, in comparison to the early stage 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 [161]. 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 miR396 was lower, and the levels of miR167 and miR166 were higher, as well (Table 6).
It has been proven for sugarcane that miR156, miR319-3p, miR396 and miR408 take part in leaf abscission [160]. In their study, both mature (5p) and passenger (3p) miRNAs from MIR167 family were up-regulated in leaf abscission sugarcane plants (LASP), comparing to leaf packaging sugarcane plants (LPSP) (which corresponds to the FPAB vs. FPNAB comparison in our study). A total of four pairs of mature and passenger miRNAs were differentially expressed in this comparison, and they were members of the MIR156, MIR167 and MIR393 families [160]. Apart from the above-mentioned sugarcane miRNAs, miR396, miR395 and miR171 were up-regulated in this library comparison, as well. Four miRNAs, namely miR164, 166, 159, 319-3, showed an opposite trend [160]. In our study, both mature and passenger members of the MIR167 family were leaders among DEMs, 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.
The situation is different for miR396. In contrast to sugarcane, passenger miR396-3p in yellow lupine is mainly accumulated in the pedicels with an inactive AZ ( Table 6). It is also abundant in flowers from the lower whorls (Table 5), which suggests that in yellow lupine miR396-3p and miR396-5p (present in LF3) are engaged in maintaining flowers on the plant. This is probably related to the fact that this miRNA* may have different target genes in different organisms due to the higher rates of change characteristic of passenger miRNAs [162]. Thus, it is possible that in L. luteus miR396-3p targets other genes than it does in sugarcane.

Auxin-related sRNAs Involved in Regulating Flower Development and Abscission
Recent studies have shown that sRNA activity is associated with hormonal regulation of plant development through influencing the spatio-temporal localization of the hormone response pathway [163].
The auxin signal transduction pathway mainly consists of the three elements discussed below. Auxin is perceived by members of the TAAR family, namely TIR1 and its homologues AUXIN SIGNALING F-BOX PROTEIN (AFB), F-box proteins comprised in the SCF complex (SKP1-Cullin -F-box-RBX1).
Downstream these receptors there are AUX/IAA repressor proteins and ARF transcription factors [151,164,165]. The expression of TAAR receptors is regulated by miR393 and secondary ta-siRNA derived from their own transcripts [13]. miR167 and miR160, by affecting the ARF6, ARF8 [88] and ARF10, ARF16, ARF17 [166] transcript accumulation, respectively, are able to modulate the expression of GH3-like [167] and, by this, influence auxin-regulated plant development. It has been proven that the expression of ARF2, together with ARF3 and ARF4, is regulated by the ta-siRNA/miR390 module [168].
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 [169].
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 Ll-miR224, a member of the MIR393 family, which guides cleavage of TIR1 mRNA , was identified as differential and unique for the upper part of lupine inflorescence, which makes it a molecular marker of flowers fated to fall off.
Three members of the MIR160 family displayed differential expression, namely Ll-miR332 and Ll-miR329 in UF2 vs LF2, and Ll-miR332 and Ll-miR333 in UF3 vs LF3, all being up-regulated. As already mentioned, these miRNAs may promote flower abscission.
Strikingly, our sRNA-seq analyses showed differential expression of three miRNAs identical with miR390-3p deposited in miRbase for other plant species. These miRNAs (Ll-miR99, Ll-miR100 and Ll-miR102) were down-regulated in all of the following comparisons: between upper vs lower flowers (UF1 vs LF1, UF4 vs LF4), between pedicels (FPAB vs FPNAB), and between flowers in different stages of development (LF2 vs LF1, UF3 vs UF2). There was one exception, namely that of Ll-miR9, which was identical to, inter alia, aly-miR390a-5p or ath-miR390a-5p (miRBase), which was up-regulated in the LF3 vs LF2 and UF4 vs UF3 comparisons. The differential expression and functioning of passenger miRNAs has already been described. The research carried out by Xie and Zhang 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. These miRNA*s were up-regulated in seedlings, but down-regulated in other growth stages [170]. Therefore, miRNA*s can be specifically expressed in various tissues to maintain the steady state of the organism. Our degradome analysis showed for yellow lupine 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, and further research is required in order to identify its accumulation and function in the organs concerned.
ARF2, ARF3 and ARF4 are possibly down-regulated in the processing that is guided by Ll-siR249 and Ll-siR308, which were shown to be differential in the LF2 vs LF1 and UF2 vs UF1 comparisons (Table   4), and which are identical to tasiR-ARFs in many plant species according to the tasiRNAdb database [171] (data not shown). These tasi-ARFs probably originate from TAS3 transcript (TRINITY_DN55534_c4_g1) containing two binding sites for miR390 (data not shown). Ll-miR9/miR390, differentially expressed in LF3 vs LF2 and UF4 vs UF3, and, surprisingly, also Ll-siR240, guide cleavage of another TAS3 mRNA (TRINITY_DN54998_c6_g5_i2) ( Figure S8) which contains only one target site for miR390 (data not shown). 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 (Additional File 2: Figure S8). 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 (Fig. 7). Ll-siR240 accumulated proportionally to TAS3, which means that it was least expressed in the pedicels, while in flowers its expression increased with time (Additional File 2: Table S21). The identified target transcripts belonging to the ARF2, ARF3 and ARF4 families showed organ-specific differential expression, which indicated their strongly local regulation by siRNA (Additional File 2: Table S21). 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.

Function of the Potential miRNA Targets
Our functional analysis of the putative targets identified for miRNAs in lupine indicated their engagement in regulating a number of metabolic -especially 'carbohydrate metabolism' and 'nucleotide metabolism' -pathways (Additional File 5: Figure S4). 'Carbohydrate metabolism' was also one of the most enriched KEGG pathways in our previous Lupinus luteus transcriptome analysis [10], and its activation may be an indication of the rebuilding of cell walls or changes in nutrient supply. The next most numerous group of miRNA targets was categorized into the 'Genetic The MAPK pathway is involved in regulating several processes, such as biotic and abiotic stress response (reviewed in [172,173]), and associated with the functioning of hormones such as ethylene [174] and abscisic acid, engaged in organ abscission and other processes (reviewed in [175,176]).
The MAPK cascade is also an element of the positive feedback loop amplifying the abscission signal [177]. Auxin is not the only hormone whose signal transduction pathway is regulated by sRNAs in yellow lupine. Our 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 yellow lupine, as well (Fig. 9). These data show how important sRNAs are for growth and development regulation.
We also conducted a GO enrichment analysis (Fig. 8, Additional File 2: Table S13). Within the 'Cellular component', 83.3 % of the target genes were categorized as 'cell', with only 4.7 % categorized as 'extracellular region' and 3.9 % as 'cell junction'. This means that miRNAs impacted processes mainly occurring within the cell. As for the 'Molecular function', most of the targets fell into the 'binding' (67.5 %) and 'catalytic activity' (52.2 %) categories. Thus, in yellow lupine flowers, miRNAs were found to most frequently influence molecular interactions and biochemical reaction rates. Within the 'Biological Process' group, most of the targets were categorized as 'cellular' (72.1 %) and 'metabolic process' (63.6 %). What is most interesting is that quite a considerable number of the target genes fell within the 'response to stimulus' and 'signaling' categories, which means that miRNAs modulated the way the plant adapted to environmental stimuli. Surprisingly few genes were categorized into the 'growth' and 'cell proliferation' categories, which indicated that these processes were probably regulated by other factors.
An in-depth analysis of KEGG and GO terms concerning plant hormones (Fig. 8b) showed that most of the miRNAs identified in yellow lupine modulated more than one hormone signaling pathway. For example, in a GO analysis performed exclusively for the targets exhibiting the 'Flower development' term ( Fig. 8c), 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 [178], but also Cullin-3A, an element of the ubiquitination complex [179]. 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 [180].
Similarly, a KEGG analysis for the MIR166 family showed that it was involved in the auxin, cytokinin, and brassinosteroid signal transduction pathways (Fig. 9).
Our GO analysis additionally showed for yellow lupine flowers that miRNAs were responsible for guiding the processing of target genes simultaneously involved in multiple processes associated with flower development (Fig. 8c). For example, AP2 is involved in the specification of floral organ identit [181], and ovule [182] and seed development [183,184], and in our study it was targeted by ten lupine miRNAs: 8 belonging to the MIR172 family (Ll-miR25, Ll-miR26, Ll-miR27, Ll-miR28, Ll-miR29, Ll-miR30, Ll-miR140, Ll-miR171), one to the MIR166 family (Ll-miR249), and Ll-miR24. On the other hand, seven of these miRNAs additionally targeted a negative flower development regulator, chromodomain-containing protein (LHP1), which is a structural component of heterochromatin involved in repressing floral homeotic genes and FLT that regulates flowering time [185,186]. 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 dependent on a complex network of interconnected factors.

Conclusions
In this paper we show our results of next generation sequencing of sRNA, transcriptome and degradome libraries, which allowed us to identify miRNAs, siRNAs and their target genes probably Declarations

Ethics approval and consent to participate
Not applicable. This is to confirm that no specific permits were needed for the described experiments, and this study did not involve any endangered or protected species.

Consent for publication
Not applicable.

Availability of data and material
The data supporting the conclusions of this article are within the paper and its additional files. All sequencing reads are deposited in the National Center for Biotechnology Information under the BioProject number PRJNA419564 with the Sequence Read Archive (SRA) study accession SUB3230840.