Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Identification of MicroRNAs in the Coral Stylophora pistillata

Abstract

Coral reefs are major contributors to marine biodiversity. However, they are in rapid decline due to global environmental changes such as rising sea surface temperatures, ocean acidification, and pollution. Genomic and transcriptomic analyses have broadened our understanding of coral biology, but a study of the microRNA (miRNA) repertoire of corals is missing. miRNAs constitute a class of small non-coding RNAs of ∼22 nt in size that play crucial roles in development, metabolism, and stress response in plants and animals alike. In this study, we examined the coral Stylophora pistillata for the presence of miRNAs and the corresponding core protein machinery required for their processing and function. Based on small RNA sequencing, we present evidence for 31 bona fide microRNAs, 5 of which (miR-100, miR-2022, miR-2023, miR-2030, and miR-2036) are conserved in other metazoans. Homologues of Argonaute, Piwi, Dicer, Drosha, Pasha, and HEN1 were identified in the transcriptome of S. pistillata based on strong sequence conservation with known RNAi proteins, with additional support derived from phylogenetic trees. Examination of putative miRNA gene targets indicates potential roles in development, metabolism, immunity, and biomineralisation for several of the microRNAs. Here, we present first evidence of a functional RNAi machinery and five conserved miRNAs in S. pistillata, implying that miRNAs play a role in organismal biology of scleractinian corals. Analysis of predicted miRNA target genes in S. pistillata suggests potential roles of miRNAs in symbiosis and coral calcification. Given the importance of miRNAs in regulating gene expression in other metazoans, further expression analyses of small non-coding RNAs in transcriptional studies of corals should be informative about miRNA-affected processes and pathways.

Citation: Liew YJ, Aranda M, Carr A, Baumgarten S, Zoccola D, Tambutté S, et al. (2014) Identification of MicroRNAs in the Coral Stylophora pistillata. PLoS ONE 9(3): e91101. https://doi.org/10.1371/journal.pone.0091101

Editor: Mónica Medina, Pennsylvania State University, United States of America

Received: November 6, 2013; Accepted: February 6, 2014; Published: March 21, 2014

Copyright: © 2014 Liew et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Funding: This project was partially funded by bursaries from Cambridge Commonwealth Trust and Trinity Hall to YJL, and an Academic Excellence Alliance (AEA) Award (Award Number 1000000533) to CRV and GM. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Scientific curiosity about corals has recently intensified, following observations of the deterioration of coral reefs at an unprecedented rate worldwide — for instance, in the Caribbean, Hughes [1] reported that coral cover has declined from over 50% in the 1970s to less than 5% in the 1990s; in the Indo-Pacific region, home to 75% of the world's coral reefs, Bruno and Selig [2] estimated that coral cover declined ∼1% annually in the past 20 years, and ∼2% annually between 1997–2003. This trend is worrying, as coral reefs are important ecosystems, supporting more marine biodiversity per unit area than any other marine habitat [3]. There are many reasons behind the global decline of coral reefs, which include, but are not limited to, accelerated warming and acidification of oceans [4], [5], overfishing [1], pollution [6], [7], and disease [8].

In recent years, the increasing use of genomics has broadened our understanding of basic coral biology. The genome sequence of the coral Acropora digitifera [9] revealed a potential dependency of some coral species on their symbiont population for synthesis of an essential amino acid, and highlighted an unexpectedly diverse repertoire of immune-response genes [9]. Furthermore, microarray and RNA sequencing studies on several coral species have shed light on their responses to environmental cues at the transcriptional level. Shifts in transcriptional landscapes have been noted, based on the composition of symbionts in the coral cell [10], [11], or as a response to stressors such as increased temperatures [12]–[15]; long-term darkness [16]; elevated CO2 levels [17], [18], and ultraviolet radiation [19]. Despite the increasing accumulation of genomic data, some aspects of the molecular machinery potentially involved in these processes, such as microRNAs (miRNAs), have yet to be studied in corals.

miRNAs are a class of small non-coding RNAs of ∼22 nucleotides (nt) in length, which regulate gene expression through posttranscriptional degradation or translational repression via the RNA interference pathway (RNAi) [20]–[22]. Recent studies in plants and metazoans have discovered pivotal roles for miRNAs in regulating developmental timing [23]–[25]; cell cycle progression [26], [27]; immune response [28], [29]; metabolism [30]; response to stress [31]–[33]; and potentially biomineralisation [34]–[36]. miRNAs have been identified in more than 200 species that span major kingdoms of life: animals, plants, and protists (based on miRBase v20, June 2013) [37]–[40]. miRNAs have also been identified in the genome and transcriptome of the coral symbiont Symbiodinium microadriaticum [41] as well as in the genomes of two other cnidarians: Chapman et al. [42] reported 17 miRNAs for Hydra magnipapillata, while Grimson et al. [43] reported 40 miRNAs in the sea anemone Nematostella vectensis. The large evolutionary distance from Hydra and Nematostella to corals (∼500 million years [9]) warranted a search for the presence of miRNAs and the corresponding RNAi machinery in scleractinian corals. Here we present a first assessment of the miRNA repertoire, the RNAi machinery, and putative gene targets in the scleractinian coral S. pistillata from the Red Sea.

Materials and Methods

Ethics statement

Corals were kept in accordance with recommendations by the Centre Scientifique de Monaco and appropriate guidelines for the Care and Use of Laboratory Animals (Permit Number DCI/89/32).

Growth conditions of S. pistillata

S. pistillata was maintained in aquaria at the Centre Scientifique de Monaco, Principality of Monaco, in controlled culture conditions: semi-open circuit, Mediterranean seawater heated to 25±0.5°C, salinity of 38.2 psu, illuminated with HQI-10000K; BLV-Nepturion at a constant irradiance of 175 µmol photons m−2 s−1 on a 12 h∶12 h day∶night light cycle. Corals were fed three times a week with a mix of Artemia salina nauplii and A. salina frozen adults and frozen krill.

mRNA sequencing and transcriptome

Total RNA extraction was performed as described previously [44]. Briefly, nubbins of coral (sampled at noon) were snap-frozen in liquid nitrogen and ground into powder in a cryogrinder (Freezer/Mill 6770, Spex Sample Prep) and then extracted with TRIzol Reagent (Invitrogen, Carlsbad, CA) according to manufacturer's instructions. Total RNA was quality-checked using a Bioanalyzer 2100 (Agilent, Santa Clara, CA) and a Nanodrop 2000c (ThermoScientific, Wilmington, DE) prior to library creation and sequencing by the KAUST Bioscience Core lab. For mRNA sequencing, paired-end reads for Illumina sequencing were generated from oligo-dT selected total RNA using the Illumina TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA) according to manufacturer's instructions. A total of 152,552,099 paired-end read pairs (read length: 101 bp, insert size: 175 bp) were sequenced on the HiSeq 2000 platform (Illumina, San Diego, CA).

