Sequencing of small RNAs of the fern Pleopeltis minima (Polypodiaceae) offers insight into the evolution of the microrna repertoire in land plants

MicroRNAs (miRNAs) are short, single stranded RNA molecules that regulate the stability and translation of messenger RNAs in diverse eukaryotic groups. Several miRNA genes are of ancient origin and have been maintained in the genomes of animal and plant taxa for hundreds of millions of years, playing key roles in development and physiology. In the last decade, genome and small RNA (sRNA) sequencing of several plant species have helped unveil the evolutionary history of land plants. Among these, the fern group (monilophytes) occupies a key phylogenetic position, as it represents the closest extant cousin taxon of seed plants, i.e. gymno- and angiosperms. However, in spite of their evolutionary, economic and ecological importance, no fern genome has been sequenced yet and few genomic resources are available for this group. Here, we sequenced the small RNA fraction of an epiphytic South American fern, Pleopeltis minima (Polypodiaceae), and compared it to plant miRNA databases, allowing for the identification of miRNA families that are shared by all land plants, shared by all vascular plants (tracheophytes) or shared by euphyllophytes (ferns and seed plants) only. Using the recently described transcriptome of another fern, Lygodium japonicum, we also estimated the degree of conservation of fern miRNA targets in relation to other plant groups. Our results pinpoint the origin of several miRNA families in the land plant evolutionary tree with more precision and are a resource for future genomic and functional studies of fern miRNAs.

Introduction exindusiate sori (cluster of sporangia) of round shape on the abaxial side (Fig 1B-1D). One particular characteristic of P. minima is its dissecation tolerance (poikilohydry) [28,29]. We have sequenced the small RNA fraction of P. minima and searched for miRNAs phylogenetically conserved in other plants. We present a group of conserved fern miRNAs and their predicted targets, thereby advancing our understanding of miRNA evolution in the fern and land plant lineage.

Sequencing and analysis of P. minimas RNAs
Small RNAs (sRNAs) were extracted from two samples of fertile fronds of P. minima, converted to cDNA and subjected to Illumina sequencing, generating a total of 25,947,729 reads, of which 20,989,537 were high-quality reads with lengths between 15 and 45 nt. Of these, 6,268,973 reads corresponding to mRNA fragments, tRNAs, rRNAs, snoRNAs, snRNAs and repetitive elements were discarded, leaving a total of 14,720,564 sRNA reads to be analysed (S1 Table).
The length of functional sRNAs are thought to lay in the 20-24 nucleotide (nt) range. In most angiosperms, sRNA length distribution typically displays two peaks at 21 and 24 nt, with miRNAs centred around 21 nt, while 24 nt small-interfering RNAs (siRNAs) are associated with silencing of retroposons by heterochromatinization [12]. In non-angiosperms, the presence of the 24-nt sRNA peak is not always evident, but it has recently been shown that the  [3,17,24]. Note that some molecular phylogenies recover a clade composed of liverworts and mosses as sister to vascular plants [25,26]. The fern ("monilophytes") clade is the sister group to seed plants, gymnosperms and angiosperms. Within ferns, P. minima belongs to the order Polypodiales and is more closely related to Ceratopteris thalictroides, while Marsilea quadrifolia belongs to the order Salviniales and Lygodium japonicum to the order Schizaeales. heterochromatic siRNA pathway is a primitive feature of land plants, being present in mosses [30]. In P. minima, we observed a prominent 21 nt peak and a smaller, but noticeable, peak at 24 nt (Fig 2). A similar distribution is found in the aquatic fern M. quadrifolia [23], showing that the formation of 24-nt sRNAs, probably associated with heterochromatin formation, is a common characteristic of ferns.

