Differential Gene Expression during Larval Metamorphic Development in the Pearl Oyster,<i> Pinctada fucata</i>, Based on Transcriptome Analysis

Li, Haimei; Zhang, Bo; Huang, Guiju; Liu, Baosuo; Fan, Sigang; Zhang, Dongling; Yu, Dahui

doi:https://doi.org/10.1155/2016/2895303

International Journal of Genomics

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusions Additional Points Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Special Issue

Recent Advances in High Throughput Sequencing Analysis

View this Special Issue

Research Article | Open Access

Volume 2016 | Article ID 2895303 | https://doi.org/10.1155/2016/2895303

Differential Gene Expression during Larval Metamorphic Development in the Pearl Oyster, Pinctada fucata, Based on Transcriptome Analysis

Haimei Li,^1,2Bo Zhang,¹Guiju Huang,¹Baosuo Liu,¹Sigang Fan,¹Dongling Zhang,³and Dahui Yu¹

Academic Editor: Quanhu Sheng

Received28 Jun 2016

Revised26 Aug 2016

Accepted20 Sept 2016

Published23 Oct 2016

Abstract

P. fucata experiences a series of transformations in appearance, from swimming larvae to sessile juveniles, during which significant changes in gene expression likely occur. Thus, P. fucata could be an ideal model in which to study the molecular mechanisms of larval metamorphosis during development in invertebrates. To study the molecular driving force behind metamorphic development in larvae of P. fucata, transcriptomes of five larval stages (trochophore, D-shape, umbonal, eyespots, and spats) were sequenced using an Illumina HiSeq™ 2000 system and assembled and characterized with the transcripts of six tissues. As a result, a total of 174,126 unique transcripts were assembled and 60,999 were annotated. The number of unigenes varied among the five larval stages. Expression profiles were distinctly different between trochophore, D-shape, umbonal, eyespots, and spats larvae. As a result, 29 expression trends were sorted, of which eight were significant. Among others, 80 development-related, differentially expressed unigenes (DEGs) were identified, of which the majority were homeobox-containing genes. Most DEGs occurred among trochophore, D-shaped, and UES (umbonal, eyespots, and spats) larvae as verified by qPCR. Principal component analysis (PCA) also revealed significant differences in expression among trochophore, D-shaped, and UES larvae with ten transcripts identified but no matching annotations.

1. Introduction

Metamorphosis is a series of key steps in the process of larval development, the success of which affect the survival of the organism. Metamorphosis is prevalent in insects, amphibians, some fishes, and many marine invertebrates, such as barnacles, sponges, shellfishes, shrimps, and echinoderms. Similar to most benthic marine invertebrates, the pearl oyster (Pinctada fucata) has a microscopic, free-swimming larval phase in their complex life cycle [1]. Oyster larvae spend several weeks in the water column before attaining competency to attach and metamorphose, commencing their sessile life. The developmental processes of P. fucata, from swimming larvae to sessile spats, have been classified into six stages: fertilized egg, trochophore, D-shaped, umbonal, juvenile, and adult stages [2]. In oysters, the transition from free-swimming larvae to the attached juvenile form often requires morphological, physiological, structural, and functional changes, which are under genetic regulatory control [3]. Therefore, the identification of key developmental genes involved in the metamorphosis of P. fucata larvae, as well as characterizing their expression patterns, is important to understand the molecular mechanism of metamorphic development of this economically important species.

Numerous studies have been conducted to explore the mechanisms of hormones, neurotransmitters, genes, and signaling pathways that regulate larval metamorphic development. Some studies have demonstrated that eight superfamily genes showed differential expression during the metamorphosis of Ciona intestinalis [4–7]. Several homeobox-containing genes were found to be responsible for larval metamorphic development in Haliotis rufescens [8–10]. In addition, abnormal dopamine and adrenaline were observed in the larval attaching stage of the Pacific oyster, Crassostrea gigas [11], while a different study observed increased expression of a molluscan growth and differentiation factor (mGDF) in the metamorphosing stage of the same organism [12]. These findings indicate the diversity of genes involved in the transitions of larval forms.

However, previous studies have focused on changes in a small number of genes and have provided a fragmented view of the genetic modulation of larval metamorphosis. Recent developments in sequencing technology have allowed for the development of new genomic tools, which can provide a more global view of changes in gene expression over the course of larval developmental stages [13–16]. In terms of genome-wide studies, transcriptome analyses are considered to be an ideal choice for obtaining comprehensive information regarding animal development and growth [17, 18]. For P. fucata, the draft genome [19] and tissue transcriptomes [20–23] have been recently reported. Based on the transcriptomic sequences from a mixture of nine developmental stages of P. fucata [19], biomineralization-related gene expression profiles during larval development have been investigated [24, 25] and genes involved in body patterning [26], transcription factors [27], and homeobox genes [28] have been identified. Nonetheless, developmentally important genes and their expression patterns during the larval stages of developing P. fucata have not been systematically studied at the transcript level to date. In the present study, the transcriptomes of five larval stages (trochophore, D-shape, umbonal, eyespots, and spats) and six tissues (gill, adductor muscle, hepatopancreas, mantle, hemocytes, and pearl sac) from P. fucata were sequenced using Illumina HiSeq 2000, with an emphasis on the molecular mechanisms underlying larval metamorphic development. This study aims to provide a valuable insight into the mechanisms of genetic modulation over the course of larval metamorphic development for P. fucata as well as for other molluscan species.

2. Materials and Methods

2.1. Larval Culture and Sample Collection

Larvae of P. fucata were bred (using several females and males of a selectively bred F3 generation as parents) through artificial insemination on March 10, 2013, in Sanya, Hainan Island, China, as described by Fujimura et al. [2]. Fertilized eggs were incubated in a 1000 L tank at 24°C. After removing nondeveloping embryos and dead larvae, trochophore, D-shaped, umbonal, eyespots, and spats larvae stages were harvested with filtering net at 12 h, 36 h, 11.5 d, 18.5 d, and 23.5 d after fertilization, respectively, and immediately preserved in RNA later (TaKaRa Bio Inc) until RNA extraction. Meanwhile, RNAs of six tissues (gill, adductor muscle, hepatopancreas, mantle, hemocytes, and pearl sac) of three other adult animals were sequenced for a more robust assembly.

2.2. RNA Extraction and cDNA Library Preparation

Following the manufacturer’s instructions, total RNA was extracted from five developmental stages (each stage with thousands of larvae) and six tissues using Trizol and RNAs of each type of tissues of the three individuals were mixed by equal weight. RNA integrity and quantity were confirmed by lab-on-chip analysis using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and visualized on a 1% agarose gel. Then cDNA was synthesized using the mRNA fragments as templates as usual and was sequenced by the BGI (Shenzhen, China) using the Illumina HiSeq 2000 system (San Diego, CA, USA) (PE100).

2.3. Sequence Assembly and Annotation

