MicroRNAs in Daphnia magna identified and characterized by deep sequencing, genome mapping and manual curation

MicroRNAs (miRNAs) are small non-coding RNAs that function in RNA silencing and post-transcriptional regulation of gene expression in most organisms. The water flea, Daphnia magna is a key model to study phenotypic, physiological and genomic responses to environmental cues and miRNAs can potentially mediate these responses. By using deep sequencing, genome mapping and manual curations, we have characterised the miRNAome of D. magna. We identified 66 conserved miRNAs and 13 novel miRNAs; all of these were found in the three studied life stages of D. magna (juveniles, subadults, adults), but with variation in expression levels between stages. Forty-one of the miRNAs were clustered into 13 genome clusters also present in the D. pulex genome. Most miRNAs contained sequence variants (isomiRs). The highest expressed isomiRs were 3′ template variants with one nucleotide deletion or 3′ non-template variants with addition of A or U at the 3′ end. We also identified offset RNAs (moRs) and loop RNAs (loRs). Our work extends the base for further work on all species (miRNA, isomiRs, moRNAs, loRNAs) of the miRNAome of Daphnia as biomarkers in response to chemical substances and environment cues, and underline age dependency.

targeting properties 12 . Previous deep sequencing projects have also revealed reads from miRNA offset RNAs (called moRs) of the pri-miRNA transcript adjacent to the 5′ end of mature miRNA of the 5p arm and to the 3′ end of mature miRNAs of the 3p arm 13,14 . Reads have also been identified from the loop of pre-miRNAs 15,16 and these loop-derived sequences are named loRs 17 . The functions of moRs and loRs are still unclear.
MiRNAs are prevalent in insects 18 and are involved in regulation of several biological processes 19 , as for example insect physiology 20 like metamorphosis 21 or immune response 22 and even targeting viruses 23 .
Daphnia magna and D. pulex are the aquatic key invertebrate model organisms in ecotoxicology, ecology and evolution 24 . Responses to different triggers and stressors, including environmental pollutants, pesticides and GMOs, have been studied on its life-history and physiology as well as on a genetic and epigenetic levels [25][26][27][28][29] . Recently, some attention has been given to describe miRNAs in the genus Daphnia. In D. pulex, 229 microRNAs were identified in silico and verified by qPCR 30 and microarrays 31 while in D. magna 205 putative mature miRNAs were identified based on small RNA sequencing and mapping against databases 32 . Eighteen of these were confirmed by secondary structure analysis of miRNA hairpins. Recently, Hearn et al. 33 reported 33 miRNAs in D. magna by deep sequencing and genome mapping. The authors also investigated whether transgenerational effect could be mediated via miR-NAs and showed that miRNAs played a role in maternal provisioning rather than longer-term transgenerational responses. Indeed, the expression profiles of miRNA are a promising step to further examine organism responses to stressors, in particular on specific genes, genetic pathways or networks 34 . Thus, miRNAs may extend their use as important biomarkers 31,33,35 . However, it has recently been reported that annotating miRNAs without mapping sequenced reads to the same species genome may largely overestimate the presence of miRNAs 36 . For instance, Yalle et al. 36 found that insects contained 65 conserved families with very low variation.
In this study, we characterize the miRNAome from three different life stages (juveniles, subadults, adults) in D. magna by a four step process: 1) deep sequencing of small RNAs, 2) alignment of the trimmed, clean reads to the D. magna genome, 3) identification of putative miRNA by aligning the mapped reads to miRBase and 4) manual inspection of sequence and reads in the genome at the position of the putative miRNA. With this approach we aimed at improving the coverage and verification of authentic presence of this important class of gene regulating nucleic acids, building on previous miRNA studies in D. magna 32,33 . We identify conserved and novel miRNAs, miRNA clusters, isomiRs, offset RNAs (moRs) and loop RNAs (loRs), and compare these data between juveniles, subadults and adults. Our work presents novel miRNA discoveries and characterisation on age-specific miRNAs in Daphnia which can be used for further verification and investigation.

Results and Discussion
Analyses of the small RNA sequencing data. The raw reads from deep sequencing of three small RNA libraries, one each from the three life stages of pooled D. magna juveniles, subadults and adults, were processed into trimmed clean reads of 17-30 nucleotides ( Table 1).
The trimmed clean reads were mapped to the genome of D. magna and about 96% of the clean reads of 17-30 nucleotides from all three libraries aligned to the genome sequence ( Table 1).
As a first step to identify miRNAs, reads mapping to the D. magna genome were aligned to miRNAs in miR-Base release 21 37 and dpu-miRNAs from D. pulex 30 . The percentage of mapped reads that aligned to miRNA sequences was similar for juveniles (35.4%) and adults (36.7%) but only 15.2% for subadults ( Table 1). The mapped read counts were normalized and distributed on read lengths (Fig. 1). For all three life-stages, reads of 22 nucleotides were most abundant, but the number of reads of 22 nucleotides from subadults was less than one third of the reads of 22 nucleotides from juveniles and adults. The next highest counts were found for reads of 18 nucleotides (juveniles, adults) and 19 nucleotides (subadults). Overall, juveniles and adults had more read counts of lengths below 24 nucleotides, while subadults had more reads with lengths between 25-30 nucleotides (Fig. 1).
Identification of conserved miRNAs in D. magna. The reads that mapped to the D. magna genome aligned to several hundred putative miRNA species in miR Base and dpu-miRNAs. To verify that the putative miRNAs identified in miRBase and dpu-miRNAs were real D. magna miRNAs, we mapped the mature miRNA sequence of each putative miRNAs back to the D. magna genome. The positions in the genome, where each of these mature miRNA sequences aligned, were then manually scanned for sequence and mapped reads in the same regions. Through these analyses we identified the presence of reads that could be annotated to 66 conserved microRNAs in miRBase and dpu-miRNAs, representing 51 miRNA-families (Table 2). This number of miRNA families is similar to a recent re-evaluation of the number of miRNAs in insects, which identified 65 conserved miRNA families 36 . Hearn et al. 33 Table S1 for details on the overlap of miRNAs in previous studies 32,33 and this study). In addition 37 miRNAs of the 229 found in D. pulex 30,31 were identified in D. magna in this work (Supplementary2 Table S1). The number of normalized reads aligned to the miRNA hairpins (pre-miRNA + adjacent sequences) ranged from a few reads (miR-153) to more than 3 million (miR-1) ( Table 2). When comparing the miRNA expressions in the three life stages, the adults had highest expression of 28 miRNA hairpins, followed by juveniles with highest expression of 24. Subadults had highest expression of only 14 miRNA hairpins (Table 2).
Notably, the sequenced reads from the 66 miRNA hairpins were from whole animals, thus lacking further details on how different body parts, organs, etc., differ in expression. However, the three age groups expressed all the same 66 species of miRNA, but with some proportional differences between the three different life stages (Fig. 3, Table 2). The heat map (Fig. 3) also depicts some larger groups of miRNAs where the expression pattern of miRNAs in each group was similar across the three life stages. Cluster analysis of the total miRNA read counts of each life stages indicated that juveniles and subadults were the two most similar whereas the adults deviated from these two (Fig. 3).
Conserved mature miRNAs. Mature guiding strands of miRNAs can be expressed from the 5p arm, 3p arm or both in the pre-miRNAs by changing from one to the other arm in different tissue in the same organism 7,38 . Usually the pre-miRNA arm containing the mature miRNA sequence with the highest read counts is treated as the mature guide miRNA, leaving the sequence of the other strand as passenger (or star) miRNA 5,39 . Sometimes miRNA sequences form both pre-miRNA arms are expressed as mature guide miRNAs and may co-mature due to different expression in different cells (so called arm switching) 38 . In this study, with the primary aim to discover new, and strengthen the description and verification of previously described miRNA in D. magna, we did a single sequencing of small RNAs from each of three life stages (relatively large, pooled batches of whole-body juveniles, subadults and adults).
Based on the counts of read sequences and folding of pre-miRNAs, we could identify the mature miRNA sequences from both arms of the vast majority of the conserved miRNA hairpins (Supplementary1 Table S2, Supplementary1 Fig. S1). It should be noted that for the identical mature miRNAs, the reads were equally divided between miR-2a-1-3p and miR-2a-2-3p, and between miR-87-1-3p and miR-87-2-3p due to the impossibility of knowing which miRNA the reads mapped to.
The majority of the mature guide miRNA sequences in D. magna were expressed from the 3p arm (44 miRNA) of the pre-miRNA (Supplementary1 Table S2). Almost all of the mature miRNA sequences had more than 50% of the total read counts mapped to its arm of the miRNA hairpin (Supplementary1 Table S3), but with variation from 0% (miR-279b-5p) to 94.6% (miR-2944-5p) (Supplementary2 Table S1). In addition to the read counts of each mature miRNA, additional read counts to each arm were mainly due to miRNA isoforms (isomiRs) (Supplementary2 Table S1. See also below). Although for most of the miRNAs the mature guide miRNA had much more reads than the passenger miRNA, several of the passengers also exhibited rather high read counts (Supplementary1 Table S2). This could indicate that some of these passenger miRNAs may be active regulatory species. However, for a few miRNAs the read counts of mature miRNA sequences from both 5p and 3p arms were similar. To indicate co-expression of mature miRNA from both arms, we used the ratio of <2 between the read counts of the most expressed arm divided by the number of reads from the less expressed arm 36 . We found that some of the mature miRNAs displayed co-expression of both arms in different combinations of juveniles, subadults and adults as shown in Table 3. Furthermore, some mature miRNAs had undergone arm-switching of the most expressed miRNA from one arm to the other in the different life stages of D magna (Table 3, Supplementary1  Table S4). The co-expression and arm-switching for miR-2a-1 and miR-2a-2 are highly uncertain since the expression levels of mature miRNAs were only estimated to be equal (see above). We also found that some miR-NAs displayed shift in expression between the mature miRNA and its isomiR in the three life stages (see below).
The mature guide miRNA sequences of D. magna were identical or similar to most of the conserved mature guide miRNAs in miR Base, release 21 37 . The mature guide miRNAs (highest expression) are shown in Supplementary1 Table S2, and their putative differential expression between the three developmental stages is displayed in the heat map in Supplementary1 Fig. S2. The clustering of the mature miRNA profiles (Supplementary1 Fig. S2) indicate a closer relationship between juveniles and subadults, as was seen for the miRNA hairpin profiles (Fig. 3). Again, we observed some groups of miRNAs where the pattern of expression in each group was similar across juveniles, subadults and adults. It is interesting to observe that miR-1, let-7, miR-184, miR-275 and miR-276 which are among the highest expressed miRNAs in D. magna, are also among the highest expressed miRNAs in the closely related, white shrimp Litopenaeus vannamei 40 . In contrast, only miR-1 and miR-276 (and then miR-276-5p and not miR-276-3p which is the highest in our result) were among the five highest expressed D. magna miRNAs analysed by Hearn et al. 33 . Many of the conserved miRNAs found in D. magna have been revealed to have diverse and important biological functions in insects 18,19,21,36 . For instance, let-7 is the most known miRNA that regulates developmental timing in animals. Mir-1, predominantly expressed in muscle tissue, in insects, is responsible for maintaining muscles integrity and mir-184 is involved in regulation of oogenesis 19 .    41,42 . In mammalian species, most of the miRNA genes (about 2/3) are clustered within 50 kb 43 . Based on this knowledge we first performed a query for putative novel miRNAs in D. magna by manual scanning of scaffolds and contigs that contained the conserved miRNAs. Using the criteria for annotation of novel miRNA 37 described in the Methods section we identified 13 putative novel miRNAs (Supplementary1 Table S5) of which six were identical with non-annotated precursor miRNAs in D.magna reported by Hearn et al. 33 (Supplementary2 Table S1). Five (miR-nov-1, miR-nov-3, miR-nov-4, miR-nov-10 and miR-nov-12) fulfilled the criteria of novel miRNAs in D. magna ( Table 4). The secondary structures of these five miRNA hairpins with their folding free energies (∆G) are shown in Fig. 5.
In addition, eight other pre-miRNAs fulfilled several of the criteria for the annotation as novel miRNAs 37 . However, for most of them mature miRNA reads were only detected from the putative miRNA sequence on one of the hairpin arms (Supplementary1 Table S5, Supplementary1 Fig. S3).
Mapping the D. pulex genome with the mature sequence of the novel D. magna miRNAs showed that several of these miRNA sequences were also present in D. pulex (Supplementary1 Table S6).
Ikeda et al. 44 reported that 21 of the novel miRNAs in Triops cacriformis (tadpole shrimp) showed more than 80% sequence similarities with the genomic sequence of these miRNAs in D. pulex. In our data we found only four of the T. cancriformis novel miRNAs as expressed miRNAs in D. magna. It is interesting that the three novel D. magna miRNAs: miR-nov-7-3p, miR-nov-8-3p and miR-nov-9-3p have sequences identical to tca-miR-n504-3p, tca-mir-n512-3p and tca-miR-n503-3p from T. cancriformis. Moreover, miR-nov-13 and tca-miR-3477 have identical seed sequences, but with five mismatches outside the seed region.
Genomic clustering of miRNAs. MiRNAs are often found as clusters in metazoan genomes. In flies and worms about 30% of all miRNAs are clustered as two or more miRNAs within 10 kb, while ~45% of conserved miRNAs are found in clusters in these species (and mammals) 42 . Moreover, many of these miRNAs clustered www.nature.com/scientificreports www.nature.com/scientificreports/ within 10 kb, and especially when clustered within 1 kb, are co-expressed [45][46][47][48] . In D. melanogaster, miRNAs separated by less than 1 kb seem to be highly co-expressed in different tissues 49 .
Identification of clustered miRNAs (both conserved and novel) within 10 kb of D. magna genome and comparison with clustering of the same miRNAs in the D. pulex genome, is shown in Supplementary1 Table S7. Both species contain the same 13 clusters with the same order of the miRNAs in each cluster. Of the 66 conserved miR-NAs in D. magna, 31 (47%) are grouped into 11 clusters. Furthermore, miR-nov-7, miR-8, and miR-nov-9, which A d u l t J u v e n i l e S u b a d u lt www.nature.com/scientificreports www.nature.com/scientificreports/ are clustered in D. magna, have their identical miRNAs tca-miR-n504, tca-miR-n512 and tca-miR-504 clustered in the same orientation in T. cancriformis 44 . We also note that miR-283, miR-nov-13, miR-12 build a cluster in D. pulex, but these miRNAs are on different scaffolds in D. magna. However, scaffold sc1551 (miRNA-283) may be connected to scaffold sc2703 (miR-nov-13, miR-12) since the Thyroid receptor-interacting protein (http:// arthropods.eugenes.org:8091/gbrowse/cgi-bin/gbrowse/daphnia_magna7/) is encoded across these two scaffolds. Thus, these three miRNAs seem to be within a 10 kb distance also in D. magna. Another support for this cluster is that similar clusters are documented in other organisms. For example, the tca-miR-12/miR-3477 (same seed as miR-nov-13)/miR-283 cluster is within 1000 bp in T. cancriformis 44 , the miR12/miR-1889/miR-283 in Anopheles gambiae 50 , and similar clusters in several other arthropods (see miRBase).
Several miRNA clusters are conserved among animals 42 . In mammals, clustered miRNAs may function as a unit in regulation of biological processes 51   www.nature.com/scientificreports www.nature.com/scientificreports/ clusters in other arthropods (miRBase release 21 37 ). Particularly many of the D. magna clusters are conserved in T. cancriformis (10 clusters) 44 , Tribolium castaneum (9 clusters) 52 and Apis mellifera (9 clusters) (miRBase release 21 37 ). Among the D. magna clusters with the highest expressed miRNAs are the miR-100/let-7/miR-125 cluster which have been found conserved in all metazoans 53 and is involved in regulating the larvae-to-adult transition as well as neural development 54 . The miR-2/miR-13/miR-71 cluster is conserved in many invertebrates 55 and was recently reported to be involved in resistance to the insecticide deltamethrin in mosquitos 56 .

miRNA isoforms (isomiRs) in D. magna.
Deep sequencing of small RNAs has shown that most miRNAs are expressed both as mature sequences and as variants (isomiRs) 8,9,57 . We have analysed the read sequences that mapped with more than 10 read counts to miRNA hairpins. The read sequences have been divided into seven groups: i) mature miRNAs; ii) 3′ template isomiRs containing reads with the mature 5′ end and deletions or additions at the 3′ end complementary to the genome sequence; iii) 3′ non-template isomiRs containing reads with mature 5′ end and 3′ additions not complementary to the genome sequence; iv) isomiR reads not part of groups ii) or iii) and with the mature 5′ end but with single nucleotide polymorphisms (SNPs) not altering the mature seed sequence; v) 5′ template isomiRs with reads containing the mature 3′ end and deletions or additions at the 5′ end complementary to the genome sequence; vi) 5′ non-template isomiRs containing reads with intact mature 3′ end and additions at the 5′ end not complementary to the genome sequence; and vii) isomiR reads not part of group v) or vi) and with SNPs altering the seed sequence 10 . Although the majority of D. magna read counts were mature miRNA sequences (conserved and novel), expressed in all the three developmental stages (juveniles, subadults, adults), all miRNA hairpin arms (5p or 3p) mapped with more than 10 reads, contained isomiRs (Supplementary2 Table S1). For each miRNA calculation of the ratio of read counts for each of the seven isomiR groups, divided on the total read count for the miRNA arm, showed large variation both for each miRNA and each isomiR group (Supplementary2 Table S1). Moreover, determining the mean ratio for each    www.nature.com/scientificreports www.nature.com/scientificreports/ isomiR group (across all miRNAs) revealed large variation also in the sizes of each isomiR group. Generally, the highest expressed isomiRs were the 3′-templated variants (mean 20.7% of the total read counts for 5p+3p arms). The mean level of expression of the other variants was as following: 3′ non-template isomiRs, 3.7%; 5′ template isomiRs, 4.2%; isomiR reads with SNPs changing seed sequence, 3.5%; isomiR reads with SNPs not altering seed, 2.2% and 5′ non-template isomiRs, 0.08%. (Fig. 6, Supplementary1 Table S8). Supplementary1 Figure S4a & b) shows mean level of expression for each arm. Although some of these variants displaying one or a few reads may be sequencing errors, most isomiR sequences have several read counts and are present in all three independent sequencing samples (juvenile, subadult, adult stages) supporting that they are real isomiR sequences. www.nature.com/scientificreports www.nature.com/scientificreports/ Interestingly, the proportional relationship between the different isomiR variants in D. magna is similar to that reported for the feline (cat family) miRNAome 58 , suggesting that the roles and functions of mature miRNAs and their isomiRs may be conserved by evolution, even across distantly related organisms. A recent review suggests that isomiRs may function as an additional repression tool beside the mature miRNA 59 . For some miRNAs we found 3′ template variants that had a similar level of expression as their mature miRNAs: miR-1175-5p (two variants), miR-12-5p (mainly one variant), miR-279b-3p (juveniles, mainly one variant), miR-71-5p (mainly three variants), miR-993-5p (mainly three variants), miR-998-3p (juveniles and subadults, one variant; adults, two variants), miR-iab-8-5p (juveniles and adults, two variants; subadults, three variants), miR-9b-5p (two variants), miR-31-5p (one variant), miR-34-5p (mainly three variants), miR-210-5p (one variant), miR-285-3p (mainly one variant), miR-2944-3p (one variant), miR-317-3p (mainly two variants), miR-92b-3p (one variant), miR-981-3p (one variant) and 5′ template variants: miR-210-3p, miR-305-3p, miR-9a-3p (summarised in Table 5). In some cases the level of expression between mature miRNA and one of its 3′ template or 5′ template isomiR variants were noticeably different between the three developmental stages, e.g miR-12-5p, miR-210, miR-285-3p, miR-31-5p miR-998-3p, miR-305-3p, miR-9a-3p (Table 5). MiR-71-5p was special in this regard since the 3′ template variants had from 55% to 120% more read counts than the mature sequence (Supplementary2 Table S1, Table 5). However, the 3′ template isomiRs reads are here mainly divided on three different sequences with one to three deletions at the end (Table 5).
Although 3′ template isomiR variants usually are the most common, very little is known about their functions. In a recent report Yu et al. 60 found that 3′ template isomiRs of miR-222 could play functional roles in human cell lines. We found that the highest expressed 3′ template variants were usually those with one deletion at the end (Supplementary 2 Table S1). However, some exceptions were also found where the 3′ end additions were non-template A or U (see below).
We observed a mean expression of 3.5% of reads with SNPs that alter the seed sequence and a mean expression of 2.2% of reads with SNPs that did not affect the seed sequence (Fig. 6, Supplementary 1 Table S8). Although some of these SNP variants may be sequencing errors, we find that several variants have the same mismatch in identical sequences with several read counts across the three independent sequencing samples (juvenile, subadult, adult) supporting that these are real isomiRs (Suplementary2 Table S1). The SNP variants with altered seed sequences were the highest expressed isomiRs of miR-190-3p, miR-279c-5p, and miR-750-5p (Suplementary2 Table S9b). Editing of miRNA may result in SNPs both altering the seed sequence and the sequence outside the seed region 65,66 .
Taken together our sequencing results show expression of a variety of isomiRs of miRNAs in D. magna. The observed proportions for all three life stages shows that the 3′ variants and variants with SNPs that do not alter the seed, have a mean percentage of 26.5% while the 5′ variants and variants with SNPs changing the seed have a mean percentage of 7.8% (Supplementary1 Table S8). The 3′ variants show the highest expression for most of the miRNAs while the 5′ variants are expressed to a lower degree. IsomiRs have been suggested as promising, additional biomarkers in diagnostic assays based on mature miRNA for studying illness such as for example cancers 67,68 . In D. magna isomiRs may be developed into important biomarkers in addition to the role of mature miRNAs for analysing biological changes.
Offset-and loop miRNAs. Analysis of reads adjacent and in the loop of pre-miRNAs revealed one or more reads which were recognised as moRs and/or loRs (Supplementary2 Table S10). Nine miRNAs display moRs and/ or loRS with more than 10 reads in at least one the three life stages of D. magna (Table 6). MiR-7 is interesting  Table 5. IsomiR variants with similar read counts as the mature miRNA ( * )putative mature passenger strand.
since it contains 5′ moR reads that add up to relatively high percentages of the mature miR-7 reads (juveniles: 27%, subadults: 14% and adults: 96%) (Supplementary2 Table S10). As reported previously 17 , we also found that moRs at the 5′ site (5p) of mature miRNAs are more frequent and with higher read counts than on the 3′ site (3p) (Supplementary2 Table S10). In the three investigated life stages the total expressions of loRs were at similar level as for 5′ moRs. However, the expression level of the sum of 5′moR + 3′moR was somewhat higher than the expression of loRs (Supplementary2 Table S10). Both moRs and loRs may play functional roles 15,69 .

