Elucidating the Small Regulatory RNA Repertoire of the Sea Anemone Anemonia viridis Based on Whole Genome and Small RNA Sequencing

Elucidating the Small Regulatory RNA Repertoire of the Sea Anemone Anemonia viridis Based on... Cnidarians harbor a variety of small regulatory RNAs that include microRNAs (miRNAs) and PIWI-interacting RNAs (piRNAs), but detailed information is limited. Here, we report the identification and expression of novel miRNAs and putative piRNAs, as well as their genomic loci, in the symbiotic sea anemone Anemonia viridis. We generated a draft assembly of the A. viridis genome with putative size of 313 Mb that appeared to be composed of about 36% repeats, including known transposable elements. We detected approximately equal fractions of DNA transposons and retrotransposons. Deep sequencing of small RNA libraries constructed from A. viridis adults sampled at a natural CO gradient off Vulcano Island, Italy, identified 70 distinct miRNAs. Eight were homologous to previously reported miRNAs in cnidarians, whereas 62 appeared novel. Nine miRNAs were recog- nized as differentially expressed along the natural seawater pH gradient. We found a highly abundant and diverse population of piRNAs, with a substantial fraction showing ping–pong signatures. We identified nearly 22% putative piRNAs potentially targeting transposable elements within the A. viridis genome. The A. viridis genome appeared similar in size to that of other hexacorals with a very high divergence of transposable elements resembling that of the sea anemone genus Exaiptasia. The genome encodes and expresses a high number of small regulatory RNAs, which include novel miRNAs and piRNAs. Differentially expressed small RNAs along the seawater pH gradient indicated regulatory gene responses to environmental stressors. Key words: coastal ecology, CO seep, ocean acidification, miRNA, piRNA, transposable elements. Introduction biological function, and origin (Ghildiyal and Zamore 2009). Two major classes of small regulatory RNAs in eumetazoans Compared with Bilateria, knowledge about cnidarian small are microRNAs (miRNAs) and PIWI-interacting RNAs (piRNAs). RNAs remains scarce. To date, small RNAs have only been These classes are distinct in terms of their sizes, biogenesis, reported in four cnidarians; the non-symbiotic sea anemone The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com 410 Genome Biol. Evol. 10(2):410–426. doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Nematostella vectensis, the stony corals Stylophora pistillata (Han et al. 2015; Mohn et al. 2015). A secondary piRNA path- and Acropora digitifera,and the hydroid Hydra magnipapillata way serves for piRNA amplification by a “ping–pong loop” (Grimson et al. 2008; Wheeler et al. 2009; Chapman et al. mechanism (Brenneckeetal. 2007; Gunawardane et al. 2010; Krishna et al. 2013; Liew et al. 2014; Moran et al. 2014; 2007). The two distinct piRNA populations (primary and sec- Gajigan and Conaco 2017). ondary) show opposite orientation and complementarity in miRNAs represent a well-studied class of small RNAs that theirfirst10ntpositions (Brennecke et al. 2007; usually range in size from 20 to approximately 24 nt Gunawardane et al. 2007). piRNAs lack universal sequence (miRBase, http://mirbase.org; last accessed January 14, conservation, except that the primary and secondary piRNAs 2018). In animals, miRNAs are initially transcribed as RNA show a preference for an uracil residue at the 5 end (1 U) and polymerase II transcripts (pri-miRNAs), which are further an adenine residue at the 10th position (10 A), respectively. In processed by the RNases Drosha and Dicer into stem–loop all cnidarians investigated so far, piRNAs appear to be highly precursor miRNAs (pre-miRNAs) and mature miRNAs, re- abundant compared with miRNAs and short-interfering RNAs spectively. One strand of the mature miRNA duplex is usually (siRNAs) (Grimson et al. 2008; Krishna et al. 2013; Juliano et al. incorporated into the RNA-induced silencing complex (RISC) 2014; Moran et al. 2014; Gajigan and Conaco 2017; Praher (Schwarz et al. 2003; Gregory et al. 2005), and it guides et al. 2017). The biological role of piRNAs is not well under- the whole complex to complementary mRNA for post- stood, but their most important function seems to be guiding transcriptional gene silencing. In plants, the miRNA biogen- PIWI proteins to suppress transposon activity in animal germ esis pathway involves a Dicer-like 1 (DCL-1) protein that is cells (Brennecke et al. 2007; Gunawardane et al. 2007). piRNA responsible for both cropping and slicing miRNA precursors profiling in cnidarians (Nematostella and Hydra) suggested a in the nucleus (Voinnet 2009). The miRNA silencing mecha- similar role in transposon silencing, but proposed broader si- nism is fundamentally different in animals and plants. lencing functionalities as well (Grimson et al. 2008; Juliano et al. Although animal miRNAs usually perform translational re- 2014; Praher et al. 2017). pression through partial base-pairing to target mRNAs, plant The sea anemone Anemonia viridis exposed to natural miRNAs mostly bind with full or nearly full complementarity ocean acidification conditions appears to be physiologically leading to targeted mRNA cleavage (Bartel 2009). Cnidarian acclimatized to low pH and optimizes its energy utilization miRNAs appear to contain plant-like features in their biogen- under elevated pCO through an increased autotrophic input esis and the post-transcriptional gene silencing follows a (Suggett et al. 2012, Horwitz et al. 2015). Our recent tran- plant-like regulatory pathway (Moran et al. 2013, 2014). scriptome sequencing from the same sampling site indicates Among cnidarians, Nematostella was reported to express increased expression of stress-related transcripts, repression of 87 distinct miRNAs, compared with 26, 31, and 126 global synthesis and boost in certain retrotransposon ele- miRNAs in Acropora, Stylophora,and Hydra, respectively ments at low pH in A. viridis (Urbarova et al., unpublished (Grimson et al. 2008; Wheeler et al. 2009; Chapman et al. results). In plants, it is known that small RNAs can regulate 2010; Krishna et al. 2013; Liew et al. 2014; Moran et al. species tolerance to stress via post-transcriptional gene silenc- 2014; Gajigan and Conaco 2017). Interestingly, only one ing (Sunkar et al. 2007). We therefore wanted to elucidate if miRNA (miR-100) was found conserved between the bilat- acclimatization responses of A. viridis that we observe at the erian and the cnidarian species. miRNAs have several impor- transcriptome level could be caused by small RNA-mediated tant regulatory roles in plants and animals (Bartel 2004; post-transcriptional regulation. Here, we report whole ge- Ghildiyal and Zamore 2009; Vashisht and Nodine 2014). nome and small RNA library sequencing of the symbiotic Expression profiling indicated the cnidarian miRNAs to be sea anemone A. viridis, sampled at a natural seawater pH involved in developmental regulation, regeneration and gradient off Vulcano Island, Sicily—Italy. We mainly aimed thermal stress resilience (Krishna et al. 2013; Moran et al. to identify novel small RNA species in A. viridis, but also elu- 2014; Gajigan and Conaco 2017). However, their roles in cidate their possible involvement in the acclimatization other biological processes, including other environmental responses to low pH conditions. We detected 70 distinct stress responses, have not been investigated in detail. miRNA species, and assessed differentially expressed small piRNAs are usually between 23 and 30 nt in size (Krishna RNAs. Most of the putative piRNAs contained features typical et al. 2013; Liew et al. 2014; Moran et al. 2014). The single- of primary piRNAs and a large fraction showed ping–pong stranded piRNA precursors are either derived from transposable signatures. Our study indicates possible regulatory gene elements or from specific piRNA genomic clusters, and they do responses of small RNAs to low pH. not require Dicer nuclease activity for their processing (Vagin et al. 2006; Houwing et al. 2007; Das et al. 2008). piRNAs Materials and Methods represent a highly diverse class of small regulatory RNAs, reach- Sampling ing several thousand distinct members within a single organism (Aravin et al. 2006; Kawamura et al. 2008). The uniqueness of The temperate symbiotic sea anemone A. viridis (the piRNAs arises from phased production of primary piRNAs Snakelocks Anemone) was collected at Levante Bay, North Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 411 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE Vulcano Island, Sicily—Italy. Acidification conditions are cre- database of proteins encoded by transposable elements ated here by the release of CO into the seawater from a (TEs; NCBI keywords: retrotransposon, transposase, reverse natural vent site at 1 m depth (Boatta et al. 2013; transcriptase, gypsy, and copia). These two databases were Johnson et al. 2013). Samplingwas performedonMay 13 separately queried against our reference genome assembly and 14, 2013, at the depth of 1–2 m at>350 m from the vent and the best annotation was chosen based on alignment cov- site along a gradient of decreasing pH ( pH 7.6 and 7.9), and erage and score. A combined tabular output from the at a control location at 800 m from the vent site with pH searches was further run through two Perl scripts, corresponding to ambient seawater levels (pH 8.2). For sim- “blast92gff3.pl” with additional options -lowscore 0.0001 - plicity, we are referring to average pH values throughout this alignmax 9999 -exonType exon (http://arthropods.eugenes. work as reported in Johnson et al. (2013).A total of nine org/EvidentialGene/evigene/scripts/blast92gff3.pl; last accessed individuals of A. viridis were sampled in 2 days (three from January 14, 2018) and the “overbestgene2.pl” (http://iubio. each location). Small pieces of tissue (0.5 cm) from body bio.indiana.edu/gmod/tandy/perls/overbestgene2.perl; last wall, tentacles and oral disc of each individual were collected accessed January 14, 2018) to create a gff file from blast results and stored separately at 4 C in RNAlater (ThermoFisher and to remove overlapping blast hits, respectively. The results Scientific, Waltham, MA, USA) during transport from the were imported into IBM SPSS Statistics software (version 23), sampling site to laboratory. Then, RNAlater solution was re- where counting of transposable elements was performed. moved and all samples were frozen at 80 C before further Sequence regions corresponding to transposable elements in processing steps. our reference genome assembly were then extracted from our scaffolds using BLAST fastacmd tool (Altschul et al. 1990, 1997) and used as a reference for piRNA analyses. Reference Genome Assembly DNA from one individual of A. viridis (pH 8.2) was extracted RNA Extraction using Wizard(R) Genomic DNA Purification Kit (Promega, Each tissue sample (without excess RNAlater solution) was Madison, Wisconsin, USA). Two whole genome paired-end immediately transferred from 80 Cto 1 ml cold TRIzol re- libraries (2  150 bp) were constructed and sequenced on agent (ThermoFisher Scientific, Waltham, MA, USA). The tis- Illumina HiSeq2500 at Eurofins MWG Operon (Germany). sue was then crushed using Precellys tissue homogenizer at The paired-end reads were processed using Trimmomatic 6,000 rpm for 30 s (Stretton Scientific, Stretton, UK) to mini- (Bolger et al. 2014). Adapters were removed and reads mize degradation of RNA. RNA was twice extracted by chlo- were trimmed and quality filtered using sliding window roform, and subsequently precipitated in isopropanol at 4 C with Phred score> 20. A bias at the first nine nucleotides overnight, washed with 70% ethanol, and rehydrated in was removed by trimming these bases, and reads with Nuclease-Free Water (ThermoFisher Scientific, Waltham, length<40 bp were discarded. SGA preqc tool was run pool- MA, USA). The RNA quality was examined using the Agilent ing the forward and reverse reads from the two libraries 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA, together (Simpson 2014). Platanus, a de novo genome assem- USA) and quantity of the samples was measured using bler for highly heterozygous diploidic organisms, was then Qubit 2.0 fluorometer (ThermoFisher Scientific, Waltham, used to assemble reads with k-mer length of 51 (Kajitani MA, USA). Only high quality samples with RNA integrity num- et al. 2014). To be able to assess our A. viridis genome as- ber (RIN) equal to 7 or higher were used in library sembly for repeat-enriched regions using RepeatMasker (Smit constructions. et al. 2013–2015), we first filtered out short reads from our reference assembly using N75 statistics, resulting in 210,233 Small RNA Sequencing sequences with sizes larger than 173 bp. The filtered assembly was then assessed for repeat-enriched regions using Nine individuals of A. viridis representing three different pH RepeatMasker (Smit et al. 2013–2015), with a custom library conditions (8.2, 7.9, and 7.6) were included in small RNA created by RepeatModeler, which integrates RECON, sequencing. Total RNA from three different tissue samples RepeatScout, and Tandem Repeats Finder (TRF) de novo re- of an individual was pooled at equal amounts. The small peat finding tools to build a repeat library for an assembly RNA fraction was enriched using PureLink miRNA Isolation (Smit and Hubley 2008–2015). In addition to RepeatMasker Kit (ThermoFisher Scientific, Waltham, MA, USA). Libraries annotation, the repeat-enriched regions were extracted from were prepared only from high quality RNA samples the assembly and transposable element annotation was per- (RIN 7) following the SOLiD Total RNA-Seq Kit protocol formed as described previously (Chapman et al. 2010; (ThermoFisher Scientific, Waltham, MA, USA). Different Baumgarten et al. 2015). The annotation pipeline included A. viridis small RNA libraries were barcoded and sequencing TM then also a TBLASTX run using RepBase database (Bao et al. was performed on three lanes of a SOLiD 6-Lane FlowChip 2015), version 22.09 (e-value< 10 ), and a BLASTX search using SOLiD 5500xl sequencer at the Nord University (Bodø, (e-value< 10 ) against a custom-made non-redundant Norway). 412 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Discovery of Novel miRNA (ORFs). Our scaffolds were searched for ORFs using OrfPredictor (http://bioinformatics.ysu.edu/tools/OrfPredictor. After removing low quality (quality score< 18) and less com- html; last accessed January 14, 2018) (Min et al. 2005). plex sequences from our raw small RNA sequencing data set, adapter sequences were trimmed away using trimSOLiDAdaptor.pl Perl script keeping only sequences equal miRNA Analyses to or longer than 18 nt. The filtered reads were mapped in Expression of selected miRNAs was confirmed by quantitative color space using Bowtie (Langmead et al. 2009)to our PCR (qPCR). Six Locked Nucleid Acid (LNA) probes targeting A. viridis genome reference with parameters –integer-quals the predicted miRNAs were designed using online miRNA -l 18 -M 20 –best –strata -e 150 –nomaqround –maxbts 800 – qPCR designer tool from Exiqon (Vedbaek, Denmark). cDNA tryhard -a –col-cqual –col-keepends –mapq 20 –threads 14 – was synthesized from three individuals (10 ng of total RNA chunkmbs 200. These options select the best alignments TM each) per condition using miRCURY LNA Universal RT based on the seed mismatches only and mismatches outside microRNA PCR, Polyadenylation, and cDNA synthesis kit II the seed region are ignored. Therefore, we needed to per- (Exiqon, Vedbaek, Denmark) following the instruction man- form additional filtering using processBowtieAlignments.pl ual. Small RNA for the qPCR analysis was isolated from the Perl script to select alignments with the minimum mismatches same samples that were used for preparation of small RNA along the whole reads. Both Perl scripts are available at https:// libraries. qPCR was performed in duplicates using miRCURY github.com/patelhardip/bitx.git (last accessed January 14, LNA microRNA PCR, ExiLENT SYBR Green master mix (Exiqon, 2018). Mapped reads from each condition were then prepro- Vedbaek, Denmark) in 10 ml. miRNAs were assessed for dif- cessed by bwa_sam_converter.pl Perl script (Friedl€anderetal. ferential expression among the sampling sites with differing 2012), outputting two files essential for running miRDeep2 pH by edgeR (FDR< 0.05) (Robinson et al. 2010). Mature software tool for the novel miRNA predictions (Friedlander miRNAs were aligned to their precursor sequences. miRNAs et al. 2012). The two files were used as input into with <20 counts in less than three conditions were not con- miRDeep2.pl Perl script (Friedlander et al. 2012)that was sidered. We searched for putative animal-like miRNA targets run for identification of novel miRNAs in A. viridis. In addition, by Probability of Interaction by Target Accessibility (PITA) soft- known miRNAs of three related species, N. vectensis, S. pis- ware, based on target complementarity and site accessibility tillata,and H. magnipapillata (Krishna et al. 2013; Liew et al. (Kertesz et al. 2007). Coding regions were predicted from the 2014; Moran et al. 2014) thatwereavailable at thetimeof A. viridis transcriptome (Urbarova et al., unpublished results) the analysis, have been used in the prediction pipeline. These by TransDecoder, version 2.0.1 (Haas et al. 2013) and used as miRNA sequences were downloaded from miRBase, release input into PITA software. The results were filtered based on 21 (Kozomara and Griffiths-Jones 2014). Output from the change in free energy (ddG) of miRNAs binding to its miRDeep2 software was then inspected manually, keeping targets (ddG<10 kcal/mol), seed length of 8 nt with no only predicted miRNAs with miRDeep2 score larger than ten mismatches and no wobble pairs. Only targets fulfilling these and with significant randfold value (p value< 0.05). Small criteria were considered further. We then checked the pre- RNA sequencing data from each individual were assessed sep- dicted targets from PITA for more extensive complementarity arately for the presence of novel miRNAs. Only miRNAs iden- to the miRNAs usingFASTA v36(Pearson andLipman1988) tifiedinat least twoindividuals were considered further. as previously described (Moran et al. 2014), and scored the Possible tRNA contamination was examined by running alignments accordingly. Blast hits were obtained using NCBI tRNAscan on the reference genome (Lowe and Eddy 1997). nr database (e-value< 10 ) and GO terms were assigned to Further, presence of rRNA sequences in predicted hairpin the potential mRNA targets using B2G4Pipe Blast2GO pipe- structures was tested by querying against a custom database line (Go ¨ tz et al. 2008). combining known rRNA sequences from N. vectensis and A. viridis. No contamination was found in either case. In addition, Search for Putative piRNAs to ensure that our miRNA candidates come from the host, assembled scaffolds of A. viridis were screened for their pos- Raw reads from SOLiD sequencing were quality filtered and sible contamination by symbiont DNA using genomes of adapter sequences were trimmed, as described previously for Symbiodinium minutum, Symbiodinium microadriaticum, miRNA discovery. However, only sequences equal to or longer and Symbiodinium kawaguti (Shoguchi et al. 2013; Lin than 23 nt were kept for the piRNA analyses. Sequences were et al. 2015; Aranda et al. 2016). 1,642 scaffolds (mainly short further filtered for reads mapping to miRNA precursors iden- ones with length 100 bp) that were highly similar (e-val- tified in the present study and reads mapping to rRNAs using ue< 10 )toSymbiodinium genomes in the A.viridis genome rRNA databases according to Praher et al. (2017). Sequences assembly were filtered out prior to the small RNA alignment. in color space were aligned using Bowtie with the same Genomic setting of aligned putative miRNAs was inspected for parameters as for the miRNA alignment, but allowing for overlapping regions corresponding to open reading frames maximum three mismatches with the “seed length” of Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 413 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE 23 nt (-l 23 -n 3). We refrained from mapping our reads to reference was found to contain repetitive regions. About unique locations in the A. viridis reference genome due to 27.5% of the repetitive sequences could be assigned to pre- presence of many sequence stretches in assembled scaffolds viously known repetitive elements, but most of these sequen- that most probably correspond to the same genomic loca- ces (25% of the genome) could not be classified into any tions. The aligned sequences were filtered for the best assigned category (supplementary table S1, Supplementary matches using the same custom Perl script as for filtering Material online). Repeat annotation identified only about the miRNA alignments and filtered for reads with 1 U or 8.3% of genome to be comprised of transposable elements 10 A sequence signatures. TE-targeting potential of putative (TEs). These are similar observations as made previously for piRNAs was then assessed by including only putative piRNAs symbiotic sea anemone Exaiptasia sp. (formerly known as mapping antisense to the transposable elements. Overlap Aiptasia sp.; Grajales and Rodrıguez 2014; Baumgarten probabilities of the putative piRNA sequences with opposite et al. 2015). From the identified TE fraction, about half orientation were analyzed using signature.py script (44.4%) were retrotransposons, and amongst them, non-long (Antoniewski 2014). The script computes the probability of terminal repeat (non-LTR) retrotransposons were predominating an antisense read overlapping a sense read with defined (supplementary table S1, Supplementary Material online). length and assigns each overlap length a z-score. The over- lapping signatures of the putative piRNAs were inspected in miRNA Discovery more detail by running PingPongPro v1.0 software (http:// sourceforge.net/projects/pingpongpro/; last accessed January Small RNA libraries from nine individual polyps of A. viridis, 14, 2018). This software also served for inspection of trans- sampled at three different seawater pH conditions (at normal poson silencing by putative piRNAs. Silencing of transposable seawater pH 8.2, as well as at pH 7.9 and 7.6) were prepared elements by the putative piRNAs was only considered if FDR (q and subjected to sequencing on the SOLiD 5500xl platform value)< 0.01 and if it was supported by at least ten putative (fig. 1). The sequencing generated approximately 116 million piRNA reads with ping–pong signatures normalized to the reads (18 nt) after adapter trimming and quality filtering transposon length. To inspect data for presence of piRNA (table 2 and fig. 2). Despite the fragmented nature of our clusters, the putative piRNA reads were mapped onto the genome draft assembly, a high proportion of the small RNA masked genome reference produced by RepeatMasker reads mapped to the genomic reference (between 88% and (Smit et al. 2013–2015) as described previously in color space 92%, table 2). This indicates that even a very preliminary ge- using Bowtie, reporting at maximum five valid alignments. nome assembly is sufficient for the discovery of small RNAs. Finally, our sorted alignment files were submitted to piClust A. viridis miRNAs were identified by the miRDeep2 software software (Jung et al. 2014). Here, the Eps parameter was set tool (Friedlander et al. 2012), and a representative analysis to 1000 and MinReads to 50. result is shown in figure 3A. We predicted in total 70 high- confidence miRNA candidates (20 to 25 nt) for A. viridis, including 61, 60, and 65 distinct miRNA species at pH 8.2, 7.9, and 7.6, respectively (table 3). Most miRNAs were detected Results in all pH conditions studied (n¼ 51), 14 miRNAs in two differ- Genome Reference Assembly and Search for Repeat- ent pH conditions, and five miRNAs were detected in one pH Enriched Regions condition only (table 3). Eight candidate miRNAs in A. viridis Total DNA from a single A. viridis polyp (normal seawater were apparently homologous to those reported in other cni- conditions, pH 8.2) was extracted and subjected to whole darian species (avi-miR-temp-100, 2022, 2023, 2025, 2028, genome sequencing on the Illumina HiSeq2500 platform 2030, 2036, and 2037) (fig. 3B). The predicted miRNA with the (fig. 1). Sequencing generated 43 billion nucleotides (nt) of highest miRDeep2 score (avi-miR-temp-100), and which was genomic data, which corresponded to 144 million paired-end highly expressed in all pH conditions, was identical in sequence reads (table 1). The basic genome characteristics were deter- to miR-100 in N. vectensis and S. pistillata (fig. 3B)(Grimson mined using the SGA preqc software tool (Simpson 2014), et al. 2008; Liew et al. 2014; Moran et al. 2014). Similarly, avi- showing an estimated genome size of 313 Mb (140x cov- miR-temp-2022, 2023, and 2025 were identical to the corre- erage). Adapters were trimmed and the reads were quality sponding miRNAs in N. vectensis (Moran et al. 2014). Other A. filtered before de novo genome assembly, which created viridis miRNAs (avi-miR-temp-2028, 2030, 2036, and 2037) about 1.1 million short scaffolds with N50¼ 2,087. This ge- have one or two nucleotides substitution compared with that nome assembly is highly fragmented, but it was sufficient for of the N. vectensis homolog (fig. 3B). Multiple precursor the mapping of small RNAs and the identification of transpos- sequences for some of the predicted miRNAs were found by able elements (fig. 1). miRDeep2 (table 3). All mature, star, and precursor sequences The A. viridis genome assembly was inspected for repeat are presented in supplementary table S2, Supplementary and low complexity regions using the RepeatMasker software Material online, with read counts for each pH condition in tool (Smit et al. 2013–2015), and about 36% of the genome supplementary table S3, Supplementary Material online. 414 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE 9 individuals RNA isolation (3 individuals/condition) 1 individual (control condition) DNA isolation body wall tentacles mouth piece Illumina HiSeq 2500 sequencing Pooling equal RNA amounts Small RNA enrichment Trimming and quality filtering Preparation of barcoded libraries Genome assembly with Platanus SOLiD 5500xl sequencing putative piRNA pipeline Identification of repetitive DNA by RepeatMasker Trimming and quality filtering Size filtering (≥ 23 nt) miRNA pipeline Annotation of repeat-enriched Size filtering (≥ 18 nt) Mapping to genome and fraction, including transposable elements transposable elements in color space Mapping to genome in color space Assessment of Novel miRNA predictions ping-pong signatures by miRDeep2 Differential expression with edgeR qPCR verification FIG.1.—Data analysis overview. DNA and RNA were isolated from A. viridis adult polyps sampled from a natural seawater pH gradient (at normal seawater pH 8.2, and at low seawater pH 7.9 and 7.6) off Vulcano Island, Sicily—Italy. Only one polyp (pH 8.2) was used for DNA extraction and was subjected to paired- end sequencing on the Illumina HiSeq2500 platform. Sequencing reads were assembled into a draft genome reference. Subsequently, repeat-enriched regions, including transposable elements were identified and annotated in this assembly. Nine polyps were used for small RNA library preparation and sequencing on the SOLiD 5500xl platform. Sequencing reads were further used for novel miRNA discovery and description of putative piRNA reads. Table 1 of putative miRNA clusters. We identified four clusters, three The Amount of Reads Gained from Genome Sequencing of Anemonia contained two miRNA sequences (avi-miR-temp-11 and 66; viridis avi-miR-temp-28 and 27; and avi-miR-temp-64 and 67) and Sequencing No. of No. of %Paired-End one cluster contained three miRNAs (avi-miR-temp-2, 13, and Index Paired-End Trimmed and Reads Kept After 39) (supplementary fig. S1, Supplementary Material online). Raw Reads Quality Filtered Filtering No open reading frames (ORFs) spanning cluster regions could Reads be predicted, implying that the miRNAs were transcribed as CTTGTA 82,428,617 71,676,503 87.0 independent transcription units. Further clustering of miRNA GCCAAT 61,246,666 52,649,432 86.0 loci could not be assessed due to the fragmented nature of Sequencing of barcoded genome libraries was performed in one lane of the genome assembly. Illumina HiSeq2500 sequencing machine in 2 150 bp mode. The amount of sequen- ces presented here is the number of raw paired-end sequences obtained after the run. We then asked if any of the expressed miRNAs were colocal- ized with predicted transposable elements in the A. viridis ge- nome reference. One miRNA (avi-miR-temp-58) was found Genome Context of Identified miRNAs encoded within a DNA transposon (fig. 4A). avi-miR-temp-58 In sea anemones, little is known about the genomic miRNA was detected only at pH 7.9 in the small RNA sequencing exper- clusters that generate pri-miRNAs. Therefore, we searched iment (but in all three individuals inspected). However, we the A. viridis reference genome sequence for the presence detected avi-miR-temp-58 in all conditions studied by a Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 415 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE Table 2 The Amount of Reads Gained from Small RNA Sequencing of Anemonia viridis Individuals Raw Filtered Small % Reads Aligned Filtered Small Reads Aligned Putative piRNA %Putative piRNAs Reads RNA Reads to the Genome RNA Reads to Genome Reads Aligned Aligned to b b c (18 nt) (18 nt) (23 nt) (23 nt) to Genome Genome pH 7.6–1 16,111,517 12,117,933 88.5 9,831,771 7,681,304 6,643,549 86.5 pH 7.6–2 12,492,090 11,175,766 89.5 9,465,543 7,721,730 6,795,032 88.0 pH 7.6–3 13,794,761 10,609,870 89.1 8,921,147 7,238,217 6,310,624 87.2 pH 7.9–1 18,815,344 15,837,619 91.5 13,235,131 11,175,365 10,001,436 89.5 pH 7.9–2 20,094,253 17,369,794 90.0 16,247,704 13,837,596 11,960,594 86.4 pH 7.9–3 16,794,089 15,567,040 89.3 14,958,985 12,687,556 11,353,938 89.5 pH 8.2–1 13,066,319 12,109,753 88.0 11,317,833 9,184,500 8,163,612 88.9 pH 8.2–2 17,379,705 10,309,153 90.3 8,124,195 6,661,625 5,543,870 83.2 pH 8.2–3 13,492,380 11,225,436 90.3 8,769,865 7,189,261 6,271,324 87.2 a TM Small RNA libraries from each individual were barcoded, pooled, and sequencing was performed on three lanes of a SOLiD 6-Lane FlowChip using the SOLiD 5500xl sequencer. The amount of sequences presented here is the sum of raw reads from the three lanes. These reads had adapter removed and were quality filtered. They differ only according to the size filtering. Reads (23 nt) aligned to the reference genome and filtered for piRNA sequence signatures (1 U and 10 A). 6e+06 (Robinson et al. 2010). Here, avi-miR-temp-37, 52, 56, 58, and 59 appeared up-regulated at pH 7.9, whereas avi-miR- 5e+06 temp-13, 29, 48, and 60 appeared down-regulated at pH 7.9 4e+06 (fig. 5). Six miRNAs were then selected for verification analysis by qPCR (avi-miR-temp-37, 58, 60, 100, 2023, and 2028), 3e+06 where three miRNAs homologous to Nematostella with ap- parently unaffected expression levels in the different pH con- 2e+06 ditions served as controls (supplementary fig. S2, 1e+06 Supplementary Material online). The control miRNAs (avi- miR-temp-100, 2023, and 2028) were detected by qPCR in all pH conditions at similar expression levels, and thus are in 18 19 20 20 21 22 23 24 25 26 27 28 29 30 30 31 32 33 34 35 good agreement with results generated from small RNA se- quencing. The miRNA avi-miR-temp-37 was detected only at sequence length in nt pH 7.9 and only in one individual at pH 7.6, and avi-miR- temp-60 was detected only at pH 8.2. However, in contrast FIG.2.—Sequence length representation of small RNAs in A. viridis. Distribution of small RNA reads after adapter trimming and quality filtering with the observation from small RNA sequencing, we in one individual sampled from pH 8.2. Two distinct peaks could be ob- detected the presence of avi-miR-temp-58 in all conditions served; first around 22 nt representing both miRNA and siRNA reads, and studied, though in higher abundance at pH 7.9 (supplemen- second around 28 nt representing putative piRNA reads. tary fig. S2, Supplementary Material online). We then searched for putative mRNA targets of 13 se- quantitative PCR (qPCR) approach (see below), though in higher lected miRNAs that were differentially expressed along the abundance at pH 7.9 (supplementary fig. S2, Supplementary pH gradient (avi-miR-temp-13, 29, 37, 48, 52, 56, 58, 59, Material online). avi-miR-temp-58 and its precursor was pre- and 60) (fig. 5), detected only in one pH condition (avir- dicted to create a 1 nt 3 overhang (fig. 4B). The latter feature miR-temp-48, 58, 59, 60, and 65), or detected only at low suggests a group II pre-miRNA that requires a 3 -end mono- pH, that is, at pH 7.6 and/or pH 7.9 (avi-miR-temp-37, 48, 50, uridylation for further Dicer processing (Heo et al. 2012). 57, 58, 59, 60, 64, and 65) (table 3). Differentially expressed mRNAs previously identified in A. viridis in the same individu- als from the same sampling experiment (Urbarova et al., Differential miRNA Expression upon Seawater pH Gradient unpublished results) were assessed as potential targets. All high-confidence miRNAs were included in differential ex- After stringent filtering criteria, including full seed matching pression analyses. Despite that 19 candidate miRNAs could with extended pairing, we identified 9 out of the 13 selected not be detected in all conditions studied, only nine miRNAs miRNAs that could potentially target 13 of the differentially were recognized as differentially expressed between expressed mRNAs along the low pH gradient (supplementary conditions by edgeR (FDR< 0.05) (fig. 5; supplementary table S6, Supplementary Material online). Although we could tables S4 and S5, Supplementary Material online) not consistently detect miRNA upregulation and its mRNA 416 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 number of reads Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE A B Provisional ID : scaffold25818_len18526_cov78_2565 Score total : 428917.1 hsa-miR-100 aacccguagauccgaacuugug- Score for star read(s) : 3.9 nve-miR-100 -acccguagauccgaacuugugg Score for read counts : 428907.3 spi-miR-100 -acccguagauccgaacuugugg Score for mfe : 1.3 Score for randfold : 1.6 adi-miR-100 -ucccguagauccgaacuugugg Score for cons. seed : 3 5' avi-miR-temp-100 -acccguagauccgaacuugugg Total read count : 841295 u u c u Mature read count : 799715 ccguagauc gaacuug gg g c ******************** Loop read count : 0 3' g g c a u c u a g u u u g a a c c c c Star read count : 41580 c u u a u c nve-miR-2022 uuugcuaguugcuuuugucccgc- avi-miR-temp-2022 uuugcuaguugcuuuugucccgc- freq. spi-miR-2022 uuugcuaguugcuuuugucccguu adi-miR-2022 uuugcuaguugcuuuugucccgu hma-miR-2022 uuugcuaguugcuuuuguccccu- 0.75 ********************* 0.5 nve-miR-2023 aaagaaguacaagugguaggg spi-miR-2023 aaagaaguacaagugguaggg 0.25 adi-miR-2023 aaagaaguacaagugguaggg avi-miR-temp-2023 aaagaaguacaagugguaggg 1 23 35 57 length ********************* nve-miR-2025 uuuuuuagcccgcggaaguugu Mature Star adi-miR-2025 auuuuuagcccgcggaaguugc 5'- guaggcaguguuguuugucgacccguagauccgaacuugugguugucuaccuccccaaguuuugaucuacggaacuaaaaucaacgcuaauaugaccaggaacauccaagga -3' obs avi-miR-temp-2025 uuuuuuagcccgcggaaguugu guaggcaguguuguuugucgacccguagauccgaacuugugguugucuaccuccccaaguuuugaucuacggaacuaaaaucaacgcuaauaugaccaggaacauccaagga exp ******************** ......(((((((.((......(((((((((.(((((((.((..(....)..)).))))))).)))))))))......)).))))))).......((.(((...)))..)). reads mm sample nve-miR-2028 uaauguuccugcuuguuccua ....................acccguagauccgaacuugugg...................................................................... 737600 0 seq avi-miR-temp-2028 aaauguuccugcuuguuccug ....................acccguagauccgaacuugu........................................................................ 26766 0 seq ******************* ....................acccguagauccgaacuugug....................................................................... 21383 0 seq hma-miR-2030 uagcauaacauuguaagaaaca avi-miR-temp-2030 uagcauaacauaguaagagauu ....................acccguagauccgaacuug......................................................................... 4646 0 seq nve-miR-2030 uagcauaacauuguaagagauu spi-miR-2030 uagcauaacauuguaagagauc ....................acccguagauccgaacuugugU...................................................................... 1627 1 seq adi-miR-2030 uagcauaacauuguaagagaucu *********** ****** * ......................................................ccaaguuuugaucuacggaacu.................................... 30488 0 seq spi-miR-2036 uauauuguacgacucucaucgugu ......................................................ccaaguuuugaucuacggaac..................................... 8550 0 seq adi-miR-2036 uauauuguacgacucucaucgug nve-miR-2036 uauauuguacgacucucaucgua- ......................................................ccaaguuuugaucuacggaa...................................... 1460 0 seq avi-miR-temp-2036 uauauuguacgacucucaucguag ********************** nve-miR-2037 ugugauuggagacuuuuaccgu avi-miR-temp-2037 ugugauuggagacuuuuaucgu ****************** *** FIG.3.—Identified miRNAs in A.viridis with similarity to known miRNAs. (A) A typical prediction result from miRDeep2 software tool showing the miRNA precursor, mature and star sequence and their abundances in the sample. Shown is avi-miR-temp-100 precursor with top sequence alignments in the sample for each strand. (B) Alignments of novel miRNAs from A. viridis to known miRNAs from other species. Sequences of our predicted miRNAs from A. viridis (denoted as avi-miR-temp) were aligned to known miRNA sequences from H. sapiens (hsa-miR), N. vectensis (nve-miR), H. magnipapillata (hma-miR), S. pistillata (spi-miR), and A. digitifera (adi-miR). (temp¼ temporary; miRNAs that are not yet registered in the miRBase). target downregulation, we made one interesting observation. plot (fig. 2), and most likely represent PIWI-interacting RNAs We detected avi-miR-temp-50, present only at pH 7.6 and pH (piRNAs). Based on an earlier report on piRNA signatures in 7.9, to target an RNase HI domain of a DIRS1 retrotransposon. cnidarians (Moran et al. 2014), we explored the trimmed and The corresponding transcript was found downregulated both quality filtered small RNA reads with minimum length of 23 nt at pH 7.6 and 7.9 compared with pH 8.2, which could mean and with the typical base preference signatures (hereafter that this domain is inactivated and reverse transcription is called putative piRNAs) in our data set. The putative piRNAs therefore inhibited. were aligned to two different reference data sets; the A. viridis reference genome, and the transposable elements identified in the genome. In total, about 83–90% of reads23 nt Search for Putative piRNAs and Their Characteristics aligned to the reference genome were putative piRNAs Most small RNA sequences present in our data set showed a (table 2), with about 14–42% mapping to transposable ele- distinct peak at 27–29 nt in the small RNA size distribution ments (supplementary table S7, Supplementary Material Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 417 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE Table 3 The Listof70 Predicted miRNAs in Anemonia viridis from Various pH Conditions Temporary miRNA Name Mature Sequence Length Stem–Loop Length Present in Condition Similarity to Known miRNAs avi-miR-temp-100 acccguagauccgaacuugugg 22 56 All nve-miR-100-5p avi-miR-temp-2022 uuugcuaguugcuuuugucccgc 23 52; 53 All nve-miR-2022-3p avi-miR-temp-2023 aaagaaguacaagugguaggg 21 53 All nve-miR-2023-3p avi-miR-temp-2025 uuuuuuagcccgcggaaguugu 22 53 All nve-miR-2025-3p avi-miR-temp-2028 aaauguuccugcuuguuccug 21 48 All nve-miR-2028-5p avi-miR-temp-2030 uagcauaacauaguaagagauu 22 52 All nve-miR-2030-5p avi-miR-temp-2036 uauauuguacgacucucaucguag 24 54 All nve-miR-2036-3p avi-miR-temp-2037 ugugauuggagacuuuuaucgu 22 54 All nve-miR-2037-3p avi-miR-temp-1 gaucaagucaaauacaucucu 21 48 All avi-miR-temp-2 uaucaaggcagucuuaccauau 22 54 All avi-miR-temp-3 uacaaauguuacgcagcagaac 22 55 All avi-miR-temp-4 ugacauugcugcccgaaucucc 22 88 All avi-miR-temp-5 uuuaauguuacugcucguucc 21 53 All avi-miR-temp-6 aauuucaaauauccacugauuga 23 55 All avi-miR-temp-7 uugagcaucuguugcaugucua 22 53 All avi-miR-temp-8 aucaucgccacuagcaucguca 22 55 All avi-miR-temp-9 aagggcaagacaauagaauuuca 23 59 All avi-miR-temp-10 cuugauaguacuuuugccuugc 22 52 pH 7.9, pH 8.2 avi-miR-temp-11 uaguagguucuuauaagcuauu 22 55 All avi-miR-temp-12 uauaagucuaggcugguuaaga 22 56 All avi-miR-temp-13 auacugaacuugaaagaagugau 23 55 All avi-miR-temp-14 aaacgcuguucuugguaguca 21 55 All avi-miR-temp-15 uaacaaagcaguuuggcuguau 22 55 All avi-miR-temp-16 ucuggcugauuugaagaaaga 21 51 All avi-miR-temp-17 acaucaaacaaagcaguuug 20 53 All avi-miR-temp-18 auuacccguaaauaaauucaau 22 54 All avi-miR-temp-19 aaccccaacgcgggccucugg 21 51 All avi-miR-temp-20 uuaguuugcacucauuugcugg 22 55 All avi-miR-temp-21 auuacccagaauggggccuuu 21 55 All avi-miR-temp-22 uauucuccaaaaauucacaagg 22 52 All avi-miR-temp-23 uaaacuaguugauaggauugu 21 51 All avi-miR-temp-24 acagauugcggcaaccgugcag 22 72; 88 All avi-miR-temp-25 ucaaauguugcgcagcagaac 21 55 All avi-miR-temp-26 ugcugcaguuuagacugaccuc 22 53 All avi-miR-temp-27 uccucaaguuuugauuguaauac 23 51; 52 All avi-miR-temp-28 uucuuaaguuuugauuguaauac 23 52 All avi-miR-temp-29 aucuacugauacuaaguauccg 22 54 All avi-miR-temp-30 uuucuguaguacuuuauccuggc 23 54 All avi-miR-temp-31 uauucaaucagucuggcuguua 22 52 All avi-miR-temp-32 ucuuuugauaaauaccaccaaca 23 56 All avi-miR-temp-33 uacucugaaguguacuuagugu 22 53; 54 All avi-miR-temp-34 gauaugauauaauauguaugug 22 57 All avi-miR-temp-35 uauacauauuuaguaucgauaucag 25 57; 58 All avi-miR-temp-36 uaaauacacaauaucuauagcagu 24 55 All avi-miR-temp-37 uaugguagugauguuuagaaa 21 49 pH 7.6, pH 7.9 avi-miR-temp-38 ccggacaaugagaauagcuga 21 56 All avi-miR-temp-39 ugaucaauaaaagaaacaucguu 23 54 All avi-miR-temp-42 uaucacauuuaaaacacucaug 22 53 All avi-miR-temp-43 ucauacgauauuuuucacuagu 22 55; 56 All avi-miR-temp-44 aaccucaugucagagaucaaa 21 53 pH 7.6, pH 8.2 avi-miR-temp-45 acagagccuccuuuaaccuccu 22 60 pH 7.6, pH 8.2 avi-miR-temp-47 ugguagaacaaguaacuugcugc 23 55 pH 7.6, pH 8.2 avi-miR-temp-48 gaaaaagacauuuagagacuug 22 56 pH 7.6 (continued) 418 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Table 3 Continued Temporary miRNA Name Mature Sequence Length Stem–Loop Length Present in Condition Similarity to Known miRNAs avi-miR-temp-49 aaugucaccaaguuucgacca 21 50 pH 7.6, pH 8.2 avi-miR-temp-50 aggcccuggggaaacaaugga 21 54 pH 7.6, pH 7.9 avi-miR-temp-52 uggaugcucaauuugccaauugc 23 75 pH 7.9, pH 8.2 avi-miR-temp-53 aacuuaaaacaaaaaucucccu 22 53 pH 7.6, pH 8.2 avi-miR-temp-54-1 aucuauucacugugggcguccagu 24 54 All avi-miR-temp-54-2 aucuauucauugugggcguccagu 24 55 All avi-miR-temp-55 uacuacuuugacaaugugaugg 22 53; 77 pH 7.6, pH 8.2 avi-miR-temp-56 aggucagucuaaacugcagca 21 54; 55 pH 7.6, pH 8.2 avi-miR-temp-57 gcuuugaaaauguaaagaaca 21 50 pH 7.6, pH 7.9 avi-miR-temp-58 ugcaguauucaguaugcacua 21 65 pH 7.9 avi-miR-temp-59 ucggcgccggucacgcgauaga 22 52 pH 7.9 avi-miR-temp-60 caagcuauaaauuccaacuga 21 50 pH 7.6 avi-miR-temp-61 ucgaguaaaauauuacagaaaug 23 54 pH 7.6, pH 8.2 avi-miR-temp-64 ucaucucuuguggcuugacauu 22 51 pH 7.6, pH 7.9 avi-miR-temp-65 uggugcaguuuagacugacccuu 23 54 pH 7.9 avi-miR-temp-66 cuagauuaugagagcuuaugu 21 53 All avi-miR-temp-67 ugugugaaaacaugacaagaucu 23 50 All Presence of two numbers in this column indicates that two different miRNA precursors (pre-miRNAs) of the same mature miRNAs have been detected in our genome assembly. All nucleotide sequences of mature and star miRNAs and pre-miRNA precursors are listed in supplementary table S2, Supplementary Material online. DNA transposon miRNA precursor star mature GCACCATTCAGTATGCAGTATTGAGTATGCAGTATTCAGTGTGCAGTATGCAGTGTGCAGTGTGCAGTGTGC AGTGTGCAGTGTGCAGTGTGCAGTGTGCAGTATGCAGTAAGCAGTATGCAGTGTGCAGTGTACAGTGTGCAG TATGCAGTATTCAGTGTGCAGTGTGCAGTGTGCAGTATGCAGTATGCAGTATGCAGTATGCAGTATTCAGTG TGCAGTATGCAGTATGCAGTGTGCAGTATACAGTATTCAGTATGCACTATTGAGTATGCAGTATTCAGTATG CACTATTGAGTATGCAGTATTCAGTGTGCAGTATTCAGTATGCACTATTCAGTATGCAGTATTCAGTGTGCA GTATTCAGTGTGCACTATTCAGTATGCACTATTGAGTATGCAGTATTCAGTATGCAGTATTCAGTATGCAGT ATGTAGTATGCAGTATGCAGTATGCACTATTGAGTATGCAGTATTCAGTGTGCACTATTCAGTATGCACTAT TCAGTATGCACTATTGAGTATGCAGTATTGAGTATGCA miRNA precursor hairpin structure u star sequence 5' a g c ugc a u a u u a u g u a ug c a c a g u u g a g u a 3' a c g u a u g a u u a u a c g u g c u g a g c u c u a c u g u a mature sequence FIG.4.—Precursor of avi-miR-temp-58 and its DNA transposon localization. (A) The miRNA precursor of avi-miR-temp-58 was found localized in a DNA transposon. A schematic representation of the scaffold region is depicted above the DNA sequence. Only the part of the scaffold with similarity to the DNA transposon is shown. The DNA transposon has homology to a transposable element from Crassostrea gigas. The miRNA precursor is marked in color and the whole sequence is underlined. Mature miRNA is indicated in red and star sequence in violet. (B) Hairpin structure of avi-miR-temp-58 precursor with 1 nt 3 overhang. Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 419 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 individual Urbarova et al. GBE Color Key −1.5 −1 −0.5 0 0.5 1 1.5 Row Z−Score avi-miR-temp_56 avi-miR-temp_29 avi-miR-temp_13 avi-miR-temp_52 avi-miR-temp_37 avi-miR-temp_58 avi-miR-temp_59 avi-miR-temp_48 avi-miR-temp_60 1 3 2 2 1 3 1 3 2 pH 7.9 pH 8.2 pH 7.6 pH 8.2 pH 7.6 condition FIG.5.—Differentially expressed miRNAs under low pH conditions. Nine miRNAs were found differentially expressed among the sampling sites (edgeR, FDR< 0.05). This included two miRNAs detected in all pH conditions studied (avi-miR-temp-13 and 29) and seven miRNAs that could be detected in only one or two different pH conditions. Five miRNAs were differentially expressed between pH 7.6 and 7.9, three downregulated (avi-miR-temp-52, 58, and 59) and two upregulated (avi-miR-temp-48 and 60) at pH 7.6 compared with pH 7.9. Only one miRNA was detected differentially expressed between pH 7.6 and 8.2 (avi-miR-temp-37), and it was upregulated at pH 7.6 compared with pH 8.2. Eight miRNAs were found differentially expressed between pH 7.9 and 8.2, three downregulated (avi-miR-temp-13, 29, and 48) and five upregulated (avi-miR-temp-37, 52, 56, 58, and 59) at pH 7.9 compared with pH 8.2. Differentially expressed miRNAs were hierarchically clustered into heatmap based on counts per million (cpm) and scaled by row. online), including unclassified fraction of repeats. Most of the Ping–Pong piRNA Amplification Signature in A. viridis putative piRNA reads mapped to the reference genome and We further investigated if ping–pong signatures, that is, transposable elements showed strong preference for 1 U 10 nt overlaps of putative piRNAs with opposite direction, (fig. 6), a feature consistent with the primary piRNA popula- were common in our data set. The probability of overlap by tion. Most of the putative piRNAs are found in genomic clus- 1–30 nt of putative piRNAs with opposite orientation was ters and the majority of piRNA cluster loci appear unistranded assessed. The data exhibited strong ping–pong signatures, (61–68%), where piRNAs are transcribed from one strand of since most reads showed preference for 10 nt 5 overlaps of the piRNA locus (supplementary fig. S3, Supplementary putative piRNAs with opposite orientation in all conditions, in Material online). all individuals, and for all reference data sets (z-score> 5). About one third of the genome scaffolds contained Probability of other overlaps was much lower (z-score< 1) expressed piRNA loci (supplementary table S8, (fig. 7A). We could detect 14–15% putative piRNAs map- Supplementary Material online). We observed more ping to identified transposable elements, and up to 42% (sup- expressed piRNA loci at pH 7.9 than at pH 8.2 or 7.6. plementary table S7, Supplementary Material online) when Putative piRNAs were found to map to about 24–30% of including the unclassified fraction of identified repeats (sup- the identified transposable elements in all conditions and in plementary table S1, Supplementary Material online). Only all individuals (supplementary table S9, Supplementary around 10% of putative piRNAs mapped to transposable ele- Material online). These included both DNA transposons and ments showed ping–pong signatures. A slightly higher pro- retrotransposons (supplementary fig. S4, Supplementary portion of putative piRNA reads with ping–pong signatures Material online). Interestingly, retrotransposons appeared was found to map to the genome reference outside the more frequently targeted by piRNAs than DNA transposons repeat-enriched regions (17%). Majority of these most (supplementary table S9, Supplementary Material online). 420 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Sense putative piRNA reads aligned to genome 2.0 1.0 C A U U C C C A UACC U UA A A UAC UCAUACU UAA AUA U AA U AUGCGUUCC GA UCU CGCU G A UUC GU AU CA G GGU CACGGACGAGG CG C CG G 0.0 5 10 15 20 25 30 WebLogo 3.3 2.0 Antisense putative piRNA reads aligned to genome 1.0 A A U U U GC G A G A U AUU U C C C C C 0.0 5 10 15 20 25 WebLogo 3.3 Sense putative piRNA reads aligned to TEs 2.0 1.0 GACA CCAGUUU AU U AU UU U U C AU U 0.0 5 10 15 20 25 WebLogo 3.3 Antisense putative piRNA reads aligned to TEs 2.0 1.0 AG C A G U A U UU U U C CA C G 0.0 5 10 15 20 25 WebLogo 3.3 FIG.6.—Base preferences of putative piRNA reads mapped to various data sets. Shown are base preferences of putative piRNA reads mapping to sense (A, C)and antisense (B, D) strand of genome (A, B) and transposable elements (TEs; C, D). Base preferences did not significantly differ at various pH conditions. Depicted is always one sequence set of specific length from one condition. The Y-axis represents the entropy score for the base bias. probably represent protein-coding genes. When assessing the sense and antisense orientation (figs. 7B and C). Shorter an- base-preferences of piRNAs with ping–pong signatures, we tisense reads (27 nt) mapped to transposable elements found that a substantial amount of the putative piRNAs had showed minor preference for 10 A (secondary piRNAs) 1 U preference (primary piRNAs) (fig. 7B). Sequences mapping (fig. 7B). Only around 8–9% of the putative piRNAs appeared to transposable elements showed preference for 1 U in both to be targeting transposable elements in all conditions studied Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 421 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 bits bits bits bits Urbarova et al. GBE FIG.7.—Ping–pong pathway signature. (A) Overlap probabilities of sense and antisense reads mapping to the genome and transposable elements (TEs). Overlap probabilities in all the different pH conditions are shown. (B) Depicted are base preferences of putative piRNA reads with ping–pong signatures mapping to the genome and TEs. The Y-axis represents the entropy score for the base bias. (C) Shown are putative piRNA reads aligned to a TE. Three regions with identified ping–pong signatures are highlighted. Green reads correspond to sense strand, and red reads to antisense strand. (supplementary table S7, Supplementary Material online). This genome (supplementary table S7, Supplementary Material fraction further increased nearly up to 22% when including online). These might potentially represent very divergent the unclassified fraction of the identified repeats in the transposable elements, as reported also for Exaiptasia sp. 422 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE (Baumgarten et al. 2015). Higher fraction of putative piRNAs non-symbiotic Nematostella was reported to carry about targeting transposable elements is also expected to be found four times more DNA transposons than retrotransposons during development or when extracting specifically germline (Putnam et al. 2007), which contrasts that of the symbiotic cells from A. viridis, as observed in N. vectensis (Praher et al. Exaiptasia with slightly more retrotransposons than DNA 2017). transposons (Baumgarten et al. 2015). Our data from the Ping–pong activity is mainly linked to transposon silencing symbiotic A. viridis does not appear to resemble any previously (Brennecke et al. 2007; Gunawardane et al. 2007). Therefore, sequenced cnidarian in terms of the transposable element we investigated if ping–pong activity changes could be ob- distribution, even though it contains approximately equal frac- served at different pH conditions. Transposable elements tions of DNA transposons and retrotransposons. It is interest- were only considered silenced if a significant ping–pong ac- ingtonote that the Gypsy element is the most frequent LTR tivity feature could be detected within the transposon region, retrotransposon in all the cnidarian species, including A. viridis. with FDR (q value)< 0.01. We found that a possible ping– We identified70distinct miRNAs in A. viridis, and 61 of pong dependent suppression varied among individuals in these were detected in normal seawater conditions at pH 8.2. each condition, and examples from the BEL and Gypsy LTR Only eight miRNAs were similar to previously known miRNAs retrotransposons are shown in supplementary figures S5 and in Nematostella (Grimson et al. 2008; Moran et al. 2014), six S6, Supplementary Material online. However, only a small in Acropora (Gajigan and Conaco 2017), five in Stylophora fraction of the identified transposable elements appeared si- (Liew et al. 2014)and two in Hydra (Krishna et al. 2013). lenced by the ping–pong pathway in all conditions studied These results support that taxonomically restricted miRNAs (< 1%; supplementary table S10, Supplementary Material are common to cnidarians, including A. viridis—an observa- online). tion seen mainly in plants, and which could be explained by high sequence turnover rates of miRNAs, as suggested by Moran et al. (2017). In agreement with other reports in cni- Discussion darians (Grimson et al. 2008; Wheeler et al. 2009; Liew et al. Here, we report a preliminary draft genome reference se- 2014; Moran et al. 2014, Gajigan and Conaco 2017), we quencing of the symbiotic sea anemone A. viridis,with an detected only one miRNA in A. viridis (avi-miR-temp-100) to estimated genome size of approximately 313 Mb. The par- be conserved with miRNAs in bilaterians. This miRNA belongs tially assembled reference genome was used to assess trans- to the miR-100 family, and it was found identical in sequence posable element and small RNA loci. We also performed small to nve-miR-100 and spi-miR-100 in Nematostella and RNA sequencing along a natural seawater pH gradient and Stylophora, respectively (Grimson et al. 2008; Wheeler et al. identified differentially expressed RNA candidates. In A. viridis, 2009; Liew et al. 2014; Moran et al. 2014). In bilaterians, we found 70 distinct miRNAs and thousands of putative including nematodes and humans, miR-100 makes a cluster piRNAs, suggesting that small RNAs are widespread regula- in the genome together with let-7 and miR-125, and regu- tors in the control of gene expression and transposable ele- lates transcripts involved in multiple cellular and developmen- ment silencing in this species. tal processes, as well as cancer progression (Christodoulou Theestimated genomesizeof A. viridis appears interme- et al. 2010; Sokol 2012; Li et al. 2015). The absence of diate compared with the sea anemones Exaiptasia sp. miR-51/miR-100 family was first reported in nematodes to (260 Mb) and Nematostella vectensis (329/450 Mb), slightly result in lethality during development (Shaw et al. 2010). less than the stony coral Acropora digitifera (420 Mb), and In A. viridis, as well as in other cnidarians, miR-100 appears substantially smaller than the freshwater hydroid Hydra mag- to be transcribed from an individual gene locus, but its bio- nipapillata (1.3 Gb) (Putnam et al. 2007; Chapman et al. logical role in gene repression is not well established. In the 2010; Shinzato et al. 2011; Baumgarten et al. 2015). Thus, coral Stylophora, Liew and coworkers speculated that miR- a general trend is that hexacorals harbor relatively small 100 could be involved in the calcification process (Liew genomes. We found that about 36% of the A. viridis genome et al. 2014). However, since sea anemones lack any sort of contains repeated sequences, which is a higher fraction than calcified skeleton, other processes have to be regulated by Exaiptasia and Nematostella (both 26%) and Acropora (13%), miR-100 in sea anemones. but less than Hydra (57%) (Putnam et al. 2007; Chapman We identified and described four miRNA clusters within the et al. 2010; Shinzato et al. 2011; Baumgarten et al. 2015). A. viridis genome. However, a more detailed analysis in regard There is a significant heterogeneity in the distribution of clas- to clustering of individual miRNAs was not possible due to ses and subclasses of transposable elements among the inves- insufficient contiguity of our draft assembly. Therefore, we tigated cnidarians. Whereas Hydra contains approximately cannot exclude that some miRNAs predicted in our study equal fractions of DNA transposons and retrotransposons, form additional miRNA clusters, where miRNA pairs are lo- Acropora harbors four times as many retrotransposons than cated further apart. Interestingly, one of the identified miRNA DNA transposons. There are also significant differences be- locus localized inside a DNA transposon (fig. 4), and this tween the sea anemones Nematostella and Exaiptasia.The miRNA (avi-miR-temp-58) appeared expressed mostly at Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 423 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE seawater pH 7.9. The formations of small RNAs from trans- detect any significantly differentially expressed piRNAs or posable element loci are not unusual among animals, and piRNA clusters at this point, we noted an increase in putative dozens of publications inspecting various species have piRNA expression at pH 7.9 compared with pH 7.6 and 8.2. reported miRNAs originating from transposable elements One possible biological implication could be less restricted (reviewed in Roberts et al. 2014). However, to our knowledge transposon activities at seawater pH 7.6 compared with avi-miR-temp-58 is the first example of a TE-encoded miRNA pH 7.9. reported in any cnidarian. We found nine miRNAs to be dif- ferentially expressed in A. viridis along the seawater pH gra- Conclusion dient, indicating that miRNA-based gene repression might be The A. viridis genome appears similar in size and in transpos- involved in compensating environmental stressors. Here, we able element divergence to that of Exaiptasia sp., a related identified few potential mRNA targets, including stress- sea anemone with a symbiotic lifestyle resembling that of related and mobile element proteins. Anemonia spp. The A. viridis genome encodes and expresses PIWI-interacting RNAs (piRNAs) have previously been a high number of small regulatory RNAs, and when compared reported in Nematostella (Grimson et al. 2008; Praher et al. with the sea anemone Nematostella, a large fraction (89%) of 2017). This sea anemone contains two piRNA classes, where miRNAs appears taxonomically restricted. A. viridis expresses a class I possesses an unknown function during germline devel- high amount of candidate piRNA sequences with putative opment and class II is involved in gene silencing, including functions in transposable element silencing and in other still transposons, by the ping–pong mechanism. In A. viridis,we unknown cellular functions. Some small RNAs appeared dif- found a high number of expressed piRNA candidates, appar- ferentially expressed along a seawater pH gradient, suggesting ently representing both piRNA classes, even though we did a regulatory role in the response to environmental stressors. not specifically extract and analyze germline cells in our study. However, we were not able to characterize piRNA gene loci at high resolution in A. viridis due to presence of many short Supplementary Material scaffolds in the genome assembly. Supplementary data are available at Genome Biology and A relatively high proportion of our piRNA reads showed a Evolution online. strong enrichment for uridine at 5 ends (1 U) and a higher probability to carry an adenine at the nucleotide 10 (10 A). In addition, the majority of sense and antisense putative Acknowledgments piRNA reads showed an overlap by exactly ten nucleotides. This bidirectional production of piRNA reads with 10 nt offset We thank to Sebastian Uhrig, Johannes Gutenberg University indicated a ping–pong dependent piRNA biogenesis. of Mainz, for advices using PingPongPro software tool, Inuk However, only a small fraction of the putative piRNA reads Jung, Seoul National University, for assistance in running that mapped to transposable elements showed ping–pong piClust software tool, and Professor Don Gilbert, Indiana signatures. Although we could detect nearly 22% putative University, for advices on filtering of the transposable element piRNAs potentially targeting transposons, only around 10% search output. We also thank members of the RAMP research showed ping–pong signatures in all pH conditions studied. group at UiT and the Genomics group at Nord University for This might be caused by very high divergence of transposable practical support and discussions. We also thank two anony- elements in A. viridis, an observation made previously in mous reviewers for their valuable comments and suggestions Exaiptasia sp. (Baumgarten et al. 2015). Another possible ex- that helped us to improve the manuscript. This work was planation is that piRNAs may fulfil various functions mainly supported by grants from the Research Council of Norway during development or in female adults. Here, TE-targeting (CoralSeq; to S.D.J.); and Tromsø Research Foundation (to piRNAs could be connected to the process of oogenesis and S.D.J.). The publication charges for this article have been serve the maintenance of the germline genome, as recently funded by a grant from the publication fund of UiT The reported by Praher et al. (2017). However, it indicates that Arctic University of Norway. piRNAs in cnidarians may also have additional function to that of transposable element silencing, a notion supported by Author Contributions observations in Hydra (Krishna et al. 2013; Juliano et al. 2014). More detailed characteristics of the putative piRNA I.U. and S.D.J. designed the study; I.U. collected, analyzed and population remain to be elucidated once better genome as- interpreted the data; S.F., H.P., B.O.K., and I.U. designed sembly and gene predictions are available for A. viridis. One workflows for the data analyses; H.P. wrote two Perl scripts important aim of our study was to identify and assess differ- for processing of miRNA sequencing data; I.U. and J.M.H.-S. entially expressed small RNAs along a natural seawater pH organized and performed the fieldwork; T.E.J. sequenced the gradient. We found high amounts of putative piRNA reads in small RNA libraries; all authors helped I.U. and S.D.J. prepare all the different pH conditions. Although it was difficult to the manuscript for publication. All authors (except S.F.) 424 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Haas BJ, et al. 2013. De novo transcript sequence reconstruction from reviewed, commented and approved the final manuscript for RNA-seq using the Trinity platform for reference generation and anal- publication. S.F. recently passed away; S.F. reviewed, com- ysis. Nat Protoc. 8(8):1494–1512. mented and approved an earlier version of the final Han BW, Wang W, Li C, Weng Z, Zamore PD. 2015. Noncoding RNA. manuscript. piRNA-guided transposon cleavage initiates Zucchini-dependent, phased piRNA production. Science 348(6236):817–821. Heo I, et al. 2012. Mono-uridylation of pre-microRNA as a key step in the Literature Cited biogenesis of group II let-7 microRNAs. Cell 151(3):521–532. Horwitz R, Borell EM, Yam R, Shemesh A, Fine M. 2015. Natural high Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. 1990. Basic local alignment search tool. J Mol Biol. 215(3):403–410. pCO increases autotrophy in Anemonia viridis (Anthozoa) as revealed Altschul SF, et al. 1997. Gapped BLAST and PSI-BLAST: a new generation from stable isotope (C, N) analysis. Sci Rep. 5:8779. Houwing S, et al. 2007. A role for Piwi and piRNAs in germ cell mainte- of protein database search programs. Nucleic Acids Res. 25(17):3389–3402. nance and transposon silencing in Zebrafish. Cell 129(1):69–82. Johnson VR, et al. 2013. Responses of marine benthic microalgae to ele- Antoniewski C. 2014. Computing siRNA and piRNA overlap signatures. vated CO . Mar Biol. 160(8):1813–1824. Methods Mol Biol. 1173:135–146. Aranda M, et al. 2016. Genomes of coral dinoflagellate symbionts high- Juliano CE, et al. 2014. PIWI proteins and PIWI-interacting RNAs function in Hydra somatic stem cells. Proc Natl Acad Sci U S A. 111(1):337–342. light evolutionary adaptations conducive to a symbiotic lifestyle. Sci Rep. 6:39734. Jung I, Park JC, Kim S. 2014. piClust: a density based piRNA clustering algorithm. Comput Biol Chem. 50:60–67. Aravin A, et al. 2006. A novel class of small RNAs bind to MILI protein in Kajitani R, et al. 2014. Efficient de novo assembly of highly heterozygous mouse testes. Nature 442(7099):203–207. Bao W, Kojima KK, Kohany O. 2015. Repbase update, a database of genomes from whole-genome shotgun short reads. Genome Res. 24(8):1384–1395. repetitive elements in eukaryotic genomes. Mob DNA. 6:11. Bartel DP. 2004. MicroRNAs: genomics, biogenesis, mechanism, and func- Kawamura Y, et al. 2008. Drosophila endogenous small RNAs bind to Argonaute 2 in somatic cells. Nature 453(7196):793–797. tion. Cell 116(2):281–297. Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E. 2007. The role of site Bartel DP. 2009. MicroRNAs: target recognition and regulatory functions. Cell 136(2):215–233. accessibility in microRNA target recognition. Nat Genet. 39(10):1278–1284. Baumgarten S, et al. 2015. The genome of Aiptasia, a sea anemone model for coral symbiosis. Proc Natl Acad Sci U S A. 112(38):11893–11898. Kozomara A, Griffiths-Jones S. 2014. miRBase: annotating high confi- Boatta F, et al. 2013. Geochemical survey of Levante Bay, Vulcano Island dence microRNAs using deep sequencing data. Nucleic Acids Res. 42(Database issue):D68–D73. (Italy), a natural laboratory for the study of ocean acidification. Mar Pollut Bull. 73(2):485–494. Krishna S, et al. 2013. Deep sequencing reveals unique small RNA reper- toire that is regulated during head regeneration in Hydra magnipapil- Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30(15):2114–2120. lata. Nucleic Acids Res. 41(1):599–616. Brennecke J, et al. 2007. Discrete small RNA-generating loci as master Langmead B, Trapnell C, Pop M, Salzberg SL. 2009. Ultrafast and memory- efficient alignment of short DNA sequences to the human genome. regulators of transposon activity in Drosophila.Cell 128(6):1089–1103. Genome Biol. 10(3):R25. Li C, et al. 2015. Multiple roles of microRNA-100 in human cancer and its Chapman JA, et al. 2010. The dynamic genome of Hydra.Nature therapeutic potential. Cell Physiol Biochem. 37(6):2143–2159. 464(7288):592–596. Christodoulou F, et al. 2010. Ancient animal microRNAs and the evolution Liew YJ, et al. 2014. Identification of microRNAs in the coral Stylophora pistillata. PLoS One 9(3):e91101. of tissue identity. Nature 463(7284):1084–1088. Das PP, et al. 2008. Piwi and piRNAs act upstream of an endogenous Lin S, et al. 2015. The Symbiodinium kawagutii genome illuminates dino- siRNA pathway to suppress Tc3 transposon mobility in the flagellate gene expression and coral symbiosis. Science 350(6261):691–694. Caenorhabditis elegans germline. Mol Cell. 31(1):79–90. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. 2012. Lowe TM, Eddy SR. 1997. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 25(5):955–964. 40(1):37–52. Min XJ, Butler G, Storms R, Tsang A. 2005. OrfPredictor: predicting protein-coding regions in EST-derived sequences. Nucleic Acids Res. Gajigan AP, Conaco C. 2017. A microRNA regulates the response of corals to thermal stress. Mol Ecol. 26(13):3472–3483. 33(Web Server issue):W677–W680. Mohn F, Handler D, Brennecke J. 2015. Noncoding RNA. piRNA-guided Ghildiyal M, Zamore PD. 2009. Small silencing RNAs: an expanding uni- verse. Nat Rev Genet. 10(2):94–108. slicing specifies transcripts for Zucchini-dependent, phased piRNA bio- Go ¨ tz S, et al. 2008. High-throughput functional annotation and data min- genesis. Science 348(6236):812–817. Moran Y, Praher D, Fredman D, Technau U. 2013. The evolution of miRNA ing with the Blast2GO suite. Nucleic Acids Res. 36(10):3420–3435. Grajales A, Rodrıguez E. 2014. Morphological revision of the genus pathway protein components in Cnidaria. Mol Biol Evol. 30(12):2541–2552. Aiptasia and the family Aiptasiidae (Cnidaria, Actiniaria, Moran Y, et al. 2014. Cnidarian microRNAs frequently regulate targets by Metridioidea). Zootaxa 3826(1):55–100. Gregory RI, Chendrimada TP, Cooch N, Shiekhattar R. 2005. Human RISC cleavage. Genome Res. 24(4):651–663. Moran Y, Agron M, Praher D, Technau U. 2017. The evolutionary origin of couples microRNA biogenesis and posttranscriptional gene silencing. Cell 123(4):631–640. plant and animal microRNAs. Nat Ecol Evol. 1(3):27. Pearson WR, Lipman DJ. 1988. Improved tools for biological sequence Grimson A, et al. 2008. Early origins and evolution of microRNAs and Piwi- comparison. Proc Natl Acad Sci U S A. 85(8):2444–2448. interacting RNAs in animals. Nature 455(7217):1193–1197. Gunawardane LS, et al. 2007. A slicer-mediated mechanism for repeat- Praher D, et al. 2017. Characterization of the piRNA pathway during de- velopment of the sea anemone Nematostella vectensis. RNA Biol. associated siRNA 5 end formation in Drosophila.Science 315(5818):1587–1590. 7:1–15. Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 425 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE Putnam NH, et al. 2007. Sea anemone genome reveals ancestral eume- Smit AFA, Hubley R. 2008–2015. RepeatModeler Open-1.0. Available tazoan gene repertoire and genomic organization. Science from: http://www.repeatmasker.org (last accessed January 14, 2018). 317(5834):86–94. Smit AFA, Hubley R, Green P. 2013–2015. RepeatMasker Open-4.0. Roberts JT, Cardin SE, Borchert GM. 2014. Burgeoning evidence indicates Available from: http://www.repeatmasker.org (last accessed January that microRNAs were initially formed from transposable element 14, 2018). sequences. Mob Genet Elements. 4:e29255. Sokol NS. 2012. Small temporal RNAs in animal development. Curr Opin Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor Genet Dev. 22(4):368–373. package for differential expression analysis of digital gene expression Suggett DJ, et al. 2012. Sea anemones may thrive in a high CO world. data. Bioinformatics 26(1):139–140. Glob Chang Biol. 18(10):3015–3025. Schwarz DS, et al. 2003. Asymmetry in the assembly of the RNAi enzyme Sunkar R, Chinnusamy V, Zhu J, Zhu JK. 2007. Small RNAs as big players in complex. Cell 115(2):199–208. plant abiotic stress responses and nutrient deprivation. Trends Plant Shaw WR, Armisen J, Lehrbach NJ, Miska EA. 2010. The conserved miR-51 Sci. 12(7):301–309. microRNA family is redundantly required for embryonic development and Vagin VV, et al. 2006. A distinct small RNA pathway silences selfish genetic pharynx attachment in Caenorhabditis elegans. Genetics 185(3):897–905. elements in the germline. Science 313(5785):320–324. Shinzato C, et al. 2011. Using the Acropora digitifera genome to under- Vashisht D, Nodine MD. 2014. MicroRNA functions in plant embryos. stand coral responses to environmental change. Nature Biochem Soc Trans. 42(2):352–357. 476(7360):320–323. Voinnet O. 2009. Origin, biogenesis, and activity of plant microRNAs. Cell Shoguchi E, et al. 2013. Draft assembly of the Symbiodinium minutum 136(4):669–687. nuclear genome reveals dinoflagellate gene structure. Curr Biol. Wheeler BM, et al. 2009. The deep evolution of metazoan microRNAs. 23(15):1399–1408. Evol Dev. 11(1):50–68. Simpson JT. 2014. Exploring genome characteristics and sequence quality without a reference. Bioinformatics 30(9):1228–1235. Associate editor: B. Venkatesh 426 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology and Evolution Oxford University Press

Elucidating the Small Regulatory RNA Repertoire of the Sea Anemone Anemonia viridis Based on Whole Genome and Small RNA Sequencing

Free
17 pages

Loading next page...
 
/lp/ou_press/elucidating-the-small-regulatory-rna-repertoire-of-the-sea-anemone-3UigfwLrJ7
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution.
ISSN
1759-6653
eISSN
1759-6653
D.O.I.
10.1093/gbe/evy003
Publisher site
See Article on Publisher Site

Abstract

Cnidarians harbor a variety of small regulatory RNAs that include microRNAs (miRNAs) and PIWI-interacting RNAs (piRNAs), but detailed information is limited. Here, we report the identification and expression of novel miRNAs and putative piRNAs, as well as their genomic loci, in the symbiotic sea anemone Anemonia viridis. We generated a draft assembly of the A. viridis genome with putative size of 313 Mb that appeared to be composed of about 36% repeats, including known transposable elements. We detected approximately equal fractions of DNA transposons and retrotransposons. Deep sequencing of small RNA libraries constructed from A. viridis adults sampled at a natural CO gradient off Vulcano Island, Italy, identified 70 distinct miRNAs. Eight were homologous to previously reported miRNAs in cnidarians, whereas 62 appeared novel. Nine miRNAs were recog- nized as differentially expressed along the natural seawater pH gradient. We found a highly abundant and diverse population of piRNAs, with a substantial fraction showing ping–pong signatures. We identified nearly 22% putative piRNAs potentially targeting transposable elements within the A. viridis genome. The A. viridis genome appeared similar in size to that of other hexacorals with a very high divergence of transposable elements resembling that of the sea anemone genus Exaiptasia. The genome encodes and expresses a high number of small regulatory RNAs, which include novel miRNAs and piRNAs. Differentially expressed small RNAs along the seawater pH gradient indicated regulatory gene responses to environmental stressors. Key words: coastal ecology, CO seep, ocean acidification, miRNA, piRNA, transposable elements. Introduction biological function, and origin (Ghildiyal and Zamore 2009). Two major classes of small regulatory RNAs in eumetazoans Compared with Bilateria, knowledge about cnidarian small are microRNAs (miRNAs) and PIWI-interacting RNAs (piRNAs). RNAs remains scarce. To date, small RNAs have only been These classes are distinct in terms of their sizes, biogenesis, reported in four cnidarians; the non-symbiotic sea anemone The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com 410 Genome Biol. Evol. 10(2):410–426. doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Nematostella vectensis, the stony corals Stylophora pistillata (Han et al. 2015; Mohn et al. 2015). A secondary piRNA path- and Acropora digitifera,and the hydroid Hydra magnipapillata way serves for piRNA amplification by a “ping–pong loop” (Grimson et al. 2008; Wheeler et al. 2009; Chapman et al. mechanism (Brenneckeetal. 2007; Gunawardane et al. 2010; Krishna et al. 2013; Liew et al. 2014; Moran et al. 2014; 2007). The two distinct piRNA populations (primary and sec- Gajigan and Conaco 2017). ondary) show opposite orientation and complementarity in miRNAs represent a well-studied class of small RNAs that theirfirst10ntpositions (Brennecke et al. 2007; usually range in size from 20 to approximately 24 nt Gunawardane et al. 2007). piRNAs lack universal sequence (miRBase, http://mirbase.org; last accessed January 14, conservation, except that the primary and secondary piRNAs 2018). In animals, miRNAs are initially transcribed as RNA show a preference for an uracil residue at the 5 end (1 U) and polymerase II transcripts (pri-miRNAs), which are further an adenine residue at the 10th position (10 A), respectively. In processed by the RNases Drosha and Dicer into stem–loop all cnidarians investigated so far, piRNAs appear to be highly precursor miRNAs (pre-miRNAs) and mature miRNAs, re- abundant compared with miRNAs and short-interfering RNAs spectively. One strand of the mature miRNA duplex is usually (siRNAs) (Grimson et al. 2008; Krishna et al. 2013; Juliano et al. incorporated into the RNA-induced silencing complex (RISC) 2014; Moran et al. 2014; Gajigan and Conaco 2017; Praher (Schwarz et al. 2003; Gregory et al. 2005), and it guides et al. 2017). The biological role of piRNAs is not well under- the whole complex to complementary mRNA for post- stood, but their most important function seems to be guiding transcriptional gene silencing. In plants, the miRNA biogen- PIWI proteins to suppress transposon activity in animal germ esis pathway involves a Dicer-like 1 (DCL-1) protein that is cells (Brennecke et al. 2007; Gunawardane et al. 2007). piRNA responsible for both cropping and slicing miRNA precursors profiling in cnidarians (Nematostella and Hydra) suggested a in the nucleus (Voinnet 2009). The miRNA silencing mecha- similar role in transposon silencing, but proposed broader si- nism is fundamentally different in animals and plants. lencing functionalities as well (Grimson et al. 2008; Juliano et al. Although animal miRNAs usually perform translational re- 2014; Praher et al. 2017). pression through partial base-pairing to target mRNAs, plant The sea anemone Anemonia viridis exposed to natural miRNAs mostly bind with full or nearly full complementarity ocean acidification conditions appears to be physiologically leading to targeted mRNA cleavage (Bartel 2009). Cnidarian acclimatized to low pH and optimizes its energy utilization miRNAs appear to contain plant-like features in their biogen- under elevated pCO through an increased autotrophic input esis and the post-transcriptional gene silencing follows a (Suggett et al. 2012, Horwitz et al. 2015). Our recent tran- plant-like regulatory pathway (Moran et al. 2013, 2014). scriptome sequencing from the same sampling site indicates Among cnidarians, Nematostella was reported to express increased expression of stress-related transcripts, repression of 87 distinct miRNAs, compared with 26, 31, and 126 global synthesis and boost in certain retrotransposon ele- miRNAs in Acropora, Stylophora,and Hydra, respectively ments at low pH in A. viridis (Urbarova et al., unpublished (Grimson et al. 2008; Wheeler et al. 2009; Chapman et al. results). In plants, it is known that small RNAs can regulate 2010; Krishna et al. 2013; Liew et al. 2014; Moran et al. species tolerance to stress via post-transcriptional gene silenc- 2014; Gajigan and Conaco 2017). Interestingly, only one ing (Sunkar et al. 2007). We therefore wanted to elucidate if miRNA (miR-100) was found conserved between the bilat- acclimatization responses of A. viridis that we observe at the erian and the cnidarian species. miRNAs have several impor- transcriptome level could be caused by small RNA-mediated tant regulatory roles in plants and animals (Bartel 2004; post-transcriptional regulation. Here, we report whole ge- Ghildiyal and Zamore 2009; Vashisht and Nodine 2014). nome and small RNA library sequencing of the symbiotic Expression profiling indicated the cnidarian miRNAs to be sea anemone A. viridis, sampled at a natural seawater pH involved in developmental regulation, regeneration and gradient off Vulcano Island, Sicily—Italy. We mainly aimed thermal stress resilience (Krishna et al. 2013; Moran et al. to identify novel small RNA species in A. viridis, but also elu- 2014; Gajigan and Conaco 2017). However, their roles in cidate their possible involvement in the acclimatization other biological processes, including other environmental responses to low pH conditions. We detected 70 distinct stress responses, have not been investigated in detail. miRNA species, and assessed differentially expressed small piRNAs are usually between 23 and 30 nt in size (Krishna RNAs. Most of the putative piRNAs contained features typical et al. 2013; Liew et al. 2014; Moran et al. 2014). The single- of primary piRNAs and a large fraction showed ping–pong stranded piRNA precursors are either derived from transposable signatures. Our study indicates possible regulatory gene elements or from specific piRNA genomic clusters, and they do responses of small RNAs to low pH. not require Dicer nuclease activity for their processing (Vagin et al. 2006; Houwing et al. 2007; Das et al. 2008). piRNAs Materials and Methods represent a highly diverse class of small regulatory RNAs, reach- Sampling ing several thousand distinct members within a single organism (Aravin et al. 2006; Kawamura et al. 2008). The uniqueness of The temperate symbiotic sea anemone A. viridis (the piRNAs arises from phased production of primary piRNAs Snakelocks Anemone) was collected at Levante Bay, North Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 411 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE Vulcano Island, Sicily—Italy. Acidification conditions are cre- database of proteins encoded by transposable elements ated here by the release of CO into the seawater from a (TEs; NCBI keywords: retrotransposon, transposase, reverse natural vent site at 1 m depth (Boatta et al. 2013; transcriptase, gypsy, and copia). These two databases were Johnson et al. 2013). Samplingwas performedonMay 13 separately queried against our reference genome assembly and 14, 2013, at the depth of 1–2 m at>350 m from the vent and the best annotation was chosen based on alignment cov- site along a gradient of decreasing pH ( pH 7.6 and 7.9), and erage and score. A combined tabular output from the at a control location at 800 m from the vent site with pH searches was further run through two Perl scripts, corresponding to ambient seawater levels (pH 8.2). For sim- “blast92gff3.pl” with additional options -lowscore 0.0001 - plicity, we are referring to average pH values throughout this alignmax 9999 -exonType exon (http://arthropods.eugenes. work as reported in Johnson et al. (2013).A total of nine org/EvidentialGene/evigene/scripts/blast92gff3.pl; last accessed individuals of A. viridis were sampled in 2 days (three from January 14, 2018) and the “overbestgene2.pl” (http://iubio. each location). Small pieces of tissue (0.5 cm) from body bio.indiana.edu/gmod/tandy/perls/overbestgene2.perl; last wall, tentacles and oral disc of each individual were collected accessed January 14, 2018) to create a gff file from blast results and stored separately at 4 C in RNAlater (ThermoFisher and to remove overlapping blast hits, respectively. The results Scientific, Waltham, MA, USA) during transport from the were imported into IBM SPSS Statistics software (version 23), sampling site to laboratory. Then, RNAlater solution was re- where counting of transposable elements was performed. moved and all samples were frozen at 80 C before further Sequence regions corresponding to transposable elements in processing steps. our reference genome assembly were then extracted from our scaffolds using BLAST fastacmd tool (Altschul et al. 1990, 1997) and used as a reference for piRNA analyses. Reference Genome Assembly DNA from one individual of A. viridis (pH 8.2) was extracted RNA Extraction using Wizard(R) Genomic DNA Purification Kit (Promega, Each tissue sample (without excess RNAlater solution) was Madison, Wisconsin, USA). Two whole genome paired-end immediately transferred from 80 Cto 1 ml cold TRIzol re- libraries (2  150 bp) were constructed and sequenced on agent (ThermoFisher Scientific, Waltham, MA, USA). The tis- Illumina HiSeq2500 at Eurofins MWG Operon (Germany). sue was then crushed using Precellys tissue homogenizer at The paired-end reads were processed using Trimmomatic 6,000 rpm for 30 s (Stretton Scientific, Stretton, UK) to mini- (Bolger et al. 2014). Adapters were removed and reads mize degradation of RNA. RNA was twice extracted by chlo- were trimmed and quality filtered using sliding window roform, and subsequently precipitated in isopropanol at 4 C with Phred score> 20. A bias at the first nine nucleotides overnight, washed with 70% ethanol, and rehydrated in was removed by trimming these bases, and reads with Nuclease-Free Water (ThermoFisher Scientific, Waltham, length<40 bp were discarded. SGA preqc tool was run pool- MA, USA). The RNA quality was examined using the Agilent ing the forward and reverse reads from the two libraries 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA, together (Simpson 2014). Platanus, a de novo genome assem- USA) and quantity of the samples was measured using bler for highly heterozygous diploidic organisms, was then Qubit 2.0 fluorometer (ThermoFisher Scientific, Waltham, used to assemble reads with k-mer length of 51 (Kajitani MA, USA). Only high quality samples with RNA integrity num- et al. 2014). To be able to assess our A. viridis genome as- ber (RIN) equal to 7 or higher were used in library sembly for repeat-enriched regions using RepeatMasker (Smit constructions. et al. 2013–2015), we first filtered out short reads from our reference assembly using N75 statistics, resulting in 210,233 Small RNA Sequencing sequences with sizes larger than 173 bp. The filtered assembly was then assessed for repeat-enriched regions using Nine individuals of A. viridis representing three different pH RepeatMasker (Smit et al. 2013–2015), with a custom library conditions (8.2, 7.9, and 7.6) were included in small RNA created by RepeatModeler, which integrates RECON, sequencing. Total RNA from three different tissue samples RepeatScout, and Tandem Repeats Finder (TRF) de novo re- of an individual was pooled at equal amounts. The small peat finding tools to build a repeat library for an assembly RNA fraction was enriched using PureLink miRNA Isolation (Smit and Hubley 2008–2015). In addition to RepeatMasker Kit (ThermoFisher Scientific, Waltham, MA, USA). Libraries annotation, the repeat-enriched regions were extracted from were prepared only from high quality RNA samples the assembly and transposable element annotation was per- (RIN 7) following the SOLiD Total RNA-Seq Kit protocol formed as described previously (Chapman et al. 2010; (ThermoFisher Scientific, Waltham, MA, USA). Different Baumgarten et al. 2015). The annotation pipeline included A. viridis small RNA libraries were barcoded and sequencing TM then also a TBLASTX run using RepBase database (Bao et al. was performed on three lanes of a SOLiD 6-Lane FlowChip 2015), version 22.09 (e-value< 10 ), and a BLASTX search using SOLiD 5500xl sequencer at the Nord University (Bodø, (e-value< 10 ) against a custom-made non-redundant Norway). 412 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Discovery of Novel miRNA (ORFs). Our scaffolds were searched for ORFs using OrfPredictor (http://bioinformatics.ysu.edu/tools/OrfPredictor. After removing low quality (quality score< 18) and less com- html; last accessed January 14, 2018) (Min et al. 2005). plex sequences from our raw small RNA sequencing data set, adapter sequences were trimmed away using trimSOLiDAdaptor.pl Perl script keeping only sequences equal miRNA Analyses to or longer than 18 nt. The filtered reads were mapped in Expression of selected miRNAs was confirmed by quantitative color space using Bowtie (Langmead et al. 2009)to our PCR (qPCR). Six Locked Nucleid Acid (LNA) probes targeting A. viridis genome reference with parameters –integer-quals the predicted miRNAs were designed using online miRNA -l 18 -M 20 –best –strata -e 150 –nomaqround –maxbts 800 – qPCR designer tool from Exiqon (Vedbaek, Denmark). cDNA tryhard -a –col-cqual –col-keepends –mapq 20 –threads 14 – was synthesized from three individuals (10 ng of total RNA chunkmbs 200. These options select the best alignments TM each) per condition using miRCURY LNA Universal RT based on the seed mismatches only and mismatches outside microRNA PCR, Polyadenylation, and cDNA synthesis kit II the seed region are ignored. Therefore, we needed to per- (Exiqon, Vedbaek, Denmark) following the instruction man- form additional filtering using processBowtieAlignments.pl ual. Small RNA for the qPCR analysis was isolated from the Perl script to select alignments with the minimum mismatches same samples that were used for preparation of small RNA along the whole reads. Both Perl scripts are available at https:// libraries. qPCR was performed in duplicates using miRCURY github.com/patelhardip/bitx.git (last accessed January 14, LNA microRNA PCR, ExiLENT SYBR Green master mix (Exiqon, 2018). Mapped reads from each condition were then prepro- Vedbaek, Denmark) in 10 ml. miRNAs were assessed for dif- cessed by bwa_sam_converter.pl Perl script (Friedl€anderetal. ferential expression among the sampling sites with differing 2012), outputting two files essential for running miRDeep2 pH by edgeR (FDR< 0.05) (Robinson et al. 2010). Mature software tool for the novel miRNA predictions (Friedlander miRNAs were aligned to their precursor sequences. miRNAs et al. 2012). The two files were used as input into with <20 counts in less than three conditions were not con- miRDeep2.pl Perl script (Friedlander et al. 2012)that was sidered. We searched for putative animal-like miRNA targets run for identification of novel miRNAs in A. viridis. In addition, by Probability of Interaction by Target Accessibility (PITA) soft- known miRNAs of three related species, N. vectensis, S. pis- ware, based on target complementarity and site accessibility tillata,and H. magnipapillata (Krishna et al. 2013; Liew et al. (Kertesz et al. 2007). Coding regions were predicted from the 2014; Moran et al. 2014) thatwereavailable at thetimeof A. viridis transcriptome (Urbarova et al., unpublished results) the analysis, have been used in the prediction pipeline. These by TransDecoder, version 2.0.1 (Haas et al. 2013) and used as miRNA sequences were downloaded from miRBase, release input into PITA software. The results were filtered based on 21 (Kozomara and Griffiths-Jones 2014). Output from the change in free energy (ddG) of miRNAs binding to its miRDeep2 software was then inspected manually, keeping targets (ddG<10 kcal/mol), seed length of 8 nt with no only predicted miRNAs with miRDeep2 score larger than ten mismatches and no wobble pairs. Only targets fulfilling these and with significant randfold value (p value< 0.05). Small criteria were considered further. We then checked the pre- RNA sequencing data from each individual were assessed sep- dicted targets from PITA for more extensive complementarity arately for the presence of novel miRNAs. Only miRNAs iden- to the miRNAs usingFASTA v36(Pearson andLipman1988) tifiedinat least twoindividuals were considered further. as previously described (Moran et al. 2014), and scored the Possible tRNA contamination was examined by running alignments accordingly. Blast hits were obtained using NCBI tRNAscan on the reference genome (Lowe and Eddy 1997). nr database (e-value< 10 ) and GO terms were assigned to Further, presence of rRNA sequences in predicted hairpin the potential mRNA targets using B2G4Pipe Blast2GO pipe- structures was tested by querying against a custom database line (Go ¨ tz et al. 2008). combining known rRNA sequences from N. vectensis and A. viridis. No contamination was found in either case. In addition, Search for Putative piRNAs to ensure that our miRNA candidates come from the host, assembled scaffolds of A. viridis were screened for their pos- Raw reads from SOLiD sequencing were quality filtered and sible contamination by symbiont DNA using genomes of adapter sequences were trimmed, as described previously for Symbiodinium minutum, Symbiodinium microadriaticum, miRNA discovery. However, only sequences equal to or longer and Symbiodinium kawaguti (Shoguchi et al. 2013; Lin than 23 nt were kept for the piRNA analyses. Sequences were et al. 2015; Aranda et al. 2016). 1,642 scaffolds (mainly short further filtered for reads mapping to miRNA precursors iden- ones with length 100 bp) that were highly similar (e-val- tified in the present study and reads mapping to rRNAs using ue< 10 )toSymbiodinium genomes in the A.viridis genome rRNA databases according to Praher et al. (2017). Sequences assembly were filtered out prior to the small RNA alignment. in color space were aligned using Bowtie with the same Genomic setting of aligned putative miRNAs was inspected for parameters as for the miRNA alignment, but allowing for overlapping regions corresponding to open reading frames maximum three mismatches with the “seed length” of Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 413 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE 23 nt (-l 23 -n 3). We refrained from mapping our reads to reference was found to contain repetitive regions. About unique locations in the A. viridis reference genome due to 27.5% of the repetitive sequences could be assigned to pre- presence of many sequence stretches in assembled scaffolds viously known repetitive elements, but most of these sequen- that most probably correspond to the same genomic loca- ces (25% of the genome) could not be classified into any tions. The aligned sequences were filtered for the best assigned category (supplementary table S1, Supplementary matches using the same custom Perl script as for filtering Material online). Repeat annotation identified only about the miRNA alignments and filtered for reads with 1 U or 8.3% of genome to be comprised of transposable elements 10 A sequence signatures. TE-targeting potential of putative (TEs). These are similar observations as made previously for piRNAs was then assessed by including only putative piRNAs symbiotic sea anemone Exaiptasia sp. (formerly known as mapping antisense to the transposable elements. Overlap Aiptasia sp.; Grajales and Rodrıguez 2014; Baumgarten probabilities of the putative piRNA sequences with opposite et al. 2015). From the identified TE fraction, about half orientation were analyzed using signature.py script (44.4%) were retrotransposons, and amongst them, non-long (Antoniewski 2014). The script computes the probability of terminal repeat (non-LTR) retrotransposons were predominating an antisense read overlapping a sense read with defined (supplementary table S1, Supplementary Material online). length and assigns each overlap length a z-score. The over- lapping signatures of the putative piRNAs were inspected in miRNA Discovery more detail by running PingPongPro v1.0 software (http:// sourceforge.net/projects/pingpongpro/; last accessed January Small RNA libraries from nine individual polyps of A. viridis, 14, 2018). This software also served for inspection of trans- sampled at three different seawater pH conditions (at normal poson silencing by putative piRNAs. Silencing of transposable seawater pH 8.2, as well as at pH 7.9 and 7.6) were prepared elements by the putative piRNAs was only considered if FDR (q and subjected to sequencing on the SOLiD 5500xl platform value)< 0.01 and if it was supported by at least ten putative (fig. 1). The sequencing generated approximately 116 million piRNA reads with ping–pong signatures normalized to the reads (18 nt) after adapter trimming and quality filtering transposon length. To inspect data for presence of piRNA (table 2 and fig. 2). Despite the fragmented nature of our clusters, the putative piRNA reads were mapped onto the genome draft assembly, a high proportion of the small RNA masked genome reference produced by RepeatMasker reads mapped to the genomic reference (between 88% and (Smit et al. 2013–2015) as described previously in color space 92%, table 2). This indicates that even a very preliminary ge- using Bowtie, reporting at maximum five valid alignments. nome assembly is sufficient for the discovery of small RNAs. Finally, our sorted alignment files were submitted to piClust A. viridis miRNAs were identified by the miRDeep2 software software (Jung et al. 2014). Here, the Eps parameter was set tool (Friedlander et al. 2012), and a representative analysis to 1000 and MinReads to 50. result is shown in figure 3A. We predicted in total 70 high- confidence miRNA candidates (20 to 25 nt) for A. viridis, including 61, 60, and 65 distinct miRNA species at pH 8.2, 7.9, and 7.6, respectively (table 3). Most miRNAs were detected Results in all pH conditions studied (n¼ 51), 14 miRNAs in two differ- Genome Reference Assembly and Search for Repeat- ent pH conditions, and five miRNAs were detected in one pH Enriched Regions condition only (table 3). Eight candidate miRNAs in A. viridis Total DNA from a single A. viridis polyp (normal seawater were apparently homologous to those reported in other cni- conditions, pH 8.2) was extracted and subjected to whole darian species (avi-miR-temp-100, 2022, 2023, 2025, 2028, genome sequencing on the Illumina HiSeq2500 platform 2030, 2036, and 2037) (fig. 3B). The predicted miRNA with the (fig. 1). Sequencing generated 43 billion nucleotides (nt) of highest miRDeep2 score (avi-miR-temp-100), and which was genomic data, which corresponded to 144 million paired-end highly expressed in all pH conditions, was identical in sequence reads (table 1). The basic genome characteristics were deter- to miR-100 in N. vectensis and S. pistillata (fig. 3B)(Grimson mined using the SGA preqc software tool (Simpson 2014), et al. 2008; Liew et al. 2014; Moran et al. 2014). Similarly, avi- showing an estimated genome size of 313 Mb (140x cov- miR-temp-2022, 2023, and 2025 were identical to the corre- erage). Adapters were trimmed and the reads were quality sponding miRNAs in N. vectensis (Moran et al. 2014). Other A. filtered before de novo genome assembly, which created viridis miRNAs (avi-miR-temp-2028, 2030, 2036, and 2037) about 1.1 million short scaffolds with N50¼ 2,087. This ge- have one or two nucleotides substitution compared with that nome assembly is highly fragmented, but it was sufficient for of the N. vectensis homolog (fig. 3B). Multiple precursor the mapping of small RNAs and the identification of transpos- sequences for some of the predicted miRNAs were found by able elements (fig. 1). miRDeep2 (table 3). All mature, star, and precursor sequences The A. viridis genome assembly was inspected for repeat are presented in supplementary table S2, Supplementary and low complexity regions using the RepeatMasker software Material online, with read counts for each pH condition in tool (Smit et al. 2013–2015), and about 36% of the genome supplementary table S3, Supplementary Material online. 414 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE 9 individuals RNA isolation (3 individuals/condition) 1 individual (control condition) DNA isolation body wall tentacles mouth piece Illumina HiSeq 2500 sequencing Pooling equal RNA amounts Small RNA enrichment Trimming and quality filtering Preparation of barcoded libraries Genome assembly with Platanus SOLiD 5500xl sequencing putative piRNA pipeline Identification of repetitive DNA by RepeatMasker Trimming and quality filtering Size filtering (≥ 23 nt) miRNA pipeline Annotation of repeat-enriched Size filtering (≥ 18 nt) Mapping to genome and fraction, including transposable elements transposable elements in color space Mapping to genome in color space Assessment of Novel miRNA predictions ping-pong signatures by miRDeep2 Differential expression with edgeR qPCR verification FIG.1.—Data analysis overview. DNA and RNA were isolated from A. viridis adult polyps sampled from a natural seawater pH gradient (at normal seawater pH 8.2, and at low seawater pH 7.9 and 7.6) off Vulcano Island, Sicily—Italy. Only one polyp (pH 8.2) was used for DNA extraction and was subjected to paired- end sequencing on the Illumina HiSeq2500 platform. Sequencing reads were assembled into a draft genome reference. Subsequently, repeat-enriched regions, including transposable elements were identified and annotated in this assembly. Nine polyps were used for small RNA library preparation and sequencing on the SOLiD 5500xl platform. Sequencing reads were further used for novel miRNA discovery and description of putative piRNA reads. Table 1 of putative miRNA clusters. We identified four clusters, three The Amount of Reads Gained from Genome Sequencing of Anemonia contained two miRNA sequences (avi-miR-temp-11 and 66; viridis avi-miR-temp-28 and 27; and avi-miR-temp-64 and 67) and Sequencing No. of No. of %Paired-End one cluster contained three miRNAs (avi-miR-temp-2, 13, and Index Paired-End Trimmed and Reads Kept After 39) (supplementary fig. S1, Supplementary Material online). Raw Reads Quality Filtered Filtering No open reading frames (ORFs) spanning cluster regions could Reads be predicted, implying that the miRNAs were transcribed as CTTGTA 82,428,617 71,676,503 87.0 independent transcription units. Further clustering of miRNA GCCAAT 61,246,666 52,649,432 86.0 loci could not be assessed due to the fragmented nature of Sequencing of barcoded genome libraries was performed in one lane of the genome assembly. Illumina HiSeq2500 sequencing machine in 2 150 bp mode. The amount of sequen- ces presented here is the number of raw paired-end sequences obtained after the run. We then asked if any of the expressed miRNAs were colocal- ized with predicted transposable elements in the A. viridis ge- nome reference. One miRNA (avi-miR-temp-58) was found Genome Context of Identified miRNAs encoded within a DNA transposon (fig. 4A). avi-miR-temp-58 In sea anemones, little is known about the genomic miRNA was detected only at pH 7.9 in the small RNA sequencing exper- clusters that generate pri-miRNAs. Therefore, we searched iment (but in all three individuals inspected). However, we the A. viridis reference genome sequence for the presence detected avi-miR-temp-58 in all conditions studied by a Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 415 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE Table 2 The Amount of Reads Gained from Small RNA Sequencing of Anemonia viridis Individuals Raw Filtered Small % Reads Aligned Filtered Small Reads Aligned Putative piRNA %Putative piRNAs Reads RNA Reads to the Genome RNA Reads to Genome Reads Aligned Aligned to b b c (18 nt) (18 nt) (23 nt) (23 nt) to Genome Genome pH 7.6–1 16,111,517 12,117,933 88.5 9,831,771 7,681,304 6,643,549 86.5 pH 7.6–2 12,492,090 11,175,766 89.5 9,465,543 7,721,730 6,795,032 88.0 pH 7.6–3 13,794,761 10,609,870 89.1 8,921,147 7,238,217 6,310,624 87.2 pH 7.9–1 18,815,344 15,837,619 91.5 13,235,131 11,175,365 10,001,436 89.5 pH 7.9–2 20,094,253 17,369,794 90.0 16,247,704 13,837,596 11,960,594 86.4 pH 7.9–3 16,794,089 15,567,040 89.3 14,958,985 12,687,556 11,353,938 89.5 pH 8.2–1 13,066,319 12,109,753 88.0 11,317,833 9,184,500 8,163,612 88.9 pH 8.2–2 17,379,705 10,309,153 90.3 8,124,195 6,661,625 5,543,870 83.2 pH 8.2–3 13,492,380 11,225,436 90.3 8,769,865 7,189,261 6,271,324 87.2 a TM Small RNA libraries from each individual were barcoded, pooled, and sequencing was performed on three lanes of a SOLiD 6-Lane FlowChip using the SOLiD 5500xl sequencer. The amount of sequences presented here is the sum of raw reads from the three lanes. These reads had adapter removed and were quality filtered. They differ only according to the size filtering. Reads (23 nt) aligned to the reference genome and filtered for piRNA sequence signatures (1 U and 10 A). 6e+06 (Robinson et al. 2010). Here, avi-miR-temp-37, 52, 56, 58, and 59 appeared up-regulated at pH 7.9, whereas avi-miR- 5e+06 temp-13, 29, 48, and 60 appeared down-regulated at pH 7.9 4e+06 (fig. 5). Six miRNAs were then selected for verification analysis by qPCR (avi-miR-temp-37, 58, 60, 100, 2023, and 2028), 3e+06 where three miRNAs homologous to Nematostella with ap- parently unaffected expression levels in the different pH con- 2e+06 ditions served as controls (supplementary fig. S2, 1e+06 Supplementary Material online). The control miRNAs (avi- miR-temp-100, 2023, and 2028) were detected by qPCR in all pH conditions at similar expression levels, and thus are in 18 19 20 20 21 22 23 24 25 26 27 28 29 30 30 31 32 33 34 35 good agreement with results generated from small RNA se- quencing. The miRNA avi-miR-temp-37 was detected only at sequence length in nt pH 7.9 and only in one individual at pH 7.6, and avi-miR- temp-60 was detected only at pH 8.2. However, in contrast FIG.2.—Sequence length representation of small RNAs in A. viridis. Distribution of small RNA reads after adapter trimming and quality filtering with the observation from small RNA sequencing, we in one individual sampled from pH 8.2. Two distinct peaks could be ob- detected the presence of avi-miR-temp-58 in all conditions served; first around 22 nt representing both miRNA and siRNA reads, and studied, though in higher abundance at pH 7.9 (supplemen- second around 28 nt representing putative piRNA reads. tary fig. S2, Supplementary Material online). We then searched for putative mRNA targets of 13 se- quantitative PCR (qPCR) approach (see below), though in higher lected miRNAs that were differentially expressed along the abundance at pH 7.9 (supplementary fig. S2, Supplementary pH gradient (avi-miR-temp-13, 29, 37, 48, 52, 56, 58, 59, Material online). avi-miR-temp-58 and its precursor was pre- and 60) (fig. 5), detected only in one pH condition (avir- dicted to create a 1 nt 3 overhang (fig. 4B). The latter feature miR-temp-48, 58, 59, 60, and 65), or detected only at low suggests a group II pre-miRNA that requires a 3 -end mono- pH, that is, at pH 7.6 and/or pH 7.9 (avi-miR-temp-37, 48, 50, uridylation for further Dicer processing (Heo et al. 2012). 57, 58, 59, 60, 64, and 65) (table 3). Differentially expressed mRNAs previously identified in A. viridis in the same individu- als from the same sampling experiment (Urbarova et al., Differential miRNA Expression upon Seawater pH Gradient unpublished results) were assessed as potential targets. All high-confidence miRNAs were included in differential ex- After stringent filtering criteria, including full seed matching pression analyses. Despite that 19 candidate miRNAs could with extended pairing, we identified 9 out of the 13 selected not be detected in all conditions studied, only nine miRNAs miRNAs that could potentially target 13 of the differentially were recognized as differentially expressed between expressed mRNAs along the low pH gradient (supplementary conditions by edgeR (FDR< 0.05) (fig. 5; supplementary table S6, Supplementary Material online). Although we could tables S4 and S5, Supplementary Material online) not consistently detect miRNA upregulation and its mRNA 416 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 number of reads Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE A B Provisional ID : scaffold25818_len18526_cov78_2565 Score total : 428917.1 hsa-miR-100 aacccguagauccgaacuugug- Score for star read(s) : 3.9 nve-miR-100 -acccguagauccgaacuugugg Score for read counts : 428907.3 spi-miR-100 -acccguagauccgaacuugugg Score for mfe : 1.3 Score for randfold : 1.6 adi-miR-100 -ucccguagauccgaacuugugg Score for cons. seed : 3 5' avi-miR-temp-100 -acccguagauccgaacuugugg Total read count : 841295 u u c u Mature read count : 799715 ccguagauc gaacuug gg g c ******************** Loop read count : 0 3' g g c a u c u a g u u u g a a c c c c Star read count : 41580 c u u a u c nve-miR-2022 uuugcuaguugcuuuugucccgc- avi-miR-temp-2022 uuugcuaguugcuuuugucccgc- freq. spi-miR-2022 uuugcuaguugcuuuugucccguu adi-miR-2022 uuugcuaguugcuuuugucccgu hma-miR-2022 uuugcuaguugcuuuuguccccu- 0.75 ********************* 0.5 nve-miR-2023 aaagaaguacaagugguaggg spi-miR-2023 aaagaaguacaagugguaggg 0.25 adi-miR-2023 aaagaaguacaagugguaggg avi-miR-temp-2023 aaagaaguacaagugguaggg 1 23 35 57 length ********************* nve-miR-2025 uuuuuuagcccgcggaaguugu Mature Star adi-miR-2025 auuuuuagcccgcggaaguugc 5'- guaggcaguguuguuugucgacccguagauccgaacuugugguugucuaccuccccaaguuuugaucuacggaacuaaaaucaacgcuaauaugaccaggaacauccaagga -3' obs avi-miR-temp-2025 uuuuuuagcccgcggaaguugu guaggcaguguuguuugucgacccguagauccgaacuugugguugucuaccuccccaaguuuugaucuacggaacuaaaaucaacgcuaauaugaccaggaacauccaagga exp ******************** ......(((((((.((......(((((((((.(((((((.((..(....)..)).))))))).)))))))))......)).))))))).......((.(((...)))..)). reads mm sample nve-miR-2028 uaauguuccugcuuguuccua ....................acccguagauccgaacuugugg...................................................................... 737600 0 seq avi-miR-temp-2028 aaauguuccugcuuguuccug ....................acccguagauccgaacuugu........................................................................ 26766 0 seq ******************* ....................acccguagauccgaacuugug....................................................................... 21383 0 seq hma-miR-2030 uagcauaacauuguaagaaaca avi-miR-temp-2030 uagcauaacauaguaagagauu ....................acccguagauccgaacuug......................................................................... 4646 0 seq nve-miR-2030 uagcauaacauuguaagagauu spi-miR-2030 uagcauaacauuguaagagauc ....................acccguagauccgaacuugugU...................................................................... 1627 1 seq adi-miR-2030 uagcauaacauuguaagagaucu *********** ****** * ......................................................ccaaguuuugaucuacggaacu.................................... 30488 0 seq spi-miR-2036 uauauuguacgacucucaucgugu ......................................................ccaaguuuugaucuacggaac..................................... 8550 0 seq adi-miR-2036 uauauuguacgacucucaucgug nve-miR-2036 uauauuguacgacucucaucgua- ......................................................ccaaguuuugaucuacggaa...................................... 1460 0 seq avi-miR-temp-2036 uauauuguacgacucucaucguag ********************** nve-miR-2037 ugugauuggagacuuuuaccgu avi-miR-temp-2037 ugugauuggagacuuuuaucgu ****************** *** FIG.3.—Identified miRNAs in A.viridis with similarity to known miRNAs. (A) A typical prediction result from miRDeep2 software tool showing the miRNA precursor, mature and star sequence and their abundances in the sample. Shown is avi-miR-temp-100 precursor with top sequence alignments in the sample for each strand. (B) Alignments of novel miRNAs from A. viridis to known miRNAs from other species. Sequences of our predicted miRNAs from A. viridis (denoted as avi-miR-temp) were aligned to known miRNA sequences from H. sapiens (hsa-miR), N. vectensis (nve-miR), H. magnipapillata (hma-miR), S. pistillata (spi-miR), and A. digitifera (adi-miR). (temp¼ temporary; miRNAs that are not yet registered in the miRBase). target downregulation, we made one interesting observation. plot (fig. 2), and most likely represent PIWI-interacting RNAs We detected avi-miR-temp-50, present only at pH 7.6 and pH (piRNAs). Based on an earlier report on piRNA signatures in 7.9, to target an RNase HI domain of a DIRS1 retrotransposon. cnidarians (Moran et al. 2014), we explored the trimmed and The corresponding transcript was found downregulated both quality filtered small RNA reads with minimum length of 23 nt at pH 7.6 and 7.9 compared with pH 8.2, which could mean and with the typical base preference signatures (hereafter that this domain is inactivated and reverse transcription is called putative piRNAs) in our data set. The putative piRNAs therefore inhibited. were aligned to two different reference data sets; the A. viridis reference genome, and the transposable elements identified in the genome. In total, about 83–90% of reads23 nt Search for Putative piRNAs and Their Characteristics aligned to the reference genome were putative piRNAs Most small RNA sequences present in our data set showed a (table 2), with about 14–42% mapping to transposable ele- distinct peak at 27–29 nt in the small RNA size distribution ments (supplementary table S7, Supplementary Material Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 417 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE Table 3 The Listof70 Predicted miRNAs in Anemonia viridis from Various pH Conditions Temporary miRNA Name Mature Sequence Length Stem–Loop Length Present in Condition Similarity to Known miRNAs avi-miR-temp-100 acccguagauccgaacuugugg 22 56 All nve-miR-100-5p avi-miR-temp-2022 uuugcuaguugcuuuugucccgc 23 52; 53 All nve-miR-2022-3p avi-miR-temp-2023 aaagaaguacaagugguaggg 21 53 All nve-miR-2023-3p avi-miR-temp-2025 uuuuuuagcccgcggaaguugu 22 53 All nve-miR-2025-3p avi-miR-temp-2028 aaauguuccugcuuguuccug 21 48 All nve-miR-2028-5p avi-miR-temp-2030 uagcauaacauaguaagagauu 22 52 All nve-miR-2030-5p avi-miR-temp-2036 uauauuguacgacucucaucguag 24 54 All nve-miR-2036-3p avi-miR-temp-2037 ugugauuggagacuuuuaucgu 22 54 All nve-miR-2037-3p avi-miR-temp-1 gaucaagucaaauacaucucu 21 48 All avi-miR-temp-2 uaucaaggcagucuuaccauau 22 54 All avi-miR-temp-3 uacaaauguuacgcagcagaac 22 55 All avi-miR-temp-4 ugacauugcugcccgaaucucc 22 88 All avi-miR-temp-5 uuuaauguuacugcucguucc 21 53 All avi-miR-temp-6 aauuucaaauauccacugauuga 23 55 All avi-miR-temp-7 uugagcaucuguugcaugucua 22 53 All avi-miR-temp-8 aucaucgccacuagcaucguca 22 55 All avi-miR-temp-9 aagggcaagacaauagaauuuca 23 59 All avi-miR-temp-10 cuugauaguacuuuugccuugc 22 52 pH 7.9, pH 8.2 avi-miR-temp-11 uaguagguucuuauaagcuauu 22 55 All avi-miR-temp-12 uauaagucuaggcugguuaaga 22 56 All avi-miR-temp-13 auacugaacuugaaagaagugau 23 55 All avi-miR-temp-14 aaacgcuguucuugguaguca 21 55 All avi-miR-temp-15 uaacaaagcaguuuggcuguau 22 55 All avi-miR-temp-16 ucuggcugauuugaagaaaga 21 51 All avi-miR-temp-17 acaucaaacaaagcaguuug 20 53 All avi-miR-temp-18 auuacccguaaauaaauucaau 22 54 All avi-miR-temp-19 aaccccaacgcgggccucugg 21 51 All avi-miR-temp-20 uuaguuugcacucauuugcugg 22 55 All avi-miR-temp-21 auuacccagaauggggccuuu 21 55 All avi-miR-temp-22 uauucuccaaaaauucacaagg 22 52 All avi-miR-temp-23 uaaacuaguugauaggauugu 21 51 All avi-miR-temp-24 acagauugcggcaaccgugcag 22 72; 88 All avi-miR-temp-25 ucaaauguugcgcagcagaac 21 55 All avi-miR-temp-26 ugcugcaguuuagacugaccuc 22 53 All avi-miR-temp-27 uccucaaguuuugauuguaauac 23 51; 52 All avi-miR-temp-28 uucuuaaguuuugauuguaauac 23 52 All avi-miR-temp-29 aucuacugauacuaaguauccg 22 54 All avi-miR-temp-30 uuucuguaguacuuuauccuggc 23 54 All avi-miR-temp-31 uauucaaucagucuggcuguua 22 52 All avi-miR-temp-32 ucuuuugauaaauaccaccaaca 23 56 All avi-miR-temp-33 uacucugaaguguacuuagugu 22 53; 54 All avi-miR-temp-34 gauaugauauaauauguaugug 22 57 All avi-miR-temp-35 uauacauauuuaguaucgauaucag 25 57; 58 All avi-miR-temp-36 uaaauacacaauaucuauagcagu 24 55 All avi-miR-temp-37 uaugguagugauguuuagaaa 21 49 pH 7.6, pH 7.9 avi-miR-temp-38 ccggacaaugagaauagcuga 21 56 All avi-miR-temp-39 ugaucaauaaaagaaacaucguu 23 54 All avi-miR-temp-42 uaucacauuuaaaacacucaug 22 53 All avi-miR-temp-43 ucauacgauauuuuucacuagu 22 55; 56 All avi-miR-temp-44 aaccucaugucagagaucaaa 21 53 pH 7.6, pH 8.2 avi-miR-temp-45 acagagccuccuuuaaccuccu 22 60 pH 7.6, pH 8.2 avi-miR-temp-47 ugguagaacaaguaacuugcugc 23 55 pH 7.6, pH 8.2 avi-miR-temp-48 gaaaaagacauuuagagacuug 22 56 pH 7.6 (continued) 418 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Table 3 Continued Temporary miRNA Name Mature Sequence Length Stem–Loop Length Present in Condition Similarity to Known miRNAs avi-miR-temp-49 aaugucaccaaguuucgacca 21 50 pH 7.6, pH 8.2 avi-miR-temp-50 aggcccuggggaaacaaugga 21 54 pH 7.6, pH 7.9 avi-miR-temp-52 uggaugcucaauuugccaauugc 23 75 pH 7.9, pH 8.2 avi-miR-temp-53 aacuuaaaacaaaaaucucccu 22 53 pH 7.6, pH 8.2 avi-miR-temp-54-1 aucuauucacugugggcguccagu 24 54 All avi-miR-temp-54-2 aucuauucauugugggcguccagu 24 55 All avi-miR-temp-55 uacuacuuugacaaugugaugg 22 53; 77 pH 7.6, pH 8.2 avi-miR-temp-56 aggucagucuaaacugcagca 21 54; 55 pH 7.6, pH 8.2 avi-miR-temp-57 gcuuugaaaauguaaagaaca 21 50 pH 7.6, pH 7.9 avi-miR-temp-58 ugcaguauucaguaugcacua 21 65 pH 7.9 avi-miR-temp-59 ucggcgccggucacgcgauaga 22 52 pH 7.9 avi-miR-temp-60 caagcuauaaauuccaacuga 21 50 pH 7.6 avi-miR-temp-61 ucgaguaaaauauuacagaaaug 23 54 pH 7.6, pH 8.2 avi-miR-temp-64 ucaucucuuguggcuugacauu 22 51 pH 7.6, pH 7.9 avi-miR-temp-65 uggugcaguuuagacugacccuu 23 54 pH 7.9 avi-miR-temp-66 cuagauuaugagagcuuaugu 21 53 All avi-miR-temp-67 ugugugaaaacaugacaagaucu 23 50 All Presence of two numbers in this column indicates that two different miRNA precursors (pre-miRNAs) of the same mature miRNAs have been detected in our genome assembly. All nucleotide sequences of mature and star miRNAs and pre-miRNA precursors are listed in supplementary table S2, Supplementary Material online. DNA transposon miRNA precursor star mature GCACCATTCAGTATGCAGTATTGAGTATGCAGTATTCAGTGTGCAGTATGCAGTGTGCAGTGTGCAGTGTGC AGTGTGCAGTGTGCAGTGTGCAGTGTGCAGTATGCAGTAAGCAGTATGCAGTGTGCAGTGTACAGTGTGCAG TATGCAGTATTCAGTGTGCAGTGTGCAGTGTGCAGTATGCAGTATGCAGTATGCAGTATGCAGTATTCAGTG TGCAGTATGCAGTATGCAGTGTGCAGTATACAGTATTCAGTATGCACTATTGAGTATGCAGTATTCAGTATG CACTATTGAGTATGCAGTATTCAGTGTGCAGTATTCAGTATGCACTATTCAGTATGCAGTATTCAGTGTGCA GTATTCAGTGTGCACTATTCAGTATGCACTATTGAGTATGCAGTATTCAGTATGCAGTATTCAGTATGCAGT ATGTAGTATGCAGTATGCAGTATGCACTATTGAGTATGCAGTATTCAGTGTGCACTATTCAGTATGCACTAT TCAGTATGCACTATTGAGTATGCAGTATTGAGTATGCA miRNA precursor hairpin structure u star sequence 5' a g c ugc a u a u u a u g u a ug c a c a g u u g a g u a 3' a c g u a u g a u u a u a c g u g c u g a g c u c u a c u g u a mature sequence FIG.4.—Precursor of avi-miR-temp-58 and its DNA transposon localization. (A) The miRNA precursor of avi-miR-temp-58 was found localized in a DNA transposon. A schematic representation of the scaffold region is depicted above the DNA sequence. Only the part of the scaffold with similarity to the DNA transposon is shown. The DNA transposon has homology to a transposable element from Crassostrea gigas. The miRNA precursor is marked in color and the whole sequence is underlined. Mature miRNA is indicated in red and star sequence in violet. (B) Hairpin structure of avi-miR-temp-58 precursor with 1 nt 3 overhang. Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 419 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 individual Urbarova et al. GBE Color Key −1.5 −1 −0.5 0 0.5 1 1.5 Row Z−Score avi-miR-temp_56 avi-miR-temp_29 avi-miR-temp_13 avi-miR-temp_52 avi-miR-temp_37 avi-miR-temp_58 avi-miR-temp_59 avi-miR-temp_48 avi-miR-temp_60 1 3 2 2 1 3 1 3 2 pH 7.9 pH 8.2 pH 7.6 pH 8.2 pH 7.6 condition FIG.5.—Differentially expressed miRNAs under low pH conditions. Nine miRNAs were found differentially expressed among the sampling sites (edgeR, FDR< 0.05). This included two miRNAs detected in all pH conditions studied (avi-miR-temp-13 and 29) and seven miRNAs that could be detected in only one or two different pH conditions. Five miRNAs were differentially expressed between pH 7.6 and 7.9, three downregulated (avi-miR-temp-52, 58, and 59) and two upregulated (avi-miR-temp-48 and 60) at pH 7.6 compared with pH 7.9. Only one miRNA was detected differentially expressed between pH 7.6 and 8.2 (avi-miR-temp-37), and it was upregulated at pH 7.6 compared with pH 8.2. Eight miRNAs were found differentially expressed between pH 7.9 and 8.2, three downregulated (avi-miR-temp-13, 29, and 48) and five upregulated (avi-miR-temp-37, 52, 56, 58, and 59) at pH 7.9 compared with pH 8.2. Differentially expressed miRNAs were hierarchically clustered into heatmap based on counts per million (cpm) and scaled by row. online), including unclassified fraction of repeats. Most of the Ping–Pong piRNA Amplification Signature in A. viridis putative piRNA reads mapped to the reference genome and We further investigated if ping–pong signatures, that is, transposable elements showed strong preference for 1 U 10 nt overlaps of putative piRNAs with opposite direction, (fig. 6), a feature consistent with the primary piRNA popula- were common in our data set. The probability of overlap by tion. Most of the putative piRNAs are found in genomic clus- 1–30 nt of putative piRNAs with opposite orientation was ters and the majority of piRNA cluster loci appear unistranded assessed. The data exhibited strong ping–pong signatures, (61–68%), where piRNAs are transcribed from one strand of since most reads showed preference for 10 nt 5 overlaps of the piRNA locus (supplementary fig. S3, Supplementary putative piRNAs with opposite orientation in all conditions, in Material online). all individuals, and for all reference data sets (z-score> 5). About one third of the genome scaffolds contained Probability of other overlaps was much lower (z-score< 1) expressed piRNA loci (supplementary table S8, (fig. 7A). We could detect 14–15% putative piRNAs map- Supplementary Material online). We observed more ping to identified transposable elements, and up to 42% (sup- expressed piRNA loci at pH 7.9 than at pH 8.2 or 7.6. plementary table S7, Supplementary Material online) when Putative piRNAs were found to map to about 24–30% of including the unclassified fraction of identified repeats (sup- the identified transposable elements in all conditions and in plementary table S1, Supplementary Material online). Only all individuals (supplementary table S9, Supplementary around 10% of putative piRNAs mapped to transposable ele- Material online). These included both DNA transposons and ments showed ping–pong signatures. A slightly higher pro- retrotransposons (supplementary fig. S4, Supplementary portion of putative piRNA reads with ping–pong signatures Material online). Interestingly, retrotransposons appeared was found to map to the genome reference outside the more frequently targeted by piRNAs than DNA transposons repeat-enriched regions (17%). Majority of these most (supplementary table S9, Supplementary Material online). 420 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Sense putative piRNA reads aligned to genome 2.0 1.0 C A U U C C C A UACC U UA A A UAC UCAUACU UAA AUA U AA U AUGCGUUCC GA UCU CGCU G A UUC GU AU CA G GGU CACGGACGAGG CG C CG G 0.0 5 10 15 20 25 30 WebLogo 3.3 2.0 Antisense putative piRNA reads aligned to genome 1.0 A A U U U GC G A G A U AUU U C C C C C 0.0 5 10 15 20 25 WebLogo 3.3 Sense putative piRNA reads aligned to TEs 2.0 1.0 GACA CCAGUUU AU U AU UU U U C AU U 0.0 5 10 15 20 25 WebLogo 3.3 Antisense putative piRNA reads aligned to TEs 2.0 1.0 AG C A G U A U UU U U C CA C G 0.0 5 10 15 20 25 WebLogo 3.3 FIG.6.—Base preferences of putative piRNA reads mapped to various data sets. Shown are base preferences of putative piRNA reads mapping to sense (A, C)and antisense (B, D) strand of genome (A, B) and transposable elements (TEs; C, D). Base preferences did not significantly differ at various pH conditions. Depicted is always one sequence set of specific length from one condition. The Y-axis represents the entropy score for the base bias. probably represent protein-coding genes. When assessing the sense and antisense orientation (figs. 7B and C). Shorter an- base-preferences of piRNAs with ping–pong signatures, we tisense reads (27 nt) mapped to transposable elements found that a substantial amount of the putative piRNAs had showed minor preference for 10 A (secondary piRNAs) 1 U preference (primary piRNAs) (fig. 7B). Sequences mapping (fig. 7B). Only around 8–9% of the putative piRNAs appeared to transposable elements showed preference for 1 U in both to be targeting transposable elements in all conditions studied Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 421 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 bits bits bits bits Urbarova et al. GBE FIG.7.—Ping–pong pathway signature. (A) Overlap probabilities of sense and antisense reads mapping to the genome and transposable elements (TEs). Overlap probabilities in all the different pH conditions are shown. (B) Depicted are base preferences of putative piRNA reads with ping–pong signatures mapping to the genome and TEs. The Y-axis represents the entropy score for the base bias. (C) Shown are putative piRNA reads aligned to a TE. Three regions with identified ping–pong signatures are highlighted. Green reads correspond to sense strand, and red reads to antisense strand. (supplementary table S7, Supplementary Material online). This genome (supplementary table S7, Supplementary Material fraction further increased nearly up to 22% when including online). These might potentially represent very divergent the unclassified fraction of the identified repeats in the transposable elements, as reported also for Exaiptasia sp. 422 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE (Baumgarten et al. 2015). Higher fraction of putative piRNAs non-symbiotic Nematostella was reported to carry about targeting transposable elements is also expected to be found four times more DNA transposons than retrotransposons during development or when extracting specifically germline (Putnam et al. 2007), which contrasts that of the symbiotic cells from A. viridis, as observed in N. vectensis (Praher et al. Exaiptasia with slightly more retrotransposons than DNA 2017). transposons (Baumgarten et al. 2015). Our data from the Ping–pong activity is mainly linked to transposon silencing symbiotic A. viridis does not appear to resemble any previously (Brennecke et al. 2007; Gunawardane et al. 2007). Therefore, sequenced cnidarian in terms of the transposable element we investigated if ping–pong activity changes could be ob- distribution, even though it contains approximately equal frac- served at different pH conditions. Transposable elements tions of DNA transposons and retrotransposons. It is interest- were only considered silenced if a significant ping–pong ac- ingtonote that the Gypsy element is the most frequent LTR tivity feature could be detected within the transposon region, retrotransposon in all the cnidarian species, including A. viridis. with FDR (q value)< 0.01. We found that a possible ping– We identified70distinct miRNAs in A. viridis, and 61 of pong dependent suppression varied among individuals in these were detected in normal seawater conditions at pH 8.2. each condition, and examples from the BEL and Gypsy LTR Only eight miRNAs were similar to previously known miRNAs retrotransposons are shown in supplementary figures S5 and in Nematostella (Grimson et al. 2008; Moran et al. 2014), six S6, Supplementary Material online. However, only a small in Acropora (Gajigan and Conaco 2017), five in Stylophora fraction of the identified transposable elements appeared si- (Liew et al. 2014)and two in Hydra (Krishna et al. 2013). lenced by the ping–pong pathway in all conditions studied These results support that taxonomically restricted miRNAs (< 1%; supplementary table S10, Supplementary Material are common to cnidarians, including A. viridis—an observa- online). tion seen mainly in plants, and which could be explained by high sequence turnover rates of miRNAs, as suggested by Moran et al. (2017). In agreement with other reports in cni- Discussion darians (Grimson et al. 2008; Wheeler et al. 2009; Liew et al. Here, we report a preliminary draft genome reference se- 2014; Moran et al. 2014, Gajigan and Conaco 2017), we quencing of the symbiotic sea anemone A. viridis,with an detected only one miRNA in A. viridis (avi-miR-temp-100) to estimated genome size of approximately 313 Mb. The par- be conserved with miRNAs in bilaterians. This miRNA belongs tially assembled reference genome was used to assess trans- to the miR-100 family, and it was found identical in sequence posable element and small RNA loci. We also performed small to nve-miR-100 and spi-miR-100 in Nematostella and RNA sequencing along a natural seawater pH gradient and Stylophora, respectively (Grimson et al. 2008; Wheeler et al. identified differentially expressed RNA candidates. In A. viridis, 2009; Liew et al. 2014; Moran et al. 2014). In bilaterians, we found 70 distinct miRNAs and thousands of putative including nematodes and humans, miR-100 makes a cluster piRNAs, suggesting that small RNAs are widespread regula- in the genome together with let-7 and miR-125, and regu- tors in the control of gene expression and transposable ele- lates transcripts involved in multiple cellular and developmen- ment silencing in this species. tal processes, as well as cancer progression (Christodoulou Theestimated genomesizeof A. viridis appears interme- et al. 2010; Sokol 2012; Li et al. 2015). The absence of diate compared with the sea anemones Exaiptasia sp. miR-51/miR-100 family was first reported in nematodes to (260 Mb) and Nematostella vectensis (329/450 Mb), slightly result in lethality during development (Shaw et al. 2010). less than the stony coral Acropora digitifera (420 Mb), and In A. viridis, as well as in other cnidarians, miR-100 appears substantially smaller than the freshwater hydroid Hydra mag- to be transcribed from an individual gene locus, but its bio- nipapillata (1.3 Gb) (Putnam et al. 2007; Chapman et al. logical role in gene repression is not well established. In the 2010; Shinzato et al. 2011; Baumgarten et al. 2015). Thus, coral Stylophora, Liew and coworkers speculated that miR- a general trend is that hexacorals harbor relatively small 100 could be involved in the calcification process (Liew genomes. We found that about 36% of the A. viridis genome et al. 2014). However, since sea anemones lack any sort of contains repeated sequences, which is a higher fraction than calcified skeleton, other processes have to be regulated by Exaiptasia and Nematostella (both 26%) and Acropora (13%), miR-100 in sea anemones. but less than Hydra (57%) (Putnam et al. 2007; Chapman We identified and described four miRNA clusters within the et al. 2010; Shinzato et al. 2011; Baumgarten et al. 2015). A. viridis genome. However, a more detailed analysis in regard There is a significant heterogeneity in the distribution of clas- to clustering of individual miRNAs was not possible due to ses and subclasses of transposable elements among the inves- insufficient contiguity of our draft assembly. Therefore, we tigated cnidarians. Whereas Hydra contains approximately cannot exclude that some miRNAs predicted in our study equal fractions of DNA transposons and retrotransposons, form additional miRNA clusters, where miRNA pairs are lo- Acropora harbors four times as many retrotransposons than cated further apart. Interestingly, one of the identified miRNA DNA transposons. There are also significant differences be- locus localized inside a DNA transposon (fig. 4), and this tween the sea anemones Nematostella and Exaiptasia.The miRNA (avi-miR-temp-58) appeared expressed mostly at Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 423 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE seawater pH 7.9. The formations of small RNAs from trans- detect any significantly differentially expressed piRNAs or posable element loci are not unusual among animals, and piRNA clusters at this point, we noted an increase in putative dozens of publications inspecting various species have piRNA expression at pH 7.9 compared with pH 7.6 and 8.2. reported miRNAs originating from transposable elements One possible biological implication could be less restricted (reviewed in Roberts et al. 2014). However, to our knowledge transposon activities at seawater pH 7.6 compared with avi-miR-temp-58 is the first example of a TE-encoded miRNA pH 7.9. reported in any cnidarian. We found nine miRNAs to be dif- ferentially expressed in A. viridis along the seawater pH gra- Conclusion dient, indicating that miRNA-based gene repression might be The A. viridis genome appears similar in size and in transpos- involved in compensating environmental stressors. Here, we able element divergence to that of Exaiptasia sp., a related identified few potential mRNA targets, including stress- sea anemone with a symbiotic lifestyle resembling that of related and mobile element proteins. Anemonia spp. The A. viridis genome encodes and expresses PIWI-interacting RNAs (piRNAs) have previously been a high number of small regulatory RNAs, and when compared reported in Nematostella (Grimson et al. 2008; Praher et al. with the sea anemone Nematostella, a large fraction (89%) of 2017). This sea anemone contains two piRNA classes, where miRNAs appears taxonomically restricted. A. viridis expresses a class I possesses an unknown function during germline devel- high amount of candidate piRNA sequences with putative opment and class II is involved in gene silencing, including functions in transposable element silencing and in other still transposons, by the ping–pong mechanism. In A. viridis,we unknown cellular functions. Some small RNAs appeared dif- found a high number of expressed piRNA candidates, appar- ferentially expressed along a seawater pH gradient, suggesting ently representing both piRNA classes, even though we did a regulatory role in the response to environmental stressors. not specifically extract and analyze germline cells in our study. However, we were not able to characterize piRNA gene loci at high resolution in A. viridis due to presence of many short Supplementary Material scaffolds in the genome assembly. Supplementary data are available at Genome Biology and A relatively high proportion of our piRNA reads showed a Evolution online. strong enrichment for uridine at 5 ends (1 U) and a higher probability to carry an adenine at the nucleotide 10 (10 A). In addition, the majority of sense and antisense putative Acknowledgments piRNA reads showed an overlap by exactly ten nucleotides. This bidirectional production of piRNA reads with 10 nt offset We thank to Sebastian Uhrig, Johannes Gutenberg University indicated a ping–pong dependent piRNA biogenesis. of Mainz, for advices using PingPongPro software tool, Inuk However, only a small fraction of the putative piRNA reads Jung, Seoul National University, for assistance in running that mapped to transposable elements showed ping–pong piClust software tool, and Professor Don Gilbert, Indiana signatures. Although we could detect nearly 22% putative University, for advices on filtering of the transposable element piRNAs potentially targeting transposons, only around 10% search output. We also thank members of the RAMP research showed ping–pong signatures in all pH conditions studied. group at UiT and the Genomics group at Nord University for This might be caused by very high divergence of transposable practical support and discussions. We also thank two anony- elements in A. viridis, an observation made previously in mous reviewers for their valuable comments and suggestions Exaiptasia sp. (Baumgarten et al. 2015). Another possible ex- that helped us to improve the manuscript. This work was planation is that piRNAs may fulfil various functions mainly supported by grants from the Research Council of Norway during development or in female adults. Here, TE-targeting (CoralSeq; to S.D.J.); and Tromsø Research Foundation (to piRNAs could be connected to the process of oogenesis and S.D.J.). The publication charges for this article have been serve the maintenance of the germline genome, as recently funded by a grant from the publication fund of UiT The reported by Praher et al. (2017). However, it indicates that Arctic University of Norway. piRNAs in cnidarians may also have additional function to that of transposable element silencing, a notion supported by Author Contributions observations in Hydra (Krishna et al. 2013; Juliano et al. 2014). More detailed characteristics of the putative piRNA I.U. and S.D.J. designed the study; I.U. collected, analyzed and population remain to be elucidated once better genome as- interpreted the data; S.F., H.P., B.O.K., and I.U. designed sembly and gene predictions are available for A. viridis. One workflows for the data analyses; H.P. wrote two Perl scripts important aim of our study was to identify and assess differ- for processing of miRNA sequencing data; I.U. and J.M.H.-S. entially expressed small RNAs along a natural seawater pH organized and performed the fieldwork; T.E.J. sequenced the gradient. We found high amounts of putative piRNA reads in small RNA libraries; all authors helped I.U. and S.D.J. prepare all the different pH conditions. Although it was difficult to the manuscript for publication. All authors (except S.F.) 424 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Elucidating the Small Regulatory RNA Repertoire of A. viridis GBE Haas BJ, et al. 2013. De novo transcript sequence reconstruction from reviewed, commented and approved the final manuscript for RNA-seq using the Trinity platform for reference generation and anal- publication. S.F. recently passed away; S.F. reviewed, com- ysis. Nat Protoc. 8(8):1494–1512. mented and approved an earlier version of the final Han BW, Wang W, Li C, Weng Z, Zamore PD. 2015. Noncoding RNA. manuscript. piRNA-guided transposon cleavage initiates Zucchini-dependent, phased piRNA production. Science 348(6236):817–821. Heo I, et al. 2012. Mono-uridylation of pre-microRNA as a key step in the Literature Cited biogenesis of group II let-7 microRNAs. Cell 151(3):521–532. Horwitz R, Borell EM, Yam R, Shemesh A, Fine M. 2015. Natural high Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. 1990. Basic local alignment search tool. J Mol Biol. 215(3):403–410. pCO increases autotrophy in Anemonia viridis (Anthozoa) as revealed Altschul SF, et al. 1997. Gapped BLAST and PSI-BLAST: a new generation from stable isotope (C, N) analysis. Sci Rep. 5:8779. Houwing S, et al. 2007. A role for Piwi and piRNAs in germ cell mainte- of protein database search programs. Nucleic Acids Res. 25(17):3389–3402. nance and transposon silencing in Zebrafish. Cell 129(1):69–82. Johnson VR, et al. 2013. Responses of marine benthic microalgae to ele- Antoniewski C. 2014. Computing siRNA and piRNA overlap signatures. vated CO . Mar Biol. 160(8):1813–1824. Methods Mol Biol. 1173:135–146. Aranda M, et al. 2016. Genomes of coral dinoflagellate symbionts high- Juliano CE, et al. 2014. PIWI proteins and PIWI-interacting RNAs function in Hydra somatic stem cells. Proc Natl Acad Sci U S A. 111(1):337–342. light evolutionary adaptations conducive to a symbiotic lifestyle. Sci Rep. 6:39734. Jung I, Park JC, Kim S. 2014. piClust: a density based piRNA clustering algorithm. Comput Biol Chem. 50:60–67. Aravin A, et al. 2006. A novel class of small RNAs bind to MILI protein in Kajitani R, et al. 2014. Efficient de novo assembly of highly heterozygous mouse testes. Nature 442(7099):203–207. Bao W, Kojima KK, Kohany O. 2015. Repbase update, a database of genomes from whole-genome shotgun short reads. Genome Res. 24(8):1384–1395. repetitive elements in eukaryotic genomes. Mob DNA. 6:11. Bartel DP. 2004. MicroRNAs: genomics, biogenesis, mechanism, and func- Kawamura Y, et al. 2008. Drosophila endogenous small RNAs bind to Argonaute 2 in somatic cells. Nature 453(7196):793–797. tion. Cell 116(2):281–297. Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E. 2007. The role of site Bartel DP. 2009. MicroRNAs: target recognition and regulatory functions. Cell 136(2):215–233. accessibility in microRNA target recognition. Nat Genet. 39(10):1278–1284. Baumgarten S, et al. 2015. The genome of Aiptasia, a sea anemone model for coral symbiosis. Proc Natl Acad Sci U S A. 112(38):11893–11898. Kozomara A, Griffiths-Jones S. 2014. miRBase: annotating high confi- Boatta F, et al. 2013. Geochemical survey of Levante Bay, Vulcano Island dence microRNAs using deep sequencing data. Nucleic Acids Res. 42(Database issue):D68–D73. (Italy), a natural laboratory for the study of ocean acidification. Mar Pollut Bull. 73(2):485–494. Krishna S, et al. 2013. Deep sequencing reveals unique small RNA reper- toire that is regulated during head regeneration in Hydra magnipapil- Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30(15):2114–2120. lata. Nucleic Acids Res. 41(1):599–616. Brennecke J, et al. 2007. Discrete small RNA-generating loci as master Langmead B, Trapnell C, Pop M, Salzberg SL. 2009. Ultrafast and memory- efficient alignment of short DNA sequences to the human genome. regulators of transposon activity in Drosophila.Cell 128(6):1089–1103. Genome Biol. 10(3):R25. Li C, et al. 2015. Multiple roles of microRNA-100 in human cancer and its Chapman JA, et al. 2010. The dynamic genome of Hydra.Nature therapeutic potential. Cell Physiol Biochem. 37(6):2143–2159. 464(7288):592–596. Christodoulou F, et al. 2010. Ancient animal microRNAs and the evolution Liew YJ, et al. 2014. Identification of microRNAs in the coral Stylophora pistillata. PLoS One 9(3):e91101. of tissue identity. Nature 463(7284):1084–1088. Das PP, et al. 2008. Piwi and piRNAs act upstream of an endogenous Lin S, et al. 2015. The Symbiodinium kawagutii genome illuminates dino- siRNA pathway to suppress Tc3 transposon mobility in the flagellate gene expression and coral symbiosis. Science 350(6261):691–694. Caenorhabditis elegans germline. Mol Cell. 31(1):79–90. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. 2012. Lowe TM, Eddy SR. 1997. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 25(5):955–964. 40(1):37–52. Min XJ, Butler G, Storms R, Tsang A. 2005. OrfPredictor: predicting protein-coding regions in EST-derived sequences. Nucleic Acids Res. Gajigan AP, Conaco C. 2017. A microRNA regulates the response of corals to thermal stress. Mol Ecol. 26(13):3472–3483. 33(Web Server issue):W677–W680. Mohn F, Handler D, Brennecke J. 2015. Noncoding RNA. piRNA-guided Ghildiyal M, Zamore PD. 2009. Small silencing RNAs: an expanding uni- verse. Nat Rev Genet. 10(2):94–108. slicing specifies transcripts for Zucchini-dependent, phased piRNA bio- Go ¨ tz S, et al. 2008. High-throughput functional annotation and data min- genesis. Science 348(6236):812–817. Moran Y, Praher D, Fredman D, Technau U. 2013. The evolution of miRNA ing with the Blast2GO suite. Nucleic Acids Res. 36(10):3420–3435. Grajales A, Rodrıguez E. 2014. Morphological revision of the genus pathway protein components in Cnidaria. Mol Biol Evol. 30(12):2541–2552. Aiptasia and the family Aiptasiidae (Cnidaria, Actiniaria, Moran Y, et al. 2014. Cnidarian microRNAs frequently regulate targets by Metridioidea). Zootaxa 3826(1):55–100. Gregory RI, Chendrimada TP, Cooch N, Shiekhattar R. 2005. Human RISC cleavage. Genome Res. 24(4):651–663. Moran Y, Agron M, Praher D, Technau U. 2017. The evolutionary origin of couples microRNA biogenesis and posttranscriptional gene silencing. Cell 123(4):631–640. plant and animal microRNAs. Nat Ecol Evol. 1(3):27. Pearson WR, Lipman DJ. 1988. Improved tools for biological sequence Grimson A, et al. 2008. Early origins and evolution of microRNAs and Piwi- comparison. Proc Natl Acad Sci U S A. 85(8):2444–2448. interacting RNAs in animals. Nature 455(7217):1193–1197. Gunawardane LS, et al. 2007. A slicer-mediated mechanism for repeat- Praher D, et al. 2017. Characterization of the piRNA pathway during de- velopment of the sea anemone Nematostella vectensis. RNA Biol. associated siRNA 5 end formation in Drosophila.Science 315(5818):1587–1590. 7:1–15. Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 425 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Urbarova et al. GBE Putnam NH, et al. 2007. Sea anemone genome reveals ancestral eume- Smit AFA, Hubley R. 2008–2015. RepeatModeler Open-1.0. Available tazoan gene repertoire and genomic organization. Science from: http://www.repeatmasker.org (last accessed January 14, 2018). 317(5834):86–94. Smit AFA, Hubley R, Green P. 2013–2015. RepeatMasker Open-4.0. Roberts JT, Cardin SE, Borchert GM. 2014. Burgeoning evidence indicates Available from: http://www.repeatmasker.org (last accessed January that microRNAs were initially formed from transposable element 14, 2018). sequences. Mob Genet Elements. 4:e29255. Sokol NS. 2012. Small temporal RNAs in animal development. Curr Opin Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor Genet Dev. 22(4):368–373. package for differential expression analysis of digital gene expression Suggett DJ, et al. 2012. Sea anemones may thrive in a high CO world. data. Bioinformatics 26(1):139–140. Glob Chang Biol. 18(10):3015–3025. Schwarz DS, et al. 2003. Asymmetry in the assembly of the RNAi enzyme Sunkar R, Chinnusamy V, Zhu J, Zhu JK. 2007. Small RNAs as big players in complex. Cell 115(2):199–208. plant abiotic stress responses and nutrient deprivation. Trends Plant Shaw WR, Armisen J, Lehrbach NJ, Miska EA. 2010. The conserved miR-51 Sci. 12(7):301–309. microRNA family is redundantly required for embryonic development and Vagin VV, et al. 2006. A distinct small RNA pathway silences selfish genetic pharynx attachment in Caenorhabditis elegans. Genetics 185(3):897–905. elements in the germline. Science 313(5785):320–324. Shinzato C, et al. 2011. Using the Acropora digitifera genome to under- Vashisht D, Nodine MD. 2014. MicroRNA functions in plant embryos. stand coral responses to environmental change. Nature Biochem Soc Trans. 42(2):352–357. 476(7360):320–323. Voinnet O. 2009. Origin, biogenesis, and activity of plant microRNAs. Cell Shoguchi E, et al. 2013. Draft assembly of the Symbiodinium minutum 136(4):669–687. nuclear genome reveals dinoflagellate gene structure. Curr Biol. Wheeler BM, et al. 2009. The deep evolution of metazoan microRNAs. 23(15):1399–1408. Evol Dev. 11(1):50–68. Simpson JT. 2014. Exploring genome characteristics and sequence quality without a reference. Bioinformatics 30(9):1228–1235. Associate editor: B. Venkatesh 426 Genome Biol. Evol. 10(2):410–426 doi:10.1093/gbe/evy003 Advance Access publication January 27, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/410/4827693 by Ed 'DeepDyve' Gillespie user on 16 March 2018

Journal

Genome Biology and EvolutionOxford University Press

Published: Feb 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial