mRNA/microRNA Profile at the Metamorphic Stage of Olive Flounder (Paralichthys olivaceus)

Flatfish is famous for the asymmetric transformation during metamorphosis. The molecular mechanism behind the asymmetric development has been speculated over a century and is still not well understood. To date, none of the metamorphosis-related genes has been identified in flatfish. As the first step to screen metamorphosis-related gene, we constructed a whole-body cDNA library and a whole-body miRNA library in this study and identified 1051 unique ESTs, 23 unique miRNAs, and 4 snoRNAs in premetamorphosing and prometamorphosing Paralichthys olivaceus. 1005 of the ESTs were novel, suggesting that there was a special gene expression profile at metamorphic stage. Four miRNAs (pol-miR-20c, pol-miR-23c, pol-miR-130d, and pol-miR-181e) were novel to P. olivaceus; they were characterized as highly preserved homologies of published miRNAs but with at least one nucleotide differed. Representative 24 mRNAs and 23 miRNAs were quantified during metamorphosis of P. olivaceus by using quantitative RT PCR or stem-loop qRT PCR. Our results showed that 20 of mRNAs might be associated with early metamorphic events, 10 of mRNAs might be related with later metamorphic events, and 16 of miRNAs might be involved in the regulation of metamorphosis. The data provided in this study would be helpful for further identifying metamorphosis-related gene in P. olivaceus.