Conclusions
Daphnia magna is an important model organism for ecotoxicology, ecology and evolution. Its remarkable ability to respond to changing environmental conditions, i.e. its plasticity, could potentially be mediated via age-specific miRNA expression. As such, miRNAs could constitute an important epigenetic mechanism in clonal organisms.
Age or life-stage in an animal represent an important organismal context and needs to be considered for future experimental designs whenever miRNA is of interest. This is also in agreement with previous work 33 .
The main purpose of this study was to characterise the miRNAome and thereby identify and validate the number of miRNAs in D. magna by deep sequencing, mapping of reads to its genome combined with manual curation of putative miRNAs. By this procedure we identified and further characterized 66 conserved and 13 novel miRNAs from three life stages of D. magna (the pre-miRNA sequences are shown in Supplementary3 Table 11). Of the 66 miRNAs we found, 42 were already reported by Hearn et al. 33 and 17 miRNA hairpins and 53 miRNA were already identified by Ünlü et al. 31 All the miRNAs we found were expressed in all the three life stages and displayed similar miRNA profiles. However, clustering of the three miRNA profiles indicated a closer relationship between juveniles and subadults. The mature guide miRNAs were generally expressed with the highest read counts (>50% of total read counts). The additional reads were sequence variants (isomiRs) which were expressed by almost all miRNAs. We identified isomiRs which were 3′-, 5′-and SNP-variants. Several miRNAs contained 3′ template-or 5′ template variants of isomiRs that had a similar expression levels as their mature miRNAs. Furthermore, we found that miRNAs with the highest expression levels could shift between the mature miRNA and one of its 3′ template -or 5′ template isomiR variants over the three life stages. Several miRNAs displayed reads adjacent (offset RNAs (moRs)) and in the loop (loRs) of the pre-miRNAs.
Our work extends the background for further work on using both mature miRNAs and their sequence variants (isomiRs) as biomarkers of stress in Daphnia and should also help to exclude age-related confounding in follow-up studies. For future studies we recommend: i) independent validation of miRNAs in other laboratories, where at least 20 million reads is necessary (since we see some differences in our data to previous reports), ii) linked analyses of mRNA and miRNA expression to reveal more of the functional roles of miRNA in general, iii) as complete and detailed as possible reporting of the atlas of miRNA expression in different organs and sub-systems in the organism, also at different stages of development, and iv) test individual miRNAs for their ability and sensitivity to function as biomarkers of different types.

