Physiol Genomics 47: 158–169, 2015. First published March 3, 2015; doi:10.1152/physiolgenomics.00001.2015.

Dynamics of gene expression patterns during early development of the European seabass (Dicentrarchus labrax) E. Kaitetzidou,1,2 J. Xiang,3 E. Antonopoulou,2 C. S. Tsigenopoulos,1 and E. Sarropoulou1 1

Institute of Marine Biology, Biotechnology and Aquaculture, Hellenic Centre for Marine Research, Greece; 2School of Biology, Faculty of Science, Aristotle University of Thessaloniki, Greece; and 3Genomics Resources Core Facility, Weill Cornell Medical College, New York, New York Submitted 9 January 2015; accepted in final form 23 February 2015

Kaitetzidou E, Xiang J, Antonopoulou E, Tsigenopoulos CS, Sarropoulou E. Dynamics of gene expression patterns during early development of the European seabass (Dicentrarchus labrax). Physiol Genomics 47: 158 –169, 2015. First published March 3, 2015; doi:10.1152/physiolgenomics.00001.2015.—Larval and embryonic stages are the most critical period in the life cycle of marine fish. Key developmental events occur early in development and are influenced by external parameters like stress, temperature, salinity, and photoperiodism. Any failure may cause malformations, developmental delays, poor growth, and massive mortalities. Advanced understanding of molecular processes underlying marine larval development may lead to superior larval rearing conditions. Today, the new sequencing and bioinformatic methods allow transcriptome screens comprising messenger (mRNA) and microRNA (miRNA) with the scope of detecting differential expression for any species of interest. In the present study, we applied Illumina technology to investigate the transcriptome of early developmental stages of the European seabass (Dicentrarchus labrax). The European seabass, in its natural environment, is a euryhaline species and has shown high adaptation processes in early life phases. During its embryonic and larval phases the European seabass lives in a marine environment and as a juvenile it migrates to coastal zones, estuaries, and lagoons. Investigating the dynamics of gene expression in its early development may shed light on factors promoting phenotypic plasticity and may also contribute to the improvement and advancement of rearing methods of the European seabass, a species of high economic importance in European and Mediterranean aquaculture. We present the identification, characterization, and expression of mRNA and miRNA, comprising paralogous genes and differentially spliced transcripts from early developmental stages of the European seabass. We further investigated the detection of possible interactions of miRNA with mRNA. RNA-Seq; transcriptomics; aquaculture; European seabass; Dicentrarchus labrax; miRNA; mRNA; development FISH DEVELOPMENT, LIKE IN all vertebrates, is determined during the development of an embryo (embryogenesis) and can be considered as a continuum of morphological changes. To enable the analysis of embryogenesis this continuum has been broken down into distinct developmental stages, and much effort has been invested to define those developmental stages along with their critical aspects. It is known that embryonic and larval development in vertebrates share main developmental processes, like the determination of the overall body plan and organogenesis (18, 82). On the other hand, studies also have shown high phenotypic plasticity during development due to species particularities (36, 43). Phenotypic plasticity in teleost

Address for reprint requests and other correspondence: E. Sarropoulou, Inst. of Marine Biology, Biotechnology and Aquaculture, Hellenic Centre for Marine Research, Thalassocosmos, Gournes Pediados, PO Box 2214, 71003 Heraklion, Crete, Greece (e-mail: [email protected]). 158

may be explained by fish-specific genome duplication (FSGD) showing different retention rates among species as well as different expression levels. This is important especially concerning transcription factors and signaling genes responsible for the regulatory control of development. Furthermore, regulatory processes by microRNA (miRNA), chromatin remodeling (DNA methylation and histone modification), transcription factors, or mRNA processing are important factors promoting phenotypic plasticity (65). Transcriptomic studies in embryogenesis may shed light on the underlining mechanisms of morphological developmental changes and larval plasticity. Here, we have applied Illumina HiSeq sequencing technology to investigate the transcriptome of the early embryogenesis of the European seabass (Dicentrarchus labrax). The European seabass has shown high adaptation processes in early life phases (15), and it is of high commercial value in the Mediterranean (FAO, 2012, http://www.fao.org/fishery/species/ 2291/en). During its embryonic and larval phases the European seabass lives in a marine environment, and as a juvenile it migrates to coastal zones, estuaries, and lagoons. The genomic toolbox of the European seabass has been significantly enriched in the recent years, facilitating the identification and characterization of genes (9, 13, 29, 81). In addition, some high-throughput gene expression studies have been performed in the European seabass to better understand the role of genes implicated in several physiological mechanisms such as growth rate and immune response (e.g., Ref. 71), nutrition (27), salinity tolerance (9), as well as late larval development and metamorphosis (16). Several high-throughput transcriptomic studies have also be performed in other species of commercial interest, like, e.g., the gilt-head sea bream (Sparus aurata), Atlantic cod (Gadus morhua), and catfish (Ictalurus punctatus) (10, 46, 55). However, little or no attention has been paid to paralogous genes despite the FSGD event that also gave rise to changes in regulatory control (22, 40, 80, 84). It is also known that paralogous genes may have altered regulation (subfunctionalization) and/or novel function (neofunctionalization). Interestingly, studies also have shown that for some genes the expression of a paralog gradually decreases, whereas that of the other paralog increases (70, 73, 74). In addition, Garcia de la Serrana et al. (26) recently highlighted the systematic differences of paralogous gene retainment between Acanthopterygii and Ostariophysi. Most of the studies looking at the role of various paralogous genes have been performed in the model fish zebrafish (Danio rerio) (40, 91, 98). Zebrafish belongs to the superorder Ostariophysi, and hence a high amount of data is available, whereas less information has been produced for the superorder Acanthopterygii as it comprises mainly nonmodel teleost species.

1094-8341/15 Copyright © 2015 the American Physiological Society

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

To better understand the mechanisms underlying developmental processes, besides the exploration of differentially expressed transcripts during early development, we must also address the effect of regulative elements such as miRNAs or small RNAs. It has been shown that miRNAs are largely expressed in animals, plants, and viruses (62, 79) and that they are involved in a wide range of gene regulation of various biological processes. In humans miRNA are extensively studied, and it is known that altered miRNA expression plays an important role in human diseases and especially in cancer (2, 38). Next-generation sequencing technologies have enabled the high-throughput research of miRNA in model as well as in nonmodel species. Thus, during the last decade some studies were undertaken to investigate tissue-specific miRNAs as well as different developmental stages (2, 6, 7, 12, 89). Most studies of teleosts have been performed in zebrafish (12, 21, 86, 88, 97). In other teleosts, e.g., in Atlantic salmon (5), rainbow trout (69), medaka (Oryzias latipes) (52), and stickleback (Gasterosteus aculeatus) (44), miRNAs have been described independently of tissue. Tissue-specific, i.e., muscle-specific, miRNA expression is shown for, e.g., tilapia and zebra finch (57, 93). The present study uses Illumina technology to investigate the transcriptome, comprising mRNA and miRNA transcripts, of the European seabass during seven early developmental stages (from morula to 24 h after hatching). We assessed paralogous genes as well as differentially spliced transcripts. Expression profiles of mRNA in each developmental stage of the European seabass are depicted. The study of miRNA focuses on the identification of miRNAs in the European seabass that have already been described in other species (new miRNA) and their expression. Finally first approaches of regulatory control during development via miRNAs are carried out. MATERIALS AND METHODS

All procedures involving the handling and treatment of fish used during this study were approved by the Hellenic Centre for Marine Research 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. The larvae were anesthetized with 100 –200 mg/l MS222 (tricaine methanesulfonate, Sigma-Aldrich) depending on fish size. Sampling. Eggs and larvae of seabass were collected in tanks of the Institute of Marine Biology, Biotechnology and Aquaculture at the Hellenic Centre for Marine Research in Crete, Greece. In total we collected samples from seven developmental stages, comprising the morula, 50% epiboly, late gastrula-organogenesis, organogenesis-embryo (⫽ 2⁄3 of meridian), late organogenesis, hatching, and 24 h after hatching. Staging was performed following gilt-head sea bream’s developmental stages (17). All samples were immersed in the RNAlater (Ambion, Austin, TX) and transferred to a ⫺80°C ultralow freezer until preparation of RNA. mRNA and miRNA extraction. mRNA and miRNA were extracted from all developmental stages sampled, with a Nucleospin miRNA kit for isolation of small and mRNA (Macherey-Nagel, Duren, Germany) according to the manufacturer’s instructions. We disrupted pools of eggs and larvae with mortar and pestle in liquid nitrogen and homogenized them in lysis buffer by passing lysate through a 23-gauge (0.64 mm) needle five times. Quantity of RNA was estimated with NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE), and the quality was further evaluated by agarose (1%)