For the transcriptome assembly, sequence adaptors were trimmed from the raw sequences and low quality ends were cut with trimmomatic [45]. The remaining read pairs were subjected to digital normalization with Diginorm at k = 20 and C = 20 [46], reducing the dataset to 51,023,864 read pairs. Further, in order to remove contaminating sequence information from endosymbiotic dinoflagellates, remaining read pairs were mapped to the transcriptome of Symbiodinium microadriaticum [41] using Bowtie 2 [47]. This resulted in 38% of the remaining read pairs being mapped to the S. microadriaticum transcriptome, and a significant reduction in potential chimeric locus assemblies for the remaining 16,555,086 read pairs.

The transcriptome was assembled with Oases [48] using k-mer values ranging from 29 to 69. To reduce redundancy within single k-mer assemblies, only contigs with a minimum coverage of 7 were reported. Based on contig lengths, number of distinct loci, and number of transcripts, single k-mer assemblies from k = 45, 47, 49 were reassembled at k = 27, resulting in a final transcriptome assembly of 43,493 unique loci/genes ≥250 bp (Supporting Information S1). For gene annotation, the longest transcript per loci was subjected to a BLASTX search (minimum e-value threshold of 10−5) against three protein databases: UniProtKB/Swiss-Prot, UniProtKB/TrEMBL, and the non-redundant GenBank nr in a subsequent manner. Hits were selected preferentially from Swiss-Prot as Swiss-Prot is manually curated, followed by TrEMBL if no matches were found in Swiss-Prot, and lastly from nr if neither Swiss-Prot nor TrEMBL yielded hits. Out of the 20,332 transcripts with an annotation, 15,177 (77.6%) were from Swiss-Prot; 4,964 (24.4%) were from TrEMBL; and 193 (0.93%) were from nr (Supporting Information S2). For transcripts with annotations from Swiss-Prot or TrEMBL, a script was written to assign GO (Gene Ontology) terms (and their parent GO terms) from UniProt-GOA [49]. 14,558 (95.9%) of the Swiss-Prot hits and 1,955 (39.4%) of the TrEMBL hits had at least 1 GO term assigned to it (Supporting Information S2).

Identification of core RNAi proteins

In order to identify homologues of the RNAi machinery in S. pistillata, sequences from six families of proteins (Argonaute, Dicer, Piwi, Drosha, Pasha, and HEN1) were drawn from five organisms (H. sapiens, D. melanogaster, C. elegans, S. pombe, and A. thaliana). These sequences were obtained from the UniRef100 database [50], and clustered into groups with 90% sequence identity using CD-HIT [51] to remove near-identical sequences. The clustered sequences were used in a TBLASTN search against the S. pistillata transcriptome to identify candidate RNAi-related transcripts. Identified homologues (TBLASTN e-value<10−10) of known RNAi proteins were then searched for domains that are required for the function of the protein using InterProScan [52]–[54]. The domains that were determined to be essential for function were: a Paz and Piwi domain for Argonaute and Piwi; a pair of RNase III domains for Dicer and Drosha; a double-stranded RNA binding domain for Pasha; and a methyltransferase (MTase) domain for HEN1. Candidate homologues were not considered further in the absence of any of these domains.

Additional support for the inferred function of candidate homologues was obtained by carrying out a reciprocal BLASTP search of these translated candidates against all proteins in the Swiss-Prot database [55] (Supporting Information S3). The candidate homologues were aligned against known RNAi proteins on a per-family basis using Clustal Omega [56], and the alignments were visualised using Jalview [57], [58]. Key residues were derived from literature [59]–[64]. In addition, for each of the six protein families, phylogenies were constructed by aligning our candidate homologues with selected sequences from Grimson et al. [43] and Moran et al. [65] with MUSCLE [66]. Aligned regions of low quality were removed with TrimAl, using the in-built “gappyout” parameter [67] (Supporting Information S4). ProtTest3 [68] was used to determine the best model for amino acid substitution, and MEGA (version 6) [69] was used to construct maximum-likelihood phylogenetic trees (support values were computed from 1,000 bootstrap replicates).

Small RNA sequencing and processing

The small RNA (smRNA) fraction from S. pistillata was selectively enriched from isolated total RNA (see above) using the mirVana miRNA isolation kit (Ambion, Austin, TX) according to manufacturer's instructions. The small RNA fraction was quality-checked using a Bioanalyzer 2100 (Agilent, Santa Clara, CA) and a Nanodrop 2000c (ThermoScientific, Wilmington, DE). The small RNA library was created using the Illumina Small RNA Sample Prep Kit (Illumina, San Diego, CA) according to manufacturer's instructions, and sequenced on 1 lane on an Illumina Genome Analyzer IIx (GA2x) machine. A total of 30.5 million small RNA reads of <40 bp in length were produced. The reads, along with associated Phred quality scores for each sequenced base, were saved in a FASTQ file.

The raw FASTQ file was processed using several scripts to remove low-quality reads resulting in a more compact FASTQ file that contained high-quality reads for downstream analyses. First, low quality 3′ ends were trimmed from the reads. The 3′ end of the resulting reads had a Phred score >20, while the average Phred score of the entire read was >20 as well. After trimming, the overall error rate of the reads was calculated from the Phred scores of individual bases. Reads were discarded if the error rate exceeded 2%. Subsequently, the Illumina 5′ and 3′ adapter sequences used in library generation were trimmed off using Cutadapt v1.0 [70]. Last, in order to remove fragments of rRNA, tRNA, and mRNA sequences, Velvet [71] was used to assemble the short reads into contigs (at k = 25), which were then compared to the GenBank nt database (nucleotide collection at NCBI). In addition, we compared the assembled contigs to the S. pistillata transcriptome assembly using BLASTN, in order to remove short reads that matched known mRNA sequences.

miRNA prediction and filtering

We used miRDeep2 [72] to identify miRNAs. Briefly, miRDeep2 mapped smRNAs to a preliminary draft genome of S. pistillata using Bowtie, discarding reads that occurred more than five times in the genome to avoid mapping to repetitive elements. Potential pre-miRNA precursor sequences were identified based on the pattern of the mappings, and subsequently folded using RNAfold to ascertain whether they had the canonical hairpin secondary structure [72]. Predicted pre-miRNAs that had a miRDeep2 score of 10 or above were retained for further analyses and inspected manually. A script was written to produce additional information not found in the miRDeep2 output (i.e. length of 3′ overhang, proportion of reads with consistent 5′ end, number of mismatched bases in stem) to further select a set of bona fide miRNAs. Conserved miRNAs were identified using BLASTN against all previously identified pre-miRNA sequences in miRBase (ver 20) [37]–[40].

Functional analysis of miRNA targets

ORFs were identified in the transcripts using TransDecoder (part of the Trinity software pipeline [73]). Sequences downstream of the longest ORF identified in the transcripts were treated as the 3′ UTR of the transcript. 3′ UTRs under 100 bp were filtered out to remove transcripts associated with short UTR sequences. Out of the 43,493 genes, 14,125 transcripts (32.4%) had a predicted 3′ UTR longer than 100 bp.

