Inventory of European Sea Bass (Dicentrarchus labrax) sncRNAs Vital During Early Teleost Development

During early animal ontogenesis, a plethora of small non-coding RNAs (sncRNAs) are greatly expressed and have been shown to be involved in several regulatory pathways vital to proper development. The rapid advancements in sequencing and computing methodologies in the last decade have paved the way for the production of sequencing data in a broad range of organisms, including teleost species. Consequently, this has led to the discovery of sncRNAs as well as the potentially novel roles of sncRNA in gene regulation. Among the several classes of sncRNAs, microRNAs (miRNAs) have, in particular, been shown to play a key role in development. The present work aims to identify the miRNAs that play important roles during early European sea bass (Dicentrarchus labrax) development. The European sea bass is a species of high commercial impact in European and especially Mediterranean aquaculture. This study reports, for the first time, the identification and characterization of small RNAs that play a part in the 10 developmental stages (from morula to all fins) of the European sea bass. From 10 developmental stages, more than 135 million reads, generated by next-generation sequencing, were retrieved from publicly available databases as well as newly generated. The analysis resulted in about 2,000 sample grouped reads, and their subsequently annotation revealed that the majority of transcripts belonged to the class of miRNAs followed by small nuclear RNAs and small nucleolar RNAs. The analysis of small RNA expression among the developmental stages under study revealed that miRNAs are active throughout development, with the main activity occurring after the earlier stages (morula and 50% epiboly) and at the later stages (first feeding, flexion, and all fins). Furthermore, investigating miRNAs exclusively expressed in one of the stages unraveled five miRNAs with a higher abundance only in the morula stage (miR-155, miR-430a, d1, d2, and miR-458), indicating possible important key roles of those miRNAs in further embryonic development. An additional target search showed putative miRNA-mRNA interactions with possible direct and indirect regulatory functions of the identified miRNAs.

Introduction

Non-coding RNAs (ncRNAs) are involved in several different regulatory pathways. In the last decade, their importance has been demonstrated in a broad range of organisms, including teleost species (Herkenhoff et al., 2018). Several types of non-coding regulatory RNAs have been identified, chief among them being long ncRNAs and small ncRNAs (sncRNAs). The major classes of sncRNAs are short interfering RNAs, piwi acting RNAs, small nuclear RNAs (snRNA), small nucleolar RNAs (snoRNA), and microRNAs (miRNAs; Labbé et al., 2017). The majority of studies investigated the role of miRNA, which is a highly conserved class of small regulatory RNAs known to act at the translational level mainly by repressing protein production (Gavery and Roberts, 2017). The rapid advancements in sequencing and computing methodologies have significantly enhanced the discovery of new miRNAs not only in humans and model species, such as the zebrafish (Danio rerio; Mishima, 2012), but also in non-model species, such as the Atlantic halibut (Hippoglossus hippoglossus) (Bizuayehu et al., 2012a) and the catfish (Ictalurus punctatus) (Barozai, 2012). Besides their functional importance in diverse biological processes, it has also been shown that miRNAs are essential for vertebrate development (Wienholds and Plasterk, 2005). In teleosts, the importance of miRNA in modulating gene expression during development has been reported in various studies evaluating both model and non-model fish species. The first miRNAs with a regulative function in fish development have been described in zebrafish (Kloosterman et al., 2004). The dynamics of miRNA expression during early ontogeny in non-model fish species, such as the Atlantic halibut (Bizuayehu et al., 2012b), turbot (Scophthalmus maximus); (Robledo et al., 2017), and Senegalese sole (Solea senegalensis; Campos et al., 2014), have also been reported.

The European sea bass (Dicentrarchus labrax) is a species of high commercial impact whose industrial production is steadily growing. Efforts have also been made to intensify fish cultivation and target faster growth rates and better feed conversion ratios. Consequently, several objectives such as the larval survival rate, alternative feeds, and disease resistance are of significant importance to the industry. The embryonic and larval stages are part of the most important periods to ensure high performance and superior quality in the following developmental phases of the life cycle, particularly for fish in captivity (Varsamos et al., 2006; Pittman et al., 2013; Tsalafouta et al., 2015). Over the last decade, these economic interests have directed increased research efforts in the rearing of European sea bass, including the significant enrichment of the molecular toolbox for the study of this species. Today, besides several transcriptome datasets (Sarropoulou et al., 2009; Sarropoulou et al., 2012; Pinto et al., 2017), the whole genome sequence of the European sea bass (Kuhl et al., 2010; Tine et al., 2014) as well as single nucleotide polymorphism markers, genetic linkage maps (Chistiakov et al., 2008; Guyon et al., 2010), and radiation hybrid maps (Guyon et al., 2010) are available. However, scarce information has been published concerning ncRNAs (Kaitetzidou et al., 2015).

The present work aims to identify sncRNAs and their targets that play an important role in the development of the European sea bass. Teleost development can be reflected as a sequence of ongoing morphological changes, whereby the embryonic and larval stages are considered to be the most significant time points in the life cycle of marine fish. In the natural environment and during its embryonic and larval phases, the European sea bass lives in the marine environment, whereas as a juvenile it migrates to coastal zones, estuaries, and lagoons. The European sea bass is therefore considered a euryhaline species with reportedly high adaptation processes during its early life phases. To detect most of the sncRNAs that are important to the European sea bass development and to obtain a list of unique miRNAs, this study analyzed the small RNA transcripts of early development (as described by Kaitetzidou et al., 2015) along with newly generated sequencing reads from three additional later stages. It reported for the first time the identification and characterization of the European sea bass small RNAs during development (from morula to all fins); using target search, it also revealed the putative miRNA–mRNA interactions resulting in the possible direct and indirect regulatory functions of the selected miRNAs, which show differential expression (DE) during development.

