Analysis of microRNA expression reveals convergent evolution of the molecular control of diapause in annual �sh

Background: Diapause is a condition of developmental arrest in anticipation of adverse environmental conditions present in a large number of diverse taxa. Diapause is a key adaptation that enabled the colonization of ephemeral habitats subject to the alternation of dry and wet seasons by annual killiﬁshes. Upon desiccation of the ponds, killiﬁsh embryos remain vital but quiescent in the clay, where they can survive months or even years. Diapause can occur at three different developmental stages, but Diapause II (DII), which occurs at mid-somitogenesis, is the primary point of developmental arrest. Physiologically, Diapause II is associated with the arrest of the cell cycle in G1 and deeply reduced oxygen consumption and protein synthesis. However, diapause is not obligatory, and some embryos can proceed to direct development, skipping one or more diapauses. The precise molecular mechanisms that regulate entry and exit from diapause are beginning to be investigated, but this knowledge


Introduction
Life-history is a vital force in the evolutionary shaping of embryonic development.Extreme examples are seasonal species where embryos undergo diapause, a suspension of development or dormancy state, as it is the case for the larvae of many insect species from a temperate climate.The seasonal life cycle has also evolved in a clade of teleost fishes (suborder Aplocheiloidei) known as annual killifishes.Annual fish inhabit ephemeral bodies of water and are adapted to the alternation of wet-and dry-season in Africa and South America.All adult fish die when their habitat dries out.The conservation of the species is ensured by desiccation-resistant eggs that enter in diapause and remain encased in the dry mud until the next rainy season.Annual life history is present in both South American and African killifishes and is the result of a series of independent events.Early studies suggested four events of loss and re-gain of an ancestral annual trait 1,2 , while more recent studies suggest that it repeatedly evolved, at least three times in Africa and three times in South America 3 .In either of the two scenarios, an annual clade always has a sister non-annual clade that is phylogenetically closer than any other annual clade.The early development of annual fish is conserved in the different lineages [3][4][5] and is characterized by three points of developmental arrest.The primary point of developmental arrest occurs after the formation of the embryonic axis and organogenesis, at mid-somitogenesis, and is called diapause II (DII).DII is a facultative stage: when embryos are incubated at high temperatures it can be skipped 3,[6][7][8][9] , but lower temperature, darkness or dehydration (all conditions occurring in natural habitats) induce DII [10][11][12] .The duration of DII is highly variable, and the embryos can remain in this stage for several months 6,13 or even years (A.C., unpublished observation).The physiological and molecular mechanisms of diapause were studied in detail in the South American species,

Fish raising
Fish were raised, fed, and bred as described in Dolfi et al. 2014 5 .The protocols of fish maintenance were carried out in accordance with all animal use practices approved by the Italian Ministry of Health (Number 96/2003a) and the local animal welfare committee of the University of Pisa.All experimental procedures were managed, following the prescription of the

RNA isolation
For each species, 20 eggs were pooled for each preparation.Total RNA was isolated using the miRNeasy Mini Kit QIAGEN according to the manufacturer's protocol.

Small RNA sequencing
For small RNA sequencing, small RNA cDNA libraries were prepared as follows: for each developmental stage, equal quantities (1 µg) of total RNA were submitted independently to Illuminas TruSeq small RNA sample preparation protocol (Illumina Inc., San Diego, USA).The library preparation was performed as follows: RNA was ligated with proprietary adapters to the 5' and 3' termini of the RNA.The adapter-ligated samples were used as templates for cDNA synthesis.The cDNA was amplified with 13 PCR cycles to produce sequencing libraries, introducing specific nucleotide codes for indexing libraries to allow multiplexed sequencing.The cDNAs were purified by 10 % Novex TBE polyacrylamide gel electrophoresis (Invitrogen) and eluted into 300 µl elution buffer (Illumina) for at least 2 hours at room temperature, to enrich for molecules containing inserts in the range of 18 --33 nt.The resulting gel slurry was passed through a Spin-X filter (IST Engeneering Inc. Milpitas, CA, USA) and precipitated by the addition of 20 µg glycogen, 30 µl of 3M NaOAc, pH5.2, and 975 µl of pre-chilled (-20 °C) ethanol.After washing with 70 % ethanol, the pellet was dried at 37 °C for 5-10 min and dissolved in 10 µl resuspension buffer (Illumina).The purified libraries were quantified on the Agilent DNA 1000 chip, diluted to 10 nM, and subjected to sequencing-by-synthesis on Illumina HiSeq 2000.