For each miRNA, we ran PITA (Probability of Interaction by Target Accessibility) [74] on the 3′ UTRs at default settings to produce a set of putative genes targeted by the miRNA. In the absence of genomic data from other closely related organisms, PITA achieves higher sensitivity and specificity than other target prediction software (e.g. miRanda, TargetScan) as the latter algorithms rely on a filter based on evolutionary conservation to reduce the false positive rate. PITA works by calculating the difference in Gibbs' free energy (ΔΔG) between the energy that is required to unfold the secondary structure of the target site (ΔGopen), and the energy of the mature miRNA binding its target (ΔGduplex) [74]. Only miRNA targets with a ΔΔG of ≤−10 kcal mol−1 were retained.

For GO enrichment of target genes, we used topGO (version 2.12.0), an R script that is available through Bioconductor 2.0. topGO is a scoring algorithm that improves GO scoring by eliminating local dependencies between related GO terms [75]. The threshold for significance was set at P<0.01, using otherwise default topGO “weight01” settings, which produced GO terms that were significantly enriched in the set of transcripts targeted by each miRNA. The resulting P values were not corrected for multiple testing, as non-independent tests are carried out on each GO term by topGO [75].

Results

Identification of core RNAi proteins

The miRNA machinery that processes and mediates the function of miRNAs encompasses several key components that appear to be conserved across the animal kingdom [76]. In order to establish the presence of a functional miRNA machinery in S. pistillata we conducted a BLAST-based search for key proteins known to be essential for miRNA processing and function.

We identified seven candidate genes that are homologues to known RNAi proteins: one Argonaute, two Piwi, one Dicer, one Drosha, one Pasha, and one HEN1 in S. pistillata. We employed several key metrics (i.e. matches to known RNAi families, presence of protein domains crucial for catalytic activity, and a reciprocal BLAST search against manually curated proteins in Swiss-Prot) to identify candidate RNAi proteins (Supporting Information S3).

The per-family alignments of candidate homologues against known sequences revealed a striking conservation of functionally important amino acid residues located within the key protein domains. Examples include the strong conservation of the DDX triad in the Piwi domain of the Argonaute and Piwi homologues; the aspartate and glutamate residues essential for Dicer activity; and the pair of alanine/alanine and alanine/serine dipeptides involved in the binding of dsRNA in Pasha (Supporting Information S5, S6, S7, S8, S9, S10). Maximum-likelihood phylogenetic trees that were constructed for all six protein families (Figures 1A–1F) placed all of the candidate S. pistillata homologues with those from other cnidarians. Judging from the presence of the key RNAi proteins in S. pistillata in comparison to other organisms, the RNAi machinery in S. pistillata is similar in composition to those from sea anemone, worm, fruit fly, and humans (Table 1).

thumbnail
Download:
Figure 1. Maximum-likelihood phylogenies of (A) Argonaute, (B) Piwi, (C) Dicer, (D) Drosha, (E) Pasha, and (F) HEN1 proteins.

(A), (C), and (D) were constructed using the amino acid substitution model LG+G, while (B), (E), and (F) were constructed using LG+I+G. Bootstrap support values are indicated above the branches. Abbreviations used in these trees are ‘Adi’: Acropora digitifera; ‘Ami’: Acropora millepora; ‘Aqu’: Amphimedon queenslandica (a sponge); ‘Cel’: C. elegans; ‘Dme’: D. melanogaster; ‘Hsa’: H. sapiens; ‘Hvu’: Hydra vulgaris; ‘Nve’: N. vectensis; and ‘Spi’: S. pistillata.

https://doi.org/10.1371/journal.pone.0091101.g001

thumbnail
Download:
Table 1. Presence of RNAi proteins in S. pistillata in comparison to Cnidaria and other model organisms (‘+’: presence, ‘−’: absence, ‘?’: not determined).

https://doi.org/10.1371/journal.pone.0091101.t001

Besides the core RNAi proteins, we have also discovered transcripts that are candidate homologues of HYL1 (one), GW182 (two), and RdRP (RNA-dependent RNA polymerase, eight) (data not shown). HYL1 is thought to be a plant-specific partner to Dicer [77], whereas GW182 helps Argonaute repress its targets [78]. Both proteins have recently been discovered in four cnidarians (Acropora digitifera, A. millepora, Hydra vulgaris, Nematostella vectensis). However, although we could identify the PAM2 and P-GL motif in one of our GW182 homologues, there were very few GW repeats in this homologue (1 in S. pistillata, compared to 14 in N. vectensis and 40 in humans) [65]. RdRPs, using small RNAs as templates, amplify the silencing effect by directing the production of secondary dsRNAs [79]. Functional RdRPs have been discovered in plants [25], [80] and C. elegans [81], but not in mammals nor flies [79]. Four candidate homologues of RdRPs have been found in N. vectensis [82], indicating that RdRPs might be present in cnidarians.

Small RNA sequencing and miRNA repertoire

Sequencing produced 30,543,433 reads, of which 23,830,932 reads (78.0%) were kept after adapter trimming. The additional step of removing short reads that matched known rRNA, tRNA, and transcript sequences removed a further 7.2% of reads, resulting in a dataset that contained 21,625,195 reads of at least 10 bp in length. Most of these reads were of 20–31 nt in length. Relative frequencies of the starting 5′ nucleotide showed a clear enrichment of 5′-terminal uracil among short reads of length 26–31 nt (Figure 2), which is consistent with the likely presence of functional Piwi-interacting RNAs (piRNAs) in S. pistillata [43].

thumbnail
Download:
Figure 2. Frequency distribution of small RNA reads in the S. pistillata sequencing library.

The bars are coloured to reflect the proportions of reads starting with A, U, G, and C (blue, red, green, and purple respectively).

https://doi.org/10.1371/journal.pone.0091101.g002

A total of 2,811,736 reads (of length >17 nt) were mapped to a preliminary assembly of the S. pistillata genome and 46 distinct miRNAs were predicted by miRDeep2, of which a subset of 30 miRNAs passed additional filtering criteria (see Materials and Methods). An exception was made for spi-miR-temp-25 – although the precursor had a 3 bp 3′ overhang, it was included in the bona fide set as it was a close match to two known miRNAs. The resulting 31 miRNAs were used in all downstream analyses (Table 2, Supporting Information S11).

thumbnail
Download:
Table 2. Set of 31 bona fide miRDeep2-predicted miRNAs in S. pistillata.

https://doi.org/10.1371/journal.pone.0091101.t002

While the majority of these 31 predicted miRNAs were novel, 5 of them matched conserved miRNAs. spi-miR-temp-1, the predicted miRNA with the highest miRDeep2 score, was highly similar to the known miR-100 family of sequences (∼2 mismatched bases, Figure 3A). This family is known to be conserved across Eumetazoa, including at least two other cnidarians (N. vectensis and Metridium senile) [43], [83]–[87]. spi-miR-temp-25 was similar to miR-2022 from N. vectensis [43] and H. magnipapillata [88] (Figure 3B), while the other three miRNAs – spi-miR-temp-4, spi-miR-temp-40, and spi-miR-temp-30 – were similar to miR-2023, miR-2030, and miR-2036 found in N. vectensis [43], respectively (Figures 3C–E). For reasons of clarity, these five conserved miRNAs in S. pistillata will be referred to as spi-miR-100, spi-miR-2022, spi-miR-2023, spi-miR-2030, and spi-miR-2036, in accordance with miRNA naming conventions.