After filtering out low quality sequences (containing more than 5% ambiguous “N” nucleotides or >20% reads) and the removal of adapters from raw data, clean sequence data was assembled into unigenes using Trinity software and subsequently clustered by TGICL v2.1 (-l 40 -c 10 -v 20) [29]. Phrap (-repeat_stringency 0.95 -minmatch 35 -minscore 35) (Release 23.0) was used to produce the longest sequence possible (http://www.phrap.org/). Assembled unigenes were annotated based on the Nr, Swissprot, KEGG, and COG databases. The sequence direction and amino sequence of the predicted coding region (CDS) of unannotated unigenes were determined using ESTScan with default settings [30]. Functional annotations and classifications were performed by using Blast2GO [31] and WEGO [32] ( value threshold 1 10⁻⁵), respectively.

2.4. Normalization and Quantification of Gene Expression

Sequencing reads were mapped to the assembled reference sequence using SOAP aligner/soap2 (-m 0 -x 500 -s 40 -l 35 -v 5 -r 2) [33], a tool designed specifically to assemble short sequence alignments. The coverage of reads from a given gene was used to calculate the expression level of that gene, which was measured by fragments per kilobase exon per million fragments (FPKM) [34], with the following formula:where FPKM is the expression level of a unigene, is the number of fragments that uniquely aligned to the unigene, is the total number of fragments that uniquely aligned to all unigenes, and is the number of bases in the CDS of the unigene. The FPKM method eliminates the influence of sequences of differing lengths and coverage level on the calculation of gene expression. Therefore, the calculated gene expression can be directly used for comparing the difference in gene expression between samples.

2.5. Differential Gene Expression (DEGs) across Developmental Stages

Differential gene expression among different larval stages was carried out via principal component analysis (PCA) using the R package (http://www.r-project.org/) according to the manual. The pairwise differential expression conducted by edgeR, with a threshold of the false discovery rate (FDR) ≤ 0.001 and an absolute value of log 2 Ratio ≥ 1, was used to judge the significance of differences in gene expression. Trends in the expression of all differentially expressed genes were sorted using STEM (Short Time-Series Expression Miner, v1.3.8) [35]. Functional annotation and classification of genes involved in significant trends were performed by using Blast2GO [31] and WEGO [32], respectively. The enriched metabolic pathways or signal transduction pathways of genes were identified based on the KEGG database [36].

2.6. Identification and Expression Profile of Genes Involved in the Larval Metamorphic Development of P. fucata

According to annotations by Nr and Swissprot, development-related genes were identified with those that had been previously identified as keywords in the significant trends from the prior step. If several unigenes were assigned to the same reference gene, the sequence with the lowest value (Nr and Swissprot annotation value) was selected as a representative. Then, the heatmap 2 module of the gplots package in R (https://cran.r-project.org/web/packages/gplots/index.html) was used to perform the clustering analysis of gene expression on the normalized, filtered sequences to identify genes that were significantly different among the five developmental stages.

2.7. qPCR Verification of Expression Trends of Development-Related Genes

In order to verify the integrity of the transcriptome sequences and the expression levels as revealed by RNA-Seq, eight development-related genes were selected randomly for qPCR verification. The genes and respective primers are given in Table 1. qPCR was performed using an Eppendorf real-time- (RT-) PCR system (Eppendorf, Hamburg, Germany) using a SYBR(R) Premix Ex TaqTM kit (TaKaRa) according to the manufacturer’s protocol. Transcript levels of target genes were normalized against the level of a reference gene (18S rRNA). The qRT-PCR reactions were performed under the following conditions: 94°C for 5 min (one cycle), 94°C for 20 s, 50°C to 60°C for 20 s, and 72°C for 20 s (50 cycles). The comparative CT method ( method) was used to determine the relative mRNA abundance [37].

3. Results

3.1. Sequence Assembly and Annotation

Over 55 million reads per sample were generated with a base call accuracy (Q20) of over 97%. The number of contigs varied from 118,010 to 215,808, with a median length (N50) of 352 to 582 bp (Table 2). The number of unigenes varied from 51,102 to 113,516 among samples, with the mean length ranging from 536 to 689 bp, while N50 ranged from 495 to 1,025, respectively. In total, 174,126 unigenes were assembled with a mean length of 866 bp and an N50 of 1, 569 (based on 11 samples). Most unigenes were 100–500 bp long, and 26% were greater than 1,000 bp (Figure 1).

In total, 60,999 unigenes were annotated (Table 3) and 74,912 CDSs (43.02%) were predicted (13,966 predicted by ESTScan) (Figure 1, Table 3). Different databases annotated different numbers of unigenes (Table 3), where the most unigenes (49,580) were annotated by Swissprot database and the least (13,465) by the GO database (Table 3). Numbers of specific and shared unigenes annotated by COG, KEGG, Nr, and Swissprot terms can be visualized in Figure 2. Among them, 12,582 unigenes were annotated by the four databases and 8,784 were annotated specifically by Nr (Figure 2). Both KEGG and Swissprot analyses shared the most unigenes (41,605) and COG and Nr shared the least (13,385).

The 23,754 COG-annotated unigenes can be further classified into 25 functional groups, half of which were sorted into the “general function” group (Figure 3). The GO analysis revealed that 10,165 unigenes were attributed to biological process, 8,442 unigenes to cell components, and 10,588 unigenes to molecular function (Figure 4). The top 26 KEGG pathways are summarized in Table 4. Most unigenes (5,184 out of 43,753) were involved in metabolic pathways and 1,401 unigenes were involved in calcium signaling pathways, some of which may be involved in shell formation. Finally, many unigenes in the top 26 KEGG pathways were involved in immune pathways.

3.2. Differential Gene Expression (DEGs) and Expression Trends during Developmental Stages

Principal component analysis (PCA) revealed that differences in the expression of unigenes were vast among trochophore, D-shaped, and UES (umbonal, eyespots, and spats) larvae, but small within UES stages (Figure 5(a)). Based on gene effects, measured by the first principal component value, a total of 10 transcripts with unknown functions were identified to be key factors involved in the larval development of P. fucata. Differences in the gene expression of these transcripts were the greatest in trochophore, D-shaped, and UES larvae. They were relatively highly expressed in trochophore larvae and then downregulated during the D-shaped stage, and some were subsequently upregulated during the UES stage, including Unigens23340_All, Unigene8217_All, Unigene50061_All, and Cl616_All (Figure 5(b)).

(a)

(b)

The numbers of up- and downregulated unigenes were also much greater during early stage transitions (Figure 6), consistent with the results of the PCA. From trochophore to D-shape larvae, there were 18,725 unigenes upregulated and 13,162 downregulated. In total, there were 57,228 DEGs among the five developmental stages (Figure 7). Additionally, 17,609 genes were preferentially expressed at a single developmental stage, which indicates that they play an important role in the corresponding developmental stage, while 39,619 were expressed preferentially during more than two stages. A total of 20,518 genes were differentially expressed in all five of the development stages. All differentially expressed genes were sorted into 29 expression trends (Figure 8), of which eight trends were significant, comprising over 45% of the total DEGs. Furthermore, 6,653 unigenes were expressed highly only during the trochophore stage. Across the five stages, 3,340 unigenes were expressed in an increasing pattern, while 2,631 unigenes were expressed in a decreasing pattern.

Figure 8

Expression trends of unigenes across trochophore, D-shaped, umbonal, eyespots, and spats larval stages of Pinctada fucata. The profiles were ordered based on the value of the number (at bottom-left corner) of genes assigned versus expected. Color square frames denote significant profiles (). Each graph displays the mean pattern of expression (black lines) of the profiled genes. The number of profiles in each cluster is indicated in the top left corner. The -axis represents stages and the -axis represents log 2-fold change of gene expression.

3.3. Functional Enrichment Analysis

A functional enrichment analysis of the unigenes from the eight significant trends showed that there were 104, 54, and 46 GO terms for biological processes, molecular functions, and cellular components, respectively, identified for GO function enrichment (see Supplementary Table in Supplementary Material available online at http://dx.doi.org/10.1155/2016/2895303), and 272 pathways identified for KEGG pathway enrichment (Supplementary Table ). For GO enrichment data, trends 0, 2, 3, 24, 26, 28, and 29 were involved in biological processes, where most unigenes belonged to trend 3 and were related to various metabolic processes. Trends 0, 3, 24, 28, and 29 were involved in molecular function, where most unigenes fell within trends 0, 3, and 28, and were related to binding and catalytic activity. Only trends 0, 3, 28, and 29 were involved in cellular components, and most unigenes fell within trend 3 and were involved in processes related to membranes and organelles.

Trends 0, 2, 3, 24, 26, 27, 28, and 29 were implicated in KEGG pathway enrichment. For 272 significant enriched pathways, 81 pathways were observed in trend 29, 77 in trend 28, 30 in trend 27, and 27 in trend 26 (Supplementary Table ). In trends 28 and 29, most unigenes were involved in immune responses.

In trend 0, GO enrichment showed that macromolecule metabolic processes were the dominant groups in biological process, followed by positive regulation of biological process. For pathways involved in molecular function, DNA binding was the most representative category, while for cellular component pathways, all of the genes participate in processes integral to the inner mitochondrial membrane and are intrinsic to the inner mitochondrial membrane. In the KEGG category, spliceosomes were prevalent, followed by genes involved in the cell cycle.

In trend 3, GO enrichment data showed that 6,653 DEG unigenes were further categorized into 43 functional groups; among them, macromolecule metabolic processes were the dominant groups in biological process, followed by cellular macromolecule metabolic processes. In the molecular function category, a high percentage of genes came from the binding and protein binding groups. Spliceosomes were the most representative, followed by RNA transport and regulation of the actin cytoskeleton.

In trends 27 and 28, GO enrichment reveled that there were no significant categories. However, the calcium signaling pathway, hedgehog signaling pathway, and insulin signaling pathway were significantly enriched in the KEGG database, as they are all involved in early development. In trend 28, small molecule metabolic processes were the dominant group in biological process followed by ion transport. Catalytic activity was the most prevalent in the molecular function category, followed by transporter activity and transmembrane transporter activity. In cellular component pathways, membrane was the most representative, followed by plasma membrane. In KEGG enrichment categories, we also found genes related to the calcium signaling pathway in trend 27.

In trend 29, translational elongation was the only enriched category for biological processes, while three categories were enriched in molecular function, including genes involved in oxidoreductase activity, catalytic activity, and lyase activity. In cellular component pathways, five categories were enriched, while vacuole was the dominant group, followed by lytic vacuole and lysosome groups. Genes in trend 29 were enriched in only one KEGG pathway, translational elongation, with a significant value.

3.4. Identification and Expression Profiling of Genes Involved in the Larval Metamorphic Development of P. fucata

In total, 80 development-related candidate DEGs were identified and summarized in Table 5, which can be mapped to known developmentally important genes, including several homeobox genes, and can be sorted into 10 trends: trend 0 (25 unigenes), trend 28 (16), trend 3 (15), trend 29 (9), and six other trends (1–5). Cluster analyses suggested that most development-related candidate genes were highly expressed in the early developmental stages (Figure 9), including engrailed-2-B, pax family, fox family members e1 and p1, Wnt-4, and BMP3/3B upregulated in the trochophore stage and LIM, foxg1, Hox3, bicaudal, hedgehog, EGFR, foxl2, and bmp 2b genes upregulated from the D-shaped stage until the eyespots stage. In the spats stage, wnt1 and notch-like protein 2 gene were upregulated. The qPCR showed that the trends in the expression of selected genes (Figure 10) were consistent with the expression trends indicated by the trend analysis of RNA-Seq data (Figure 9), indicating that the sequence data in our study are reliable.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

4. Discussion

Not only does the pearl oyster, P. fucata, make an ideal model organism for studies of biomineralization, but also it is a good model to study the early stage metamorphic development of invertebrates. In this study, we sequenced the transcriptomes of five developmental stages in P. fucata, with the aim of developing a better understanding of the molecular mechanisms driving the change of one larval stage to the next during early life history. In our study, the de novo assembly was performed with six tissue transcriptomes, as the draft genome of P. fucata is not complete [19]. As a result, we obtained 174,126 unigenes, with a mean length of 866 bp. A total of 60,999 unigenes (35%) were annotated, a value slightly higher than previous reports [21–23, 38]. Poor annotation efficiencies have been widely prevalent in many marine organisms, likely owing limited genomic resources from aquaculture species in public databases to date [21–23, 38]. Alternatively, poor annotation efficiencies could be the result of the short length of the assembled unigene sequences [22] and great divergence among the genomes of marine organism. Similar scenarios have been reported in other marine organisms [39, 40]. In the KEGG annotation, we observed that many pathways were related to immunity, indicating that innate protection is vital in the early developmental stages.

Differential gene expressions (DEGs) occurred mainly during early stage transitions (Figure 6). Most genes were up- or downregulated from trochophore to D-shaped and from D-shaped to umbonal stages, indicating that processes associated with these transitions are very complicated. Principal component analyses yielded consistent results, where we identified 10 unigenes attributed to the divergence among trochophore, D-shaped, and UES (umbonal, eyespots, and spats) stages in P. fucata, being highly expressed in the trochophore stage. However, no functional annotations match these functionally important sequences, indicating that further research would help to elucidate the molecular mechanism of metamorphosis in this species in the future.

The analysis of expression trends indicated that 12,009 of 13,277 unigenes are sorted into eight significant expression trend groups. Among the significant trends, there were 10,031 (trends 3, 0, and 2), 11,978 (trends 28, 29, 21, 24, 26, and 27), 9,046 (trend 28, 29, 26, and 27), 9,518 (trend 28, 29, 24, and 27), and 5,214 (trend 29, 24, and 26) unigenes displaying increased expression in trochophore, D-shaped, umbonal, eyespots, and spats stages, respectively. This conveys that more genes are expressed in the early stages, consistent with the DEG and PCA analyses in our study. Particularly, 6,653 unigenes (trend 3) were highly expressed only in the trochophore stage, 3,340 unigenes (trend 29) expressed in an increasing pattern over the course of development, and 2,631 unigenes (trend 0) expressed in a decreasing pattern. These genes are worth further investigation.

The KEGG pathway enrichment analysis indicated that most unigenes in trend 3 were involved in pathways of spliceosome or RNA transport, indicating that, in the early stage of P. fucata, RNA synthesis is more predominant. On the contrary, genes in trend 29 showed significant enrichment in translational elongation pathways, suggesting that protein synthesis is more and more prevalent during larval development. In trends 27 and 28, a large number of unigenes were involved in the calcium signaling pathway, synchronizing with the shell formation of prodissoconchs I and II in D-shaped and umbonal stages [24, 41]. In addition, immune pathways were also enriched, indicating that innate protection is important during the entire course of larval development [3, 42, 43].

In our study, 80 known development-related, differentially expressed unigenes were identified throughout the five larval stages (Table 5). Half of them were homeobox-containing genes, including genes known to be involved in the development of body patterning (engrailed, SIX3, Pax-7, LIM, and Hox family members), suggesting that these genes play important roles in the metamorphic changes of P. fucata larvae. Nearly half of the homeobox-containing genes were upregulated in trochophore and D-shaped stages (Figure 9). We identified two Hox genes, Hox5 and Hox3 (Figure 9), which are highly expressed in D-shaped veliger, indicating that they are involved in the growth of D-shaped larvae. We also found early developmentally relevant signaling molecules such as Hedghog, TGFβ, and Wnt family, which are known to play important roles in axis formation, muscle differentiation, and nervous system development [26]. Recent evidence has suggested that classic morphogens, such as Wnts, TGFβ/BMP family members, and Hedgehogs, may all serve as axon guidance cues for a variety of axons in different organisms [44]. Several studies have provided increasing evidence that Sonic hedgehog (Shh) is an important axon guidance cue throughout vertebrate neural development [45, 46]. In our study, one hedgehog gene (Unigene23139_All) was identified and highly expressed in umbonal and eyespots stages (Figure 9), suggesting that increased neural development was likely taking place during those stages.

The Wnt signal pathway has been shown to play an important role in the segmentation of the marine polychaete Capitella capitata [47, 48], while the maintenance of primitive hematopoiesis has been attributed to Wnt4 in the vertebrate embryo [49]. Both Wnt4 and Wnt1 were observed in this study; Wnt4 was highly expressed in the trochophore stage, while Wnt1 was expressed in an increasing pattern over the course of larval development, suggesting possible involvement in blood formation from the beginning of development and in body transformation during all stages. Ten classes of fox genes were also found in our study, comprising the largest number of genes identified in the DEGs and displaying different trends in expression. FoxL2, XX-dominantly expressed in the differentiating ovaries of mammals [50], birds [51], and fish [52–54], was expressed highly in the D-shaped stage, suggesting the possible beginning of sexual development. Some important growth-related genes were also identified and differentially expressed among the five development stages, including EGF-like, MAPK, and MAPKK genes, which were actively expressed during the five developmental stages, and may contribute significantly to the transitions between developmental stages in P. fucata larvae.

Nonetheless, the body form transformations that take place during larval development involve a series of morphological and physiological changes and corresponding molecular changes, which have not been systematically studied and remain unclear. Therefore, a more broad understanding of the molecular underpinnings of important biological processes still merits further investigation.

5. Conclusions

In this study, a total of 174,126 unique transcripts were assembled and 60,999 were annotated. The number of unigenes varied between the five larval stages. The expression profiles of trochophore, D-shaped, and UES (umbonal, eyespots, and spats) larvae were distinctly different. Most unigenes were up- or downregulated in early stage transitions and 29 expression trends were sorted, eight of which were significant. In total, 80 development-related, differentially expressed unigenes were identified and eight were verified by qPCR. These observations should be helpful in understanding the molecular mechanisms of the larval metamorphic development of P. fucata.

Additional Points

Highlights. (i) A large number of assembled transcripts from Pinctada fucata are reported for the first time. (ii) Large variations in expression of DEGs related to development were observed in early larval stages. (iii) Twenty-nine expression trends were identified for the first time.

Competing Interests

The authors declare that they have no competing interests.

Authors’ Contributions

Haimei Li and Dahui Yu conceived and designed the project. Haimei Li, Bo Zhang, and Baosuo Liu cultured and collected the P. fucata larval samples and the adult samples. Haimei Li, Guiju Huang, and Sigang Fan carried out the transcriptome analysis and qPCR experiments. Dongling Zhang provided technical assistance. Haimei Li and Dahui Yu wrote the manuscript. All listed authors have read, edited, and approved the final manuscript.

Acknowledgments

This work was supported by The National Natural Science Foundation of China (NSFC) (31372525), the Earmarked Fund for China Agriculture Research System (Grant no. CARS-48), and Special Fund for Marine Fisheries Research and Extension of Guangdong Province (A201301A02, A201301A08, Z2014003 Z2014006, and Z2015009). The authors are also grateful to the Guangzhou Gene denovo Biotechnology Co., Ltd., for assisting in the data analysis.

Supplementary Materials

Supplement files Supplement table 1 GO function enrichment in significant trends Supplement table 2 KEGG pathway enrichment in significant trends

Supplementary Material

References

P.-Y. Qian, “Larval settlement of polychaetes,” Hydrobiologia, vol. 402, pp. 239–253, 1999.
View at: Publisher Site | Google Scholar
T. Fujimura, K. Wada, and T. Iwaki, “Development and morphology of the pearl oyster larvae, Pinctada fucata,” Venus, vol. 54, pp. 25–48, 1995.
View at: Google Scholar
A. Heyland and L. L. Moroz, “Signaling mechanisms underlying metamorphic transitions in animals,” Integrative and Comparative Biology, vol. 46, no. 6, pp. 743–759, 2006.
View at: Publisher Site | Google Scholar
C. Devine, V. F. Hinman, and B. M. Degnan, “Evolution and developmental expression of nuclear receptor genes in the ascidian Herdmania,” The International Journal of Developmental Biology, vol. 46, no. 4, pp. 687–692, 2002.
View at: Google Scholar
J. M. Arnold, R. Eri, B. M. Degnan, and M. F. Lavin, “Novel gene containing multiple epidermal growth factor-like motifs transiently expressed in the papillae of the ascidian tadpole larvae,” Developmental Dynamics, vol. 210, no. 3, pp. 264–273, 1997.
View at: Publisher Site | Google Scholar
A. Nakayama, Y. Satou, and N. Satoh, “Isolation and characterization of genes that are expressed during Ciona intestinalis metamorphosis,” Development Genes and Evolution, vol. 211, no. 4, pp. 184–189, 2001.
View at: Publisher Site | Google Scholar
R. Eri, J. M. Arnold, and V. F. Hinman, “Patterning of dopaminergic neurotransmitter identity among Caenorhabditis elegans ray sensory neurons by a TGF family signaling pathway and a Hox gene,” Development, vol. 126, no. 24, pp. 5819–5831, 1999.
View at: Google Scholar
B. M. Degnan and D. E. Morse, “Identification of eight homeobox-containing transcripts expressed during larval development and at metamorphosis in the gastropod mollusc Haliotis rufescens,” Molecular Marine Biology and Biotechnology, vol. 2, no. 1, pp. 1–9, 1993.
View at: Google Scholar
B. M. Degnan, S. M. Degnan, G. Fentenany, and D. E. Morse, “A Mox homeobox gene in the gastropod mollusc Haliotis rufescens is differentially expressed during larval morphogenesis and metamorphosis,” FEBS Letters, vol. 411, no. 1, pp. 119–122, 1997.
View at: Publisher Site | Google Scholar
A. F. Giusti, V. F. Hinman, S. M. Degnan, B. M. Degnan, and D. E. Morse, “Expression of a Scr/Hox5 gene in the larval central nervous system of the gastropod Haliotis, a non-segmented spiralian lophotrochozoan,” Evolution and Development, vol. 2, no. 5, pp. 294–302, 2000.
View at: Publisher Site | Google Scholar
S. L. Coon, W. K. Fitt, and D. B. Bonar, “Competence and delay of metamorphosis in the Pacific oyster Crassostrea gigas,” Marine Biology, vol. 106, no. 3, pp. 379–387, 1990.
View at: Publisher Site | Google Scholar
C. Lelong, M. Mathieu, and P. Favrel, “Structure and expression of mGDF, a new member of the transforming growth factor-β superfamily in the bivalve mollusc Crassostrea gigas,” European Journal of Biochemistry, vol. 267, no. 13, pp. 3986–3993, 2000.
View at: Publisher Site | Google Scholar
T. J. Fiedler, A. Hudder, S. J. McKay et al., “The transcriptome of the early life history stages of the California sea hare Aplysia californica,” Comparative Biochemistry and Physiology Part D: Genomics and Proteomics, vol. 5, no. 2, pp. 165–170, 2010.
View at: Publisher Site | Google Scholar
J. D. Lambert, X. Y. Chan, B. Spiecker, and H. C. Sweet, “Characterizing the embryonic transcriptome of the snail Ilyanassa,” Integrative and Comparative Biology, vol. 50, no. 5, pp. 768–777, 2010.
View at: Publisher Site | Google Scholar
P. Huan, H. Wang, and B. Liu, “Transcriptomic analysis of the clam Meretrix meretrix on different larval stages,” Marine Biotechnology, vol. 14, no. 1, pp. 69–78, 2012.
View at: Publisher Site | Google Scholar
Z.-X. Huang, Z.-S. Chen, C.-H. Ke et al., “Pyrosequencing of Haliotis diversicolor transcriptomes: insights into early developmental molluscan gene expression,” PLoS ONE, vol. 7, no. 12, Article ID e51279, 2012.
View at: Publisher Site | Google Scholar
J. Qin, Z. Huang, J. Chen, Q. Zou, W. You, and C. Ke, “Sequencing and de novo analysis of Crassostrea angulata (Fujian Oyster) from 8 different developing phases using 454 GSFLx,” PLoS ONE, vol. 7, no. 8, Article ID e43653, 2012.
View at: Publisher Site | Google Scholar
S. Bassim, A. Tanguy, B. Genard, D. Moraga, and R. Tremblay, “Identification of Mytilus edulis genetic regulators during early development,” Gene, vol. 551, no. 1, pp. 65–78, 2014.
View at: Publisher Site | Google Scholar
T. Takeuchi, T. Kawashima, R. Koyanagi et al., “Draft genome of the pearl oyster Pinctada fucata: a platform for understanding bivalve biology,” DNA Research, vol. 19, no. 2, pp. 117–130, 2012.
View at: Publisher Site | Google Scholar
Y. Shi, C. Yu, Z. Gu, X. Zhan, Y. Wang, and A. Wang, “Characterization of the pearl oyster (Pinctada martensii) mantle transcriptome unravels biomineralization genes,” Marine Biotechnology, vol. 15, no. 2, pp. 175–187, 2013.
View at: Publisher Site | Google Scholar
A. Wang, Y. Wang, Z. Gu, S. Li, Y. Shi, and X. Guo, “Development of expressed sequence tags from the pearl oyster, Pinctada martensii Dunker,” Marine Biotechnology, vol. 13, no. 2, pp. 275–283, 2011.
View at: Publisher Site | Google Scholar
X. Zhao, Q. Wang, Y. Jiao et al., “Identification of genes potentially related to biomineralization and immunity by transcriptome analysis of pearl sac in pearl oyster Pinctada martensii,” Marine Biotechnology, vol. 14, no. 6, pp. 730–739, 2012.
View at: Publisher Site | Google Scholar
S. Kinoshita, N. Wang, H. Inoue et al., “Deep sequencing of ESTs from nacreous and prismatic layer producing tissues and a screen for novel shell formation-related genes in the pearl oyster,” PLoS ONE, vol. 6, no. 6, Article ID e21238, 2011.
View at: Publisher Site | Google Scholar
Y. Miyazaki, T. Nishida, H. Aoki, and T. Samata, “Expression of genes responsible for biomineralization of Pinctada fucata during development,” Comparative Biochemistry and Physiology—B: Biochemistry and Molecular Biology, vol. 155, no. 3, pp. 241–248, 2010.
View at: Publisher Site | Google Scholar
J. Liu, D. Yang, S. Liu et al., “Microarray: a global analysis of biomineralization-related gene expression profiles during larval development in the pearl oyster, Pinctada fucata,” BMC Genomics, vol. 16, no. 1, article 325, 2015.
View at: Publisher Site | Google Scholar
D. H. E. Setiamarga, K. Shimizu, J. Kuroda et al., “An in-silico genomic survey to annotate genes coding for early development-relevant signaling molecules in the pearl oyster, Pinctada fucata,” Zoological Science, vol. 30, no. 10, pp. 877–888, 2013.
View at: Publisher Site | Google Scholar
H. Koga, N. Hashimoto, D. G. Suzuki et al., “A genome-wide survey of genes encoding transcription factors in Japanese pearl oyster Pinctada fucata: II. Tbx, Fox, Ets, HMG, NFκB, bZIP, and C2H2 zinc fingers,” Zoological Science, vol. 30, no. 10, pp. 858–867, 2013.
View at: Publisher Site | Google Scholar
Y. Morino, K. Okada, M. Niikura, M. Honda, N. Satoh, and H. Wada, “A genome-wide survey of genes encoding transcription factors in the Japanese pearl oyster, Pinctada fucata: I. Homeobox genes,” Zoological Science, vol. 30, no. 10, pp. 851–857, 2013.
View at: Publisher Site | Google Scholar
G. Pertea, X. Huang, F. Liang et al., “TIGR gene indices clustering tools (TGICL): a software system for fast clustering of large EST datasets,” Bioinformatics, vol. 19, no. 5, pp. 651–652, 2003.
View at: Publisher Site | Google Scholar
C. Iseli, C. V. Jongeneel, and P. Bucher, “ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences,” in Proceedings of the International Conference on Intelligent Systems for Molecular Biology (ISMB '99), pp. 138–148, Heidelberg, Germany, 1999.
View at: Google Scholar
A. Conesa, S. Götz, J. M. García-Gómez, J. Terol, M. Talón, and M. Robles, “Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research,” Bioinformatics, vol. 21, no. 18, pp. 3674–3676, 2005.
View at: Publisher Site | Google Scholar
J. Ye, L. Fang, H. Zheng et al., “WEGO: a web tool for plotting GO annotations,” Nucleic Acids Research, vol. 34, pp. W293–W297, 2006.
View at: Publisher Site | Google Scholar
R. Li, C. Yu, Y. Li et al., “SOAP2: an improved ultrafast tool for short read alignment,” Bioinformatics, vol. 25, no. 15, pp. 1966–1967, 2009.
View at: Publisher Site | Google Scholar
A. Mortazavi, B. A. Williams, K. McCue, L. Schaeffer, and B. Wold, “Mapping and quantifying mammalian transcriptomes by RNA-Seq,” Nature Methods, vol. 5, no. 7, pp. 621–628, 2008.
View at: Publisher Site | Google Scholar
J. Ernst and Z. Bar-Joseph, “STEM: a tool for the analysis of short time series gene expression data,” BMC Bioinformatics, vol. 7, article 191, 2006.
View at: Publisher Site | Google Scholar
M. Kanehisa, M. Araki, S. Goto et al., “KEGG for linking genomes to life and the environment,” Nucleic Acids Research, vol. 36, no. 1, pp. D480–D484, 2008.
View at: Publisher Site | Google Scholar
K. J. Livak and T. D. Schmittgen, “Analysis of relative gene expression data using real-time quantitative PCR and the $2^{- Δ Δ C_{T}}$ method,” Methods, vol. 25, no. 4, pp. 402–408, 2001.
View at: Publisher Site | Google Scholar
X.-D. Huang, M. Zhao, W.-G. Liu et al., “Gigabase-scale transcriptome analysis on four species of pearl oysters,” Marine Biotechnology, vol. 15, no. 3, pp. 253–264, 2013.
View at: Publisher Site | Google Scholar
E. Meyer, G. V. Aglyamova, S. Wang et al., “Sequencing and de novo analysis of a coral larval transcriptome using 454 GSFlx,” BMC Genomics, vol. 10, article 219, pp. 1–8, 2009.
View at: Publisher Site | Google Scholar
R. Bettencourt, M. Pinheiro, C. Egas et al., “High-throughput sequencing and analysis of the gill tissue transcriptome from the deep-sea hydrothermal vent mussel Bathymodiolus azoricus,” BMC Genomics, vol. 11, article 559, 2010.
View at: Publisher Site | Google Scholar
D. Fang, G. Xu, Y. Hu, C. Pan, L. Xie, and R. Zhang, “Identification of genes directly involved in shell formation and their functions in pearl oyster, Pinctada fucata,” PLoS ONE, vol. 6, no. 7, Article ID e21860, 2011.
View at: Publisher Site | Google Scholar
E. A. Williams, B. M. Degnan, H. Gunter, D. J. Jackson, B. J. Woodcroft, and S. M. Degnan, “Widespread transcriptional changes pre-empt the critical pelagic-benthic transition in the vetigastropod Haliotis asinina,” Molecular Ecology, vol. 18, no. 5, pp. 1006–1025, 2009.
View at: Publisher Site | Google Scholar
D. E. Clapham, “Calcium signaling,” Cell, vol. 131, no. 6, pp. 1047–1058, 2007.
View at: Publisher Site | Google Scholar
F. Charron and M. Tessier-Lavigne, “Novel brain wiring functions for classical morphogens: a role as graded positional cues in axon guidance,” Development, vol. 132, no. 10, pp. 2251–2262, 2005.
View at: Publisher Site | Google Scholar
F. Charron, E. Stein, J. Jeong, A. P. McMahon, and M. Tessier-Lavigne, “The morphogen sonic hedgehog is an axonal chemoattractant that collaborates with netrin-1 in midline axon guidance,” Cell, vol. 113, no. 1, pp. 11–23, 2003.
View at: Publisher Site | Google Scholar
J. K. Chen, J. Taipale, M. K. Cooper, and P. A. Beachy, “Inhibition of Hedgehog signaling by direct binding of cyclopamine to Smoothened,” Genes & Development, vol. 16, no. 21, pp. 2743–2748, 2002.
View at: Publisher Site | Google Scholar
E. C. Seaver and L. M. Kaneshige, “Expression of ‘segmentation’ genes during larval and juvenile development in the polychaetes Capitella sp. I and H. elegans,” Developmental Biology, vol. 289, no. 1, pp. 179–194, 2006.
View at: Publisher Site | Google Scholar
J. L. Christian, B. J. Gavin, A. P. McMahon, and R. T. Moon, “Isolation of cDNAs partially encoding four Xenopus Wnt-1int-1-related proteins and characterization of their transient expression during embryonic development,” Developmental Biology, vol. 143, no. 2, pp. 230–234, 1991.
View at: Publisher Site | Google Scholar
H. T. Tran, B. Sekkali, G. Van Imschoot, S. Janssens, and K. Vleminckx, “Wnt/β-catenin signaling is involved in the induction and maintenance of primitive hematopoiesis in the vertebrate embryo,” Proceedings of the National Academy of Sciences of the United States of America, vol. 107, no. 37, pp. 16160–16165, 2010.
View at: Publisher Site | Google Scholar
D. Baron, F. Batista, J. S. Chaffaux Cocquet et al., “Foxl2 gene and the development of the ovary: a story about goat, mouse, fish and woman,” Reproduction Nutrition Development, vol. 45, no. 6, pp. 729–729, 2005.
View at: Publisher Site | Google Scholar
M. S. Govoroun, M. Pannetier, E. Pailhoux et al., “Isolation of chicken homolog of the FOXL2 gene and comparison of its expression patterns with those of aromatase during ovarian development,” Developmental Dynamics, vol. 231, no. 4, pp. 859–870, 2004.
View at: Publisher Site | Google Scholar
D. Baron, J. Cocquet, X. Xia, M. Fellous, Y. Guiguen, and R. A. Veitia, “An evolutionary and functional analysis of FoxL2 in rainbow trout gonad differentiation,” Journal of Molecular Endocrinology, vol. 33, no. 3, pp. 705–715, 2004.
View at: Publisher Site | Google Scholar
M. Nakamoto, M. Matsuda, D.-S. Wang, Y. Nagahama, and N. Shibata, “Molecular cloning and analysis of gonadal expression of Foxl2 in the medaka, Oryzias latipes,” Biochemical and Biophysical Research Communications, vol. 344, no. 1, pp. 353–361, 2006.
View at: Publisher Site | Google Scholar
D.-S. Wang, T. Kobayashi, L.-Y. Zhou et al., “Foxl2 up-regulates aromatase gene transcription in a female-specific manner by binding to the promoter as well as interacting with Ad4 binding protein/steroidogenic factor 1,” Molecular Endocrinology, vol. 21, no. 3, pp. 712–725, 2007.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2016 Haimei Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1867

Downloads

1155

Citations