Identification of conserved P. minima miRNAs
Since P. minima lacks a genome sequence and no fern genome is available for comparison, we focused on finding miRNAs in our samples that are conserved in other land plant genomes. To this end, the 14,720,564 sRNA reads of P. minima were compared to version 21 of miRBase [31], yielding 1,127,970 reads that could be aligned to miRNAs or MIRNA genes in the database, corresponding to 3,912 unique sequences (S1 Table).
It is known that the confidence on the existence of miRNA families in miRBase varies widely, and the evidence for many miRNA families is rather weak, specially for miRNAs coming from species lacking a genome sequence [11]. In view of this, to identify bona fide P. minima miRNAs we required the miRNA families 1) to be represented by at least 10 reads of 20-22 nt in our data; 2) to exhibit consistent 5' cleavage processing of the mature miRNA strand and 3) to be identical or very similar to highly confident sequences of mature miRNA or miRNA Ã strands of other plant species, as defined by [11]. For some miRNA families, namely miR162, miR395 and miR477, even though the total number of reads was under 10, we still considered them to be real miRNAs since the reads included both miRNA and miRNA Ã strands that were identical or nearly identical in sequence to the corresponding miRNA and miRNA Ã sequences from other plants (S2 Table).
The analysis revealed a total of 57 conserved miRNAs in P. minima, belonging to 23 miRNA families (Table 1). MiR156 is related in sequence and function to miR529 [32,33], leading to them being often considered to be members of the same family [11], the same happening between miR159 and miR319 [34,35]. Thus, if miR156/miR529 and miR159/miR319  are regarded as two families, the total number of conserved P. minima miRNA families would be 21 instead of 23. MiR1030 is a miRNA that has been recognised by this name in the moss, P. patens [13], as well as in the liverwort, M. polymorpha [36]. However, it is closely related in sequence to miR530 of seed plants, and we decided to consider it as an unique family, miR530/ 1030, as suggested by Taylor et al [11]. For five families (pmi-miR162, pmi-miR390, pmi-miR395, pmi-miR477 and pmi-miR529), sequences that represent the putative miRNA/miRNA Ã duplex strands (i.e., the 5p and 3p strands of the miRNAs) were identified (S2 Table). Sixty-one percent (35/57) of the conserved P. minima miRNAs, belonging to 18 of the 23 miRNA families, have a 5' terminal uridine residue (Table 1), a conserved feature of miRNAs recognized by the AGO1 protein, while miRNAs from the pmi-miR390, pmi-477 and pmi-miR529 families have mostly a 5' terminal adenine residue, a feature found in miRNAs recognized preferentially by AGO2 and AGO4 in angiosperms [37].