Materials and Methods
Experimental setup. All individuals of D. magna used in the experiments were from the same clonal population. Animals were kept in M7 medium 70 , fed Desmodesmus subspicatus, ad libitum and never crowded. The juveniles (n = 120) that were sequenced were born within 24 h from the third clutch of a set of n = 25 mothers. The subadults (n = 40) were also born within 24 hours from a set of n = 25 mothers and reared for 6 days. We checked for visible eggs in the brood chamber of each individual and confirmed that eggs were not present. For the adults, we selected animals with visible eggs at day 9 (n = 30). During the course of the experiment the subadult and adult groups were fed Desmodesmus subspicatus green algae with 0.15 mg Carbon per animal per day.
Every third day, we transferred each animal to a new glass with fresh M7 medium. Temperature was held constant at 22 + /− 1.5 degrees C. Light regime was 16 hours of light and 8 hours of darkness. Throughout the experiment, the average pH of the medium was 7.8, oxygen saturation always >95% and average conductivity 575 µs cm −1 .  Table 6. moRs and loRs of miRNAs with more than ten reads in at least one life stage (not normalized read counts). (2019) 9:15945 | https://doi.org/10.1038/s41598-019-52387-z www.nature.com/scientificreports www.nature.com/scientificreports/ (Bertin instruments). The Trizol extraction of total RNA was according to the manufacturer's protocol with modifications involving prolonged precipitation and centrifugation steps to ensure recovery of the small RNA fraction. Purified total RNAs were dissolved in nuclease free water and stored at −70 °C. RNA quantity was determined with a NanoDrop spectrophotometer (NanoDrop Technologies). The quality of RNA was analysed with an Agilent 2100 Bioanalyzer using an Agilent RNA 600 Nano kit (Agilent Technologies). Analyses of sequenced reads. Three datasets containing de-multiplexed, clean reads (reads with adaptor clipped off and removed non-clipped reads, adaptor-only reads, N reads, reads < 10 nucleotides) were received from NSC as fastaq files. These clean reads were analysed on a CLC Genomics Workbench (Qiagen) using the following pipeline: 1) The clean reads were trimmed to 17-30 nucleotides ("Trimmed Clean reads"). 2) The trimmed reads were mapped to the D. magna genome (dmagna-v2.4-20100422-assembly (http://wfleabase.org/) 71 using "Map Reads to Reference" module and following mapping parameters: mismatch costs 2, insertion cost and deletion cost each 3, length fraction 0.85, similarity fraction 0.8, map randomly. Reads mapped to D. magna genome were then aligned to miRNAs in miRBase release 21 (http://www.mirbase.org) 37 and in-silico determined dpu-miRNAs 30 using the various "Mature length variants" and "Alignment settings" parameters. The DNA sequence of potential D. magna miRNAs recognised from miRBase and dpu-miRNAs alignments (matches) were mapped to the D. magna genome and manually inspection of sequence and reads in the genome at the position of the putative miRNA was used to identify and annotate miRNAs. Pre-miRNA sequences were manually identified and secondary structures and their folding free energy (ΔG) of miRNA hairpins were determined by using the "Predict Secondary Structure RNA" parameter in the CLC software.