thumbnail
Download:
Figure 3. Alignments of predicted S. pistillata miRNAs against (A) members of the miR-100 family; (B) nve- and hma-miR-2022; (C) nve-miR-2023; (D) nve-miR-2030; and (E) nve-miR-2036.

The mature sequences are shown on the left, while star sequences are on the right. Sequences were obtained from miRBase (version 20). The mature hma-miR-2030 aligned best with miR-2030* sequences from N. vectensis and S. pistillata. Sequences marked with a tilde (nve-miR-2022*, hma-miR-2022*, and hma-miR-2030) are miRNAs that we derived based on the alignment of the respective pre-miRNA sequences obtained from miRBase against S. pistillata miRNAs. Bases were coloured to provide visual indication of conservation (dark blue: >80%; blue: >60%; light blue: >40%; uncoloured otherwise). Abbreviations used are ‘dme’: D. melanogaster; ‘hma’: H. magnipapillata; ‘hsa’: H. sapiens; ‘nve’: N. vectensis; and ‘spi’: S. pistillata.

https://doi.org/10.1371/journal.pone.0091101.g003

Although nve-miR-100 has been identified in two separate studies (both of which utilised next-generation sequencing of short reads to identify miRNAs in basal metazoans), Grimson et al. [43] and Wheeler et al. [88] predicted mature miR-100 sequences which are offset by a single nucleotide. From our read data, we had >30,000 reads that exactly match the nve-miR-100 from Grimson et al. [43], but none that matched the alternative version from Wheeler et al. [88]. A negligible minority of reads (∼10) did start one nucleotide upstream, but unlike the Wheeler et al. [88] version, this residue was an adenine, making this form identical to the hsa-, dme- and xtr-miR-100s.

Functional analysis of putative miRNA targets

In order to identify processes in S. pistillata that are potentially regulated by miRNAs, we conducted a GO term enrichment analysis. Briefly, for the set of 31 miRNAs, we searched for animal-like target sites in the 3′ UTRs of those 16,513 S. pistillata genes that had available 3′ UTR and GO annotation. This analysis was performed for all miRNAs individually, and indicated that miRNAs are likely to play roles in diverse processes (Table 3, Supporting Information S12).

thumbnail
Download:
Table 3. Enriched immunity- and biomineralisation-related GO terms (default topGO settings, P<0.01) associated with predicted miRNA target genes from S. pistillata.

https://doi.org/10.1371/journal.pone.0091101.t003

We categorised the resulting enriched GO terms under several high-level groups based on known miRNA function – “immunity”, “biomineralisation”, “transcription”, “cell cycle”, “cytoskeleton”, “metabolism”, “transport/signalling”, and “differentiation/development”. Other GO terms that did not fall under one of these umbrella terms were categorised under “miscellaneous” (Supporting Information S12).

We paid particular attention to the “immunity” and “biomineralisation” high-level groups, as we considered these terms to be specifically relevant to the understanding of coral physiology (Table 3). For the former, it is likely that immunity-related transcripts are involved in the formation and retention of symbiotic relationships between the coral host and its Symbiodinium symbionts. Five miRNAs were predicted to be involved in this regulation, with a large fraction of the predictions involving spi-miR-temp-15. For the latter, 10 miRNAs were predicted to be involved in the formation of the calcium carbonate-based coral skeleton, with none of the 10 miRNAs being predominant in the predictions.

Discussion

Identification of core RNAi proteins

Advances in our understanding of miRNAs have shown that these small molecules have a big impact on the regulation of gene expression. While the biogenesis and downstream functions of miRNAs have been fairly well studied in the primary model organisms, little is known about the presence or function of these miRNAs in corals. In this study, we identified the presence of core RNAi proteins encoded by the S. pistillata transcriptome. The alignment of our candidate homologues against known RNAi proteins revealed the conservation of key protein domains and residues crucial for protein function. Additionally, maximum-likelihood phylogenetic trees of our candidate homologues with RNAi proteins from other cnidarian and model organisms showed broad agreement with those from other studies [43], [65] – all of our seven candidate homologues clustered best with those from other cnidarians, as expected from its closer phylogenetic distance to other cnidarians than to bilaterians or poriferans. This also signifies that our homologues originate from the coral host, not from its dinoflagellate symbionts. Interestingly, the S. pistillata Dicer homologue clustered better with N. vectensis Dcr2, which is thought to be involved in processing of long dsRNA into siRNAs, and not associated with the biogenesis of miRNAs [65], [89]. A reverse search of our candidate Dicer homologue against the S. pistillata draft genome revealed an open reading frame that encodes for another Dicer-like protein, which appears to be a good match (e-value of <1×10−10) of N. vectensis Dcr1 (data not shown). However, the absence of transcriptomic support for that open reading frame excluded it from being a candidate Dicer in S. pistillata in this study. Nonetheless, both observations serve to indicate the presence of a functional miRNA-processing machinery in S. pistillata. This, to our knowledge, has not been demonstrated previously for any other coral.

Small RNA sequencing and miRNA repertoire

Besides a functional RNAi machinery, and based on our analysis of short reads, we also predicted the presence of 31 bona fide miRNAs (out of a total of 46), of which 5 were conserved: the miR-100 family found in many other metazoans; miR-2022, which is conserved in N. vectensis and H. magnipapillata; miR-2023, miR-2030, and miR-2036, which are conserved in N. vectensis only. The dearth of conserved Hydra miRNAs in S. pistillata echoes the findings of Chapman et al. [42], [43], who found only one conserved N. vectensis miRNA among the H. magnipapillata miRNAs. This might be due to the evolutionary distance separating the anthozoans and hydrozoans, or, more likely, due to the incomplete coverage of short reads used in the identification of miRNAs in H. magnipapillata – only 9,654 reads were used to identify potential miRNA genes in H. magnipapillata [42]. In contrast, we (and Grimson et al. [43]) identified miRNAs from a much larger pool of short reads. We believe that the repertoire of miRNAs that are conserved across both cnidarian classes (i.e. Anthozoa and Hydrozoa) could be expanded if miRNA predictions were ran on a larger pool of small RNA reads.

The conservation of miRNA families across and within different bilaterian phyla have been fairly well-covered, with the general consensus that the continuous acquisition of miRNA families with minimal secondary losses rapidly expanded the bilaterian miRNA repertoire relative to cnidarians, which contributes to the increased morphological complexity of bilaterians [83], [88], [90]–[92]. As one of the few cnidarians with its small RNA fraction extensively sequenced, S. pistillata has demonstrated that conservation of miRNA families does occur within cnidarians too, as five of its miRNAs are conserved in N. vectensis despite the ∼500 mya evolutionary distance that separate both species. However, due to the dearth of sequenced small RNA reads from other cnidarians, we are unable to make further conclusions regarding the rate at which cnidarians acquire their own phylum-specific miRNA families. Also, recent evidence has surfaced that demonstrated the gradual loss of conserved (up to 50% in more derived species) and gain of novel miRNA families in Platyhelminthes, the first that was reported for a major lineage within Bilateria, and might be related to morphological simplifications in some of the studied flatworms [93]. Similar observations could apply to specific classes of cnidarians, but this type of study would need to include more than just a few species of cnidarians in order to elucidate the true rate underlying the gains and losses of miRNA families.