Ancient land plant miRNA families
MiRBase contains miRNAs from well-characterized plant species for which genome information is available, like the moss P. patens, the lycopod S. moellendorfii and the angiosperms Arabidopsis thaliana and Oryza sativa [11,13], but miRNAs coming from many recent sequencing surveys are not present in the database. To expand the scope of our phylogenetic analysis, we compared the conserved P. minima miRNAs to those of Marsilea quadrifolia [23], a water fern that diverged from polypod ferns like P. minima around 220 MYA [17,24]. The miRNAs of M. quadrifolia and seed plants identified by Chávez-Montes et al. [23] were reanalysed to identify bona fide conserved miRNAs (see Materials and Methods).
We also compared P. minima sequences with the miRNAs of various gymnosperms and basal angiosperms [23,[38][39][40][41], as well as conserved miRNAs identified recently in the liverworts Pellia endiviifolia [42] and Marchantia polymorpha [36,43]. This allowed us to trace the history of miRNA origin and conservation in land plants, with emphasis on the fern miRNA repertoire.
Particularly noteworthy are miRNA families miR156/529, miR159/319, miR160, miR166, miR171 and miR408, which are found in all land plant groups, including liverworts and mosses, which represent the earliest-branching groups in land plant evolution (Figs 1A and 3). Some ancestral miRNAs have a more patchy phylogenetic distribution and are missing in a few basal groups. For instance, miR396 is absent from moss, miR167 is missing from liverworts and lycopod, while miR530/1030, miR535 and miR390 are specifically missing in the lycopod. In any case, the P. minima miRNA repertoire is conserved compared to that of basal land plant groups, since all 14 ancient miRNAs present in liverworts (P. endiviifolia and/or M. polymorpha) and all 15 ancient miRNAs of the moss, are present in P. minima. Similarly, all 10 ancient miRNAs found in the lycopod are also present in P. minima (Fig 3). For gymnosperms we took into account miRNAs of P. abies [6,41], C. lanceolata [38] and T. mairei [39] and species reported in [23]; for angiosperms we included miRNAs of species reported in [23] and in miRBase release 21.0 [31]. Pairs of miRNAs that can be considered from the same family are indicated by coloured dots. Note that miRNAs are not necessarily present in all members of a group. Among liverworts, miR530/1030, miR529 and miR536 are only found in M. polymorpha, while miR156, miR395 and miR535 are only found in P. endiviifolia [36,42,43]. MiR403 and miR529 are absent from most non-dicot angiosperms [44] and eudicots [32,33], respectively. Some miRNA families are only present in a few basal groups. One is miR1024, a miRNA from the moss P. patens [13] which is found in P. minima and in no other species, not even M. quadrifolia (Fig 3). This miRNA is present in two copies in the P. patens genome and is considered to be a high-confidence miRNA, based on RNA-seq analyses and precursor structure [13,30]. Although this miRNA has only 21 reads in our samples, it displays 100% identity over 20 nt to P. patens miR1024, indicating that it is a real miRNA. Thus, it might represent an ancient miRNA that has been independently lost from several lineages, including seed plants. Another case is that of miR1083, which is found in P. minima (but not in M. quadrifolia, see Materials and Methods) with over 1,700 reads (Table 1). This miRNA is present in the lycopod, S. moellendorffii [13], is abundant in gymnosperms [23,39,41], but cannot be identified in angiosperms (Fig 3; see Materials and Methods). Thus, miR1083 seems to be an ancient tracheophyte miRNA that has been specifically lost from angiosperms.

Euphyllophyte miRNAs
We found that some miRNAs families are only common to ferns (P. minima, M. quadrifolia or both) and seed plants (gymnosperms and angiosperms), suggesting that these miRNAs originated in the early evolution of euphyllophytes. These are miR162, miR168, miR169, miR172 and miR403 (Fig 3). All of these miRNAs, except for miR403, are also found in M. quadrifolia [23], while miR172 is also found in C. thalictroides [22].
Although miR403 has been reported as being present in M. quadrifolia [23], we were unable to find this miRNA in the miRNAome of this fern and did not include it in Fig 3 (see Materials and Methods). MiR403 is an unusual miRNA in that it seems to have a very patchy phylogenetic distribution. It is usually considered to be missing in non-dicot angiosperms [44], but it has been detected in some monocots [23]. MiR403 has also been detected in some gymnosperms like the Chinese fir, Cunninghamia lanceolata [38], the spruce, Picea abies [41] and the Maire yew, Taxus mairei [39], although not in Ginkgo biloba [40]. All of this suggests a complex evolutionary history for this miRNA family, and its presence in ferns should be confirmed in the future with further RNA-seq data from other fern species and the sequencing of fern genomes.

Conserved fern miRNA targets
Extensive experimental evidence from several angiosperms and the moss P. patens, have unveiled that the mRNAs targeted by ancestral miRNAs are often also deeply conserved [10, 13,14,45]. Recently, degradome data has shown that the targets of liverwort ancestral miRNAs are also partially conserved with that of other lineages [41][42][43]. As for ferns, the only available information consists of putative mRNA targets for C. thalictroides miR160, miR171 and miR172, which were also found to be conserved with those of other plant lineages [22].
To confirm and expand the evidence for conservation of ancestral miRNA targets in ferns, we computationally predicted the mRNA targets of conserved P. minima miRNAs. Since no P. minima transcriptome is available, we made use of the extensive transcriptome data obtained for the fern Lygodium japonicum [20]. L. japonicum belongs to the Schizaeales order, which places it as a distant cousin to polypod ferns like P. minima [17,24]. We searched the L. japonicum transcriptome for potential targets of the conserved P. minima miRNA repertoire using the psRNATarget platform [46] in order to identify putative fern miRNA targets. Since functional data on fern miRNA targets is extremely limited, we restricted our analyses to targets which have been experimentally determined as targets in other plant groups.
The analysis indicates that 11 out of the 23 P. minima miRNA families have the potential to target L. japonicum RNAs that are also targeted in other species (Table 2). Thus, pmi-miR156 and the related pmi-miR529 target a mRNA that encodes a Squamosa Promoter Binding (SPB) transcription factor; pmi-miR159 and the related pmi-miR319 target mRNAs that encode transcription factors of the MYB family; pmi-miR160 targets mRNAs encoding transcription factors of the Auxin Response Factor (ARF) family; pmi-miR166 targets a mRNA encoding a transcription factor of the class III homeodomain-leucine zipper (Class III HD-Zip) family; pmi-miR171 targets a GRAS-domain transcription factor mRNA; miR172 targets the mRNA of a Apetala-2 transcription factor ( Table 2). Targets homologous to these have been observed in experimental work done in angiosperms, moss, or both [13,14,45]. Importantly, the L. japonicum mRNA regions targeted by pmi-miRNAs are located in the same positions as in the targets of other species, strongly suggesting that these are ancient, conserved targets of land plant miRNAs (Fig 4). In addition, sequence alignments show that the RNA regions targeted by miRNAs are strongly conserved between mRNAs of L. japonicum and homologous mRNAs from other species; indeed the miRNA-targeted elements are generally more conserved than mRNA regions that encode conserved protein domains (S1-S10 Figs). The P. minima miR160, miR171 and miR172 predicted targets confirm those of Axtell and Bartel [22] for C. thalictroides. Recently, a thorough survey found that ferns possess three clades of Class III HD-Zip genes [47], and we found that pmi-miR166 can potentially target all Class III HD-Zip paralogues of the fern Psilotum nudum identified in that study (S5 Fig and  data not shown).
The miR390 family does not control protein-coding mRNAs, but rather targets a group of non-coding transcripts called TAS3 RNAs (trans-acting siRNAs). In the moss P. patens and in angiosperms, two miR390 molecules hybridize to different segments of a TAS3 RNA, which is then processed by a complex mechanism to produce one or more 21-nt ta-siRNA molecules that are incorporated into a RISC complex that targets the mRNAs of ARF3 and ARF4 transcription factors [48]. In addition, P. patens TAS3 can also produce ta-siRNAs that target the transcription factor AP2 [13], and there is evidence that the ta-siRNAs generated from the TAS3 of the liverwort M. polymorpha can target AP2 as well [36,49]. As for the P. minima miR390 family, we found that they can potentially target two L. japonicum non-coding transcripts which exhibit characteristics of TAS3 RNAs (Table 2, Fig 5, S9 Fig). These L. japonicum TAS3 transcripts carry two sites predicted to hybridize to pmi-miR390; in addition, between the miR390 sites one potential ta-siRNA, derived from the (+) strand of the RNA, can be produced that could target L. japonicum ARF3/4 mRNAs (Fig 5 and S9 Fig), indicating that the control of ARF genes by miR390 is conserved in ferns.
Of particular interest are some predicted P. minima miRNA targets that are conserved with liverwort miRNAs but not other plants. Specifically, apart from transcription factors of the MYB family, pmi-miR159 and pmi-miR319 potentially target a mRNA that encode a RWP-RK domain-containing protein (RKD, Table 2, Fig 4, S3 Fig). Degradome studies in the liverworts P. endiviifolia and M. polymorpha identified mRNAs encoding RWP-RK proteins as bona fide targets of miR159 and miR319 homologues of these species [42,43]. Interestingly, the pairing between pmi-miR159/319 and its potential L. japonicum RWP-RK target occurs in the 3'-UTR of the mRNA, the same region that is targeted in the homologous liverwort M. polymorpha mRNA by miR319 (Fig 4, S3 Fig). In addition, although a RWP-RK mRNA has not yet been experimentally identified as a target for miR159/319 in moss [13,14,45], we found in the P. patens transcriptome a RWP-RK mRNA that harbours a potential miR159/319 target site near its 3' end, with a similar position and sequence to the RKD targets in liverworts and fern (S3 Fig). As for seed plants, mRNAs encoding RWP-RK proteins have not been identified as targets of miR159/319 in gymnosperms or angiosperms [14,38,39,41,50] and we could not find potential miR159/319 target elements in the RWP-RK mRNAs of these groups, suggesting that the miR159/319-mediated regulation of RWP-RK proteins is an ancient land plant feature that may have been lost from seed plants. Another target possibly conserved between ferns and liverworts is the mRNA encoding the enzyme polyphenol oxydase (PPO), which was identified as a potential target of pmi-miR408 (Table 2). In the liverwort P. endiviifolia, PPO was also identified as an experimental target of miR408 by degradome analysis [42]. Significantly, miR408 targets the 3'-UTR of PPO mRNAs in both ferns (Fig 4, S10 Fig) and liverworts [42].
In seed plants, some miRNAs target genes involved in the biogenesis and function of miRNA themselves. Thus, in angiosperms, miR162 targets Dicer-like 1 (DCL1), miR168 targets AGO1 and miR403 targets AGO2 [14]. In our analyses, P. minima miR162 and miR403 are not predicted to target genes similar to DCL1 or AGO2, but pmi-miR168 is predicted to target L. japonicum AGO1 (Table 2 and Fig 4). The sequence predicted to be targeted by pmi-miR168 is located within the 5' half of the AGO1 coding region, outside the segments that microRNAs of the fern Pleopeltis minima encode the conserved domains of the protein, and corresponds to the homologous AGO1 region targeted by miR168 in A. thaliana and other seed plants (Fig 4 and S6 Fig) [51]. Thus, the negative feedback loop involving miR168 and AGO1 may be an ancestral euphyllophyte feature.

Discussion
In this work, we identify the conserved miRNAs of a polypod fern, P. minima, by comparing its miRNA repertoire to those of other land plant groups. Fig 6 shows a revised scenario for the emergence of different miRNA families during early land plant evolution. Concerning the fern miRNAome, our analyses allowed us to reach a series of conclusions. First, we found that 15 ancient miRNA families that originated in the beginning of land plant evolution, as evidenced by their presence in liverworts and/or mosses, are all present in P. minima, showing that the fern miRNA repertoire is very conserved in this regard (Fig 6). This is not a trivial observation, as illustrated by the lycopod S. moellendorffii, which lacks half of the ancestral miRNAs that are present in bryophytes (Fig 3). Second, by comparing P. minima miRNAs and those of the aquatic fern, M. quadrifolia, with the rest of land plants, we identified five miRNAs that seem to have originated at the beginning of euphyllophyte evolution, before the split of the lineages of ferns and seed plants (gymnosperms and angiosperms; Fig 6). The miRNAome of P. minima  4). Note that, according to recent molecular phylogenies [25,26], liverworts and mosses might belong together in a clade that is sister to vascular plants, in which case nodes 1 and 2 would be the same. Also note that five miRNA families seem to be restricted to euphyllophytes and that only miR1083 seems to have originated in the branch leading to vascular plants. seems to be more conservative than that of M. quadrifolia, since six ancestral miRNAs found in the former could not be identified in the latter (Fig 3). Third, by searching putative targets of P. minima miRNAs using the transcriptome of another fern, L. japonicum, we provide evidence that the targets of 11 fern miRNA families are deeply conserved in land plant evolution.
Molecular studies on basal land plant groups can give crucial clues on the evolution of plant form [3]. As gene expression regulators, miRNAs are essential for proper plant development and may have played important roles in the early evolutionary diversification of land plants, making comparative functional studies of ancestral miRNAs in different plant groups very interesting in this regard [10]. Some ancestral miRNAs present in P. minima exhibit an intriguing phylogenetic distribution and are interesting candidates for functional studies. One of these is miR1083, a tracheophyte-specific miRNA that is present in lycopods and is abundantly expressed in P. minima and gymnosperms, but is apparently absent from angiosperms (Figs 3 and 6). The other is miR1024, which is present in moss and P. minima, but has not been detected in any other species investigated so far (Fig 3) [13]. Neither the expression pattern nor the targets of miR1083 and miR1024 have been studied in any plant, but their antiquity and conservation make them interesting subjects for future study.
Our bioinformatic search for potential targets of fern miRNAs indicate that 11 of the ancestral pmi-miRNAs might target genes that are homologous to those of other land plants ( Table 2, Figs 4 and 5). Although our evidence is based on an in silico approach, the conservation of potential miRNA/target pairs in species that are separated by hundreds of millions of years of independent evolution is strong evidence that these fern transcripts are real targets. Among these, we found that pmi-miR390 could target TAS3 RNAs with the potential to generate tasi-RNAs that, in their turn, target genes encoding ARF3/4 transcription factors from L. japonicum (Fig 5 and S9 Fig). This represents the first full description of TAS3 transcripts in ferns and an indication that miR390 function is conserved in this group, as suggested by a previous PCR-based survey [49].
Another interesting fern conserved target was found for pmi-miR166, which can target transcription factors of the Class III HD-Zip (C3HDZ) family. In angiosperms, C3HDZ genes are expressed in the adaxial side of developing leaves and are necessary for the proper dorsoventral patterning of these structures [52]. In the abaxial side, on the other hand, C3HDZ expression is repressed by miR166, which is also necessary for proper leaf patterning in angiosperms [53,54]. Even though it is often considered that the leaves of ferns and seed plants evolved independently [3,55], Vasco et al. [47] have recently observed that C3HDZ genes are specifically expressed in the adaxial side of developing fern leaves, just like in angiosperms, giving credence to the hypothesis that euphyllophyte megaphylls are homologous structures. Considering our observation that pmi-miR166 could target fern C3HDZ genes, future studies on the expression and function of fern miR166 may help establish the extent of molecular homology in the gene regulatory network that controls leaf development in euphyllophytes.
A recurrent feature in the evolution of plant gene regulation is that certain mRNAs encoding proteins involved in miRNA processing and function are themselves controlled by miR-NAs. In angiosperms, AGO1, DCL1 and AGO2 are targets of miR168, miR162 and miR403, respectively [14], and AGO1 is also a miR168 target in the Chinese fir, a gymnosperm [38]. Interestingly, liverworts and mosses also target AGO1 and DCL1 with a different set of miR-NAs. Thus, miR902 and miR1047 target AGO1 and DCL1 in the moss P. patens [13,56], while miR11707 targets AGO1 in the liverwort, M. polymorpha [43]. As for ferns, we found that pmi-miR168 has the potential to target L. japonicum AGO1 in a site that, given its position in the mRNA, seems to be homologous to the AGO1 site targeted by miR168 in angiosperms (Fig 4,  S6 Fig). Since miR168 originated in the lineage leading to ferns and seed plants (Fig 6), it can be hypothesized that the negative feedback loop involving AGO1 and miR168 has been present throughout the evolutionary history of euphyllophytes. The confirmation of the miR168/ AGO1 regulatory relationship, as well as the study of miR162 and miR403 targets, constitutes an interesting avenue of research in fern regulatory RNA research.
Target prediction of P. minima miRNAs yielded two candidate targets that seem to be conserved only between ferns and liverworts. One is the gene encoding polyphenol oxydase (PPO), which originated in the early stages of the colonisation of land by plants and has roles in pathogen defence [57]. In both L. japonicum and the liverwort, P. endiviifolia, PPO is targeted at the 3'-UTR by miR408 (Fig 4 and S10 Fig). Another common potential target is RKD, a gene encoding a transcription factor of the RWP-RK family. In the liverworts, P. endiviifolia and M. polymorpha, RKD is targeted at the 3'-UTR by miR319 as revealed by degradome analyses [42,43]. Similarly, we observed that a L. japonicum homologue of RKD is predicted to be targeted at the 3'-UTR by pmi-miR319 and pmi-miR159 (Fig 4, S3 Fig). Recently, it has been found that RKD is necessary for proper germ cell development in liverworts, as RKD-deficient M. polymorpha plants exhibit egg cells with aberrant cellular differentiation and proliferation properties, along with defects in gemma cup formation, which are vegetative reproductive organs [58,59]. Work in A. thaliana indicates that RKD homologues might have related functions in egg development of seed plants [60]. Interestingly, the overexpression of miR319 in M. polymorpha causes defects in thalli growth and gemma cup formation in gametophytes [36]. Thus, given the important functions of RKD, its potential regulation by miR319 in liverworts and ferns might represent an ancient gene regulatory interaction that plays a role in gametophyte growth and germ cell differentiation in basal land plants.
Despite their ecological importance and their key phylogenetic position as sister group to seed plants, the genetic and molecular mechanisms underlying fern development and physiology are poorly known. We hope that our work will add to the growing genomic resources available to researchers dedicated to the study of this group of land plants.