Materials and Methods

All experiments/methods in the present study were performed in accordance with the approved guidelines and regulations from the Hellenic Center for Marine Research (HCMR) Institutional Animal Care and Use Committee following the three Rs (Replacement, Reduction, and Refinement) guiding principles for more ethical use of animals in testing, which was first described by Russell and Burch in 1959 (EU Directive 2010/63). These principles are now followed in many testing establishments worldwide before the initiation of experiments.

Sampling and RNA Extraction

For the generation of new data from three later developmental stages, samples were collected from i) the first feeding stage (FF), when the mouth opens, the pectoral fins have been developed, and the yolk sac has been completely absorbed; ii) the flexion (FLX) stage, when the notochord flexion has been completed; and iii) the all-fins stage (FINS), when all fins have been established (Pavlidis et al., 2011). All samplings were carried out at the installations of the Institute of Marine Biology, Biotechnology, and Aquaculture (IMBBC), HCMR (Heraklion, Crete, Greece). For each of the newly sampled stages, i.e., FF, FLX, and FINS, three biological replicates were sampled, flash frozen in liquid nitrogen, and transferred to a −80°C ultra-low freezer until miRNA library preparation. For miRNA library construction and sequencing, the same method as described by Kaitetzidou et al. (2015) was followed. In brief, total RNA was extracted from all developmental stages using Nucleospin miRNA Kit for the isolation of small RNA (sncRNA) and mRNA (Macherey-Nagel GmbH & Co. KG, Duren, Germany) according to the manufacturer’s instructions. Larvae were disrupted with mortar and pestle in liquid nitrogen and homogenized in lysis buffer by passing lysate through a 23-gauge (0.64 mm) needle five times. The quantity of RNA was estimated with NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Inc., Wilmington, DE, USA) and the quality was further evaluated by agarose (1%) gel electrophoresis and Agilent 2100 Bioanalyzer using RNA Nano Bioanalysis chip.

miRNA Library Preparation and Sequencing

For all samples, miRNA libraries were prepared according to the manufacturer’s instructions and single-end sequencing was carried out using Illumina sequencing technology platform (Illumina HiSeq 2000) at the Cornell University Core Laboratories Center. The use of multiplex identifier tags for each library allowed the pooling of the samples to be run in one HiSeq 2000 lane. The quality control of all reads was assessed by running Fastqc version 0.10.0 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc).

Sequence Reads Analysis

Raw sequencing data from seven early developmental stages, i.e., morula (M), 50% epiboly (E), late gastrulation-organogenesis (GO), organogenesis (O), late organogenesis (LO), hatching (HA), and 24 h post-hatching (hph), were retrieved from Kaitetzidou et al. (2015). Sequence reads from a total of 10 developmental stages, i.e., reads obtained from Kaitetzidou et al. as well as the newly generated reads from the FF, FLX, and FINS stages, were quality and adapter trimmed using Trimmomatic software 0.30 (Bolger et al., 2014) and imported into CLC Genomics Workbench version 10.1 (CLC Bio, Aarhus, Denmark). The annotation of unique sncRNAs was performed against Gasterosteus_aculeatus.BROADS1.ncrna as well as against available miRNA annotations of teleosts, humans, and mice in miRBase (Griffiths-Jones et al., 2008; Cyprinus carpio, D. rerio, Fugu rubripes, I. punctatus, Oryzias latipes, Petromyzon marinus, Salmo salar, Tetraodon nigroviridis, Gallus gallus, H. hippoglossus, Homo sapiens, Mus musculus, and Paralichthys olivaceus).

Differential Expression

