The miRNA biogenesis in marine bivalves

Sequence retrieval and analysis

The Mg WGS project (ID APJB000000000.1(Nguyen, Hayes & Ingram, 2014)) and the Cg genome draft (GCA_000297895 (Zhang et al., 2012)) were retrieved from GenBank, whereas the oyster genome annotations were obtained from Ensembl Metazoa release 29 (http://metazoa.ensembl.org/Crassostrea_gigas/Info/Index). A Mg reference transcriptome was produced using 18,788 ESTs of mixed tissues previously obtained by Sanger sequencing (Venier et al., 2009) and 453 million reads obtained by paired-end (2 × 100 bp) Illumina Hiseq2000 sequencing of digestive gland from North Adriatic Sea mussels (ID: PRJNA88481) (Gerdol et al., 2014), and from haemocytes, gills, mantle and muscle of Spanish mussels (ID: SRP033481) (Moreira et al., 2015). The quality of the sequencing readout was evaluated by the FastQC suite (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) discarding the reads with PHRED quality below 20 and presenting more than two ambiguous nucleotides. De-novo assembly was performed with Trinity, release 2013-08-14 (Grabherr et al., 2011), setting the minimum contig length at 200 bp and using default settings. Subsequently, protein coding sequences (cds) were predicted with Transdecoder (Grabherr et al., 2011). Transcriptomic reads of 30 bivalve spp. (Cg plus other 29 species) were retrieved from the SRA archive and assembled as described above (details in File S1). The protein predictions of further 33 organisms were directly retrieved from public repositories or extracted from the corresponding genome releases. The NCBI transcriptome shotgun assembly (TSA) database was interrogated to retrieve hits for two additional cnidarians, Porites australiensis and Anthopleura elegantissima (Table 1).

Protein domain searches

To investigate the presence of eight key proteins of miRNA biogenesis (DROSHA, DGCR8, XPO5, RAN, DICER, TARBP2, AGO and PIWI), we downloaded their predictive HMM from PFAM v.27 (listed in Table 2) and we scanned the sequence datasets with HMMer v3.1 (Eddy, 2011) applying a cut-off E-value of 0.01. To achieve a meaningful comparison of proteins from different organisms, we retained only hits presenting all diagnostic domains. Moreover, we identified several mussel transcripts related to protein interactions occurring in the miRNA biogenesis. To identify such proteins, we retrieved from PFAM the diagnostic domains of human homologs (listed in Table 3) and we scanned their presence in the Mg transcriptome as described above. Protein domain organization was reconstructed using SMART (Letunic, Doerks & Bork, 2012).

Gene structure analysis

We used the transcript sequences of DROSHA, DGCR8, XPO5, DICER and TARBP2 as blast queries against all Mg genomic contigs (blastn) in order to recover the related gene structures. Positive hits having e-value lower than 10⁻²⁰ were extracted and assembled on the corresponding transcript, used as backbone. RNA-seq read mappings with adapted parameters (CLC Genomic Workbench large gap mapping tool, with similarity and length fraction set at 0.9) allowed us to ascertain the correct gene assembly. Homolog gene structures were retrieved by interrogating genomic browsers, like Metazome v.3 (for C. intestinalis, B. floridae, D. rerio, S. kowalevskii, S. purpuratus, N. vectensis, T. castaneum, L. gigantea, O. bimaculoides, C. elegans and H. sapiens) and Ensembl Metazoa v.29 (for C. gigas, C. quinquefasciatus, D. melanogaster, N. vitripennis, A. mellifera, A. queenslandica, P. bachei, M. leidyi, T. adhaerens, N. vectensis, D. pulex, S. mansoni, S. mediterranea, A. vaga, L. anatine, H. robusta and C. telata) or by local blastn against the downloaded genomes (A. pisum and L. albipes).

Phylogenetic analysis

The inferred protein sequences were aligned using MUSCLE, release 2014-05-29 (Edgar, 2004). Subsequently, the fasta alignments were analyzed using Gblocks v.0.91 (Castresana, 2000) to extract conserved positions (positions common to 51% of the locally aligned sequences). Trees were built using neighbor joining or maximus likelihood clustering methods with 1,000 bootstrap replicates. Bayesian phylogenies were reconstructed using MrBayes v.3.2.5 (Ronquist et al., 2012), with GTR substitution evolutionary model with gamma-distributed rate variation across sites, evaluating the convergence after 1,000,000 runs (0.5 was considered as cut-off value). Trees were visualized and edited with FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

Digital expression analysis

To analyze the expression of the selected genes in Cg and Mg RNA datasets, we retrieved all available RNA-seq samples from the NCBI SRA archive. For Cg, we analyzed 123 Illumina RNA-seq samples related to adult tissues or developmental stages. For Mg, we analyzed 13 RNA samples from gills (1), digestive gland (6), haemocytes (2), mantle (2) and muscle (2). Overall, we included in the expression analysis 2,271 and 453M reads for Cg and Mg, respectively (File S2). The trimmed reads were mapped to Cg and Mg genes using the CLC Genomics Workbench v.8.0 (Qiagen, Hilden, Germany) mapping tool, with length and similarity fractions set at 0.75 and 0.95, respectively, and mismatch/insertion/deletion penalties at 3/3/3. The number of uniquely mapped reads of each dataset were counted and used to calculate digital expression values as TPM (Transcripts Per Kilobase Million mapped reads) as described by (Wagner, Kin & Lynch, 2013), considering 3 TPMs as lower detection limit.

Mussel transcripts related to the miRNA biogenesis

We identified Mytilus galloprovincialis transcripts involved in the miRNA biogenesis by systematic searches of diagnostic domains (Table 2) in a transcriptome assembly produced from 453 million Illumina reads. Thus, we recovered nine transcripts coding for DROSHA, DGCR8, XPO5, RAN, DICER, TARBP2 and for three Argonaute genes (one Ago and two Piwi-like proteins, Table 3). We also identified 21 mussel proteins expected to play a role in the miRNA maturation or involved in RNAi processes (File S3).

Figure 1 relates the general process of eukaryotic miRNA biogenesis to the mussel proteins identified in this work. MgDROSHA and MgDGCR8 are expected to start the maturation of pri-miRNAs produced by RNA polymerase II. MgDROSHA codes for a 1,377 aa length protein containing all the canonical domains (2 RIBOc domains in positions 959–1,093 and 1,139–1,271 and one DSRM domain in position 1,278–1,351) whereas MgDGCR8 is a 728 aa length protein having one WW domain in position 229–258, necessary for the interaction with DROSHA, and two DSRM domains (positions 472–536 and 578–642) necessary for pri-miRNA binding. MgXPO5 is expected to cooperate with MgRAN in the pre-miRNA cytoplasmic translocation. MgRAN encodes a 214 aa protein whereas MgXPO5 has a length of 1,201 aa and includes two 5’ conserved domains (IBN_N and Xpo1) and one conserved region necessary for the interaction with interleukin enhancer-binding factor 3 (position 525–562). In mussels, the RISC complex uploading pre-miRNAs is defined by the endoribonuclease MgDICER (1,850 aa) and MgTARBP2 (321 aa). Like in Lophotrocozoa, mussel DICER is encoded by a unique gene and contains the seven canonical domains, namely two helicase domains, one DICER-dimer domain, one PAZ, two RIBOc and a final DSRM domain. MgTARBP2 displays three DSRM domains in positions 9–73, 101–166 and 249–314. Moreover, M. galloprovincialis possess three argonaute proteins ranging from 861 to 941 aa in length and representative of one AGO (DUF1785, PAZ and PIWI domains) and two PIWI-like proteins (PAZ and PIWI domains). We considered the above mentioned gene products as the key complement of the miRNA biogenesis.

Among the possible interacting proteins, we identified MgGW182, a transcript encoding a protein shorter than the human counterparts but holding all the features considered significant for its interaction with AGOs and the CCR4-NOT complex. In fact, MgGW182 possesses 19 N-terminal GW stretches, followed by one UBA domain, a Q-rich region (M domain) and a C-terminal RNA recognition motif (RRM domain). Moreover, we recognized a C-terminal conserved site known as PAM2 (Kozlov et al., 2010), expected to interact with the poly(A) binding protein 1 (MgPABPC1) through the MLLE motif and inhibit the mRNA translation by interfering with the mRNA circularization process (Piao et al., 2010; Van Kouwenhove, Kedde & Agami, 2011). In the mussel transcriptome, we also found putative homologs for a number of CNOT complex proteins (CNOT1, 2, 3, 6, 7, 9, and 10), for the eukaryotic translation initiation factor 4 gamma, 1 eIF4G, PAB-dependent poly(A)-specific ribonuclease subunits PAN2, PAN3, the decapping complex proteins DCP1 and DCP2, and several RNA helicases demonstrated to be crucial in the miRNA maturation (DDX5) and RNAi (DDX5- 6- 20 and 42). Finally, we recognized the putative mussel homologs of protein arginine methyltransferase 5 (MgPRMT5), tudor domain containing protein (MgTDRD-11) and maelstrom spermatogenic transposon silencer (MgMAEL).

Mussel genes related to the miRNA biogenesis

Taking advantage of mussel WGS data (Nguyen, Hayes & Ingram, 2014) we investigated the organization of the main genes involved in the mussel miRNA biogenesis. Fragmentation of the genomic mussel assembly (2.3 million contigs; 700 bp on average) and considerable dimension of the analyzed genes (9.6–17.6 kbp gene size in the case of Cg) prevented the recovery of the full gene sequences. Nevertheless, we can describe the complete gene structures of DROSHA, DGCR8, EXP5, DICER and TARBP2 (i.e., five of eight searched sequences) whose length varies between 7.5 and 27 kbp, confirmed by the back-mapping of 115,377 Illumina paired reads (Fig. 2, File S4). Moreover, these mussel genes showed a remarkable conservation in terms of exon number when compared with a selection of homolog genes from deuterostome and protostome organisms (Table 4).

Transcripts related to the miRNA biogenesis in bivalve spp

To identify the miRNA biogenesis complements in marine mollusks, we used homologous genes retrieved from the genomes of C. gigas, L. gigantea and A. californica. Since the C. gigas genome includes annotations only for the cds regions, we exploited full-length transcripts obtained from a locally assembled oyster transcriptome to expand the genome annotations in this species. In particular, we updated the annotation of CgDGCR8 and CgDICER and we added new annotations for CgPIWI-1 (CGI_10008757: genomic contig JH815696, position 184178–187825) and CgTARBP2 (JH818440, 414703–419857).

Since many marine bivalve spp. do not have at present a sequenced genome, we used publicly available RNA-seq data to build 29 specie-specific transcriptome assemblies and retrieve the homologous sequences of interest. After domain searching, we carefully considered the high number of positive hits to retain only proteins including all the expected protein features. Thus, we retrieved 132 complete hits from marine mollusks: 10 DROSHAs, 9 DGCR8s, 14 XPO5s, 34 RANs, 7 DICERs, 13 TARBP2s and 45 Argonaute-like proteins, the latter classified in 13 AGO and 32 PIWI proteins by phylogenetic analysis (Table 1, File S5).

Phylogenetic analysis of the miRNA biogenesis proteins

The inferred sequences of single miRNA biogenesis proteins were aligned together with those retrieved from 34 sequenced genomes. Here, we report the phylogenetic analysis of the five proteins centrally involved in the miRNA biogenesis, namely DROSHA, DGCR8, DICER, TARBP2 and AGOs (File S5 includes all protein sequences). We back-traced the presence of a canonical DROSHA up to Cnidaria, although we found only incomplete hits in Porifera and Placozoa and the genomes of Ctenophora spp. lack of both DROSHA and DGCR8, as reported by other authors (Maxwell et al., 2012). The DROSHA sequences from Cnidaria’s appeared as general outgroup whereas those of Chordata clustered as outgroup of the other protostomes. DROSHAs from Mollusca and Arthropoda clustered consistently with the different taxa whereas those from Platyhelmintes, Rotifera, Brachiopoda and Annelida grouped together, with DROSHA from Caenorhabditis elegans (Nematoda) being the most far-related (Fig. 3A). Contrary to DROSHA, we identified a complete DGCR8 also in the Porifera Amphimedon queenslandica, suggesting that also DROSHA should be present in this taxa. Following phylogenetic analysis, we highlighted Cnidaria and Porifera proteins as outgroup, with mollusks (and Annelida) clustering with Arthropoda and more distantly Platyhelmintes and Rotifera hits. The Chordata sequences clustered as a separate group (Fig. 3B).

The finding of putative DICER sequences in Ctenophora spp. supports the presence of this gene through the whole Opisthokonta evolution (Maxwell et al., 2012). Also plants possess DICER homologues which occur in different copy number among taxa: two genes in Porifera, Placozoa, Cnidaria, Platyhelminthes and Arthropoda (with the exception of D. pulex that possess three genes); four genes in plants like A. thaliana and P. trichocarpa and one gene in Ctenophora, Rotifera, Cephalopoda Mollusca and Chordata. Moreover, the presence of DICER was reported in some Protozoa and fungi (Mukherjee, Campos & Kolaczkowski, 2013). Phylogenetic analyses, separate insect DICER-2,plant DICERs from DICER-1. DICER-1 clade shows a consistent clustering of Arthropoda, Mollusca and Chordata hits, whereas some branches of basal metazoans and Platyhelminthes are not well resolved (Fig. 3C). Likewise, the phylogenetic tree regarding TARBP2 displays a clear cut-off between the proteins of mollusks, chordates and arthropods (Fig. 3D). We back-traced the miRNA cytoplasm export complex composed by RAN and XPO5 in all analyzed metazoans. Both RAN and XPO5 represent widely expressed sequences since we found them also in transcriptome assemblies, although with suboptimal sequence coverage.

Several AGO and Piwi proteins can be present in individual organisms and, in fact, we identified a total of 235 proteins. Whereas humans possess eight proteins, we found four proteins in the majority of the analyzed insect spp. (with the exception of 15 proteins in A. pisum) and three or four different proteins in bivalve spp. Also, basal Metazoa possess Argonaute-like sequences: four in the genomes of Ctenophora and Cnidaria spp., one in the Placozoa T. adhaerens and two in A. queenslandica. The case of C. elegans is remarkable since it holds several Argonaute gene families and at least 24 proteins (Hoogstrate et al., 2014). In agreement with other phylogenetic studies (Swarts et al., 2014), the Argonaute proteins from plants and the majority of those from C. elegans formed distinct clades and, moreover, a clear separation was evident between AGO and PIWI proteins. Bivalve protein sequences clustered always separately forming one cluster for AGO-like hits and two clusters for PIWI-like proteins (Fig. 4).

Digital expression analysis of mussel and oyster miRNA biogenesis genes

We used the 13 Mg and 124 Cg RNA-seq samples to evaluate the expression levels of miRNA biogenesis genes in different tissues and conditions. Based on total mapped reads, we computed TPM values and we used elongation factor 1 α (El1α) as normalizer housekeeping gene to compare the expression level of the different genes in each sample.

For Mg, the sequence analysis indicated a scarce basal expression of the genes mentioned above in five adult tissues: gill, digestive gland, haemolymph, muscle and mantle (below 2% of El1α, except for DDX5, RAN and CNOT9). Mantle and muscle appeared the most responsive tissues whereas haemolymph was the least responsive one. In particular, the genes that we considered as the core components of miRNA biogenesis were expressed at levels below 0.5% of El1α (File S6).

For Cg, we analyzed a considerable number of RNA-seq libraries representative of adult tissues (85) and developmental stages (39) (File S6). In adult oysters we observed low basal expression, as detected in the mussel samples. In fact, none of the experimental conditions reported for the analyzed RNA-seq samples influenced substantially the expression of the core miRNA pathway genes (expression levels below 2% of El1α), with the exception of the high levels of CgPIWI-1 levels in male and female gonads (around 3.5%, Fig. 5). Conversely, most of the miRNA biogenesis genes were expressed at remarkable levels during the early stages of the oyster development: mainly from two cells to the rotary movement and, for some genes, also in the next developmental stages until D-shaped larvae, with no detectable signals afterward in spat and juveniles. Hence, these genes are particularly active in the early development, in particular one AGO (CGI_10020511) and two PIWI transcripts from the egg to trocophora (Fig. 5). In the same developmental stages we also noticed a remarkable expression of the key miRNA genes, with the co-expression of DROSHA and DGCR8 evident in all the analyzed samples.