Plant material
Specimens of Pleopeltis minima (Bory) J. Prado & R.Y. Hirai (= Polypodium squalidum Vell.; = Pleopeltis squalida (Vell.) de la Sota.) [27] were collected in a peri-urban environment at Departamento General Alvear, Corrientes Province, Argentina (29˚05'56.0"S, 56˚33'12.0"W). Annual precipitation in the area reaches 17,000 mm approximately. At the time, no special permits were needed to collect plant specimens from a peri-urban environment. The ferns were found growing on branches of large trees, and a couple of specimens with the branches attached were maintained in a greenhouse at 23˚C under a 16 h/8 h light/dark cycle. Field work did not involve endangered or protected species.

RNA preparation and high-throughput sequencing
Two independent samples of P. minima fertile fronds from two different specimens each were subjected to a protocol for total RNA extraction enriched for small RNAs using the mirVana miRNA Isolation Kit (Thermo Fisher Scientific). The quantity and quality of the total RNAs were analised with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific) and an Agilent 2100 Bioanalyzer (for concentration, 28S/18S and RIN detection; Agilent Technologies). Small RNA libraries were generated from both P. minima independent RNA samples using the Illumina Truseq SmallRNA Preparation kit(Illumina, San Diego, USA). Total sRNAs were ligated to 3p and 5p adapters (ADTs), and the corresponding cDNA was obtained by reverse-transcription PCR. The purified cDNA library was used for cluster generation on Illumina's Cluster Station and then sequenced on Illumina GAIIx (LC Sciences, Houston, USA).

Identification of conserved miRNAs
Initial analysis of reads were done with the ACGT101-miR bioinformatics program (LC Sciences, Houston, USA). Raw reads were filtered to remove adaptor sequences and reads which were too short (<15 nt) or too long (>40 nt). Sequences corresponding to mRNA fragments, rRNA, tRNA, snoRNA and snRNA mapping to Rfam (http://rfam.janelia.org) were removed, as well as reads mapping to RepeatMasker (http://www.repeatmasker.org/). The non-redundant sRNA reads were mapped to plant precursor and mature miRNA sequences stored in miRBase 21.0 (http://www.mirbase.org/) [31] using Bowtie software to identify an initial group of putative conserved miRNAs and precursors. Length variation (1-3 nt) at both 3' and 5' ends and up to two mismatches inside of the sequence were allowed in the initial alignments, and only sequences of 20-24 nt were selected for further analysis. Putative miRNAs identified by Bowtie were manually compared to miRBase (http://www.mirbase.org/search.shtml) using BLASTN and SSEARCH to confirm the identity of the miRNAs. The final set of conserved P. minima miRNAs (Table 1) are either 100% identical or display minor changes (one or two mismatches or 1-nt length variation) to annotated miRBase sequences. The identified P. minima miRNAs were compared to a recent re-evaluation of miRBase entries [11] to discard lowconfidence miRNAs from our sample.
To better characterize the phylogenetic distribution of plant miRNAs, P. minima sequences were also checked against a large plant miRNA database generated recently [23] which includes the aquatic fern, M. quadrifolia. The miRNAs that were reported for M. quadrifolia by Chávez-Montes et al [23] were reanalysed more stringently in the following way: First, miRNA sequences for each putative miRNA family were compared against miRBase (http://www.mirbase.org/) and miRNA sequences that exhibited low similarily to miRBase entries (E-value of less than 0.01) were discarded. Next, the abundance of each miRNA sequence was checked against the Comparative Sequencing of Plant miRNA database (http://smallrna.danforthcenter.org/) [23] and sequences that were expressed at low levels (less than 1 read per million) were discarded. After this conservative reassessment, several miRNAs that seemed to be present in M. quadrifolia and other plants [23] were discarded from our evolutionary analysis. The resulting bona fide conserved miRNAs were incorporated into the descriptions of miRNA distribution and evolution shown in Figs 3 and 6.

Data availability
The raw sequencing data of P. minima small RNA transcriptome was submitted to the NCBI Sequence Read Archive under the accession number SRR4294182.
Supporting information S1 Fig. Predicted targeting of a fern SPL mRNA by miR156/miR529. (A) Sequence of a transcript (Locus_14997) encoding a Squamosa-binding protein-like (SPL) homologue from the fern L. japonicum. The region predicted to be targeted by pmi-miR156 or pmi-miR529 is indicated in yellow. The starting ATG and stop codon are highlighted in blue. Sequence of a transcript (Locus_11016) encoding a Scarecrow-like(SCL)/GRAS domain transcription factor protein from the fern L. japonicum. The region predicted to be targeted by pmi-miR168 is indicated in yellow. The transcript is incomplete and lacks the starting ATG, but the stop codon is highlighted in blue. (B) Alignment of part of the SCL transcripts from fern L. japonicum (Lja), the liverwort M. polymorpha (Mpo), the lycopod S. moellendorffii (Smo), the moss P. patens (Ppa), the gymnosperms Pinus tabulliformis (Pta) and Pinus radiata (Pra); the basal angiosperm Amborella trichopoda (Atr), the dicots A. thaliana (Ath), Solanum lycopersicum (Sly) and Populus trichocarpa (Ptr), and the monocots O. sativa (Osa) and Brachypodium distachyon (Bdi). Residues displaying over 75% identity are highlighted. The region targeted by miR171 is indicated in red, and the GRAS domain in green. Note that the miRNAtargeted region is conserved in all mRNAs and species. (C) Predicted pairing between pmi-miR171 and L. japonicum Locus_11016. The E-complementarity score between miRNA and target RNA as estimated by the psRNATarget program is shown. (DOCX) S8 Fig. Predicted targeting of a fern AP2 transcription factor mRNA by miR172. (A) Sequence of a transcript (Locus_9880) encoding a Apetala-2 transcription factor protein from the fern L. japonicum. The region predicted to be targeted by pmi-miR172 is indicated in yellow. The starting ATG and stop codon are highlighted in blue. (B) Alignment of part of the AP2 transcripts from ferns L. japonicum (Lja) and Ceratopteris thalictroides (Cta), the gymnosperms Pinus tabulliformis (Pta) and Ginkgo biloba (Gbi); the basal angiosperm Amborella trichopoda (Atr), the dicot A. thaliana (Ath, and the monocots O. sativa (Osa) and Brachypodium distachyon (Bdi). Residues displaying over 75% identity are highlighted. The region targeted by miR172 is indicated in red, and the AP2 DNA-binding domain in green. Note that the miRNA-targeted region is conserved in all mRNAs and species. (C) Predicted pairing between pmi-miR172 and L. japonicum Locus_9880. The E-complementarity score between miRNA and target RNA as estimated by the psRNATarget program is shown. and Locus 39179) from L. japonicum predicted to be targeted by pmi-miR390. The regions predicted to pair with miR390 are indicated in blue, and the trans-acting small-interfering RNAs (tasi-RNAs) derived from TAS3 processing are shown in yellow. The tasi-RNAs are predicted to target mRNAs encoding ARF3/4 transcription factors. (B) TAS3 RNAs were aligned with the CLUSTAL Omega program. Regions predicted to be targeted by pmi-miR390 are indicated in blue, with the base pairing with miR390 indicated. Note that the pairing between pmi-miR390 is more extensive at the 5' site than at the 3' site, indicating that RISC-mediated cleavage of the RNAs probably occurs only at the 5' site. The tasiRNAs derived from TAS3 processing are shown in yellow. (C) Alignment of part of the TAS3 transcripts from L. japonicum (Lja) and a TAS3 from the dicot A. thaliana (Ath). Residues displaying 100% identity are highlighted. (D) Sequence of a transcript (Isotig23398) encoding a ARF transcription factor protein from the fern L. japonicum. The region predicted to be targeted by pmi-tasiR-ARF is indicated in yellow. The starting ATG and stop codon are highlighted in blue. (E) Alignment of part of ARF transcripts from fern L. japonicum (Lja), the gymnosperm Pinus pinaster (Ppi); the basal angiosperm Amborella trichopoda (Atr), the dicot A. thaliana (Ath) and the monocot Brachypodium distachyon (Bdi). Residues displaying 100% identity are highlighted. The region targeted by pmi-tasiR-ARF is indicated in red, and the B3 DNA binding and the Auxin responsive factor domains are also indicated. Note that the tasiR-ARF-targeted region is conserved in all mRNAs and species. (DOCX) S10 Fig. Predicted targeting of a fern PPO mRNA by miR408. (A) Sequence of a mRNA encoding Polyphenol oxydase (PPO, Locus 5437) from L. japonicum predicted to be targeted by pmi-miR408. The targeted region is located in the 3'-UTR of the mRNA (yellow). The starting ATG and stop codon are highlighted in blue. (B) Predicted pairing between pmi-miR408 and L. japonicum Locus_9880. The E-complementarity score between miRNA and target RNA as estimated by the psRNATarget program is shown. (DOCX) S1