Aplocheylus lineatus
No. of samples Figure 1.Overview of the sequenced killifish embryos and the corresponding amount of samples.Fish marked in red are annual species, whereas all others are non-annual.For the annual species, the numbers of samples sequenced during diapause are given in brackets.Aplocheylus lineatus originates from India.

RNA-Seq analysis
Clipping the RA3 adapter of the TruSeq small RNA preparation kit (5'-TGGAATTCTCGGGTGCCAAGG) and trimming the read ends by a quality threshold of 20 was performed with PRINSEQ (v0.20.3) 40.Read quality was monitored by FastQC (v0.11.3; http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Mapping was done using Hisat2 (v2.0.4) 41 onto the N. furzeri reference genome 23 , for details see https://osf.io/49mn2/.Read counting was done on the improved version of our published miRNA annotation of N. furzeri 42 , counting reads for each mature miRNA of each pre-miRNA annotation separately.Analysis of differentially expressed miRNAs was performed by the Bioconductor package DESeq2 (v1.10.0) 43within three defined sample comparisons (see https://osf.io/epf6z/).Only mature miRNAs having a count of more than ten normalized reads were considered to be actively transcribed and used for subsequent analyzes.The resulting p-values were adjusted using the Benjamini and Hochberg's FDR approach 44 .We considered mature miRNAs with an adjusted p-value of < 0.05 as differentially expressed.

Enrichment analysis of miRNA targets
We used the hypergeometric test to calculate enrichment scores of miRNA targets within the published set of differentially expressed aging-related Nothobranchius furzeri genes 23 : where N is the total amount of known protein-coding genes in Nothobranchius furzeri, R the number of differentially expressed genes of the given set, n the number of protein-coding genes with predicted miRNA target sites, and r the number of differentially expressed genes with predicted miRNA target sites.Enrichment of individual miRNAs, in the set of aging-related genes, was calculated similarly, with N being the total number of protein-coding genes with predicted miRNA target sites and n the number of genes, containing a target site of the respective miRNA.The resulting p-values were adjusted using Benjamini and Hochberg's FDR approach 44 and were considered significant if p was less than 0.05.

Estimating miR-430 copy number from short genomic reads
Using raw Illumina short sequencing reads of 50 killifish species from previous studies [45][46][47] , we counted the number of reads containing the following conserved motifs of eight isoforms of miR-430 identified in N. furzeri (UAAGUGCUAUU-UGUUGGGGAAG, UAAGUGCUAACUGUUGGGGUAG, UAAGUGCUAACUGUUGGGGUAU, UAAGUGCUACAU-GUUGGGGCAG, CACCUCAAAGAAUCCACUGA, ACUCUGACAAAGGCACUGACU, UAAGUGCUACAUGUUG-GAGUCA, UAAGUGCUUCUCUUUGGGGUUG), allowing for a maximum of two mismatches to account for possible divergence between species using a custom C++ program (https://github.com/melop/mirseek/),which depends on the edlib library.Genome sizes for these species were taken from previous studies.The estimated miR-430 counts per haploid genome can be computed as: h (l/g) , where h is the number of sequence hits, l the total read length, and g the genome size.

Phylogenetic sampling
Annual killifish are divided into three main clades distributed in Africa (one West-and one East-of the Dahomey gap) and South America (Fig. 1).In each of these regions, the clade contains both annual-and non-annual genera.The position of the non-annual genus Aplocheilus is debated.Earlier analyses considered it basal to all annual species 1 , but it was later suggested to be nested between African-and South American species 3 .Therefore, Aplocheloidei offers a unique opportunity to study the parallel evolution of life-history adaptations (see also 5 ).We collected eggs in diapause II from annual species or in the corresponding developmental stage (mid-somitogenesis) from non-annual species of all three clades, in addition to Aplocheilus lineatus, and analyzed miRNA expression by miRNA-Seq.We used as reference an updated version of the (N.furzeri) miRNA catalog we published previously (see, STable 1).This catalog contains 754 loci containing 874 unique miRNA sequences, of which 498 belong to known conserved miRNA families and 376 to new, potentially killifish-specific miRNA families, i.e., miRNAs that were only identified in N. furzeri so far.All killifish miRNA sequences were annotated against this reference, allowing two mismatches within the pre-miRNA sequence.From South America, we analyzed miRNA expression in the annual species Austrofundulus leohoignei and Nematolebias papiliferous and the related non-annual species Rivulus cylindraceus.From Africa, West of Dahomey gap, we analyzed the annual species Callopanchax occidentalis and the non-annual species Epiplatys dageti monroviae.From Africa, East of Dahomey gap, we analyzed the annual species Nothobranchius furzeri (5 replicates), and the non-annual species Aphyosemion striatum (4 replicates).Finally, we incubated eggs from N. furzeri (3 replicates) at a higher temperature, a condition known to promote direct development i.e., diapause skipping 3,7 .

Overall expression activity of conserved and killifish-specific miRNAs
First, we investigated the total number of expressed miRNAs in all analyzed killifish species.As depicted in Fig. 2, the largest number of expressed miRNAs, either conserved or killifish-specific, could be found in Nothobranchius furzeri and the more closely related species Aphyosemion striatum.For these species, multiple replicates were analyzed and higher sensitivity in Nothobranchius furzeri and the more closely related species was expected, because the killifish miRNA annotation was based on the Nothobranchius furzeri reference genome 42 .Nevertheless, we observed the third-highest activity in miRNA genes in the more distantly related Austrofundulus leohoignei and the lowest number of active miRNAs in Callopanchax occidentalis, an African species.Most of these genes showed only low (10 -50 reads) or moderate expression (51 -1000 reads).Combing all samples, we detected expression activity for roughly 45 % of the predicted miRNAs genes described in our previous study 42 .Interestingly, some predicted Nothobranchius furzeri miRNA genes did not show active transcription in any of the Nothobranchius furzeri samples.But were expressed in the other investigated killifish species, particularly in those not undergoing diapause.In particular, 41 miRNAs (19 conserved and 22 killifish-specific) were expressed in A. striatum corresponding to 8 % of the total number of miRNAs expressed in this species.Out of the killifish-specific miRNAs not expressed in N. furzeri, seven are expressed in E. dageti, one in A. lineatus and seven in R. cylindraceus.This result indicates that some miRNA genes present in the common ancestor of killifish are preserved, but transcriptionally inhibited in the N. furzeri genome.Phylogenetic tree of the investigated killifish species and the number of their actively transcribed miRNA genes (red -miRNAs belonging a known miRNA family; green -miRNAs specific to killifish; blue -miRNAs not being transcribed in any of the Nothobranchius furzeri samples).Mature miRNAs having a count of more than ten reads were considered to be actively transcribed.

Expression of miRNAs can separate killifish based on their diapause status
Based on the global expression patterns of all miRNAs, we have clustered the analyzed samples using principal component analysis (see Figure 3) and observed a clear separation not only between the individual species, but also between annual and non-annual embryos.Replicates of Aphyosemion striatum and Nothobranchius furzeri form defined clusters, whereas the three diapause-skipped Nothobranchius furzeri embryos are only weakly clustered.Nevertheless, these results demonstrate that samples are divided by their physiological status and not by their phylogenetic relationships, indicating the existence of a convergent transcriptional program associated with diapause across species.

Conserved expression of developmental-related miRNAs
Next, we identified the most abundant miRNAs across all investigated samples.Interestingly, four conserved miRNA families (miR-10a/b/c, miR-92a, miR-181a and miR-430a/b/c) and two killifish-specific miRNA families (currently named miR-19337 and miR-19344) were always the most highly expressed in all the species.This observation may indicate a preminent role during the early development of the investigated fishes (see, Fig. 4 and https://osf.io/epf6z/).This finding is consistent with the report of Romney and Pobrasky 48 , who also report these four conserved miRNA families to be the highest-expressed in embryos of the South American annual killifish Austrofundulus limnaeus at mid-somitogenesis.Of particular interest is the high expression of the two killifish-specific miRNAs miR-19337 and miR-19344, because this indicates that their evolution is linked to the control of some aspects of embryonic killifish development.High expression of members of the miR-10 family is not unexpected as they are known regulators of several hox genes, which are essential transcription factors during embryonic development and particularly for the anterior-posterior patterning during somitogenesis 49 .Indeed, the expression of miR-10 sharply increases during somitogenesis in the annual killifish Austrofundulus limnaeus 48 .Additionally, it is also known that the expression of miR-10 has a positive effect on the ribosomal protein synthesis machinery 50 , which is vital for an accurate progression of embryonic development [51][52][53] .Members of the teleost-specific miR-430 family are also key mediators of early embryogenesis.In the initial phase of maternal to zygotic transcription they cause degradation of maternal transcripts.In addition, they regulate other early developmental processes, such as meso-endodermal fate specification, brain morphogenesis and development of primordial germ cells [54][55][56][57] .Usually, miR-430 pre-miRNAs are produced by transcription of exceptionally large genomic clusters that contain multiple copies of the same miRNA and whose organization differs between teleost species 42 .Recent work in Austrofundulus limnaeus has shown that expression of miR-430 is down-regulated during somitogenesis and it is higher in embryos that escape diapause and thus proceed to direct development, underlining its broad relevance for the embryogenesis of teleost for both conserved and lineage-specific patters 48 .Also, miR-181 is necessary for embryonic development, and the deletion of all its paralogs induces embryonic lethality 58,59 .Besides, miR-181 targets another hox gene, namely hox-a11 in mammals, and the homeobox transcription factor Prox1, therefore being a regulator of embryonic development hinting at a potential similar function in fish embryos 60,61 .
The miR-92 is a member of the miR-17/92 cluster and an important well-studied positive regulator of the cell cycle 62 and therefore cell proliferation during embryonic development.
For both killifish-specific miRNAs, target mRNAs are only predicted from sequence analysis in Nothobranchius furzeri 42 .Among these, TSC2 and FAM83D are both putative targets of miR-19344.TSC2 is of particular interest as it is part of the TSC complex, a major upstream negative regulator of the activity of the mTORC1 complex that represents a central regulatory hub that integrates nutrient sensing and growth-factor signaling to regulate cell growth, protein synthesis, and cell proliferation 63 .mTORC1 is also a known major regulator of aging and longevity in multiple organisms 64 .FAM83D, the other potential target of miR-19344, is involved in cell proliferation and motility 65 .
Relevant potential targets of miR-19337 are CECR2, being part of a protein complex that regulates neurulation, STARD13B a  The conserved miRNAs of the miR-10 and miR-430 families, miR-92a, miR-181a, and the two killifish-specific miR-19344 and miR-19337 showed the highest relative expression in all investigated samples.All of these miRNAs are known or predicted to play critical roles in the regulation of developmental processes in embryonic cells.For details on the miRNA expression levels, see https://osf.io/epf6z/.
known inhibitor of cell growth 66 and SRF (Serum response factor), which modulates the expression of many immediate-early genes and therefore is an essential key player in embryonic development 67 .
From these examples, we can deduce that the developmental stage of the sampled embryos has a strong impact onto the expression of their respective miRNomes and highly abundant miRNA are also those that regulate developmental or associated processes.Especially the dominating abundance of miR-10 and miR-430 genes in almost all samples indicates the specialized roles of both miRNAs as well as their evolutionarily conserved regulatory importance for teleost embryogenesis.

Differentially expressed miRNAs related to diapause regulation in killifish
To identify differentially expressed mature miRNAs (DEMs) linked to diapause regulation of annual species, we set up three different comparisons of the examined small RNA-Seq samples: (I) the "large interspecies comparison" contrasting all diapause samples from annual killifish against the non-annual fish-derived samples, (II) the "intraspecies comparison" contrasting the Nothobranchius furzeri samples that underwent diapause against those that skipped it and (III) the "selected interspecies comparison" comparing the Nothobranchius furzeri diapause samples against the samples of its sister non-annual taxon, Aphyosemion striatum.Mature miRNA may arise form either the 5'-or 3'-side of the stem-loop that originate from a single miRNA gene and these forms target different mRNAs.Therefore, we have examined the mature 5'-p and 3'-p miRNAs of single miRNA genes separately in the subsequent analyses.
The large interspecies comparison (comparison I) detected 242 DEMs (122 up-regulated and 120 down-regulated in the annual fishes compared with the non-annual ones) belonging to 156 distinct miRNA loci (see https://osf.io/c53g2/).Some of the most significantly changing miRNAs in the annual comparison were members of the miR-430 family, essential regulators of developmental processes, and belonging to the above mentioned highest expressed miRNA genes.Interestingly, we found that the particular miR-430a/c-3p mature miRNAs are down-regulated in all the annual species, whereas its mature complementary form miR-430a/c-5p is up-regulated in the same species (see, Fig. 5).This observation is striking, since the mature 5p transcripts are expressed at much lower levels as the corresponding 3p mature transcripts.This finding is a prominent example of a miRNA mature isoform switch that we describe in more detail in the next subsection.All the mature miRNAs of miR-430a and miR-430c act in developmental processes but target different mRNAs [54][55][56][57] .In addition to miR-430, we identified other DEMs with essential development-related functions, such as miR-19 being part of the miR-17/92 cluster and associated with endothelial cell differentiation in embryonic stem cells 68 and more generally with the positive control of the cell cycle.Its downregulation is therefore expected given the profound depression of cell cycle in diapausing embryos 69 .Similarly, miR-200 regulating the epithelial/mesenchymal transition during embryonic development 70 and miR-221, which is involved in cell proliferation and angiogenesis 71,72 .All of these miRNAs seem to be necessary for a the proper regulation of different developmental processes and may in part be responsible for the carefully orchestrated diapause regulation and embryogenesis in killifish.
The intraspecies comparison (comparison II) between the Nothobranchius furzeri samples in diapause and those that have skipped it due to high-temperature incubation, revealed 17 DEMs (10 up-and 7 down-regulated in the samples undergoing diapause) belonging to 14 distinct miRNA loci (see https://osf.io/kqx9y/).Observing only a few differentially regulated mature miRNAs between the diapause and non-diapause embryos of Nothobranchius furzeri can be due to the fact that not all embryos incubated at high temperature necessarily skip diapause and the samples are not pure samples on direct-developing embryos, but may also indicate a more simple miRNA control layer of diapause at least within this species.Prominent is the up-regulation of miR-9-3p and miR-29a-5p during diapause.MiR-9 is an important regulator of neurogenesis 73 by balancing quiescent and activated state of neuronal progenitors 74 .The miR-29 family, on the other hand, is strongly up-regulated during development and promotes neuronal differentiation 75,76 .In addition, both miRNAs are regulated during aging of N. furzeri, making them promising key regulators for diapause 39,77 .Most interestingly, the killifish-specific miR-19344 shows a 10-fold up-regulation within the intraspecies comparison (comparison II).As discussed already above, one of its potential targets is TSC2, regulating the mTORC1 complex, which is an important modulator of cell growth, protein synthesis, and cell proliferation 63 .This finding strongly indicates that miR-19344 evolved in the killifish lineage as part of the molecular adaptations that enable diapause.Not much is known about the other DEMs of this comparison or their absolute change in expression is relatively low: One study suggests that miR-7641 is involved in endothelial differentiation in embryonic stem cells 78 (this miRNA was differentially expressed in the annual comparison too), a recent study links miR-563 to the development of the spine in vitro 79 and expression of miR-210 appears to promote angiogenesis by inhibiting efna3, a suppressor of blood vessel formation 80 .Within the selected interspecies comparison (comparison III), we observed a total of 292 DEMs (171 up-regulated and 121 down-regulated in the annual Nothobranchius furzeri samples compared with the non-annual Aphyosemion striatum) belonging to 203 distinct miRNA loci (see https://osf.io/kahsb/).We identified a large proportion of DEMs in common between comparison I and comparison III, including the same directional changes in expression (see https://osf.io/92tvr/).This overlap includes the miR-430, miR-7641-3p, miR-10, and miR-17 mature miRNAs discussed above.Still, we identified unique DEMs belonging to 25 and 68 miRNA loci in the large interspecies comparison and the selected interspecies comparison, respectively.Interestingly, we found both mature miR-19344 forms up-regulated within the annual samples compared to the non-annual samples.Those miRNAs unique to the Nothobranchius furzeri and Aphyosemion striatum comparison might be in part responsible for the difference in embryogenesis in their clade but not necessarily in the other annual/non-annual killifish clades, because annual life appears to have evolved independently several times [1][2][3] .In general, the regulatory involvement of miRNAs during the embryogenesis of killifish seems to be pervasive.In contrast, only a few miRNAs seem to be involved in the maintenance of diapause in Nothobranchius furzeri, making these miRNAs critical regulators of this dormancy state.

mature miRNA switch
Within our differential expression analysis, we observed some interesting cases of DEMs.In particular, we found some miRNA genes that showed a "switch" in the expression of their main mature miRNA: We did not observe a change in the total transcribed amount of the pre-miRNA transcripts between contrasted conditions, but in the relative abundance of the 5p and 3p mature miRNAs.The reason could be an alteration of the processing machinery responsible for the biogenesis of the mature miRNA products.In addition to the miR-430 family, we identified two such cases in comparison I (miR-200a and miR-181a) and four in comparison III (miR-128, miR-130c, miR-217 and miR-223), for an example see Fig. 6 and https://osf.io/u6vbf/.We already briefly mentioned the regulatory targets of miR-200a and its potential role during diapause in the previous subsection.But also the mature miRNA switch cases show potential regulatory roles in different maturation processes: Developmental regulation of T cells and the vascular system is modulated by miR-181a 61,81 , miR-217 shows strong expression in late stages of the fetal rat development and an expression disruption showed damaged lung tissue 82 , and miR-223 is a well-studied modulator in hematopoiesis and osteoclast differentiation [83][84][85] .The functions of miR-128 are more linked to the formation and progression of various cancers, suggesting a general role in either cell development or proliferation, and it is known to inhibit telomerase activity [86][87][88] .Not much is known about miR-130c; however, its paralog miR-130a appears to be present in hematopoietic progenitor cell but not in mature blood cells 89 .

Strong correlation between the identified differentially expressed miRNAs and expression levels of developmental-related mRNAs
Next, we examined to which extent the results from the differential expression analysis of the miRNAs are reflected on the transcript regulation level.To do so, already identified aging-related differentially expressed genes 23 were compared with the predicted targets of the identified DEMs.We found that 143 DEMs have 88 unique mRNA targets that show a significant expression change within the selected interspecies comparison (see https://osf.io/s2rjp/).If we look at targets of some of the DEMs discussed above, we see a clear overlap to the identified differentially expressed mRNAs.Almost all of them are associated with cell developmental or proliferation processes and show opposite expression patterns compared with their targeting miRNAs.This observation indicates a strong correlation between the expression of the identified DEMs and the regulation of their target mRNAs.As one example, cyclin D1 is one of the targets of the DEM miR-430a/c-3p, promoting G1 progression of the cell cycle 90 and was found differentially expressed.Additionally, for each DEM, we tested if there was a specific enrichment of its targets among the DEGs and could identify 23 miRNAs where this was the case.Some of the most significant enrichments were observed for miR-23a-3p, miR-731-5p, miR-17a/b-3p, miR-430a-5p and miR-19b-5p, all of them discussed already above (see https://osf.io/gxcse/).No significant change in expression can be observed when comparing all reads of the pre-miR-200a annotation between both conditions (p-Value=0.43).However, a clear difference can be observed for the 5p and 3p mature miRNAs separately (p-Value of 1.0e −9 and 2.1e −4 ).In contrast to the non-annual A. striatum (n=4) were both mature sites are equally expressed, the 3p mature miRNA is increased in expression, whereas the 5p mature miRNA appears to be "switched off" within the annual N. furzeri (n=5).The green bases indicate the mature sequences.For more examples, see https://osf.io/u6vbf/.

The expression of miRNAs involved in embryogenesis and diapause regulation is modulated in an agingrelated manner
Our previous work revealed a correlation between transcriptome regulation during diapause and aging 23 .Here, we extend this analysis to miRNA regulation.We found miR-29a-5p and miR-29a-3p to be similarly up-regulated during diapause and aging within the aged brain, liver, and skin samples of Nothobranchius furzeri.The miR-29a gene was already known to be up-regulated in highly aged individuals and is also considered a tumor suppressor 39 .Another up-regulated miRNA in both conditions was miR-210, which acts in the hypoxia pathway 91 .When active, this miRNA positively regulates the adaptation to a hypoxic environment, helping cells to survive with only low available oxygen 91 .Diapausing killifish embryos show severe repression of aerobic metabolism 13,14 and are incredibly resistant to hypoxia 10 ; therefore, under normal oxidative conditions, hypoxia-sensitive genes may be down-regulated.Hypoxia is also a known condition in aged tissues, indicating that miR-210 is generally used to respond to this oxidative stress 92 .We observed the expression of miR-222a to be equally and strongly down-regulated during diapause and aging.When active, miR-222a promotes the expression of proteins responsible for muscle cell development by silencing the expression of the translational repressor cpeb3 93 .However, neither miR-9 nor the killifish-specific miR-19344 could be correlated with aging and seem to be specific for diapause and early developmental regulation.When including DEMs of the selected interspecies comparison (comparison III), we observed a considerable overlap of 57 and 40 DEMs being up-and down-regulated in the same way during embryogenesis and aging.Additionally, we found 99 DEMs being modulated in opposite directions (see https://osf.io/cqvnz/).Closer inspection revealed unusual behavior of some of the already discussed miRNAs.Both 3p mature transcripts of miR-430a/c are down-regulated during diapause but up-regulated during aging in the brain.In the aged brain and liver of Nothobranchius furzeri, miR-10b-5p is down-regulated, whereas it is up-regulated in embryonic development.The same is true for miR-92b-3p in the liver and skin.Also interesting, miR-181b-3p is down-regulated in both diapause and aging, but miR-181b-5p appears to be up-regulated during diapause but stays down-regulated in the aged liver and skin samples and not significantly changed in the brain during aging.These exemplarily described behaviors of specific miRNAs already indicate that there is an overlap between aging processes and the regulatory network controlling diapause and embryogenesis in Nothobranchius furzeri.However, to elucidate the precise mechanisms or regulatory responses that act in both conditions, more functional analyses need to be performed.

Evolution of diapause is linked to a contraction of miR-430 family copy numbers
The miR-430 locus in teleost fishes is large, variable across species, and characterized by a large number of repetitions of its family members.We reasoned that differential expression of miR-430 family members in the embryos of annual and non-annual killifish could be associated with variations in the copy number of miR-430 stem-loops.To test this hypothesis, we took advantage of the recently-sequenced genomes of Nothobranchius orthonotus, Aphysosemion australe, Callopanchax, and Pachypanchax platyfairii and the low-coverage genomes of 45 annual and non-annual killifish species from Africa.Highly repetitive loci are difficult to assemble from short-reads and thus may be under-represented in genome assemblies.To obtain an unbiased estimate of copy numbers, we, therefore, quantified the fraction of reads that match to any sequence of the miR-430 family (two mismatches allowed) in these different genomes, accounting for genome size.Although this approach only provides a crude estimate of copy number in each species, its precision should allow a comparison between different genera, which contain multiple species.The apparent result is that the copy number is lower in the annual genera Nothobranchius (N=10) and Callopanchax (N=4) as opposed to the non-annual genera Aphyosemion (N=9), Scriptaphyosemion (N=4), Archiaphyosemion (N=2), and Epiplatys (N=7) (see Figure 7).Only two genomes are available for the African clade, the annual Austrofundulus limnaeus and the closely related non-annual Kryptolebias marmoratus.Similarly, in this case, the difference is in the expected direction.The semi-annual genus Fundulopanchax showed an estimated copy number distribution intermediate between the non-annual genus Aphyosemion and the annual genus Nothobranchius, although they do not differ statistically.

Conclusion
In the present paper, we provide one of the first miRNA studies of vertebrate diapause.We could observe that regulation of miRNA expression in different clades of annual fishes shows convergent evolution, and the samples of different taxa do not cluster according to their phylogenetic relationships, but according to their physiological status (i.e., being in diapause or not).This observation represents a further example of how annual fish, belonging to different clades, independently converged on similar adaptations to support embryonic development in their ephemeral habitat 3,5 .Additionally, we describe different miRNAs as key modulators during fish embryogenesis and diapause regulation, such as members of the miR-10 and miR-430 family, miR-92, and miR-29a, but also killifish-specific ones, such as miR-19344 and miR-19337.Some of these miRNAs are of particular interest because their implications in specific developmental processes are already further studied in other species.For example, miRNA-430 acts primarily by inhibiting protein synthesis 94 .However, after the onset of zygotic transition and in adult cells, mRNA degradation appears to be the predominant mechanism 95 .Diapause is characterized by a prominent depression of protein synthesis that is reduced to 10 % of the pre-diapause levels 16 .It is therefore highly likely that both proteins and mRNAs are stabilized during diapause.In this context, miRNAs should act primarily, if not exclusively, by inhibiting protein synthesis.A prominent feature of miRNA regulation during diapause II is a dramatic down-regulation of miR-430 expression.This is paralleled by a reduction of the miR-430 copy numbers in the genomes of annual killifish, and this reduction may well represent the basis for this down-regulation.Upon exit of diapause, the release of miRNAs would allow the immediate onset of translation.This concept is further supported by the observation that diapausing embryos of Nothobranchius furzeri show paradoxical up-regulation of genes related to translational elongation 23 , suggesting that these embryos are primed for a catch-up process upon exit from diapause.Studies in the worm Caenorhabditis elegans revealed miRNAs that are regulated both during diapause and aging.In particular, miR-71 does not influence diapause entry but is necessary for survival during diapause in Caenorhabditis elegans 24 .Deletion of miR-71 results in a shortened lifespan, and miR-71 abundance increases with age before dropping late in life 27 .The timing of the miR-71 expression drop can be used to predict the longevity of individual worms 25 .Our results show a similar overlap between miRNAs regulated during diapause and aging in the annual fish Nothobranchius furzeri.Several miRNAs significantly modulated in diapausing Nothobranchius furzeri were modulated during aging as well.One of these miRNAs, miR-101a, is particularly interesting.He et al. showed that miR-101 acts as a highly-connected hub in gene regulatory networks of transcription factors and epigenetic modulators involved in cell cycle progression 96 .Overexpression of miR-101 is also known to induce cell cycle arrest in different cell types 96 .So, high levels of miR-101 in diapausing embryos would contribute to the G1 block typical of diapause 15 .On the other hand, we observed down-regulation of members of the miR-17/92 cluster, also known as oncomiR-1.The organization of this cluster is conserved in Nothobranchius furzeri.Furthermore, it was observed to be down-regulated during the aging of the brain 39 and adult neuronal progenitors 97 .OncomiR-1 is a significant regulator of cell cycle and overexpressed in several tumors 98 .On the other hand, it is well-established that oncomiR-1 is down-regulated during the aging of human mitotically-active cells where it targets the cell cycle inhibitor p21 99 .An overlap in the regulation of cell cycle protein-coding genes between diapause and aging was also demonstrated by RNA-seq analysis in Nothobranchius furzeri 23 , further supporting the concept that, like in Caenorhabditis elegans, aging and diapause are controlled by overlapping genetic pathways.Our study results serve as a sound basis for further research to infer the impact and precise control of miRNA regulation of these overlapping pathways.

Figure 2 .
Figure 2. Actively transcribed conserved and killifish-specific miRNAs.Phylogenetic tree of the investigated killifish species and the number of their actively transcribed miRNA genes (red -miRNAs belonging a known miRNA family; green -miRNAs specific to killifish; blue -miRNAs not being transcribed in any of the Nothobranchius furzeri samples).Mature miRNAs having a count of more than ten reads were considered to be actively transcribed.

Figure 3 .
Figure 3. Principal component analysis of the investigated small RNA samples.PCA of the 18 RNA-Seq libraries based on the differentially expressed miRNAs identified in the large interspecies comparison (comparison I).The libraries generated from annual and non-annual species are relatively well separated and samples from Nothobranchius furzeri and Aphyosemion striatum cluster together, respectively.Samples from N. furzeri embryos that skipped diapause (marked with a white dot) show a visible shift in PC2 towards the A. striatum cluster.

Figure 4 .
Figure 4. Relative amounts of the total expressed reads within all 18 sequenced killifish samples of the six most highly expressed miRNA families.The conserved miRNAs of the miR-10 and miR-430 families, miR-92a, miR-181a, and the two killifish-specific miR-19344 and miR-19337 showed the highest relative expression in all investigated samples.All of these miRNAs are known or predicted to play critical roles in the regulation of developmental processes in embryonic cells.For details on the miRNA expression levels, see https://osf.io/epf6z/.

Figure 6 .
Figure 6.Expression profiles of the miR-200a gene from the large interspecies comparison (comparison I).No significant change in expression can be observed when comparing all reads of the pre-miR-200a annotation between both conditions (p-Value=0.43).However, a clear difference can be observed for the 5p and 3p mature miRNAs separately (p-Value of 1.0e −9 and 2.1e −4 ).In contrast to the non-annual A. striatum (n=4) were both mature sites are equally expressed, the 3p mature miRNA is increased in expression, whereas the 5p mature miRNA appears to be "switched off" within the annual N. furzeri (n=5).The green bases indicate the mature sequences.For more examples, see https://osf.io/u6vbf/.
Difference of genomic miR-430 copy number within different killifish species.Estimated number of miR-430 gene copies within different killifish species based on Illumina short read data.A significant difference between the annual genera Nothobranchius (n=10) and the non-annual genera Aphyosemion (n=9) (Wilcoxon rank sum test: W=11, p-value=0.002089),as well as, between the annual genera Callopanchax (n=4) and the non-annual genera of