159

gel electrophoresis and Agilent 2100 Bioanalyzer using RNA Nano Bioanalysis chip. All samples had an RNA integrity number value between 8.9 and 9.9. RNA sequencing. mRNA and miRNA of all developmental stages was sequenced at Cornell University Core Laboratories Center using Illumina HiSeq vs 2000. For mRNA paired-end libraries and for miRNA single-end libraries were prepared and sequenced over two Illumina lanes (⬃38 Gb). Use of multiplex identifier tags allowed us to distinguish among RNA of different stages. Quality control and de novo assembly. Raw data were sorted by developmental stage. Quality control was assessed with FastQC (version 0.10.0; http://www.bioinformatics.babraham.ac.uk/projects/ fastqc). Low-quality reads were discarded by Trimmomatic software (8). Sequences of all stages were assembled into a reference transcriptome using Trinity version 2012-06-08 (28). Differential expression analysis. Reads of each developmental stage were mapped to the reconstructed transcriptome with Bowtie (47, 48), allowing a maximum of three mismatches for each read. RSEM (version 1.2.3; Ref. 50) was used for quantification of read abundances, and transcripts represented less than once per million mappable reads were excluded from the following analysis. Differential expression between different stages was assessed with the R Bioconductor package DESeq (3), and transcripts P ⬍ 0.05 in their expression were assigned as differentially expressed. BLAST, annotation, and classification. The filtered sequences were compared using BLAST tools (version 2.2.25; Ref. 1) to a nonredundant protein database and a nonredundant nucleotide database using. The results were analyzed with Blast2GO software (14). First, sequences were mapped to Gene Ontology (GO) annotations according to the GO database. Thereafter, sequences that were successfully mapped were annotated according to the following GO terms: cellular component, molecular function, and biological process. To determine putative paralogous genes, we undertook manual evaluation of the alignment analysis to exclude the classification to different clusters due to alternative splicing. Subsequently all putative paralogous genes were mapped onto the published genome of the European seabass (81), and BLAT searches (41) were performed against five teleost species (zebrafish medaka, stickleback, tetraodon, and fugu). Cluster analysis. Significantly differentially expressed transcripts were partitioned according to their expression pattern by the k-means clustering method in the R statistical package (version 3.0.2), with the “Hartigan-Wong” algorithm (iteration number ⫽ 200). The number of centers was determined by the plot of the within-groups sum of squares by number of clusters (32). Venn diagram. A Venn diagram was constructed with the web application “jvenn” provided for free and nonprofit use from GenoTool Bioinformatics (http://bioinfo.genotoul.fr/jvenn) (4) from all the differentially expressed genes in six stages and comparing them with the first (morula) stage. Enrichment analysis. Fischer’s exact test tool in Blast2GO was used with term filter modes “FDR” and “pval,” term filter values: false discovery rate (FDR) ⬍ 0.05, FDR ⬍ 0.01, and P ⬍ 0.005, as well as two-tailed test. Transcripts grouped into the first and the fourth cluster after k-means partitioning were used as test datasets that were compared with the annotated transcriptome (reference dataset). miRNA analysis. All reads were trimmed with Trimmomatic software (8), and a quality check was performed with FastQC (version 0.10.0; http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Trimmed sequences were clustered with the cdHIT (53) cluster program, and poorly supported transcripts were filtered out for downstream analysis. Differential expression analysis was carried out with the R Bioconductor package DESeq (3), and transcripts with P value ⬍ 0.005 were used for miRNA identification by National Center for Biotechnology Information (NCBI) database search as well as miRBase search. Putative targets were identified with the online version of RNAhybrid (45) and the miRanda web tool (19). Prediction was proceeded with the default parameters and an energy threshold of mfe ⬍ ⫺20 Kcal.

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

160

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

Network analysis. A module network was created by the LeMoNe software (61) from highly significant (P value ⬍ 0.005) differentially expressed transcripts detected between all stages under study and expression data from the 42 known miRNA. The normalized log2-fold values were used as input of transcript expression, and the 42 known miRNAs were used as potential regulators. The output of the software presents a set of clusters of coexpressed transcripts and one or more assigned regulator(s) of particular weight. Cytoscape (version 2.7.0) (76) was used for the visualization of the clusters and their correlation with particular regulators and with each other.

Table 2. Summary of reads used for de novo assembly before and after filtering Assembly Statistics

Transcripts Contig N50 Total bases Min. sequence length, bp Max. sequence length, bp Average sequence length, bp

Before Filtering

After Filtering

201,130 3,015 296,031,440 201 27,326 1,471

69,794 2,769 119,954,223 201 27,326 1,718

RESULTS

Transcriptome sequencing of early developmental stages of seabass. RNA Illumina sequencing of European seabass’ seven developmental stages resulted in 149,642,793 of 100 bp pairedend reads in total. Raw sequence data have been deposited in NCBI’s Short Read Archive database under the accession numbers SRX796243 (morula), SRX796245 (50% epiboly), SRX796246 (late gastrula-organogenesis), SRX796264 (organogenesis-embryo ⫽ 2⁄3 of meridian), SRX796265 (late organogenesis), SRX796276 (hatching), and SRX796587 (24 h after hatching). Trimming of low-quality reads and adaptors yielded 129,467,381 paired-end reads of high quality, which were used for the assembly analysis (Table 1). De novo assembly of the 129,467,381 reads resulted in 201,130 transcripts of an average length of 1,472 bp. The filtering process led to a final reference transcriptome of 69,794 transcripts, which was used for downstream analysis (Table 2). Gene identification and annotation. More than half of the transcripts, 40,898 (58.6%), had a significant BLASTX similarity matches with sequences deposited in the NCBI nonredundant database (nonredundant GenBank CDS translations ⫹ RefSeq a ⫹ PDP ⫹ SwitProt ⫹ PIR ⫹ PRF). Almost 28,000 transcripts had significant similarity with sequences of two fish species belonging to Perciformes: Maylandia zebra and Oreochromis niloticus. Among the retrieved BLAST matches a set of transcripts showed similar or even the same BLAST matches indicating either different splice forms or paralogous genes. A subset was further examined by a comparative mapping approach. Table 3 shows putative paralogous genes, their genome location, as well as their homolog group in five teleost species. Blast2GO successfully mapped 20,991 sequences and gene ontology (GO) terms were assigned to 16,459 transcripts. Under the category “biological process” 3,858 transcripts (7.24%) were classified to the GO term “developmental processes”. Other GO child terms of biological process are cellular process (20.21%), metabolic process (15.93%), single-organism process (14.74%), biological regulation (9.57%), multicelTable 1. Sequencing and quality filtering summary Stage

Input Read Pairs

Both Ends Survived

Morula 50% epiboly Late gastrulation-organogenesis Organogenesis-embryo ⫽ 3/4 of the meridian Late organogenesis Hatching 24 h after hatching

19,436,879 21,404,651 21,101,942

17,416,068 (89.60%) 17,614,326 (82.29%) 17,831,063 (84.50%)

17,234,236 22,929,410 20,954,714 26,562,961

14,913,368 (86.53%) 21,179,558 (92.37%) 18,412,965 (87.87%) 22,100,033 (83.20%)

lular organismal process (7.38%), response to stimulus (6.48%), localization (5.08%), signaling (4.60%), cellular component organization or biogenesis (4.42%), locomotion (1.21%), immune system process (0.97%), growth (0.70%), biological adhesion (0.65%), reproduction (0.56%) and multiorganism process (0.27%). Assessment of differential gene expression. Transcripts with P value ⱕ 0.05 were assigned as differentially expressed. The number of differentially expressed genes comparing each pair of condition is given in Table 4. Thus after pairwise comparison, 18,140 unique transcripts were found to be significantly differentially expressed and were used for downstream analysis. Subsequently those transcripts found to be differentially expressed at least in one of the stages were subjected to k-means clustering. Cluster methods are widely used for categorizing genes according to their expression patterns. k-means clustering partitioned significantly differentially expressed transcripts into four clusters determined from the within sum of squares by number of clusters plot (Fig. 1A). Further discriminant analysis showed that 98.8% of original grouped cases are correctly classified. Gene expression profiles of each cluster are illustrated as heat maps in Fig. 1B. Each row of the heat map represents the pattern of one transcript’s expression, while columns represent the developmental stage. Cluster 1 contained 2,939 transcripts, while cluster 2, cluster 3, and cluster 4 were composed of 842; 12,968; and 1,391 transcripts, respectively. Enrichment analysis of cluster 1 and cluster 4 is shown in Fig. 2. The majority of GO terms in cluster 1 (FDR ⱕ 0.05) are related to developmental process. In cluster 4 the analysis with term mode FDR ⱕ 0.05 and FDR ⱕ 0.1 resulted in four and five GO terms, respectively. To retrieve a broader range of GO terms in this cluster we applied a term mode value of P ⱕ 0.005, resulting in GO terms mainly involved in regulation of gene expression and transcription. Additionally, differentially expressed genes of all seven stages studied were submitted to principal component analysis (PCA) shown in Fig. 3. The gene expression values of the stages as shown in the plot reflect the successiveness of the developmental stages. The first studied stage morula is noticeably separated from the remaining six stages. The last two stages, hatching and 24 h after hatching, are also grouped more closely together and further away from the other four stages. To extract stage-specific transcripts for the later stages, we constructed a Venn diagram using the morula stage as the reference stage. Therefore, differential expressed transcripts between morula and one of the remaining six stages were used, resulting in a total of 11,486 transcripts and representing 4,930 unique transcripts. The number of transcripts that

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

161

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

Table 3. Transcripts considered candidate paralogs in seabass after mapping to different chromosomes or gene groups of other fish species TranscriptID

comp48446_c1_seq3, comp48446_c1_seq5 comp50487_c0_seq1, comp50487_c0_seq2 comp44838_c0_seq2 comp50898_c0_seq3 comp53857_c0_seq1 comp51024_c0_seq1 comp57791_c14_seq5 comp50517_c0_seq3 comp55416_c0_seq10 comp55416_c0_seq3 comp55416_c0_seq7 comp55416_c0_seq5 comp50528_c0_seq2, comp50528_c0_seq4 comp50528_c0_seq3 comp57506_c7_seq1, comp57506_c7_seq4 comp57506_c7_seq2, comp57506_c7_seq7 comp26500_c0_seq1, comp49142_c0_seq2 comp47884_c0_seq2 comp47884_c2_seq1, comp47884_c2_seq2 comp56751_c9_seq2 comp51561_c0_seq5, comp51561_c0_seq9, comp51561_c0_seq20 comp54952_c0_seq3 comp56596_c1_seq8, comp56596_c1_seq15 comp26234_c0_seq1 comp38539_c0_seq1 comp47939_c0_seq2 comp52627_c0_seq6 comp54317_c2_seq1, comp54317_c2_seq3, comp54317_c2_seq7, comp54317_c2_seq9 comp60664_c0_seq1

BLAST Match

Accession No.

D. labrax

D. rerio

O. latipes

G. aculeatus

T. nigroviridis

T. rubripes

myozenin-2 like

XP_008284262

LG3

n/a

22

VII

Un

8

myozenin-2-like

XP_008284938

LG7

1

scf461

IX

Un

17

junctional adhesion molecule b-like junctional adhesion molecule b-like hepatocyte nuclear factor 4 alpha hepatocyte nuclear factor 4-beta-like isoform x2 orthodenticle 2 (otx2) protein orthodenticle 2 (otx2) protein homeobox protein engrailed-1-like homeobox protein engrailed-1-like homeobox protein engrailed ⫺2a like homeobox protein engrailed ⫺2b like homeobox protein orthopedia B-like homeobox protein orthopedia B-like transcritpion factor SOX-6like isoform Transcritpion factor SOX-6like isoform transcritpion factor SOX-6like isoform diencephalon mesencephalon homeobox protein 1-a-like diencephalon mesencephalon homeobox protein 1-a-like diencephalon mesencephalon homeobox protein 1-like serpin h1-like

XP_003452956

LG15

n/a

scf2132

XVI

2

1

XP_003440998

LG24

n/a

21

I

3

8

ACO56245

LG22–25

23

7

XII

9

HE591758

XP_004543170

LG5

n/a

7

II

n/a

n/a

XP_005916803.1

LG12

17

22

XV

10

2

XP_003453097

LG17

17

22

XVIII

14

16

XP_003456262

LG15

1

Uc 256

XVI

Un

1

XP_003440904

LG24

1

2

I

Un

8

XP_006795130

LG18–21

7

20

XXI

6

10

XP_003459215

LG10

7

17

III

Un

22

XP_005929900

LG19

21

9

XIV

4

6

XP_003451102

LG20

21

9

XIII

12

21

XP_004564551

LG6

7

6

XIX

13

9

XP_004564551

LG6

7

6

XIX

13

9

XP_004543251

LG5

7

3

II

5

13

XP_003975776

LG10

2

17

III

15

22

XP_004570628

LG10

6

17

III

15

22

XP_004553852

LG4

n/a

4

VIII

1

HE591782

XP_003441717

LG13

15

13

I

n/a

11

serpin h1-like serpin h1-like

XP_004561658 XP_003457702

LG14 LG3

10 n/a

13 18

VII VII

n/a Un

11 8

nodal homolog nodal homolog nodal homolog 2-a-like histone-lysine nmethyltransferase mll3like isoform x7 histone-lysine nmethyltransferase mll3like isoform x4

XP_004544784 XP_003965673 XP_003961216 XP_004542769

LG20 LG19 LG1B LG10

n/a 21 n/a 24

9 12 19 17

XIII XIV V III

12 4 2_random 15

21 6 1 HE591741

XP_004542766

LG10

2

17

III

15

HE591741

histone-lysine nmethyltransferase mll3-like

XP_004570816

UN

n/a

20

Un

4

10

scf, scaffold; Un, not random; n/a, not available.

are either shared between two or more particular sets of comparisons or unique in just one set of comparison is shown in the Venn diagram (Fig. 4). In total 271 differentially expressed transcripts were found to be in common in all six sets of comparisons. The first comparison (morula and 50% epi-

boly) included 576 unique transcripts, while the second (morula and late gastrulation-organogenesis), third (morula and organogenesis-embryo ⫽ ¾ of meridian), fourth (morula and late gastrulation), fifth (morula and stage of hatching), and sixth comparison (morula and 24 h after hatching stage) com-

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

162

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

Table 4. Double entry table with numbers of differentially expressed genes (P value ⬍ 0.05) comparing each pair of condition Morula

Morula 50 % epiboly Late gastrulation-organogenesis Organogenesis-embryo ⫽ 3/4 of meridian Late organogeneisis Hatching 24 h after hatching

50 % Epiboly

Late Gastrulation-Organogenesis

Organogenesis-Embryo ⫽ 3/4 of Meridian

Late Organogenesis

Hatching

24 h After Hatching

2,172 3,530

1,952 3,571 4,328

1,960 3,658 4,252 4,350

1,716 3,445 4,140 4,199

1,678 3,075 3,593 3,687

4,629

4,135 4,472

2,008 — — —

— —

—

— — —

— — —

— — —

prised 400, 235, 230, 243, and 434 unique transcripts, respectively. miRNA identification during seabass’ early developmental stages. Illumina sequencing of miRNA of European seabass’ early developmental stages produced 58,216,215 reads of an average length of 50 bp. Quality trimming and clustering decreased the number of unique reads to 745,139. Differential expression analysis further decreased the amount of miRNA reads to 1,928. In the present study we aimed to identify new miRNA (miRNAs described for other species than the European seabass) involved in European seabass development. Thus, after annotation a total of 42 miRNAs (Supplemental File 1) were identified, showing differential expression during the stages studied.1 The expression pattern of these 42 miRNAs is illustrated by the heat map in Fig. 5. 1

The online version of this article contains supplemental material.

— — —

— —

—

Network analysis and miRNA target site prediction. Fourteen out of 42 known miRNAs were correlated with 173 clusters of a total of 7,740 highly differentially expressed mRNA transcripts (Supplemental File 4, Fig. 6). LeMoNe software combined transcripts with common expression patterns in all pair-wise stage comparisons, to create clusters. Afterward, clusters were associated with the expression pattern of one or more of the 14 miRNA that were either compatible or opposite, indicating an indirect or direct regulation of mRNA expression, respectively. Six of the miRNAs that seem to have direct regulation of clusters were tested with RNAhybrid software (45) and the MiRanda v3.3a software (19) (Supplemental File 2). The software tests the possible hybridization of miRNA with a corresponding mRNA, providing a value for the energy needed. Transcripts assigned to the GO term developmental process and with a minimum free energy (mfe) of ⫺25 are shown in Fig. 7.

Fig. 1. k-means clustering of all differentially expressed transcripts. A: the optimal number of clusters were obtained by printing the cluster centroids obtained with k-means function using R. Number of clusters are given on the x-axes, and the sums of squares within groups on the y-axes. The plot determines 4 clusters as the best solution to partition the differentially expressed transcripts (red cycle). B: the 4 clusters determined in A are illustrated by heat map imaging (black arrow). Cluster 1 comprises 2,939 transcripts, cluster 2 842, cluster 3 12,968, and cluster 4 1,391. Stages are indicated under each column of clusters. Each row represents the expression of 1 transcript where shades of red represent upregulation and shades of green represents downregulation. Heat maps were drawn on the normalized counts of the transcripts of each developmental stage. Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

163

Fig. 2. Enrichment analysis. A: enrichment analysis with the set of transcripts grouped into the 1st (transcripts downregulated at the morula stage) cluster using the total of annotated transcriptome as reference set and false discovery rate (FDR) ⱕ 0.05. Gene Ontology (GO) terms related to biological process are shown with a red frame. B: enrichment analysis of the transcripts grouped into the 4th cluster (transcripts upregulated at the morula stage) using the annotated transcriptome as reference set and P ⱕ 0.005. GO terms found to be enriched with FDR ⱕ 0.1 and FDR ⱕ 0.05 are marked * or **, respectively. Red and green frames include GO terms related to biological process and gene expression. C: enrichment analysis with different term filter modes and values produced different number of GO terms. Symbol ⽧ marks the cases that are shown in A and B. DISCUSSION

The onset of next-generation sequencing technologies has enabled the research community to retrieve genomic and transcriptomic information for any organism of interest (30, 58, 63, 68, 85). The main characteristics of these platforms is their speed as well as the assessment of high amount of data at low cost. These advantages clearly have advanced clinical, environmental, as well as biotechnological research. Nevertheless, next-generation sequencing also comprises error rates as well as pitfalls during analysis and data interpretation. Several RNA sequencing data have been produced by applying GS454 FLX pyrosequencing, the first high-throughput next-generation sequencing platform. The advantages of GS454 FLX certainly are the capability to produce long reads (up to 1,000 bp with GSFLX⫹ technology) facilitating the following assembly process. However, for robust gene expression results highthroughput data are necessary, and thus a number of GSFLX runs have to be executed (25), which in turn results in high consumable cost. Alternatively, 3=-untranslated region (UTR) tagged cDNAs for library construction may decrease the amount of transcripts needed for robust differential gene expression studies (71). Today mainly the Illumina platform is applied due to its high throughput, low cost, and low sequencing error rates (54). In the present study, the transcriptome

comprising mRNA as well as miRNA of the European seabass was sequenced using paired-end reads and single-end reads, respectively. Each stage was individually tagged for downstream expression analysis. For mRNA sequencing the generation of paired-end reads resolves assembly problems due to repetitive regions. A reference transcriptome was constructed, resulting in 69,794 contigs with an average length of 1,718 bp and N50 of 2,769 bp (Table 2). The Illumina technology of sequencing paired-end reads facilitates the assembly, and the use of a combined assembly and mapping strategy results in a sound and accurate reference transcriptome that is also strongly supported by the obtained results of sequence average length and N50 length. We assessed transcript expression profiles by mapping paired-end reads against the reconstructed assembly. To obtain a robust transcript dataset we applied stringent filters, resulting in 69,794 transcripts. It has to be noted that the obtained transcripts do not reflect the gene number of an organism. Next-generation sequencing technologies provide a number of different transcript variants. Thus, BLAST searches may result in the same annotation and/or they result in matches entitled “transcript variant X1,” “transcript variant X2,” without any further information whether transcript variants are paralogous genes or alternative spliced transcript. It has been shown that paralogous genes arose due to whole genome

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

164

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

Fig. 3. Principal Component Analysis (PCA) of the expression profiles of the studied developmental stages. The 7 developmental stages are projected in the 2-dimensional plane by the expression profile of all transcripts differentially expressed. Red ellipses group together the stages that were neighboring both after PCA and in the developmental process.

duplication events, and some may have altered or acquired novel functions (sub- or neofunctionalization). However, detection of differential expression of transcript variants alone is not sufficient to identify putative paralogs as differentially spliced transcripts may also show stage-specific expression, such as the ATP-binding cassette subfamily D member 4-like

protein in the present study, where one transcript is significantly upregulated at the morula and the 50% epiboly stage, whereas the other two are upregulated at the later stage with the second having a peak at gastrulation and hatching and the third with an expression peak at the organogenesis-embryo ⫽ ¾ of the meridian and at the organogenesis stage. Consequently, to distinguish paralogous genes from simple transcript variants in the present study we pursued manual evaluation of transcript alignments as well as a comparative mapping approach. Transcripts mapping on different linkage groups of the European seabass or on different chromosomes in any of the five model species were considered as putative paralogous genes (Table 3). Linkage groups as well as chromosomes found were homologous groups as described in comparative mapping studies in seabass but also in the gilt-head sea bream compared with the five model fish species (29, 75). Examples of putative paralogs involved in developmental processes, found in the present study, are orthodenticle 2 (otx2) (59), nodal (13a), the cell surface receptors junctional adhesion molecule b-like (jamb) (66), diencephalon/ mesencephalon homeobox 1, and orthopedia (otp). The latter three genes were expressed especially during organogenesis of seabass embryos. It has been shown that both paralogs of diencephalon/mesencephalon homeobox play an important role in brain development of zebrafish embryos (91) and that both paralogs of otp are essential for the development of a subset of dopaminergic neurons (20). Jamb on the other hand has been shown to play a pivotal role along with jamc during somitogenesis in fusion of myocyte precursors (66). For further downstream analysis only transcripts significantly differentially expressed in at least one of the studied stages were used. The in total 18,140 differentially expressed transcript were clustered into four groups with specific expression patterns. Among the four clusters, clusters 1 and 4 are the most appeal-

Fig. 4. Venn diagram. Venn diagram showing the relation between the 4,930 significant differentially expressed transcripts (P ⱕ 0.05) after the comparison of morula stage with the other 6 stages (50% epiboly, late gastrula-organogenesis, organogenesis-embryo ⫽ 2⁄3 of meridian, late organogenesis, hatching, and 24 h after hatching). Total number of transcripts differentially expressed in each stage comparison is given next to the Venn diagram.

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

Fig. 5. Heat map imaging of microRNAs (miRNAs). Heat map of the 42 identified miRNAs throughout the 7 developmental stages of the European seabass. Stages are indicated under each column. Each row represents the expression of 1 miRNA, where red represents upregulation and green represents downregulation. Heat map was drawn on the normalized counts of miRNAs of each developmental stage.

ing ones, comprising transcripts being either down- (cluster 1) or up- (cluster 4) regulated at the morula stage (Fig. 1B). Morulation implements the repeated cleavage resulting in the morula stage, a bundle of cohering and sticky blastomers loosely arranged. It has been shown that embryonic stages

165

before gastrulation are more vulnerable (83a), and as early as three decades ago in medaka (O. latipes) stage sensitivity and decrease of mortality after the morula stage was shown (60). It is also known that the early embryonic processes are carried out by maternal factors and the zygotic genome is activated only at the midblastula transition. Enrichment analysis using transcripts belonging to cluster 1 against the total number of annotated transcripts showed a significant accumulation of GO terms involved in the developmental process (Fig. 2A). Transcripts in cluster 1 were mainly downregulated during the morula stage, which shows that at this stage transcripts important for the development of the organism are not yet activated. Transcripts found to be mainly upregulated during the morula stage (cluster 4) were classified in GO terms with regulative functions (Fig. 2B). In addition, GO terms “pole plasm” as well as “germ plasm” were enriched in cluster 4. Both terms belong to the ontology of cellular components and encode for differentiated cytoplasm, which is associated with a pole of an oocyte, egg, or early embryo (39). Apparently, the morula stage differentiates itself clearly form the other six stages. This is shown not only by cluster analysis as discussed before but also by a PCA plot (Fig. 3), where transcripts differentially expressed during the morula stage are isolated from the other six stages. In this line the morula stage was used as a reference stage to extract stage-specific differentially expressed transcripts (in total 11,486 transcripts, representing 4,930 unique transcripts) illustrated by a Venn diagram in Fig. 4. Out of the unique transcripts, only 5.5% were in common for all six stages, while transcripts found only at the stages late organogenesis, hatching, and l’ of meridian were between 4.6 and 4.9%. Transcripts solely detected at the stages late gastrulation and hatching after 24 h were 8.1 and 8.8%, respectively. The highest amount of transcripts (11.7%) uniquely found were at 50% epiboly (Fig. 4). The largest amount of transcripts significantly regulated only at the 50% epiboly stage was also found in zebrafish (83). The high amount of unique transcripts at 50% epiboly detected in zebrafish as well as in the present study highlights the importance of this stage and pinpoints an enhanced transcriptional activity. Besides differential expression of paralogous genes and/or developmentally specific genes, the importance of miRNAs during development has been shown in several studies. Therefore, a first approach to describe new miRNAs for the European seabass was performed in the present study. This is to our knowledge the first report about miRNAs in the early development of the European seabass. miRNAs are small noncoding RNA regulating gene expression mainly in a suppressive way at the posttranscriptional level (33). In teleosts, 1,044 of mature miRNAs have been identified (reviewed in Ref. 6), and it has been shown that they play an important role in several biological process including development. The developmental profiles of miRNAs have been reported in zebrafish (64), fugu Tetraodon (6), medaka (52), Asian seabass (92), Atlantic cod (46), Japanese flounder (23), and Atlantic halibut (7). The identification of novel miRNAs, thus not yet described in any other species, has not been aimed in the present work, instead miRNAs characterized in other species but not yet described in the European seabass were detected. In total, 42 miRNAs have

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

166

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

Fig. 6. Network analysis. Representation of the network inferred by the LeMoNe algorithm showing the interaction between differential expressed transcripts during developmental stages and the 14 of the 42 identified miRs of European seabass. Clusters of coexpressed transcripts are shown as black circles, and miRNAs are shown as red triangles.

not critical for the earliest developmental processes (88, 89). The low expression of miRNA found in the earlier stages in the present study may also suggest that miRNAs are not critical for the first developmental processes in European seabass. It has to

Transcript ID

Blast match

miR21

comp54330_c1_seq2 H+ V0 subunit b

XP_003438118

-29.6*

miR30a

comp58836_c0_seq19 major vault protein

XP_003979245

-25.6*

comp53312_c0_seq3 delta-like protein B-like

XP_003449875

-30.6*

comp59023_c0_seq23 protein strawberry notch homolog 1

XP_003455166

-26.2*

comp57660_c1_seq1 fukutin-like

XP_003965670

-25.9

comp57361_c0_seq2 transcription factor Sp8-like

XP_003454480

-25.7

comp55922__c0_seq3 thyroid hormone receptor associated protein complex component TRAP230/KIAA0192

CAC44632

-25.5

miR455b

RNAhybrid

mfe (kcal/mol)

miR

miR221-2

been described in the present work, and their stage-specific expression is illustrated in Fig. 5. In the present work clearly most of the miRNAs are highly expressed at the late organogenesis stage. It has been shown that in zebrafish miRNAs are

Fig. 7. Transcripts assigned to the GO term developmental process were tested for hybridization with corresponding miRNA. *Positive hybridization with both RNAhybrid software (minimum free energy ⱕ ⫺25.0) and miRanda software. Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

be noted that of the 42 identified miRNAs, some of the classical miRNAs like lin-4, let-7, lsy-6, and miR-273 (89) have not been detected. The two miRNA lin-4 and let-7 are the best studied ones and regulate developmental timing in Caenorhabditis elegans (49, 90), whereas let-7 and lsy-6 are known to regulate left/right asymmetry in worms (11). A possible explanation could be that those miRNAs are highly expressed even earlier in development. On the other hand, other classical miRNAs involved in important development processes have been found, such as miR-196. miR-196 is involved in developmental patterning in mouse, and its target has been identified to be Hoxb8 (96). Other important miRNAs identified later on in various teleost fish species are also found to be differentially expressed in the present work. Ramachandra et al. (67) for instance described the function of miR-21, 23a, 26a, 200b, and 455 during the early embryo, whereas for example miR-1 has been shown in zebrafish (77), common carp (Cyprinus carpio) (95), and blunt snout bream (Megalobrama amblycephala) to play a distinct role in muscles. Furthermore, Yan et al. (94) showed the involvement of miR-203b in the muscle development of the Nile tilapia (Astatotilapia burtoni). The remainder of the discussion will be focused on the interaction of the identified miRNA with the detected differentially expressed genes (Supplemental File 4). We achieved this by performing network analysis based on the expression data. Network analysis is mostly performed in human cancer research to reveal regulators, mainly miRNAs, responsible for cancer outbreak. Molecular networks comprise nodes and modules, where the nodes are in the present study miRNAs and the modules coexpressed transcripts. The module network is a probabilistic model based on probabilistic graphical models and Bayesian networks and aims to identify regulatory modules from gene expression data. The procedure identifies modules of co-regulated genes and their regulators. Here in total 14 miRNAs were identified as nodes linked to modules based on expression patterns using the LeMoNe software package (Fig. 6). One of the nodes presents miR-1 regulating modules, which were revealed to regulate transcripts involved in the proper function of muscle as well as muscle development (12, 24, 95). As network analysis does not distinguish whether regulators repress or enhance possible transcripts it has to be taken into account that not only direct interaction of miRNAs is illustrated here but also indirect involvement. Indirect involvement for example is seen in the case of miR-1, where transcripts clustered in the modules are upregulated in concordance with the expression of miR-1 (expression data not shown). Direct regulation was found for 22 modules linked to seven miRNAs. The final step to detect and evaluate interaction between miRNAs and mRNAs of the present study was to allocate possible targets of the miRNAs. Therefore, 3=-UTRs of the sequences belonging to the clusters found to be regulated by one of the miRNAs were searched by two software packages for putative target sites. Figure 7 shows the most significant ones comprising miRNA, target site, as well as structure and mfe values. All miRNAs shown in Fig. 7 have been described to be involved in developmental processes (89). One of the miRNAs, miR-30a, is linked to several modules, including one module with putative direct regulation of miR-196a. Both miRNAs have been previously described to be involved in development with the former one in early muscle development (42) and the latter one in developmental

167

patterning (96). In the present study, two putative target sites in the major vault protein (mvp) and in the tata box binding protein 1 (tbp1) (mfe ⫽ ⫺25.6 and mfe ⫽ ⫺22.7, respectively) were identified for miR-30a (Fig. 7 and Supplemental File 2). The former has been shown to be induced during regeneration (98) and is considered to belong to the highly conserved ribonucleoprotein particles. It is suggested that mvp might be involved in the regulation of important cell signaling pathways including the PI3K/Akt and the MAPK pathways (78). Interestingly, miR-30a is also known to be required for hepatobiliary development (31), and MAPK, as well as Akt, is activated in extrahepatic biliary tract carcinoma (35), a malfunction of the extrahepatic biliary tract. Nevertheless, further experimental studies have to be carried out to confirm the predicted target sites of the present work. Conclusion This study reports the characterization of the European seabass (D. labrax) transcriptome during early development (from morula to hatching). mRNAs comprising paralogous genes and differentially spliced transcripts are described, and their expression patterns during development assessed. In addition, we identified 42 new miRNA, as well as their expression during the seven early developmental stages. Putative direct and indirect regulatory function of the miRNA studied in the present work were determined by network analysis identifying putative miRNA/mRNA interactions as well as putative target sites that have not been described before. ACKNOWLEDGMENTS We thank the Department of Aquaculture for providing fish eggs and assistance in sampling and the Informatics group for computational support. We also acknowledge the Génopole-Plateforme bio-informatique Centre INRA de Toulouse-Unité de BIA for providing Venn diagram calculation for six groups after request. GRANTS 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), cofunded by the EU and the Hellenic Republic through the European Social Fund. DISCLOSURES No conflicts of interest, financial or otherwise, are declared by the author(s). AUTHOR CONTRIBUTIONS Author contributions: E.K. and J.X. performed experiments; E.K. and E.S. analyzed data; E.K. and E.S. prepared figures; E.K. and E.S. drafted manuscript; E.K., J.X., E.A., C.S.T., and E.S. approved final version of manuscript; E.A., C.S.T., and E.S. edited and revised manuscript; E.S. conception and design of research; E.S. interpreted results of experiments. REFERENCES 1. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol 215: 403–410, 1990. 2. Alvarez-Garcia I, Miska EA. MicroRNA functions in animal development and human disease. Development (Cambridge, UK) 132: 4653–4662, 2005. 3. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol 11: R106, 2010. 4. Bardou P, Mariette J, Escudié F, Djemiel C, Klopp C. jvenn: an interactive Venn diagram viewer. BMC Bioinformat 15: 293, 2014. 5. Bekaert M, Lowe NR, Bishop SC, Bron JE, Taggart JB, Houston RD. Sequencing and characterisation of an extensive Atlantic salmon (Salmo salar L.) microRNA repertoire. PLoS One 8: e70136, 2013.

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

168

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax

6. Bizuayehu TT, Babiak I. MicroRNA in teleost fish. Genome Biol Evol 6: 1911–1937, 2014. 7. Bizuayehu TT, Lanes CF, Furmanek T, Karlsen BO, Fernandes JM, Johansen SD, Babiak I. Differential expression patterns of conserved miRNAs and isomiRs during Atlantic halibut development. BMC Genomics 13: 11, 2012. 8. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30: 2114 –2120, 2014. 9. Boutet I, Long Ky CL, Bonhomme F. A transcriptomic approach of salinity response in the euryhaline teleost, Dicentrarchus labrax. Gene 379: 40 –50, 2006. 10. Calduch-Giner JA, Sitja-Bobadilla A, Davey GC, Cairns MT, Kaushik S, Perez-Sanchez J. Dietary vegetable oils do not alter the intestine transcriptome of gilthead sea bream (Sparus aurata), but modulate the transcriptomic response to infection with Enteromyxum leei. BMC Genomics 13: 470, 2013. 11. Chang S, Johnston RJ, Frøkjaer-Jensen C, Lockery S, Hobert O. MicroRNAs act sequentially and asymmetrically to control chemosensory laterality in the nematode. Nature 430: 785–789, 2004. 12. Chen PY, Manninga H, Slanchev K, Chien M, Russo JJ, Ju J, Sheridan R, John B, Marks DS, Gaidatzis D, Sander C, Zavolan M, Tuschl T. The developmental miRNA profiles of zebrafish as determined by small RNA cloning. Genes Dev 19: 1288 –1293, 2005. 13. Chistiakov DA, Hellemans B, Tsigenopoulos CS, Law AS, Bartley N, Bertotto D, Libertini A, Kotoulas G, Haley CS, Volckaert FAM. Development and linkage relationships for new microsatellite markers of the sea bass (Dicentrarchus labrax L.). Anim Genet 35: 53–57, 2004. 13a.Collignon J, Varlet I, Robertson EJ. Relationship between asymmetric nodal expression and the direction of embryonic turning. Nature 381: 155–158, 1996. 14. Conesa A, Götz S. Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int J Plant Genom 2008: 619832, 2008. 15. Cucchi P, Sucré E, Santos R, Leclère J, Charmantier G, Castille R. Embryonic development of the sea bass Dicentrarchus labrax. Helgoland Mar Res 66: 199 –209, 2012. 16. Darias MJ, Zambonino-Infante JL, Hugot K, Cahu CL, Mazurais D. Gene expression patterns during the larval development of European sea bass (Dicentrarchus labrax) by microarray analysis. Mar Biol 10: 416 – 428, 2008. 17. Divanach P. Contribution to biology and rearing of six Mediterranean Sparids: Sparus aurata, Diplodus sargus, Diplodus vulgaris, Diplodus annularis, Lithognathus mormyrus, Puntazzo puntazzo (Teleostean fishes) (PhD thesis). Montpellier II, France, University of Sciences & Techniques of Languedoc, 1985. 18. Dooley K, Zon L. Zebrafish: a model system for the study of human disease. Curr Opin Genet Dev 10: 252–256, 2000. 19. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol 5: R1, 2003. 20. Fernandes AM, Beddows E, Filippi A, Driever W. Orthopedia transcription factor otpa and otpb paralogous genes function during dopaminergic and neuroendocrine cell specification in larval zebrafish. PLoS One 8: e75002, 2013. 21. Fjose A, Zhao XF. Exploring microRNA functions in zebrafish. New Biotechnology 27: 250 –255, 2010. 22. Force A, Lynch M, Pickett FB, Amores A, Yan YL, Postlethwait J. Preservation of duplicate genes by complementary, degenerative mutations. Genetics 151: 1531–1545, 1999. 23. Fu Y, Shi Z, Wu M, Zhang J, Jia L, Chen X. Identification and differential expression of microRNAs during metamorphosis of the Japanese flounder (Paralichthys olivaceus). PLoS One 6: e22957, 2011. 24. Gagan J, Dey BK, Dutta A. MicroRNAs regulate and provide robustness to the myogenic transcriptional network. Curr Opin Pharmacol 12: 383– 388, 2012. 25. Garcia de la Serrana D, Estevez A, Andree K, Johnston IA. Fast skeletal muscle transcriptome of the gilthead sea bream (Sparus aurata) determined by next generation sequencing. BMC Genomics 13: 181, 2012. 26. Garcia de la Serrana D, Mareco EA, Johnston IA. Systematic variation in the pattern of gene paralog retention between the teleost superorders ostariophysi and acanthopterygii. Genome Biol Evol 6: 981–987, 2014. 27. Geay F, Ferraresso S, Zambonino-Infante JL, Bargelloni L, Quentel C, Vandeputte M, Kaushik S, Cahu CL, Mazurais D. Effects of the total replacement of fish meal and fish oil by plant sources on the hepatic transcriptome of European sea bass (Dicentrarchus labrax). Aquaculture Europe 11. Mediterranean Aquaculture 2020: 377–378, 2011.

28. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol 29: 644 –652, 2011. 29. Guyon R, Senger F, Rakotomanga M, Sadequi N, Volckaert FAM, Hitte C, Galibert F. 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, 2010. 30. Haas BJ, Zody MC. Advancing RNA-Seq analysis. Nat Biotechnol 28: 421–423, 2010. 31. Hand NJ, Master ZR, Eauclaire SF, Weinblatt DE, Matthews RP, Friedman JR. The microRNA-30 family is required for vertebrate hepatobiliary development. Gastroenterology 136: 1081–1090, 2009. 32. Hartigan JA, Wong MA. Algorithm AS 136: a K-means clustering algorithm. J Roy Stat Soc C 28: 100 –108, 1979. 33. He L, Hannon GJ. MicroRNAs: small RNAs with a big role in gene regulation. Nat Rev Genet 5: 522–531, 2004. 35. Hori H, Ajiki T, Mita Y, Horiuchi H, Hirata K, Matsumoto T, Morimoto H, Fujita T, Ku Y, Kuroda Y. Frequent activation of mitogen-activated protein kinase relative to Akt in extrahepatic biliary tract cancer. J Gastroenterol 42: 567–572, 2007. 36. Iwamatsu T. Stages of normal development in the medaka Oryzias latipes. Mech Dev 121: 605–618, 2004. 38. Jansson MD, Lund AH. MicroRNA and cancer. Mol Oncol 6: 590 –610, 2012. 39. Jin Z, Xie T. Germline specification: small things have a big role. Curr Biol 16: R966 –R967, 2006. 40. Kassahn KS, Dang VT, Wilkins SJ, Perkins AC, Ragan a M. Evolution of gene function and regulatory control after whole-genome duplication: comparative analyses in vertebrates. Genome Res 19: 1404 –1418, 2009. 41. Kent WJ. BLAT-The BLAST-like alignment tool. Genome Res 12: 656 –664, 2002. 42. Ketley A, Warren A, Holmes E, Gering M, Aboobaker AA, Brook JD. The miR-30 microRNA family targets smoothened to regulate hedgehog signalling in zebrafish early muscle development. PLoS One 8: e65170, 2013. 43. Kimmel CB, Ballard WW, Kimmel SR, Ullmann B, Schilling TF. Stages of embryonic development of the zebrafish. Dev Dyn 203: 253– 310, 1995. 44. Kitano J, Yoshida K, Suzuki Y. RNA sequencing reveals small RNAs differentially expressed between incipient Japanese threespine sticklebacks. BMC Genomics 14: 214, 2013. 45. Krüger J, Rehmsmeier M. RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res 34: W451–W454, 2006. 46. Lanes CFC, Bizuayehu TT, de Oliveira Fernandes JM, Kiron V, Babiak I. Transcriptome of Atlantic cod (Gadus morhua L.) early embryos from farmed and wild broodstocks. Mar Biotechnol 15: 677–694, 2013. 47. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Meth 9: 357–359, 2012. 48. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memoryefficient alignment of short DNA sequences to the human genome. Genome Biol 10: R25, 2009. 49. Lee RC, Feinbaum RL, Ambros V. The C elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell 75: 843–854, 1993. 50. Li B, Dewey CN. RSEM: accurate transcript quantification from RNASeq data with or without a reference genome. BMC Bioinformatics 12: 323, 2011. 51. Li P, Peatman E, Wang S, Feng J, He C, Baoprasertkul P, Xu P, Kucuktas H, Nandi S, Somridhivej B, Serapion J, Simmons M, Turan C, Liu L, Muir W, Dunham R, Brady Y, Grizzle J, Liu Z. Towards the ictalurid catfish transcriptome: generation and analysis of 31,215 catfish ESTs. BMC Genomics 8: 177, 2007. 52. Li SC, Chan WC, Ho MR, Tsai KW, Hu LY, Lai CH, Hsu CN, Hwang PP, Lin W. Discovery and characterization of medaka miRNA genes by next generation sequencing platform. BMC Genomics 11, Suppl 4: S8, 2010. 53. Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22: 1658 – 1659, 2006.

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org