Functional analysis of putative miRNA targets

Functional analysis of all 31 miRNAs, using target predictions for each miRNA followed by a GO enrichment analysis on the predicted target genes, revealed several putative processes and pathways that are regulated by miRNAs in corals. For the miR-100 homologue in S. pistillata, the GO terms “embryonic forelimb morphogenesis” and “bone development” were enriched (P<0.01, Supporting Information S12) in the predicted targets, which is reminiscent of its reported function: in humans, miR-100 has been shown to target genes involved in growth and development. Examples include Plk1, a key mitotic checkpoint regulatory protein [26]; RBSP3, involved in cell proliferation and myeloid cell differentiation [27]; BMPR2, involved in osteogenesis [94]; and FRAP1/mTOR, which regulates cell growth [95]. It is possible that miR-100 plays an analogous role in coral calcification, making this miRNA a potentially important piece of the puzzle in coral physiology, as well as a gene of interest when investigating coral responses to ocean acidification. However, as miRNA-mRNA target recognition depends critically on the miRNA seed sequence (bases 2–7 of the mature RNA), it is possible that the targets of bilaterian and cnidarian miR-100 will differ due to the one nucleotide offset between the two miRNA sequences. This 5′ offset has also been observed for miR-2, miR-10, miR-133, and miR-210 that are otherwise well-conserved across two phylogenetically-related taxa, and presumably able to regulate non-overlapping sets of target mRNAs [91]. Thus, further experimentation is required to confirm the bona fide function of cnidarian miR-100 in corals. Nonetheless, our spi-miR-100 adds to the existing literature documenting the strong conservation of miR-100 within metazoans.

Besides the only miRNA with documented function, we identified miRNAs whose targets are involved in high-level functions such as immunity, biomineralisation, regulation of cell cycle, cellular motility, metabolism, signalling, and development, analogous to functions that were previously ascribed to miRNAs in other organisms [23]–[36]. We were interested in the first two high-level groups, as immunity genes might regulate the relationship with the symbiotic dinoflagellate Symbiodinium, and biomineralisation genes may control the rate of coral skeleton growth, two processes that are arguably of importance to corals under conditions of environmental change.

Out of the 5 miRNAs that were predicted to regulate coral immunity genes, we speculate that spi-miR-temp-15 should warrant further investigation due to the significant enrichment of multiple immunity-related GO terms in the transcripts targeted by this miRNA. Indeed, several of the predicted target genes of spi-miR-temp-15 have homologues that are known to be regulated by other miRNAs: Nod2 is repressed by miR-122 [96]; TLR2 is regulated by miR-19 and miR-105 [97], [98]; while caspase-8 is targeted by miR-874 [99]. Interestingly, this miRNA is not conserved in N. vectensis, which does not form long-term symbiotic relationships with Symbiodinium.

In contrast to the previous category, 10 miRNAs were predicted to have roles in biomineralisation – one of which being miR-100, which regulates growth and development in humans [26], [27], [94], [95]. Further, among the targeted transcripts, we found several transcripts which are predicted homologues of genes involved in calcium and bicarbonate ion transport that are directly regulated by miRNAs (miR-506 targets human anion exchange protein 2 [100], while miR-17 targets polycystin-2 [101]). A potential involvement of miRNAs in regulating ion transport is intriguing, given the significance of these processes in relation to ocean acidification and associated consequences to coral calcification [102]. However, future experiments (e.g in-situ hybridisations, gene expression assays, or immunoprecipitation studies) are essential in unequivocally verifying these predicted interactions.

In conclusion, our study provides strong support for the presence of a functional RNAi machinery in S. pistillata as highlighted by our phylogenetic analyses, the strong conservation of key RNAi protein domains, and the presence of conserved miRNAs. miRNAs seem to affect a variety of biological processes in corals, but further studies that focus on the coordinated expression of miRNAs and associated target mRNAs under different conditions, as well as their interaction with RNAi proteins, are needed in order to identify, characterise, and understand the operational miRNAome in scleractinian corals.

Data access

All small RNA and RNASeq data are available in the NCBI SRA (Sequence Read Archive) under accessions SRR1130519 and SRR1125978 respectively. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GARY00000000. The version described in this paper is the first version, GARY01000000. Names of the miRNAs are temporary, as miRBase (the microRNA registry) requires acceptance of manuscripts prior to assigning names to newly identified miRNAs. Other data are available as Supporting Information.

Supporting Information

Supporting Information S1.

Stylophora pistillata transcriptome (43,493 genes/loci ≥250 bp, DDBJ/EMBL/GenBank accession GARY00000000).

https://doi.org/10.1371/journal.pone.0091101.s001

(ZIP)

Supporting Information S2.

Stylophora pistillata transcriptome BLASTX and GO annotation (43,493 genes/loci ≥250 bp).

https://doi.org/10.1371/journal.pone.0091101.s002

(ZIP)

Supporting Information S3.

Candidate RNAi proteins in Stylophora pistillata.

https://doi.org/10.1371/journal.pone.0091101.s003

(DOCX)

Supporting Information S4.

Alignment of sequences used to construct maximum-likelihood phylogenetic trees (FASTA format).

https://doi.org/10.1371/journal.pone.0091101.s004

(FA)

Supporting Information S5.

Graphical alignment of the PAZ domains in Argonaute and Piwi proteins. Of note are the strong conservation of glutamate (E) at position 137 (mutants produce insoluble protein) and phenylalanine (F) at position 72 (required for RNA binding). However, the phenylalanine at position 48 in D. melanogaster AGO2 (also required for RNA binding) was not conserved at all. Key residue positions are marked with red asterisks.

https://doi.org/10.1371/journal.pone.0091101.s005

(EPS)

Supporting Information S6.

Graphical alignment of the Piwi domains in Argonaute and Piwi proteins. The catalytic DDX triad, which contributes to the slicing activity of the ribonuclease (marked in red asterisks), is located at positions 46, 140 and 284 or positions 46, 140 and 155. This triad is present in one S. pistillata candidate, but not in two others, most likely due to the transcript sequences being incomplete.

https://doi.org/10.1371/journal.pone.0091101.s006

(EPS)

Supporting Information S7.

Graphical alignment of the first RNase III domain in Dicer and Drosha proteins. Remarkably, all of the key acidic aspartate (D) and glutamate (E) residues, which are involved in the coordination of a divalent metal cation, are conserved across the candidate homologues and known sequences.

https://doi.org/10.1371/journal.pone.0091101.s007

(EPS)

Supporting Information S8.

Graphical alignment of the second RNase III domain in Dicer and Drosha proteins. Similarly, most of the aspartate (D) and glutamate (E) residues involved in the coordination of a divalent metal cation are conserved - perfectly conserved for the Drosha candidate (“Locus_18820”), while the Dicer candidate (“Locus_10081”) only has the first two key residues. Both sequences however align well to known Dicer and Drosha proteins.

https://doi.org/10.1371/journal.pone.0091101.s008

(EPS)

Supporting Information S9.

Graphical alignment of the dsRNA-binding domain in Pasha. The key alanine/alanine pair (AA, positions 21 and 22) and alanine/serine pair (AS, positions 139 and 140) involved in the binding of dsRNA are also present in the S. pistillata candidate Pasha. As Pasha is an essential cofactor of Drosha, it lends support to the positive discovery of Drosha in S. pistillata.

https://doi.org/10.1371/journal.pone.0091101.s009

(EPS)

Supporting Information S10.

Graphical alignment of the methyltransferase domain in HEN1. The residues involved in Mg2+ coordination (positions 118, 121, 122 and 123) are well-conserved across the aligned sequences; residues associated with the cofactor AdoHcy and 3′ terminus (other positions marked by a red asterisk) are also well conserved.

https://doi.org/10.1371/journal.pone.0091101.s010

(EPS)

Supporting Information S11.

List of additional criteria used to select bona fide miRNAs in S. pistillata from miRDeep2 results.

https://doi.org/10.1371/journal.pone.0091101.s011

(XLSX)

Supporting Information S12.

Enriched GO terms (P<0.01) associated with the set of 31 bona fide miRNAs identified in Stylophora pistillata.

https://doi.org/10.1371/journal.pone.0091101.s012

(XLSX)

Acknowledgments

We thank the KAUST Bioscience Core Lab and S. Ali for Illumina library generation and sequencing, and Till Bayer and Mike Lyne for helpful discussions in regard to data analyses.

Author Contributions

Conceived and designed the experiments: CRV GM MA. Performed the experiments: YJL MA AC SB DZ. Analyzed the data: YJL MA CRV. Contributed reagents/materials/analysis tools: DZ ST DA. Wrote the paper: YJL MA CRV.