For DE analysis, transcript counting was carried out with CLC Genomics Workbench 10.1.1 (CLC Bio; https://www.qiagenbioinformatics.com/). For all not yet normalized transcripts, read counts greater than five were considered for DE analysis. For DE sequencing reads, which could not be annotated through the Gasterosteus_aculeatus.BROADS1.ncrna database, additional effort was made by submitting them to publicly available databases of the National Center for Biotechnology Information (NCBI; non-redundant nucleotide and expressed sequence tags). Transcripts per million normalized count values were transformed into a decimal logarithmic scale. Obtained small RNA reads were mapped to the D. labrax genome for validation purposes. DE was assessed based on the empirical analysis method (EDGE test) provided within the CLC Genomics Workbench with default parameters. Transcripts were considered significant DE if false discovery rate (FDR) values were below 0.05 and log₂fold change > 2. To find stage-specific miRNAs, transcripts found only in one of the stages but zero counts in the other were considered as stage-specific miRNA.

miRNA Target Search

For target search, the coding sequence information (mRNA) for the European sea bass was extracted from the NCBI nucleotide databases and open reading frame prediction was carried out using Transdecoder version 5.5.0 software program (Haas et al., 2013). 3′ Untranslated regions were extracted by ExUTR version 0.1.0 and RNAhybrid version 2.1.2 (Krüger and Rehmsmeier, 2006) as well as miRanda web tool (Enright et al., 2003) and used for target prediction with the default parameters and an energy threshold of mfe < −30 kcal. The complete workflow is illustrated in Figure 1.

FIGURE 1

Figure 1 Workflow overview comprising all steps and methods used in the present work. Annotation with miRBase was carried out by applying the databases of C. carpio, D. rerio, F. rubripes, I. punctatus, O. latipes, P. marinus, S. salar, T. nigroviridis, G. gallus, H. hippoglossus, H. sapiens, M. musculus, and P. olivaceus.

Gene Ontology (GO) Classification and Enrichment Analysis of Target Genes

GO terms were retrieved using the PANTHER classification system (http://www.pantherdb.org/). The chosen analysis type was “Statistical overrepresentation test” and the annotation data set the PANTHER GO-Slim Biological Process data set (PANTHER version 14.1 Released 2019-03-12). The analyzed list consisted of the zebrafish target genes for miR-430 and miR-21 retrieved from TargetScanFish (TargetScanFish, release 6.2, June 2012; Ulitsky et al., 2012). The reference gene list for the test consisted of all genes in the PANTHER database for zebrafish. The test type was Fisher’s exact test with FDR correction. The statistically overrepresented GO terms of miR-430 and miR-21 target genes were illustrated using PANTHER Overlaid Area Chart of Difference (observed vs. expected).

Results

Small RNA Libraries from Developmental Stages of the European Sea Bass

For samples FF, FLX, and FINS, new miRNA libraries were generated. RNA quality control by DNAnalyzer Nano RNA chip before library construction showed high-quality total RNA with an RIN number > 8 (Supplemental Figure S1). Only samples passing this evaluation step were used for next-generation sequencing (NGS). For each sample after trimming, an average of about 7 million reads were obtained and about 40% were annotated. Annotated read counts for each stage were grouped (“sample grouped”) and resulted in an average of 1,380 annotated reads per stage. An overview of sequencing read numbers is listed in Table 1. Sample grouped reads of all stages together amounted to 2.115 reads, whereas sequence lengths of 20 to 24 bp were more than 50% (Figure 2A). The total number of successfully mapped small RNA (sample grouped) onto the European sea bass genome amounted to 1,524 reads (∼72%; Supplemental Table S1) with 1,169 (77%) reads corresponding to miRNA (Figure 2B). Raw sequencing reads were submitted to the SRA database of the NCBI under accession numbers PRJNA269278 and PRJNA369460.

TABLE 1

Table 1 Summary of sequencing results.

FIGURE 2

Figure 2 (A) Read length distribution of all sequencing reads after trimming and sample grouping. Y-axes show the percentage of read counts at a specific read length, whereas X-axes show read lengths. (B) Successfully mapped small RNA on the available European sea bass (D. labrax) genome. The outer ring illustrates the chromosome number, whereas the inner rings present the relative amount of mapped miRNAs.

Differentially Expressed miRNAs During Sea Bass Development

Pair-wise DE analysis was carried out and transcripts with FDR values < 0.05 and log₂fold change > 2 were considered as DE, resulting in a total of 1,157 annotated reads in at least one comparison among all the 10 stages under study (Supplemental Table S2). The majority of transcripts were annotated as miRNA (59%) followed by snRNA (25%) and snoRNA (10%; Figure 3). Principal component analysis (PCA) showed that earlier stages separated themselves clearly from the later stages: the FLX and FINS stages (Figure 4A). A similar pattern was obtained by hierarchical clustering (Figure 4B). According to the PCA plot and the heatmap, the M stage was selected as the reference stage to generate a Venn diagram (Figure 5). Furthermore, the comparisons of stages GO, O, LO, HA, and 24 hph to stage M were combined in one list as well as the comparison of stage M to stages FLX and FINS. The Venn diagram showed that the majority of transcripts exclusively expressed belonged to the FLX and FINS stages (263, 26%), with 62% of them classified as snRNA and 24% as miRNA. The list comprising stages GO, O, LO, HA, and 24 hph had 40 (3.9%) unique transcripts (36 miRNAs and four snoRNAs), the FF had 13 (1.3%), and the E stage included only 4 (0.4%) miRNA.

FIGURE 3

Figure 3 Distribution of obtained small RNA types differentially expressed within the whole sample set. misc_RNA, miscellaneous RNA; Mt_rRNA, mitochondrial rRNA; rRNA gene_Symbol, rRNA gene.

FIGURE 4

Figure 4 (A) Principal component analysis (PCA) plot of normalized mean values of all stages. (B) Heatmap of normalized mean values of all significant DE transcripts.

FIGURE 5

Figure 5 Venn diagram of differentially expressed reads comparing four lists: (1) morula (M) compared to epiboly (E); (2) combined list of the comparisons of stages gastrulation-organogenesis (GO), organogenesis (O), late organogenesis (LO), hatching (HA), and 24 h post-hatching (24 hph) to M; (3) comparison of M to first feeding (FF) stage; and (4) combined list of the comparisons of M to stages flexion (FLX) and all-fins (FINS).

Stage-Specific miRNAs

Transcripts found only in one of the stages and zero counts in the other were considered as stage-specific miRNA. In the present work, no stage-specific miRNAs were detected for stages E, FLX, FINS, and HA. For the first stage, the M stage, five miRNAs were found, among which three belonged to the miR-430 family. The last studied stage, the FINS stage, also did not show the presence of any uniquely expressed miRNA. Nevertheless, one of the highly abundant small RNAs at this stage was miR-462. Among the few miRNAs regulated between FLX and FINS stages (a total of 9 miRNAs) were hs-miR-7641-1 and hs-miR-7641-2, which were found to be more abundant at the FLX stage.

Target Search

miRNA target search in European sea bass, applying positive hybridization scores with an energy threshold < −30.0, was carried out at the gene level for the 14 stage-specific miRNAs (transcript found only in one stage and zero counts in the other stages). The putative targets for the European sea bass, as identified by RNAhybrid, are listed in Supplemental Table S3. miRNA target search in zebrafish with miR-430 (stage-specific miRNA) and miR-21 (up-regulated in FF) resulted in specific GO terms illustrated in Figure 6 and listed in Supplemental Table S4.

FIGURE 6

Figure 6 Biological process GO terms of miR-430 and miR-21 target genes statistical overrepresented illustrated using PANTHER Overlaid Area Chart of Difference (observed vs. expected). The numbers of the six most represented GO terms of the miR-430 target genes correspond to (1) cellular response to endogenous stimulus, (2) cellular response to growth factor stimulus, (3) establishment or maintenance of cell polarity, (4) protein targeting, (5) response to endogenous stimulus, and (6) response to growth factor.

Discussion

Deep Sequencing of European Sea Bass Developmental Stages

Today, NGS has made it possible to assess not only gene expression but also the abundance of small RNA within a given sample. Here, we analyzed more than 135 million reads obtained from 10 developmental stages of the European sea bass comprising NGS data from a previously published work (Kaitetzidou et al., 2015) as well as from the present study. Read numbers varied across the samples, but the numbers of annotated reads, which were “sample grouped,” were alike (Table 1), pinpointing the fact that sufficient reads for each sample were obtained. The majority (>50%) of the trimmed and sample grouped reads were 20 to 24 bp long (Figure 2A). Including the transcripts read lengths of 15 to 19 bp as well as 25 to 29 bp resulted in more than 70% of the sample grouped reads, demonstrating that the libraries were significantly enriched with miRNAs.

Small RNA Analysis

Two main approaches exist for analyzing the small RNA reads obtained by NGS. Either the different types of small RNAs in the data are counted or the read counts are first mapped to an appropriately annotated reference genome and then submitted to publicly available databases for annotation purposes. The first approach does not require an annotated genome for mapping, and small RNAs not mapping to the reference genome (due to gaps for instance) can still be measured. In the present study, the first approach was implemented to detect all possible small RNAs obtained by NGS. Nevertheless, with the aim of locating the obtained small RNAs within the genome, reads were mapped to the European sea bass genome. Notably, the majority of miRNAs mapped to four linkage groups (LG5, LG14, LG15, and LG25; Figure 2B). Other studies, in about 20 species, have also investigated the detailed chromosome-specific location of miRNA genes. The authors of these studies showed that, in all species, certain chromosomes accumulated a higher number of miRNAs and that miRNAs involved in specific diseases also accumulated in a specific location (e.g., Ghorai and Ghosh, 2014). The aforementioned study by Ghorai and Gosh (Ghorai and Ghosh, 2014) also included three teleost species: the zebrafish, medaka (O. latipes), and tetraodon (T. nigroviridis). Although miRNAs clearly accumulated in specific chromosomes in zebrafish and medaka, in tetraodon, no clear preference was shown. Nevertheless, the authors listed the top four highest numbers of miRNA genes containing chromosomes, which were chromosomes 7, 4, 14, and 17 for medaka and chromosomes 2, 10, 9, and 13 for tetraodon (chromosome numbers ordered by the highest miRNA counts found). Medaka chromosome 7 is the homolog group to tetraodon 9 (Sarropoulou et al., 2008) and corresponds to the European sea bass LG25, which also showed an abundance of miRNA in the present study (Figure 2B). The European sea bass LG14 corresponds to medaka chromosome 14 (third highest number of miRNA in medaka) and LG15 to medaka chromosome 21 (fifth highest number of miRNA in medaka) and tetraodon 2 (highest number of miRNA in tetraodon). Of the four LGs with enriched miRNA counts (LG5, LG14, LG15 and LG25), only European sea bass LG5 did not match with any of the chromosomes with enriched miRNA counts in medaka and tetraodon. In humans, it has been shown that clustered miRNAs are involved in specific pathways or, in particular, cell functions and that clustered miRNAs may be maintained through evolution (Guo et al., 2014). In teleosts, this aspect has not yet been investigated and may be the objective of future studies.

Differentially Expressed miRNAs During Sea Bass Development

DE miRNAs showed that the investigated earlier stages (i.e., M and E) and the latest stages (i.e., FLX and FINS) separated themselves clearly from the other stages (Figure 4A). Therefore, a first comparison to the M stage was carried out, grouping together stages GO, O, LA, HA, and 24 hph as well as stages FLX and FINS (Figure 5).

Between Morula and 50% Epiboly

A total of 220 transcripts were found to be DE between the two earliest stages studied within the present work. Of the 220 transcripts, 165 (75%) were annotated as miRNAs, and 70 miRNAs were found in higher abundance in the E stage. The other 95 miRNAs were detected to be more abundant in the M stage. Four transcripts were exclusively DE in the M vs. E comparison (Figure 5). All of them were found in higher abundance in the M stage and appeared to belong to the miR-30a family (miR-30a-1, miR-30a-2, miR-30a-3, and miR-30a-4). For comparison, in Xenopus, studies have shown that pri-miRNAs of miR-30a are active only during maternal stages and that no zygotic transcription could be detected within the studied stages (Nepal et al., 2015). Higher miR-30a abundance during the blastula period in Xenopus was also shown to be linked to early neural crest development (Ward et al., 2018). In contrast, the highest miRNA fold change between the M and E stages was found for miR-196b with 0 and 1.111 copy numbers, respectively. In zebrafish, research has shown that miR-196 is involved in the regulation of axial patterning and pectoral appendage initiation and that the first appearance of miR-196b has been seen to be later than it was detected in the present study, i.e., at the bud stage (He et al., 2011).

Between Morula, 50% Epiboly and Gastrulation-Organogenesis, Organogenesis, Late Organogenesis, Hatching, and 24 Hours Post-Hatching

According to the PCA plot, the stages from GO to 24 hph are grouped closely together in comparison to the other stages. Therefore, in the Venn diagram, DE transcripts between those stages and the M stage were considered as one list (Figure 5). Here, a total of 40 transcripts, comprising 36 miRNAs and 4 snoRNAs, were found to be exclusively DE. The fact that 90% of the reads exclusively found in the M-GO comparison to 24 hph stages are annotated as miRNA may indicate that miRNAs are more active during early development. Similar findings were described in zebrafish, where an increase of miRNA expression was observed as early development proceeded (Yao et al., 2014).

Between Morula and First Feeding Stage

Likewise, 13 transcripts (6 miRNAs, 2 snoRNAs, and 5 snRNAs) were found to be exclusively present in the FF stage. Five of the six identified miRNAs, namely miR-4792, let-7b-1, let-7b-2, miR-155-1, and miR-155-2, showed higher abundance in the M stage than in the FF stage (Figure 5). Whereas let-7 is a well-studied miRNA belonging to the first founding miRNAs with possible roles in growth development (Kloosterman et al., 2004; Zhao et al., 2017), both miR-155 and miR-4792 have been examined less. Of the remaining miRNAs, miR-21 (ENSGACT00000028044.1) showed higher expression levels in the FF stage. It is also among the first miRNAs to have been identified (Kim, 2005) and has been found to be involved in various biological processes, including development (Kumarswamy et al., 2011). In zebrafish, miR-21 expression was found at very early developmental stages (Chen et al., 2005), whereas, in rainbow trout, miR-21 was suggested to play an important role in degrading maternally inherited mRNAs (Ramachandra et al., 2008). Enrichment analysis of target genes in zebrafish, identified via TargetScanFish (Ulitsky et al., 2012), revealed a high number of genes classified in the biological process termed as “tube development” (GO:0035295; Figure 6). The formation of tubes such as epithelial and endothelial tubes are of importance in view of gases, liquids, and cell transport. In tilapia (Oreochromis niloticus), it has been shown that miR-21 is involved in the modulation of alkalinity stress (Zhao et al., 2016). During European sea bass development, stress, as indicated by water cortisol measurement, was first detected at the FF stage (Tsalafouta et al., 2015). The present finding may indicate that miR-21 in the European sea bass is also involved in mechanisms related to the modulation of stress.

Between Morula and Flexion/All Fins

The comparison between M and FLX/FINS stages revealed the highest number (263) of exclusively found transcripts. Notably, of the 263 transcripts, only 62 (24%) were classified as miRNA, whereas 164 (62%) reads were annotated as snRNA; all of them were found in higher copy numbers in the M stage. snRNA molecules are known to be an abundant component of eukaryotic cells (Valadkhan and Gunawardane, 2013), their importance for development having already been demonstrated in Xenopus nearly four decades ago (Forbes et al., 1983). The exclusively high abundance of snRNA, as illustrated in Figure 5, may indicate the importance of studying snRNAs during development in the near future; however, it is beyond the scope of the present work.

Stage-Specific miRNAs and Their Targets

Annotated miRNAs detected only in one stage and not in the others were found for the M, GO, O, LO, 24 hph, and FF stages. Unique miRNAs found in only one stage and their targets in the European sea bass are listed in Supplemental Table S3. For the first stage of the present study (i.e., the M stage), five miRNAs were found; of these, three belonged to the miR-430 family. In zebrafish, it has been shown that the miR-430 family comprised 72 members and that they targeted a large number of maternal mRNAs (Giraldez et al., 2006). The results of miR-430 target search using TargetScanFish in the present study are illustrated in Figure 6. Among the significantly enriched biological process categories were the GO terms “establishment or maintenance of cell polarity,” “protein targeting,” “response to growth factor,” “cellular response to growth factor stimulus,” “response to endogenous stimulus,” and “cellular response to endogenous stimulus.” It has been documented in zebrafish that miR-430a directs cell division, which in turn leads to neural tube development (Takacs and Giraldez, 2016); this is also consistent with the findings of this study. Concerning miR-430d, two clusters were found to comprise the mature mirR-430d sequence. These clusters are located on different regions of the genome and appear to fold into the miRNA typical stem-loop (Supplemental Figure S2).

The FINS stage, determined by the end of metamorphosis and the start of squamation, did not show the presence of any unique DE miRNA. Nevertheless, one of the highly abundant small RNAs at this stage, compared to the others, was miR-462, which has been linked to growth and muscle development in the blunt snout bream (Megalobrama amblycephala; Yi et al., 2013). The authors have found miR-462 to have a higher count number in small-sized fish than in bigger fish. However, in both cases, a high count number was reported. Among the few miRNAs being regulated between stages FLX and FINS (9), hs-miR-7641-1 and hs-miR-7641-2 were found to be more abundant at the FLX stage. miR-7641 in humans is known as a regulator of ribosomal proteins (Reza et al., 2017). For humans, about 3,500 targets have been identified (TargetScanHuman, release 7.2, March 2018). The respective search in zebrafish did not result in any match. Target search using the European sea bass transcriptome identified the gene paralemmin-1 as a putative target for miR-7641-1/miR-7641-2. Paralemmin-1 is also among the 3,500 targets identified for humans and is known to play an important role in filopodia induction and spine maturation (Arstikaitis et al., 2008). In contrast, hhi-miR-7641 (H. hippoglossus) was detected in all stages, except in GO and O, with the highest expression of hhi-miR-7641 found at the 24 hph stage.

In conclusion, investigating the miRNA repertoire during the early development of a commercially important fish species such as the European sea bass may contribute to a better understanding of regulatory processes during early ontogenesis. In the present study, we identified 2,115 ncRNA transcripts in the European sea bass; of these, 684 were annotated as miRNA. Distinct ncRNA expression profiles of 10 developmental stages were recognized and stage-specific miRNA were identified. Putative targets were also detected to provide first insights into miRNA involvement during early development.

Data Availability

The datasets generated for this study can be found in NCBI SRA database, PRJNA269278 and PRJNA369460.

Ethics Statement

All experiments/methods in the present study were performed in accordance with the approved guidelines and regulations from the HCMR institutional animal care and use committee following the three Rs (Replacement, Reduction, Refinement) guiding principles for more ethical use of animals in testing, first described by Russell and Burch in 1959 (EU Directive 2010/63). These principles are now followed in many testing establishments worldwide prior to initiation of experiments.

Author Contributions

ES conceived the study, drafted the manuscript, and performed the bioinformatics analyses. EK carried out RNA extractions and validation and helped in the miRNA work. AT and NP conceived the sampling and collected them for further analysis. MP initiated the study, critically revised the manuscript, and approved the final version.

Funding

Financial support for this study has been provided by the Ministry of Education and Religious Affairs, under the Call “ARISTEIA I” of the National Strategic Reference Framework 2007–2013 (ANnOTATE), co-funded by the EU and the Hellenic Republic through the European Social Fund, the European Union Seventh Framework Program (FP7 2010–2014) under the grant agreement no 265957 (CopeWell).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Financial support for this study was provided by the Ministry of Education and Religious Affairs, under the Call “ARISTEIA I” of the National Strategic Reference Framework 2007–2013 (ANnOTATE), co-funded by the European Union and the Hellenic Republic through the European Social Fund as well as the European Union Seventh Framework Program (FP7 2010–2014) under grant agreement no. 265957 (CopeWell). We further, would like to thank the Cornell University Core Laboratories Center for sequencing provision and the Informatics Group of the IMBBC for computational support. We acknowledge support of this work by the project “Centre for study and sustainable exploitation of Marine Biological Resources (CMBR)” (MIS 5002670) which is implemented under the Action “Reinforcement of the Research and Innovation Infrastructure”, funded by the Operational Programme “Competitiveness, Entrepreneurship and Innovation” (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00657/full#supplementary-material

Supplemental Figure S1 | Total RNA extraction was evaluated on an RNA Nano Bioanalysis chip (Agilent) as well as by gel electrophoreses. Only samples passing this evaluation step as presented here were used for NGS. (A) Total RNA. (B) Small RNA.

Supplemental Figure S2 | miR-430d sequences along with their identified stem-loop.

Supplemental Table S1 | Identified small RNAs mapped onto the European sea bass genome with threshold e-value of < 0.00005.

Supplemental Table S2 | Normalized mean read counts of significant DE transcripts.

Supplemental Table S3 | miRNAs found in only one stage and zero counts in the other stages along with their putative target genes, annotation, and accession number.

Supplemental Table S4 | GO terms of putative target genes regulated by miR-21 and miR-430.

References

Arstikaitis, P., Gauthier-Campbell, C., Carolina Gutierrez Herrera, R., Huang, K., Levinson, J. N., Murphy, T. H., et al. (2008). Paralemmin-1, a modulator of filopodia induction is required for spine maturation. Mol. Biol. Cell 19, 2026–2038. doi:10.1091/mbc.e07-08-0802

PubMed Abstract | CrossRef Full Text | Google Scholar

Barozai, M. Y. K. (2012). The microRNAs and their targets in the channel catfish (Ictalurus punctatus). Mol. Biol. Rep. 39, 8867–8872. doi: 10.1007/s11033-012-1753-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bizuayehu, T. T., Babiak, J., Norberg, B., Fernandes, J. M. O., Johansen, S. D., Babiak, I. (2012a). Sex-biased miRNA expression in atlantic halibut (Hippoglossus hippoglossus) brain and gonads. Sex. Dev. 6, 257–266. doi: 10.1159/000341378

PubMed Abstract | CrossRef Full Text | Google Scholar

Bizuayehu, T. T., Lanes, C. F. C., Furmanek, T., Karlsen, B. O., Fernandes, J. M. O., Johansen, S. D., et al. (2012b). Differential expression patterns of conserved miRNAs and isomiRs during Atlantic halibut development. BMC Genomics. 13, 11. doi: 10.1186/1471-2164-13-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Campos, C., Sundaram, A. Y. M., Valente, L. M. P., Conceição, L. E. C., Engrola, S., Fernandes, J. M. O. (2014). Thermal plasticity of the miRNA transcriptome during Senegalese sole development. BMC Genomics. 15, 525. doi: 10.1186/1471-2164-15-525

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, P. Y., Manninga, H., Slanchev, K., Chien, M., Russo, J. J., Ju, J., et al. (2005). The developmental miRNA profiles of zebrafish as determined by small RNA cloning. Genes Dev. 19, 1288–1293. doi: 10.1101/gad.1310605

PubMed Abstract | CrossRef Full Text | Google Scholar

Chistiakov, D. A., Tsigenopoulos, C. S., Lagnel, J., Guo, Y. M., Hellemans, B., Haley, C. S., et al. (2008). A combined AFLP and microsatellite linkage map and pilot comparative genomic analysis of European sea bass Dicentrarchus labrax L. Anim. Genet. 39, 623–634. doi: 10.1111/j.1365-2052.2008.01786.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., Marks, D. S. (2003). MicroRNA targets in Drosophila. Genome Biol. 5, R1. doi: 10.1186/gb-2003-5-1-r1

PubMed Abstract | CrossRef Full Text | Google Scholar

Forbes, D. J., Kornberg, T. B., Kirschner, M. W. (1983). Small nuclear RNA transcription and ribonucleoprotein assembly in early Xenopus development. J. Cell Biol. 97, 62–72. doi: 10.1083/jcb.97.1.62

PubMed Abstract | CrossRef Full Text | Google Scholar

Gavery, M. R., Roberts, S. B. (2017). Epigenetic considerations in aquaculture. Peer J. doi: 10.7717/peerj.4147

CrossRef Full Text | Google Scholar

Ghorai, A., Ghosh, U. (2014). miRNA gene counts in chromosomes vary widely in a species and biogenesis of miRNA largely depends on transcription or post-transcriptional processing of coding genes. Front. Genet. 5, 100. doi: 10.3389/fgene.2014.00100

PubMed Abstract | CrossRef Full Text | Google Scholar

Giraldez, A. J., Mishima, Y., Rihel, J., Grocock, R. J., Van Dongen, S., Inoue, K., et al. (2006). Zebrafish MiR-430 promotes deadenylation and clearance of maternal mRNAs. Science 312, 75–79. doi: 10.1126/science.1122689

PubMed Abstract | CrossRef Full Text | Google Scholar

Griffiths-Jones, S., Saini, H. K., Van Dongen, S., Enright, A. J. (2008). miRBase: tools for microRNA genomics. Nucleic Acids Res. 36, D154–D158. doi: 10.1093/nar/gkm952

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Z., Maki, M., Ding, R., Yang, Y., Zhang, B., Xiong, L. (2014). Genome-wide survey of tissue-specific microRNA and transcription factor regulatory networks in 12 tissues. Sci. Rep. 4, 5150. doi:10.1038/srep05150.

PubMed Abstract | CrossRef Full Text | Google Scholar

Guyon, R., Senger, F., Rakotomanga, M., Sadequi, N., Volckaert, F. A. M., Hitte, C., et al. (2010). A radiation hybrid map of the European sea bass (Dicentrarchus labrax) based on 1581 markers: synteny analysis with model fish genomes. Genomics 96, 228–238. doi: 10.1016/j.ygeno.2010.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Haas, B. J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P. D., Bowden, J., et al. (2013). De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8, 1494–1512. doi: 10.1038/nprot.2013.084

PubMed Abstract | CrossRef Full Text | Google Scholar

He, X., Yan, Y. L., Eberhart, J. K., Herpin, A., Wagner, T. U., Schartl, M., et al. (2011). MiR-196 regulates axial patterning and pectoral appendage initiation. Dev. Biol. 357, 463–477. doi: 10.1016/j.ydbio.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Herkenhoff, M. E., Oliveira, A. C., Nachtigall, P. G., Costa, J. M., Campos, V. F., Hilsdorf, A. W. S., et al. (2018). Fishing into the microRNA transcriptome. Front. Genet. 9, 88. doi: 10.3389/fgene.2018.00088

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaitetzidou, E., Xiang, J., Antonopoulou, E., Tsigenopoulos, C. S., Sarropoulou, E. (2015). Dynamics of gene expression patterns during early development of the European seabass (Dicentrarchus labrax). Physiol. Genomics. 47, 158–169. doi: 10.1152/physiolgenomics.00001.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, V. N. (2005). MicroRNA biogenesis: coordinated cropping and dicing. Nat. Rev. Mol. Cell Biol. 5, 376–385. doi: 10.1038/nrm1644

CrossRef Full Text | Google Scholar

Kloosterman, W. P., Wienholds, E., Ketting, R. F., Plasterk, R. H. A. (2004). Substrate requirements for let-7 function in the developing zebrafish embryo. Nucleic Acids Res. 34, 6284–6291. doi: 10.1093/nar/gkh968

CrossRef Full Text | Google Scholar

Krüger, J., Rehmsmeier, M. (2006). RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 34, W451–W454. doi: 10.1093/nar/gkl243

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhl, H., Beck, A., Wozniak, G., Canario, A. V. M., Volckaert, F. A. M., Reinhardt, R. (2010). The European sea bass Dicentrarchus labrax genome puzzle: comparative BAC-mapping and low coverage shotgun sequencing. BMC Genomics. 11, 68.doi: 10.1186/1471-2164-11-68

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumarswamy, R., Volkmann, I., Thum, T. (2011). Regulation and function of miRNA-21 in health and disease. RNA Biol. 8, 706–713. doi: 10.4161/rna.8.5.16154

PubMed Abstract | CrossRef Full Text | Google Scholar

Labbé, C., Robles, V., Herraez, M. P. (2017). Epigenetics in fish gametes and early embryo. Aquaculture. 472, 93–106. doi: 10.1016/j.aquaculture.2016.07.026

CrossRef Full Text | Google Scholar

Mishima, Y. (2012). Widespread roles of microRNAs during zebrafish development and beyond. Dev. Growth Differ. 54, 55–65. doi: 10.1111/j.1440-169X.2011.01306.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nepal, C., Coolen, M., Hadzhiev, Y., Cussigh, D., Mydel, P., Steen, V. M., et al. (2015). Transcriptional, post-transcriptional and chromatin-associated regulation of pri-miRNAs, pre-miRNAs and moRNAs. Nucleic Acids Res. 44, 3070–3081. doi: 10.1093/nar/gkv1354

PubMed Abstract | CrossRef Full Text | Google Scholar

Pavlidis, M., Karantzali, E., Fanouraki, E., Barsakis, C., Kollias, S., Papandroulakis, N. (2011). Onset of the primary stress in European sea bass Dicentrarchus labrax, as indicated by whole body cortisol in relation to glucocorticoid receptor during early development. Aquaculture 315, 125–130. doi: 10.1016/j.aquaculture.2010.09.013

CrossRef Full Text | Google Scholar

Pinto, P. I. S., Thorne, M. A. S., Power, D. M. (2017). European sea bass (Dicentrarchus labrax) skin and scale transcriptomes. Mar. Genomics. 35, 35–37. doi: 10.1016/j.margen.2017.05.002

CrossRef Full Text | Google Scholar

Pittman, K., Yúfera, M., Pavlidis, M., Geffen, A. J., Koven, W., Ribeiro, L., et al. (2013). Fantastically plastic: fish larvae equipped for a new world. Rev. Aquac. 5, S224–S267. doi: 10.1111/raq.12034

CrossRef Full Text | Google Scholar

Ramachandra, R. K., Salem, M., Gahr, S., Rexroad, C. E., Yao, J. (2008). Cloning and characterization of microRNAs from rainbow trout (Oncorhynchus mykiss): their expression during early embryonic development. BMC Dev. Biol. 8, 41. doi: 10.1186/1471-213X-8-41

PubMed Abstract | CrossRef Full Text | Google Scholar

Reza, A. M. M. T., Choi, Y. J., Yuan, Y. G., Das, J., Yasuda, H., Kim, J. H. (2017). MicroRNA-7641 is a regulator of ribosomal proteins and a promising targeting factor to improve the efficacy of cancer therapy. Sci. Rep. 7, 8365. doi: 10.1038/s41598-017-08737-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Robledo, D., Martin, A. P., Álvarez-Dios, J. A., Bouza, C., Pardo, B. G., Martínez, P. (2017). First characterization and validation of turbot microRNAs. Aquaculture. 472, 76–83. doi: 10.1016/j.aquaculture.2016.05.002

CrossRef Full Text | Google Scholar

Sarropoulou, E., Nousdili, D., Magoulas, A., Kotoulas, G. (2008). Linking the genomes of nonmodel teleosts through comparative genomics. Mar. Biotechnol. 10, 227–233. doi: 10.1007/s10126-007-9066-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarropoulou, E., Sepulcre, P., Poisa-Beiro, L., Mulero, V., Meseguer, J., Figueras, A., et al. (2009). Profiling of infection specific mRNA transcripts of the European seabass Dicentrarchus labrax. BMC Genomics 10, 157. doi: 10.1186/1471-2164-10-157

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarropoulou, E., Galindo-Villegas, J., García-Alcázar, A., Kasapidis, P., Mulero, V. (2012). Characterization of European sea bass transcripts by RNA SEQ after oral vaccine against V. anguillarum. Mar. Biotechnol. 14, 634–642. doi: 10.1007/s10126-012-9466-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Takacs, C. M., Giraldez, A. J. (2016). miR-430 regulates oriented cell division during neural tube development in zebrafish. Dev. Biol. 409, 442–450. doi: 10.1016/j.ydbio.2015.11.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Tine, M., Kuhl, H., Gagnaire, P. A., Louro, B., Desmarais, E., Martins, R. S. T., et al. (2014). European sea bass genome and its variation provide insights into adaptation to euryhalinity and speciation. Nat. Commun. 5, 5770. doi: 10.1038/ncomms6770

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsalafouta, A., Papandroulakis, N., Pavlidis, M. (2015). Early life stress and effects at subsequent stages of development in European sea bass (D. labrax). Aquaculture. 436, 27–33. doi: 10.1016/j.aquaculture.2014.10.042

CrossRef Full Text | Google Scholar

Ulitsky, I., Shkumatava, A., Jan, C. H., Subtelny, A. O., Koppstein, D., Bell, G. W., et al. (2012). Extensive alternative polyadenylation during zebrafish development. Genome Res. 10, 2054–2066. doi: 10.1101/gr.139733.112

CrossRef Full Text | Google Scholar

Valadkhan, S., Gunawardane, L. S. (2013). Role of small nuclear RNAs in eukaryotic gene expression. Essays Biochem. 54, 79–90. doi: 10.1042/bse0540079

PubMed Abstract | CrossRef Full Text | Google Scholar

Varsamos, S., Flik, G., Pepin, J. F., Bonga, S. E. W., Breuil, G. (2006). Husbandry stress during early life stages affects the stress response and health status of juvenile sea bass, Dicentrarchus labrax. Fish Shellfish Immunol. 20, 83–96. doi: 10.1016/j.fsi.2005.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ward, N. J., Green, D., Higgins, J., Dalmay, T., Münsterberg, A., Moxon, S., et al. (2018). microRNAs associated with early neural crest development in Xenopus laevis. BMC Genomics 19, 59. doi: 10.1186/s12864-018-4436-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wienholds, E., Plasterk, R. H. A. (2005). MicroRNA function in animal development. FEBS Lett. 579, 5911–5922. doi: 10.1016/j.febslet.2005.07.070

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, Y., Ma, L., Jia, Q., Deng, W., Liu, Z., Zhang, Y., et al. (2014). Systematic characterization of small RNAome during zebrafish early developmental stages. BMC Genomics 15, 117. doi: 10.1186/1471-2164-15-117

PubMed Abstract | CrossRef Full Text | Google Scholar

Yi, S., Gao, Z. X., Zhao, H., Zeng, C., Luo, W., Chen, B., et al. (2013). Identification and characterization of microRNAs involved in growth of blunt snout bream (Megalobrama amblycephala) by Solexa sequencing. BMC Genomics 14, 754. doi: 10.1186/1471-2164-14-754

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Wu, J. W., Wang, Y., Zhao, J. L. (2016). Role of miR-21 in alkalinity stress tolerance in tilapia. Biochem. Biophys. Res. Commun. 471, 26–33. doi: 10.1016/j.bbrc.2016.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, B. W., Zhou, L. F., Liu, Y. L., Wan, S. M., Gao, Z. X. (2017). Evolution of fish let-7 microRNAs and their expression correlated to growth development in blunt snout bream. Int. J. Mol. Sci. 18, 3. doi: 10.3390/ijms18030646

CrossRef Full Text | Google Scholar