DYNAMISM IN GENE EXPRESSION OF Dicentrarchus labrax 54. Liu L, Li Y, Li S, Hu N, He Y, Pong R, Lin D, Lu L, Law M. Comparison of next-generation sequencing systems. J Biomed Biotechnol 2012: 251364, 2012. 55. Liu S, Zhang Y, Zhou Z, Waldbieser G, Sun F, Lu J, Zhang J, Jiang Y, Zhang H, Wang X, Rajendran K, Khoo L, Kucuktas H, Peatman E, Liu Z. Efficient assembly and annotation of the transcriptome of catfish by RNA-Seq analysis of a doubled haploid homozygote. BMC Genomics 13: 595, 2012. 57. Luo GZ, Hafner M, Shi Z, Brown M, Feng GH, Tuschl T, Wang XJ, Li X. Genome-wide annotation and analysis of zebra finch microRNA repertoire reveal sex-biased expression. BMC Genomics 13: 727, 2012. 58. McGettigan PA. Transcriptomics in the RNA-seq era. Curr Opin Chem Biol 17: 4 –11, 2013. 59. Mercier P, Simeone A, Cotelli F, Boncinelli E. Expression pattern of two otx genes suggests a role in specifying anterior body structures in zebrafish. Int J Dev Biol 39: 559 –573, 1995. 60. Michibata H, Nojima Y, Kojima MK. Stage sensitivity of eggs of the teleost Oryzias latipes to cadmium exposure. Environ Res 42: 321–327, 1987. 61. Michoel T, Maere S, Bonnet E, Joshi A, Saeys Y, Van den Bulcke T, Van Leemput K, van Remortel P, Kuiper M, Marchal K, Van de Peer Y. Validating module network learning algorithms using simulated data. BMC Bioinformat 8, Suppl 2: S5, 2007. 62. Mishra NS, Mukherjee SK. A peep into the plant miRNA world. Open Plant Sci J 1: 1–9, 2007. 63. Oshlack A, Robinson MD, Young MD. From RNA-seq reads to differential expression results. Genome Biol 11: 220, 2010. 64. Pauli A, Valen E, Lin MF, Garber M, Vastenhouw NL, Levin JZ, Fan L, Sandelin A, Rinn JL, Regev A, Schier AF. Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome Res 22: 577–591, 2012. 65. Pittman K, Yúfera M, Pavlidis M, Geffen AJ, Koven W, Ribeiro L, Zambonino-Infante JL, Tandler A. Fantastically plastic: fish larvae equipped for a new world. Rev Aquacult 5, Suppl S1: S224 –S267, 2013. 66. Powell GT, Wright GJ. Jamb and jamc are essential for vertebrate myocyte fusion. PLoS Biol 9: e1001216, 2011. 67. Ramachandra RK, Salem M, Gahr S, Rexroad CE, Yao J. Cloning and characterization of microRNAs from rainbow trout (Oncorhynchus mykiss): their expression during early embryonic development. BMC Dev Biol 8: 41, 2008. 68. Rapaport F, Khanin R, Liang Y, Pirun M, Krek A, Zumbo P, Mason CE, Socci ND, Betel D. Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data. Genome Biol 14: R95, 2013. 69. Salem M, Xiao C, Womack J, Rexroad CE, Yao J. A microRNA repertoire for functional genome research in rainbow trout (Oncorhynchus mykiss). Mar Biotechnol 12: 410 –429, 2010. 70. Sarropoulou E, Fernandes JMO, Mitter K, Magoulas A, Mulero V, Sepulcre MP, Figueras A, Novoa B, Kotoulas G. Evolution of a multifunctional gene: the warm temperature acclimation protein Wap65 in the European seabass Dicentrarchus labrax. Mol Phylogenet Evol 55: 640 –649, 2010. 71. Sarropoulou E, Galindo-Villegas J, Garcia-Alcazar A, Kasapidis P, Mulero V, García-Alcázar A, Kasapidis P, Mulero V, Garcia-Alcazar A. Characterization of European sea bass transcripts by RNA SEQ after oral vaccine against V. anguillarum. Mar Biotechnol (New York, NY) 14: 634 –642, 2012. 72. Sarropoulou E, Kotoulas G, Power DM, Geisler R. Gene expression profiling of gilthead sea bream during early development and detection of stress-related genes by the application of cDNA microarray technology. Physiol Genomics 23: 182–191, 2005. 73. Sarropoulou E, Moghadam HK, Papandroulakis N, De La Gándara F, Garcia AO, Makridis P. The Atlantic Bonito (Sarda sarda, Bloch 1793) transcriptome and detection of differential expression during larvae development. PLoS One 9: e87744, 2014. 74. Sarropoulou E, Nousdili D, Kotoulas G, Magoulas A. Functional divergences of GAPDH isoforms during early development in two perciform fish species. Mar Biotechnol 13: 1115–1124, 2011. 75. Sarropoulou E, Nousdili D, Magoulas A, Kotoulas G. Linking the genomes of nonmodel teleosts through comparative genomics. Mar Biotechnol (New York, NY) 10: 227–233, 2008. 76. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment

169

for integrated models of biomolecular interaction networks. Genome Res 13: 2498 –2504, 2003. 77. Soares AR, Pereira PM, Santos B, Egas C, Gomes AC, Arrais J, Oliveira JL, Moura GR, Santos MAS. Parallel DNA pyrosequencing unveils new zebrafish microRNAs. BMC Genomics 10: 195, 2009. 78. Steiner E, Holzmann K, Elbling L, Micksche M, Berger W. Cellular functions of vaults and their involvement in multidrug resistance. Curr Drug Targets 7: 923–934, 2006. 79. Stern-Ginossar N, Elefant N, Zimmermann A, Wolf DG, Saleh N, Biton M, Horwitz E, Prokocimer Z, Prichard M, Hahn G, GoldmanWohl D, Greenfield C, Yagel S, Hengel H, Altuvia Y, Margalit H, Mandelboim O. Host immune system gene targeting by a viral miRNA. Science (New York, NY) 317: 376 –381, 2007. 80. Taylor JS, Braasch I, Frickey T, Meyer A, Van de Peer Y. Genome duplication, a trait shared by 22,000 species of ray-finned fish. Genome Res 13: 382–390, 2003. 81. Tine M, Kuhl H, Gagnaire PA, Louro B, Desmarais E, Martins RS, Hecht J, Knaust F, Belkhir K, Klages S, Dieterich R, Stueber K, Piferrer F, Guinand B, Bierne N, Volckaert FA, Bargelloni L, Power DM, Bonhomme F, Canario AV, Reinhardt R. European sea bass genome and its variation provide insights into adaptation to euryhalinity and speciation. Nat Commun 5: 5770, 2014. 82. Ton C, Hwang DM, Dempsey AA, Tang HC, Yoon J, Lim M, Mably JD, Fishman MC, Liew CC. Identification, characterization, and mapping of expressed sequence tags from an embryonic zebrafish heart cDNA library. Genome Res 10: 1915–1927, 2000. 83. Vesterlund L, Jiao H, Unneberg P, Hovatta O, Kere J. The zebrafish transcriptome during early development. BMC Dev Biol 11: 30, 2011. 83a.von Westernhagen H. Sublethal effects of pollutants on fish eggs and larvae. In: Fish Physiology (vol. XI, part A), edited by Hoar WS, Randall DJ. San Diego, CA: Academic, 1988, p. 254 –258. 84. Wagner A. Does evolutionary plasticity evolve? Evolution 50: 1008 – 1023, 1996. 85. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 10: 57–63, 2009. 86. Wienholds E, Kloosterman WP, Miska E, Alvarez-Saavedra E, Berezikov E, de Bruijn E, Horvitz HR, Kauppinen S, Plasterk RHA. MicroRNA expression in zebrafish embryonic development. Science 309: 310 –311, 2005. 88. Wienholds E, Koudijs MJ, van Eeden FJM, Cuppen E, Plasterk RHA. The microRNA-producing enzyme Dicer1 is essential for zebrafish development. Nat Genet 35: 217–218, 2003. 89. Wienholds E, Plasterk RHA. MicroRNA function in animal development. FEBS Lett 579: 5911–5922, 2005. 90. Wightman B, Bürglin TR, Gatto J, Arasu P, Ruvkun G. Negative regulatory sequences in the lin-14 3=-untranslated region are necessary to generate a temporal switch during Caenorhabditis elegans development. Genes Dev 5: 1813–1824, 1991. 91. Wong L, Weadick CJ, Kuo C, Chang BSW, Tropepe V. Duplicate dmbx1 genes regulate progenitor cell cycle and differentiation during zebrafish midbrain and retinal development. BMC Dev Biol 10: 100, 2010. 92. Xia JH, He XP, Bai ZY, Yue GH. Identification and characterization of 63 microRNAs in the Asian seabass Lates calcarifer. PLoS One 6: e17537, 2011. 93. Yan B, Guo JT, Zhao LH, Zhao JL. microRNA expression signature in skeletal muscle of Nile tilapia. Aquaculture 364 –365: 240 –246, 2012. 94. Yan B, Guo JT, Zhu CD, Zhao LH, Zhao JL. miR-203b: a novel regulator of MyoD expression in tilapia skeletal muscle. J Exp Biol 216: 447–4451, 2013. 95. Yan X, Ding L, Li Y, Zhang X, Liang Y, Sun X, Teng CB. Identification and profiling of microRNAs from skeletal muscle of the common carp. PLoS One 7: e30925, 2012. 96. Yekta S, Shih IH, Bartel DP. MicroRNA-directed cleavage of HOXB8 mRNA. Science (New York, NY) 304: 594 –596, 2004. 97. Ying SY, Lin SL. MicroRNA: fine-tunes the function of genes in zebrafish. Biochem Biophys Res Commun 335: 1–4, 2005. 98. Yoshinari N, Ishida T, Kudo A, Kawakami A. Gene expression and functional analysis of zebrafish larval fin fold regeneration. Dev Biol 325: 71–81, 2009.

Physiol Genomics • doi:10.1152/physiolgenomics.00001.2015 • www.physiolgenomics.org