References

  1. 1. Hughes TP (1994) Catastrophes, phase shifts, and large-scale degradation of a Caribbean coral reef. Science 265: 1547–1551.
  2. 2. Bruno JF, Selig ER (2007) Regional decline of coral cover in the Indo-Pacific: timing, extent, and subregional comparisons. PLoS One 2: e711.
  3. 3. Knowlton N, Brainard RE, Fisher R, Moews M, Plaisance L, et al. (2010) Coral reef biodiversity. Life in the World's Oceans: Diversity Distribution and Abundance 65–74.
  4. 4. Hughes TP, Baird AH, Bellwood DR, Card M, Connolly SR, et al. (2003) Climate change, human impacts, and the resilience of coral reefs. Science 301: 929–933.
  5. 5. Hoegh-Guldberg O, Mumby PJ, Hooten AJ, Steneck RS, Greenfield P, et al. (2007) Coral reefs under rapid climate change and ocean acidification. Science 318: 1737–1742.
  6. 6. Bak R (1987) Effects of chronic oil pollution on a Caribbean coral reef. Marine Pollution Bulletin 18: 534–539.
  7. 7. Pastorok RA, Bilyard GR (1985) Effects of sewage pollution on coral-reef communities. Marine ecology progress series Oldendorf 21: 175–189.
  8. 8. Green EP, Bruckner AW (2000) The significance of coral disease epizootiology for coral reef conservation. Biological Conservation 96: 347–361.
  9. 9. Shinzato C, Shoguchi E, Kawashima T, Hamada M, Hisata K, et al. (2011) Using the Acropora digitifera genome to understand coral responses to environmental change. Nature 476: 320–323.
  10. 10. DeSalvo MK, Sunagawa S, Fisher PL, Voolstra CR, Iglesias-Prieto R, et al. (2010) Coral host transcriptomic states are correlated with Symbiodinium genotypes. Mol Ecol 19: 1174–1186.
  11. 11. Voolstra CR, Schwarz JA, Schnetzer J, Sunagawa S, Desalvo MK, et al. (2009) The host transcriptome remains unaltered during the establishment of coral-algal symbioses. Mol Ecol 18: 1823–1833.
  12. 12. DeSalvo MK, Voolstra CR, Sunagawa S, Schwarz JA, Stillman JH, et al. (2008) Differential gene expression during thermal stress and bleaching in the Caribbean coral Montastraea faveolata. Mol Ecol 17: 3952–3971.
  13. 13. Polato NR, Voolstra CR, Schnetzer J, DeSalvo MK, Randall CJ, et al. (2010) Location-specific responses to thermal stress in larvae of the reef-building coral Montastraea faveolata. PLoS One 5: e11221.
  14. 14. Voolstra CR, Schnetzer J, Peshkin L, Randall CJ, Szmant AM, et al. (2009) Effects of temperature on gene expression in embryos of the coral Montastraea faveolata. BMC Genomics 10: 627.
  15. 15. Meyer E, Aglyamova GV, Matz MV (2011) Profiling gene expression responses of coral larvae (Acropora millepora) to elevated temperature and settlement inducers using a novel RNA-Seq procedure. Mol Ecol 20: 3599–3616.
  16. 16. DeSalvo M, Estrada A, Sunagawa S, Medina M (2012) Transcriptomic responses to darkness stress point to common coral bleaching mechanisms. Coral Reefs 31: 215–228.
  17. 17. Moya A, Huisman L, Ball EE, Hayward DC, Grasso LC, et al. (2012) Whole transcriptome analysis of the coral Acropora millepora reveals complex responses to CO(2)-driven acidification during the initiation of calcification. Mol Ecol 21: 2440–2454.
  18. 18. Vidal-Dupiol J, Zoccola D, Tambutte E, Grunau C, Cosseau C, et al. (2013) Genes related to ion-transport and energy production are upregulated in response to CO2-driven pH decrease in corals: new insights from transcriptome analysis. PLoS One 8: e58652.
  19. 19. Aranda M, Banaszak AT, Bayer T, Luyten JR, Medina M, et al. (2011) Differential sensitivity of coral larvae to natural levels of ultraviolet radiation during the onset of larval competence. Mol Ecol 20: 2955–2972.
  20. 20. Lagos-Quintana M, Rauhut R, Lendeckel W, Tuschl T (2001) Identification of novel genes coding for small expressed RNAs. Science 294: 853–858.
  21. 21. Lau NC, Lim LP, Weinstein EG, Bartel DP (2001) An abundant class of tiny RNAs with probable regulatory roles in Caenorhabditis elegans. Science 294: 858–862.
  22. 22. Lee RC, Ambros V (2001) An extensive class of small RNAs in Caenorhabditis elegans. Science 294: 862–864.
  23. 23. Pasquinelli AE, Reinhart BJ, Slack F, Martindale MQ, Kuroda MI, et al. (2000) Conservation of the sequence and temporal expression of let-7 heterochronic regulatory RNA. Nature 408: 86–89.
  24. 24. Axtell MJ, Bowman JL (2008) Evolution of plant microRNAs and their targets. Trends Plant Sci 13: 343–349.
  25. 25. Chen X (2009) Small RNAs and their roles in plant development. Annu Rev Cell Dev Biol 25: 21–44.
  26. 26. Li BH, Zhou JS, Ye F, Cheng XD, Zhou CY, et al. (2011) Reduced miR-100 expression in cervical cancer and precursors and its carcinogenic effect through targeting PLK1 protein. Eur J Cancer 47: 2166–2174.
  27. 27. Zheng YS, Zhang H, Zhang XJ, Feng DD, Luo XQ, et al. (2012) MiR-100 regulates cell differentiation and survival by targeting RBSP3, a phosphatase-like tumor suppressor in acute myeloid leukemia. Oncogene 31: 80–92.
  28. 28. Chen CZ, Schaffert S, Fragoso R, Loh C (2013) Regulation of immune responses and tolerance: the microRNA perspective. Immunol Rev 253: 112–128.
  29. 29. O'Connell RM, Rao DS, Chaudhuri AA, Baltimore D (2010) Physiological and pathological roles for microRNAs in the immune system. Nat Rev Immunol 10: 111–122.
  30. 30. Horie T, Ono K, Nishi H, Iwanaga Y, Nagao K, et al. (2009) MicroRNA-133 regulates the expression of GLUT4 by targeting KLF15 and is involved in metabolic control in cardiac myocytes. Biochem Biophys Res Commun 389: 315–320.
  31. 31. Leung AK, Sharp PA (2010) MicroRNA functions in stress responses. Mol Cell 40: 205–215.
  32. 32. Sunkar R, Li YF, Jagadeeswaran G (2012) Functions of microRNAs in plant stress responses. Trends Plant Sci 17: 196–203.
  33. 33. Babenko O, Golubov A, Ilnytskyy Y, Kovalchuk I, Metz GA (2012) Genomic and epigenomic responses to chronic stress involve miRNA-mediated programming. PLoS One 7: e29441.
  34. 34. Cao H, Wang J, Li X, Florez S, Huang Z, et al. (2010) MicroRNAs play a critical role in tooth development. J Dent Res 89: 779–784.
  35. 35. Jiao Y, Zheng Z, Du X, Wang Q, Huang R, et al. (2014) Identification and Characterization of MicroRNAs in Pearl Oyster Pinctada martensii by Solexa Deep Sequencing. Mar Biotechnol (NY) 16: 54–62.
  36. 36. van Wijnen AJ, van de Peppel J, van Leeuwen JP, Lian JB, Stein GS, et al. (2013) MicroRNA functions in osteogenesis and dysfunctions in osteoporosis. Curr Osteoporos Rep 11: 72–82.
  37. 37. Griffiths-Jones S (2004) The microRNA Registry. Nucleic Acids Res 32: D109–111.
  38. 38. Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ (2006) miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res 34: D140–144.
  39. 39. Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ (2008) miRBase: tools for microRNA genomics. Nucleic Acids Res 36: D154–158.
  40. 40. Kozomara A, Griffiths-Jones S (2011) miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res 39: D152–157.
  41. 41. Baumgarten S, Bayer T, Aranda M, Liew YJ, Carr A, et al. (2013) Integrating microRNA and mRNA expression profiling in Symbiodinium microadriaticum, a dinoflagellate symbiont of reef-building corals. BMC Genomics 14: 704.
  42. 42. Chapman JA, Kirkness EF, Simakov O, Hampson SE, Mitros T, et al. (2010) The dynamic genome of Hydra. Nature 464: 592–596.
  43. 43. Grimson A, Srivastava M, Fahey B, Woodcroft BJ, Chiang HR, et al. (2008) Early origins and evolution of microRNAs and Piwi-interacting RNAs in animals. Nature 455: 1193–1197.
  44. 44. Moya A, Tambutte S, Beranger G, Gaume B, Scimeca JC, et al. (2008) Cloning and use of a coral 36B4 gene to study the differential expression of coral genes between light and dark conditions. Mar Biotechnol (NY) 10: 653–663.
  45. 45. Lohse M, Bolger AM, Nagel A, Fernie AR, Lunn JE, et al. (2012) RobiNA: a user-friendly, integrated software solution for RNA-Seq-based transcriptomics. Nucleic Acids Res 40: W622–627.
  46. 46. Brown CT, Howe A, Zhang Q, Pyrkosz AB, Brom TH (2012) A reference-free algorithm for computational normalization of shotgun sequencing data. arXiv preprint arXiv:12034802.
    • 47. Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nat Methods 9: 357–359.
    • 48. Schulz MH, Zerbino DR, Vingron M, Birney E (2012) Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics 28: 1086–1092.
    • 49. Dimmer EC, Huntley RP, Alam-Faruque Y, Sawford T, O'Donovan C, et al. (2012) The UniProt-GO Annotation database in 2011. Nucleic Acids Res 40: D565–570.
    • 50. Suzek BE, Huang H, McGarvey P, Mazumder R, Wu CH (2007) UniRef: comprehensive and non-redundant UniProt reference clusters. Bioinformatics 23: 1282–1288.
    • 51. Li W, Godzik A (2006) Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22: 1658–1659.
    • 52. Mulder N, Apweiler R (2007) InterPro and InterProScan: tools for protein sequence classification and comparison. Methods Mol Biol 396: 59–70.
    • 53. Quevillon E, Silventoinen V, Pillai S, Harte N, Mulder N, et al. (2005) InterProScan: protein domains identifier. Nucleic Acids Res 33: W116–120.
    • 54. Zdobnov EM, Apweiler R (2001) InterProScan–an integration platform for the signature-recognition methods in InterPro. Bioinformatics 17: 847–848.
    • 55. Magrane M, Consortium U (2011) UniProt Knowledgebase: a hub of integrated protein data. Database (Oxford) 2011: bar009.
    • 56. Sievers F, Wilm A, Dineen D, Gibson TJ, Karplus K, et al. (2011) Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol Syst Biol 7: 539.
    • 57. Clamp M, Cuff J, Searle SM, Barton GJ (2004) The Jalview Java alignment editor. Bioinformatics 20: 426–427.
    • 58. Waterhouse AM, Procter JB, Martin DM, Clamp M, Barton GJ (2009) Jalview Version 2–a multiple sequence alignment editor and analysis workbench. Bioinformatics 25: 1189–1191.
    • 59. Huang Y, Ji L, Huang Q, Vassylyev DG, Chen X, et al. (2009) Structural insights into mechanisms of the small RNA methyltransferase HEN1. Nature 461: 823–827.
    • 60. Lee YS, Nakahara K, Pham JW, Kim K, He Z, et al. (2004) Distinct roles for Drosophila Dicer-1 and Dicer-2 in the siRNA/miRNA silencing pathways. Cell 117: 69–81.
    • 61. Lingel A, Simon B, Izaurralde E, Sattler M (2003) Structure and nucleic-acid binding of the Drosophila Argonaute 2 PAZ domain. Nature 426: 465–469.
    • 62. Rivas FV, Tolia NH, Song JJ, Aragon JP, Liu J, et al. (2005) Purified Argonaute2 and an siRNA form recombinant human RISC. Nat Struct Mol Biol 12: 340–349.
    • 63. Song JJ, Smith SK, Hannon GJ, Joshua-Tor L (2004) Crystal structure of Argonaute and its implications for RISC slicer activity. Science 305: 1434–1437.
    • 64. Yeom KH, Lee Y, Han J, Suh MR, Kim VN (2006) Characterization of DGCR8/Pasha, the essential cofactor for Drosha in primary miRNA processing. Nucleic Acids Res 34: 4622–4629.
    • 65. Moran Y, Praher D, Fredman D, Technau U (2013) The Evolution of MicroRNA Pathway Protein Components in Cnidaria. Mol Biol Evol
    • 66. Edgar RC (2004) MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res 32: 1792–1797.
    • 67. Capella-Gutierrez S, Silla-Martinez JM, Gabaldon T (2009) trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25: 1972–1973.
    • 68. Darriba D, Taboada GL, Doallo R, Posada D (2011) ProtTest 3: fast selection of best-fit models of protein evolution. Bioinformatics 27: 1164–1165.
    • 69. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S (2013) MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol Biol Evol 30: 2725–2729.
    • 70. Martin M (2011) Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet journal 17: 10–12.
    • 71. Zerbino DR, Birney E (2008) Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res 18: 821–829.
    • 72. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N (2012) miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res 40: 37–52.
    • 73. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, et al. (2011) Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol 29: 644–652.
    • 74. Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E (2007) The role of site accessibility in microRNA target recognition. Nat Genet 39: 1278–1284.
    • 75. Alexa A, Rahnenfuhrer J, Lengauer T (2006) Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics 22: 1600–1607.
    • 76. Tarver JE, Donoghue PC, Peterson KJ (2012) Do miRNAs have a deep evolutionary history? Bioessays 34: 857–866.
    • 77. Voinnet O (2009) Origin, biogenesis, and activity of plant microRNAs. Cell 136: 669–687.
    • 78. Huntzinger E, Izaurralde E (2011) Gene silencing by microRNAs: contributions of translational repression and mRNA decay. Nat Rev Genet 12: 99–110.
    • 79. Ghildiyal M, Zamore PD (2009) Small silencing RNAs: an expanding universe. Nat Rev Genet 10: 94–108.
    • 80. Song X, Wang D, Ma L, Chen Z, Li P, et al. (2012) Rice RNA-dependent RNA polymerase 6 acts in small RNA biogenesis and spikelet development. Plant J 71: 378–389.
    • 81. Maniar JM, Fire AZ (2011) EGO-1, a C. elegans RdRP, modulates gene expression via production of mRNA-templated short antisense RNAs. Curr Biol 21: 449–459.
    • 82. Zong J, Yao X, Yin J, Zhang D, Ma H (2009) Evolution of the RNA-dependent RNA polymerase (RdRP) genes: duplications and possible losses before and after the divergence of major eukaryotic groups. Gene 447: 29–39.
    • 83. Sempere LF, Cole CN, McPeek MA, Peterson KJ (2006) The phylogenetic distribution of metazoan microRNAs: insights into evolutionary complexity and constraint. J Exp Zool B Mol Dev Evol 306: 575–588.
    • 84. Mourelatos Z, Dostie J, Paushkin S, Sharma A, Charroux B, et al. (2002) miRNPs: a novel class of ribonucleoproteins containing numerous microRNAs. Genes Dev 16: 720–728.
    • 85. Sempere LF, Sokol NS, Dubrovsky EB, Berger EM, Ambros V (2003) Temporal regulation of microRNA expression in Drosophila melanogaster mediated by hormonal signals and broad-Complex gene activity. Dev Biol 259: 9–18.
    • 86. Chen PY, Manninga H, Slanchev K, Chien M, Russo JJ, et al. (2005) The developmental miRNA profiles of zebrafish as determined by small RNA cloning. Genes Dev 19: 1288–1293.
    • 87. Michalak P, Malone JH (2008) Testis-derived microRNA profiles of African clawed frogs (Xenopus) and their sterile hybrids. Genomics 91: 158–164.
    • 88. Wheeler BM, Heimberg AM, Moy VN, Sperling EA, Holstein TW, et al. (2009) The deep evolution of metazoan microRNAs. Evol Dev 11: 50–68.
    • 89. Mukherjee K, Campos H, Kolaczkowski B (2013) Evolution of animal and plant dicers: early parallel duplications and recurrent adaptation of antiviral RNA binding in plants. Mol Biol Evol 30: 627–641.
    • 90. Campo-Paysaa F, Semon M, Cameron RA, Peterson KJ, Schubert M (2011) microRNA complements in deuterostomes: origin and evolution of microRNAs. Evol Dev 13: 15–27.
    • 91. Peterson KJ, Dietrich MR, McPeek MA (2009) MicroRNAs and metazoan macroevolution: insights into canalization, complexity, and the Cambrian explosion. Bioessays 31: 736–747.
    • 92. Prochnik SE, Rokhsar DS, Aboobaker AA (2007) Evidence for a microRNA expansion in the bilaterian ancestor. Dev Genes Evol 217: 73–77.
    • 93. Fromm B, Worren MM, Hahn C, Hovig E, Bachmann L (2013) Substantial loss of conserved and gain of novel microRNA families in flatworms. Mol Biol Evol 30: 2619–2628.
    • 94. Zeng Y, Qu X, Li H, Huang S, Wang S, et al. (2012) MicroRNA-100 regulates osteogenic differentiation of human adipose-derived mesenchymal stem cells by targeting BMPR2. FEBS Lett 586: 2375–2381.
    • 95. Nagaraja AK, Creighton CJ, Yu Z, Zhu H, Gunaratne PH, et al. (2010) A link between mir-100 and FRAP1/mTOR in clear cell ovarian cancer. Mol Endocrinol 24: 447–463.
    • 96. Chen Y, Wang C, Liu Y, Tang L, Zheng M, et al. (2013) miR-122 targets NOD2 to decrease intestinal epithelial cell injury in Crohn's disease. Biochem Biophys Res Commun 438: 133–139.
    • 97. Philippe L, Alsaleh G, Suffert G, Meyer A, Georgel P, et al. (2012) TLR2 expression is regulated by microRNA miR-19 in rheumatoid fibroblast-like synoviocytes. J Immunol 188: 454–461.
    • 98. Benakanakere MR, Li Q, Eskan MA, Singh AV, Zhao J, et al. (2009) Modulation of TLR2 protein expression by miR-105 in human oral keratinocytes. J Biol Chem 284: 23107–23115.
    • 99. Wang K, Liu F, Zhou LY, Ding SL, Long B, et al. (2013) miR-874 regulates myocardial necrosis by targeting caspase-8. Cell Death Dis 4: e709.
    • 100. Banales JM, Saez E, Uriz M, Sarvide S, Urribarri AD, et al. (2012) Up-regulation of microRNA 506 leads to decreased Cl-/HCO3- anion exchanger 2 expression in biliary epithelium of patients with primary biliary cirrhosis. Hepatology 56: 687–697.
    • 101. Sun H, Li QW, Lv XY, Ai JZ, Yang QT, et al. (2010) MicroRNA-17 post-transcriptionally regulates polycystic kidney disease-2 gene and promotes cell proliferation. Mol Biol Rep 37: 2951–2958.
    • 102. Erez J, Reynaud S, Silverman J, Schneider K, Allemand D (2011) Coral Calcification Under Ocean Acidification and Global Change. In: Dubinsky Z, Stambler N, editors. Coral Reefs: An Ecosystem in Transition: Springer Netherlands. pp. 151–176.