Construction of small RNA library and
For examination of the presence of D. magna miRNAs in D. pulex, miRNA sequences were mapped to the D. pulex genome (http://www.ftp://ftp.ensemblgenomes.org/pub/release-34/metazoa/fasta/daphnia_pulex/dna/) 72 Annotating novel D. magna miRNAs was performed mainly by the criteria described previously 17,37,73 except that we also included putative novel miRNA with reads aligning only to one of the pre-miRNA arms.
MoRNAs were defined as reads that were immediately adjacent either to the 5′ end of mature miRNA of the 5′ arm and/or adjacent to the 3′ end of the mature miRNA of the 3′ arm 14 . LoRNAs were annotated as reads between mature miRNAs of the 5′ and 3′ arm 15 . For both moR and loR we also included reads that started outside the mature miRNA and ended inside the miRNA and which could not be regarded as isomiRs.
We used DESeq2 package (R package) to normalise the raw read counts as described previously 74 . Shortly, as input the DESeq2 package expects raw count data obtained from sequencing of miRNA. The DESeq2 model internally corrects for differences in library size (estimateSizeFactor function) and delivers normalized values. Normalized counts were accessed by "counts" function on recalculated counts object. We did not perform statistics to look for significant differences in miRNA expression between the groups due to our experimental design (one independent library per age class). Hierarchical clustering, presented as heat maps, were generated by Heatmapper 75 using average linkage and Pearson's distant methods.

Data availability
The data sets supporting the results of this article are included within the article and its supplementary information files. Fastaq files of clean raw reads are uploaded to ENA under our Bioproject accession number PRJEB29358.