Introduction
Flatfish is famous for the asymmetric transformation during metamorphosis, especially one eye migrating to the other side. Other metamorphosis events include cranium deformation, asymmetric pigmentation, and 90-degree rotation in posture with a lifestyle transition from pelagic to benthic. The molecular mechanism of morphologic left/right asymmetry in Olive flounder, Paralichthys olivaceus, was thought to be different from that of interior organ asymmetry in vertebrate [1]. Thyroid hormone (TH) was proposed to regulate metamorphosis in flatfish [2][3][4][5]. As the nuclear receptor of TH, thyroid hormone receptor (TR) should be involved in the TH-inducing signal pathway. The spatial expression of the TR genes has been investigated in the metamorphosing Olive flounder [6]; however, it still cannot determine which metamorphosis events were regulated by TR in flatfish. In TH-TR signal pathway, the downstream genes will unavoidably be investigated in the future. To date, very few genes were investigated in metamorphosing flatfish [6][7][8][9]. Even though the cDNA libraries of various tissue types in P. olivaceus were constructed, especially immune-related tissues [10], the gene expression profile in metamorphosing P. olivaceus was still unavailable.
Expressed sequence tags (ESTs) analysis is an efficient approach to characterize transcriptome. Large-scale EST sequencing project as a part of genome project has been conducted for several teleost species, such as salmonid and catfish [11][12][13]. Small-scale ESTs analysis has also been carried out for some aquaculture teleosts [10,[14][15][16]. In this study, we tried to enrich ESTs data and investigated the gene expression profile via cDNA library random sequencing in the premetamorphosing and prometamorphosing P. olivaceus. In addition, Johnston and Hobert reported that one microRNA termed lsy-6 controlled neuronal left/right asymmetric expression of chemosensory receptor in Caenorhabditis elegans [17]. Accordingly, the possible roles of microRNAs in regulating metamorphosis in flatfish should not be neglected. This is the reason that we constructed a microRNA library and analyzed its expression profile in the metamorphosing P. olivaceus in this research as well.

Fish Maintenance and Sampling.
Larvae were obtained from the Central Experiment Station of Chinese Academy of Fisheries Sciences (Beidaihe, Hebei, China) and then transported to the laboratory in Shanghai Ocean University, Shanghai, China. The larvae were reared in the laboratory according to the methods provided in [8]. Larvae were fed live brine shrimp (Artemia) nauplii until the end of metamorphosis. We use the following classifications for the metamorphic stages of P. olivaceus in this study [18,19]: Premetamorphosis (17 DAH, days after hatching), the stage prior to the start of eye migration; Prometamorphosis (19 DAH), from the start of eye migration until the start of resorption of several elongated dorsal fin rays; Climax (23 DAH), from the start of resorption of the elongated dorsal fin rays until the completion of fin resorption and eye migration; Postclimax (27 DAH), after the completion of fin resorption and eye migration. All samples were frozen using liquid nitrogen and stored at −80 • C until proceeding to total RNA isolation.  [20] to be compared against the nonredundant protein database using BLASTX. The E-value cutoff was 1e−5. Novel ESTs were also identified by comparison with P. olivaceus EST sequences in dbEST at NCBI using BLASTN. All ESTs that were not identified as orthologs of known genes were designated as unknown EST clones. Sequences with BLASTX hits were mapped and annotated according to gene ontology terms (GO) in AmiGO database (http://amigo.geneontology.org/cgi-bin/amigo/go. cgi). The distribution of genes in each of the main ontology categories was examined, and the percentages of unique sequences in each of the assigned GO terms were calculated. In each of the three main categories of GO, namely, biological process, molecular function, and cellular component [21], 100% was considered as the total number of unique sequences having an assigned GO term. Thus, in each main category, the percentages of 2nd level do not add up to 100% because some deduced proteins have more than one GO category assigned to them [22].

MicroRNA Library Construction and Sequencing.
RNA with size less than 200 nt from premetamorphosing or metamorphosing larvae was isolated using mirVana miRNA Isolation Kit (Ambion, Austin, Tex, USA) following manufacturer's instructions with minor modifications. In brief,    Since a gene product could be assigned to more than one GO term, the percentages in each main category do not add up to 100%.

Quantitative Real-Time Reverse Transcription Polymerase
Chain Reaction. Total RNA from the whole larvae at different metamorphic stages was isolated using TRIzol reagent (Invitrogen) followed by DNase treatment. The abundance of mRNA or miRNA was quantified by qRT-PCR or stemloop qRT-PCR [24], respectively. 2 μg of DNase-treated RNA was converted to cDNA using MMLV reverse transcriptase (Promega). qRT-PCR primers for each mRNA are listed in Table 1 and for miRNA in Table 2. The relative expression of mRNA and miRNA was normalized using β-actin mRNA and U6 snRNA as control, respectively. For qRT-PCR or stem-loop qRT-PCR, the thermocycler was set at 95 • C for 15 s and 55 • C for 60 s per cycle for a total of 40 cycles, followed by 95 • C for 1 min and 55 • C for 1 min. Relative changes in mRNA or miRNA abundances were quantified by using the C t method; β-actin mRNA and U6 RNA were used as reference amplicons for data normalization [25,26].  (Figure 1). Gene ontology (GO) categories were assigned to 656 unique ESTs using AmiGO database. The percentage distributions of gene ontology terms (2nd-level GO terms) according to the GO consortium are shown in Figure 2. Cellular Process (87%) was the most dominant 2nd-level term out of the 455 unique sequences which were annotated to the Biological Process GO category. This was followed by Metabolic Process Metabolism at 71%. It is noted that 9% were assigned to the Negative Regulation of Biological Process. Protein Binding (60%) was the most dominant out of 467 ESTs with significant protein hits which were assigned to Molecular Function category at 2nd level. This was followed by Nucleic Acid Binding at 20%. Cell Part (97%) was the most dominant out of 474 ESTs which were annotated to the Cellular Component GO category. Intracellular and intracellular parts occupied 90% and 89%, respectively. ESTs that fell in each of the three main GO categories are given in Figure 2.

Results and Discussions
Compared with normalized cDNA library, the nonnormalized cDNA library is much more redundant. 2,193 cDNA clones from the nonnormalized cDNA library in this study only generated 1,051 unique genes. However, the nonnormalized cDNA library can provide raw information on the structure of gene expression level [27]. Among 656 identified distinct known genes in metamorphic P. olivaceus in this study, 413 known genes (63.0%) were sequenced only once, 180 genes (27.4%) were sequenced 2-5 times, and 63 genes (9.6%) were sequenced over 5 times. The vast majority of known genes were sequenced only once; however, a small number of genes accounted for a large proportion of transcripts in premetamorphosing and prometamorphosing P. olivaceus (Figure 3). The most abundantly expressed gene was parvalbumin accounting for 3.88% of the 2,193 clones sequenced ( Table 3). The expressed gene beta-actin accounted for only 0.05%. The other most abundant expressed genes included cytochrome c oxidase subunit II (1.28%), ribosomal protein S2 (1.23%), cytochrome c oxidase subunit III (1.00%), creatine kinase 1 (1.00%), myosin light chain 3 (1.00%), 40S ribosomal protein S8 (1.00%), nuclease diphosphate kinase B (0.87%), ribosomal protein L18a (0.87%), and antifreeze protein type IV (0.87%). Altogether, the ten most abundantly expressed genes occupied 19.39% of all clones.

miRNA Profile in Premetamorphosing and
Prometamorphosing P. olivaceus . MicroRNAs are small 19-23-nucleotide noncoding RNAs that bind to recognition sequences on 3untranslated regions (3 -UTRs) of mRNAs and target them   Figure 4: Alignment of the novel P. olivaceus miRNAs with highly preserved homologous miRNAs from other species.
for degradation or translational repression. MiRNAs have been found to play important roles in zebrafish development [28,29]. miRNAs resources were developed only in few teleosts such as zebrafish, puffer fish, and Oncorhynchus mykiss [30][31][32]. No miRNAs have been identified in flatfish. In this study, total 143 clones picked randomly were sequenced (Table 4). Sequence analysis identified 29 microR-NAs that showed the same as at least one published miR-NAs in the database (http://www.mirbase.org/search.shtml). Table 5. Four sequences had not been found to have the same sequences, but they showed significant similarities with published miRNAs in miRBase. In addition, there are four sequences identified as snoRNA by searching NCBI database. Overall, 23.08% of small RNAs in the library might be microRNAs and 2.80% were snoRNAs ( Table 4). Names of the P. olivaceus miRNA were assigned based on the homologies between the cloned sequence and published miRNA sequences (Table 5). 19 unique miRNAs are conserved across several species and affiliated to 15 subfamilies, of which 9 unique miRNA (pol-let-7a, pol-miR-7f, pol-miR-26a, pol-miR-125b, pol-miR-128, pol-miR-181a, pol-miR-200a, pol-miR-221, and pol-miR-429) are conserved higher across ten or more species.

Representing 19 unique miRNAs are shown in
The pol-miR-125b is conserved across 43 species. While mirRNAs, pol-let-7j, pol-miR-21a, pol-miR-181f, and pol-miR-724, are conserved across only one species (Table 5). Four miRNAs (pol-miR-20c, pol-miR-23c, pol-miR-130d, and pol-miR-181e) are novel to P. olivaceus characterized as having high homologies with published miRNAs but differed by at least one nucleotide. These 4 miRNAs only observed in P. olivaceus are of special interest because of their unique sequences and possibly unique targeting mechanisms ( Figure 4). Pol-miR-20c has a U to G mismatch with miR-20 of Fugu rubripes, Tetraodon nigroviridis, and Monodelphis domestica or miR-20a of Danio rerio, Xenopus tropicalis, Gallus gallus, Equus caballus, Canis familiaris, and Homo sapiens (Figure 4(a)). Pol-miR-23c has a U to C mismatch at positions 23 with miR-23b of Bos taurus, Pongo pygmaeus, Pan paniscus, and Pan troglodytes. However, the position 23 is absent between P. olivaceus and other fishes. In the miR-23b of E. caballus and H. sapiens, the positions 23 and 22 are both absent (Figure 4(b)). Pol-miR-130d has an A to G  (Figure 4(c)). Compared with other species, pol-miR-181e has a G to U mismatch at the position 19. It is interesting that there is a one base absent at position 23 in other fishes and Xenopus laevis, whereas the position of U is conserved in P. olivaceus and higher vertebrates (Figure 4(d)).

The Expression Pattern of Representative Genes in
Metamorphosing P. Olivaceus . Among miRNAs-targeted sequences in silico predicted by RNA22 miRNA target detection software [33], expression 24 genes of was confirmed in premetamorphosing or metamorphosing flounders by qRT-PCR. Only gene atcay was expressed stably at different metamorphic stages, indicating that it should not be associated with metamorphic events (Figure 5(a)).

Conclusion
In summary, we generated a collection of 1,051 unique ESTs, 23 unique miRNAs, and 4 unique snoRNAs in premetamorphosing and prometamorphosing P. olivaceus. Even though so far there were 3143 nucleotides and 13869 ESTs available in NCBI database, 1005 novel ESTs were identified successfully in this study, suggesting that special gene expression profile existed in metamorphic stage. Representative 24 mRNAs of 1051 unique ESTs were quantified during the metamorphosis of P. olivaceus using quantitative RT PCR, and the results showed that 20 genes might be associated with early metamorphic events and 10 genes might be related with later metamorphic events. In addition, the abundances of 23 miRNAs were quantified using stem-loop qRT PCR. 9 miRNAs might be associated with metamorphosis, and 7 miRNAs might play roles at metamorphic climax. The data provided in this research would be helpful for further identifying of metamorphosis-related genes in P. olivaceus.