Scholarly article on topic 'Transcriptome analysis elucidates key developmental components of bryozoan lophophore development'

Transcriptome analysis elucidates key developmental components of bryozoan lophophore development Academic research paper on "Biological sciences"

Share paper
Academic journal
Sci. Rep.
OECD Field of science

Academic research paper on topic "Transcriptome analysis elucidates key developmental components of bryozoan lophophore development"

Transcriptome analysis elucidates key developmental components of bryozoan lophophore development

Yue Him Wong1*, Taewoo Ryu2*, Loqmane Seridi2, Yanal Ghosheh2, Salim Bougouffa2, Pei-Yuan Qian1 & Timothy Ravasi2,3

1The division of Life Science, School of Science, the Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, 2Computational Biosciences Research Center, Division of Biological and Environmental Sciences & Engineering, King Abdullah University of Science and Technology, Thuwal, Kingdom of Saudi Arabia, 3Division of Applied Mathematics and Computer Sciences, King Abdullah University of Science and Technology, Thuwal, Kingdom of Saudi Arabia.

The most recent phylogenomic study suggested that Bryozoa (Ectoprocta), Brachiopoda, and Phoronida are monophyletic, implying that the lophophore of bryozoans, phoronids and brachiopods is a synapomorphy. Understanding the molecular mechanisms of the lophophore development of the Lophophorata clade can therefore provide us a new insight into the formation of the diverse morphological traits in metazoans. In the present study, we profiled the transcriptome of the Bryozoan (Ectoproct) Bugula neritina during the swimming larval stage (SW) and the early (4 h) and late (24 h) metamorphic stages using the Illumina HiSeq2000 platform. Various genes that function in development, the immune response and neurogenesis showed differential expression levels during metamorphosis. In situ hybridization of 23 genes that participate in the Wnt, BMP, Notch, and Hedgehog signaling pathways revealed their regulatory roles in the development of the lophophore and the ancestrula digestive tract. Our findings support the hypothesis that developmental precursors of the lophophore and the ancestrula digestive tract are pre-patterned by the differential expression of key developmental genes according to their fate. This study provides a foundation to better understand the developmental divergence and/or convergence among developmental precursors of the lophophore of bryozoans, branchiopods and phoronids.


Received 7 May 2014

Accepted 6 August 2014

Published 10 October 2014

Correspondence and requests for materials should be addressed to P.Q. (boqianpy@ust.

hk) or T.R. (timothy.

* These authors contributed equally to this work.

n the Lophotrochrozoa clade, the phylogenetic relationship between Bryozoa (Ectoprocta), Phoronida and Brachiopoda has been elusive. They are often grouped as lophophorates, based on the superficial similarities of their filtering apparatus "lophophore", a mesosomal tentacle crown with an upstream-collecting ciliary band1. Histological studies, however, suggested that the lophophore of phoronids and brachiopods may not be homologous to that of bryozoans. Nielsen et al2 pointed out that tentacles of ectoproct lophophore have multiciliate cells in lateral ciliary bands and lacking longitudinal haemal vessels, whereas phoronid and brachiopod tentacles are monociliate and have a longitudinal haemal vessel. Although various molecular analyses have suggested that lophophorates are polyphyletic3-7, the most recent phylogenomic analysis has once again united the three phyla under the Lophophorata clade8, implying that the lophophore of bryozoans, phoronids and brachiopods is a synapomorphy, despite the differences in the ontology and anatomy of the lophophore.

In addition to the possession of lophophore, bryozoans, phoronids and brachiopods are all marine benthos with biphasic life cycles. Larval metamorphoses in these phyla are catastrophic, involving drastic morphological and anatomical transformations9-11. Because metamorphosis recapitulates some key developmental processes such as the remodeling of nervous system and the morphogenesis of lophophore9-13, studying the regulatory mechanisms of metamorphosis may provide a new insight into the formation of diverse morphological traits in metazoans. However, molecular data on these phyla are scanty. To date, only two studies revealed distinctive expression patterns of developmental genes in the swimming larvae of a bryozoan and related these genes to metamorphosis14,15. There has not been any gene expression study on metamorphosis of phoronid and brachiopod larvae to date. The findings on bryozoan metamorphosis have led to the proposal of a "pre-patterning" developmental mechanism, in which the apical blastema, the developmental precursors of the lophophore and ancestrula digestive tract, in the larval body is pre-patterned according to their future fate14. Although an intriguing developmental mechanism was proposed in these studies, it was based solely on the expression patterns of a few developmental genes and only during the swimming larval stage. Therefore, it remains unclear whether \ > _

Table 1 | Statistics for the non-redundant contigs

Dataset Number of contigs Total length of assembly Largest contig length N20 N50 N80

SW 150,683 74,929,894 14,906 1,275 601 301

4 h 116,761 59,227,507 11,688 1,317 619 307

24 h 156,556 83,057,230 15,051 1,425 668 319

Pooled data 240,173 106,885,897 11,861 1,073 501 274

*Redundant contigs were removed using CD-HIT-EST with a sequence identity threshold of 0.99. *Length statistics are presented in base pairs (bp).

the described genes are related to the morphogenesis of the lopho-phore and ancestrula digestive tract or are required only for bryozoan larval development.

In the present study, we aimed to analyze transcriptomic changes during metamorphosis of the bryozoan Bugula neritina. We profiled and compared the transcriptome of the swimming larvae (SW), early (4 h) and late (24 h) preancestrulae (intermediate metamorphic stages) of B. neritina using high-throughput transcriptome sequencing. 240,137 contigs was compiled, representing the first and, to date, the most comprehensive dataset for the Bryozoa. To study the molecular mechanism of bryozoan metamorphosis, we performed enrichment analysis on the functional annotation of differentially expressed genes. Specifically, our analysis focused on axial patterning genes including transcription factors and several different well-implicated morphogens such as Wnt, BMP, Sonic Hedgehog and Notch. The spatial expression patterns of these axial patterning genes at different developmental stages were studied using in situ hybridization.


De novo assembly of the B. neritina transcriptome. Using Illumina paired-end sequencing, we obtained a total of 54,613,482 X 2, 46,157,987 X 2, and 58,852,506 X 2 raw read pairs with poly(A)-selected cDNA from SW, 4 h and 24 h, respectively. All these raw Illumina paired-end reads were submitted to NCBI short read archive (SRA) (Biosample no. SAMN02724736-SAMN2724738). The number of contigs ($200 bp) obtained by de novo transcriptome assembly using ABySS16 with a single k value (55) was 53,270, 42,673, 51,008, and 73,378 for SW, 4 h, 24 h, and the pooled data, respectively. When multiple k values (every odd k-mer between 45 and 81) were used in ABySS and Trans-ABySS de novo transcriptome assembly16,17, the total number of contigs ($200 bp) obtained increased to 208,280, 163,116, 222,555, and 309,137 for SW, 4 h, 24 h, and the pooled data, respectively (N50 values were 651 bp, 667 bp, 741 bp, and 553 bp, respectively; see Supplemental Table S1). As we would like to present a comprehensive transcriptome for the Bryozoa, we adopted the assemblies obtained from multiple k values method. The numbers of non-redundant contigs ($200 bp) after CD-HIT-EST18 (sequence identity 0.99, i.e. 10 bp mismatches per 1 kb) were: 150,683,116,761, 156,556 and 240,173 for SW, 4 h, 24 h, and the pooled data, respectively (see Table 1).

The pooled assembly (generated from the pooled data) was the most representative when compared with the assemblies from individual stages and the re-assembled 454 pyrosequencing data from previous study19 (see Table 2). Therefore, we selected the contigs from the pooled assembly (see Supplemental dataset 1 online) for the downstream analysis. To avoid potential problems in the downstream expression analysis, assemblies from individual stages, which contain a small portion of unique contigs, were not combined to the pooled assembly. The 454 pyrosequencing assembly was also excluded because 454 pyrosequencing have much higher error rates compared with the Illumina platform, which can lead to frame-shifted proteins and inaccurate contig sequences20. In fact, transcript sequences from Illumina pooled assembly covered majority of contigs assembled from 454 data (see Table 2). Considering the relatively high error rates of the 454 reads, heterogeneity between the 454 and Illumina sequencing platforms20 and high coverage of current Illumina data, we decided not to perform hybrid assembling of 454 and Illumina data in order to maintain high accuracy of the tran-scriptome dataset.

Functional annotation. Among the pooled assembly, 64,175 contigs have significant matches in the NCBI non-redundant (nr) protein database (blastp bit score: > 50; minimal match length: > 20 amino acids). Of these nr-annotated contigs, 55,840 have at least one Gene Ontology (GO) functional annotation; 30,156 were annotated with at least one Pfam domain; and 26,474 were annotated with at least one Reactome pathway ID. A summary of functional annotation of contigs is showed in Table 3. A brief overview of the GO annotations is shown in Supplemental Figure S1. An analysis using CPAT21 showed that the coding potentials between annotated and unannotated contigs were not substantially different (mean value, 0.399 and 0.324; p-value, 2.2e-16 by the Wilcoxon test; see Supplemental Figure S2 online), which raises the possibility that these are lineage-specific coding genes.

Highly expressed contigs are enriched in developmental and metabolic functions. We examined the correlation between gene expression levels and gene functions. Contigs were grouped into 10 equal-sized bins based on the mean base-level read coverage (refer as ''the coverage'' hereafter) such that transcripts in Bin 1 had the lowest and those in Bin 10 had the highest coverage. The overall coverage of contigs is shown in Figure 1A. Enriched GO terms in each bin (false discovery rate <0.05 after Fisher's exact test) are reported in

Table 2 | Comparison of contigs from the pooled assembly and contigs from assemblies of individual samples

Dataset A Dataset B Number of contigs in A/number of contigs from A with hits in B1 (percentage) Number of contigs in B/number of contigs from B with hits in A (percentage)

Pooled data Pooled data Pooled data Pooled data Swimming 4h 24 h Larva2 240,173/195,417(81.4) 240,173/169,761 (70.7) 240,173/194,967 (81.2) 240,173/27,272 (11.4) 150,683/148,991 (98.9) 116,761/115,203 (98.7) 156,556/154,978 (99) 7,238/6,612(91.4)

*Non-redundant contigs ($200 bp) were used for the comparisons. 'Contigs were assumed to have hits on the other dataset when the blastn bit score was >50 and the minimum length of the 2The 454 reads from a previous study6 were assembled using MIRA7, and non-redundant contigs $200 bp were used. match was 100 bp.

Table 3 | Annotation summary

Number of contigs 240,173

Number of contigs with at least one ortholog in the nr database 64,175

Number of contigs with at least one annotated Pfam domain 30,156

Number of contigs with at least one GO annotation 55,840

Number of contigs with at least one Reactome annotation 26,474

Supplemental data file 8. A summary of the 10 most enriched GO terms in each bin is showed in Figure 1B. Overall, GO terms for developmental functions (e.g., ''Neurogenesis''; ''Cell Differentiation''; ''Instar Larval or Pupal Morphogenesis'') and metabolic processes associated with energy reserve molecules (e.g., ''Phospholipid Metabolic Process''; ''Organophosphate Metabolic Process'') were most frequently over-represented in the high coverage bins Bin 7 to 10 (see Figure 1B). The similar GO annotations such as ''Sensory Perception'', ''Ciliary Motility'', ''Flagellar Motility'' and

Figure 1 | Relationship between gene expression level and gene function. All of the contigs were binned into 10 equal-sized bins according to their baselevel coverage, which was calculated using SAMtools and genomeCoverageBed of BEDTools for each stage and then summed for all three stages. The overall contig coverage and boundaries for Bins 5 to 10 are shown in (A). The top 10 over-represented GO terms in each bin are shown in (B). In the case of duplicated GO, only the term for the highest level of GO hierarchy is listed.

''Microtubule-based Movement'' were enriched in the low coverage bins (Bins 1 to 5, mean coverage below 80)(see Figure 1B). Low expression levels of these genes might reflect the degradation of the larval ciliary field and nervous system during metamorphosis4.

Differential expression analysis reveals changes in energy metabolism and activation of immune system during metamorphosis.

The expression level and the result of differential expression analysis (using GFOLD) from three pair-wise comparisons (SW-to-4 h, 4 h-to-24 h and SW-to-24 h) are listed in Supplemental dataset 2. In each comparison, contigs whose GFOLD values are within the highest (the most up-regulated) and the lowest (the most down-regulated) top 1% and 5% were shortlisted for GO enrichment analysis. The most enriched Biological Process terms in each DEG list are showed in Figure 2. Genes involved in phosphocreatine metabolism were dramatically down-regulated in 4 h and 24 h compared with SW. The GO terms ''Phosphagen Biosynthetic Process'' and ''Phosphocreatine Biosynthesis Process'' were over-represented in the list of top 1% down-regulated contigs in the SW-to-4 h and SW-to-24 h comparisons. Phosphagens and phosphocreatine serve as a reserve of high-energy phosphates that can be metabolized rapidly to generate energy22. The dramatic down-regulation of genes required for phosphocreatine metabolism in 4 h suggested that 1) the level of energy production decreased when the larvae ceased swimming and 2) energy metabolism might be changing from a mode of rapid production to a more gradual, long-term, sustainable mode. This result is consistent with our previous proteomic analysis, in which proteins for energy production were down-regulated in 4 h23.

Contigs associated with immune response functions were highly expressed in 4 h compared with SW and 24 h. The GO terms associated with the immune response, such as ''Inflammatory Response'', ''Phagocytosis'', and ''Positive Regulation Of Macrophage Activation'' (GO:0043032), were over-represented in the top 1% or top 5% up-regulated contig lists in the SW-to-4 h comparisons. In addition, up-regulation of contigs enriched in ''Cellular Response to Nitric Oxide'', ''Hydrogen Peroxide Catabolic Processes'' and ''Chitin Catabolic Processes'' in 4 h compared with SW suggested an increase in chemical cues for macrophage activation. Specifically, chitinase-3-like 1 (CHI3L1), a macrophage genetic marker24, was significantly up-regulated in 4 h compared with SW, suggesting an increase in macrophage population at the 4 h stage compared with SW. The over-representation of genes involved in the ''Cellular Response to Tumor Necrosis Factor'' suggesting secretion of cyto-kines from activated mature macrophages25.

Up-regulation of axial patterning transcription factors and signaling pathway genes. One of the primary objectives of this study is to understand the molecular mechanisms underlying bryozoan lophophore morphogenesis. Consistent with the restructuring process during metamorphosis, the top 1% most highly expressed genes in 4 h and 24 h compared with SW were enriched in ''Regeneration'', ''Tissue Development'', ''Developmental Growth'' and ''Neuron projection development''.

Based on their widely known functions in development, we selected 16 transcription factors (TF) that were dramatically up-regulated (top 1% up-regulated contigs) in 4 h or 24 h compared with SW. Their relative expression levels (presented in RPKM) and differential expression values (presented in GFOLD values) are showed in Table 4. One-third of these up-regulated transcription factors (6 out of 16) are homeobox TFs, and excluding BnISL-1-like, all of them displayed a peak expression level in 4 h. In addition to homeobox TFs, axial patterning fork-head and PAX TFs (BnPAX6, BnSLP1, BnFoxD) were also up-regulated. A number of multifunctional TFs that regulate cell cycle, apoptosis, cellular differentiation and proliferation were also up-regulated in 4 h or 24 h compared with SW, including an E26 transforming-specific protein (BnETS), a

helix-loop-helix protein (BnHLHP1), a basic leucine zipper protein (BnCEBP1) and an estrogen nuclear receptor (BnESR). Up-regulation of BnDSM3, BnDMRTA2 and BnPouP might be related to the development of the ancestrula nervous system because the Pou and DM domains of TFs are implicated in neuroendocrine and neuroblast differentiation26,27.

In addition to the TF analysis, we focused on enrichment of GO terms associated with Wnt, BMP and Notch signaling pathways in our differential expression analysis as these pathways 1) are highly conserved28-30 and 2) have been implicated in regulating metamorphosis of various taxa31'32. Genes that were dramatically up-regulated (top 1%) in 4 h compared with SW are enriched in the regulation of Wnt and BMP signaling pathways but not in the Notch signaling pathway, suggesting that these pathways are highly regulated in bryozoan metamorphosis and implicating them in the metamorphic processes. Genes involved in ''Regulation of the Wnt Signaling Pathway'' were over-represented only in the top 1% up-regulated contig lists in the SW-to-24 h and 4 h-to-24 h comparisons, which suggested that this pathway was very active in 24 h but not in 4 h. Genes involved in ''Regulation of the BMP Signaling Pathway'' were enriched in top 1% up-regulated contig lists in all three pairwise comparisons, suggesting a highly dynamic regulation of this pathway. Notch signaling pathway genes were enriched in the top 5% up-regulated but not in the top 1% up-regulated contig lists in SW-to-24 h comparison.

Next, we examined the expression levels of contigs with functions annotated to the Wnt, BMP, Notch and Hedgehog signaling pathways. The maps for these four signaling pathways were constructed using the Agilent Literature Search Cytoscape plugin33 and are shown in Figure 3. The majority of the genes in these four pathways were (on average) up-regulated in either 4 h or 24 h compared with SW. The expression patterns of BnsFRP, Wnt genes and Bnftcatenin were consistent with the qPCR results obtained in our previous study15 showing the highest expression level of these genes at 4 h. Hedgehog signaling genes (except growth arrest-specific [GAS]) displayed relatively low expression levels in comparison to the other three pathways with average RPKM values less than 2. The downstream targets in the Wnt and Notch pathways, HES, matrix metalloproteinase 7 (BnMMP7) and cyclinD (BnCCND), were up-regulated by more than 4-fold in 24 h compared with SW. In contrast, ID, a downstream target in the BMP pathway, had a relatively high average expression level (above 50 RPKM) with a profile similar to BMPs and peak expression level at 24 h. Overall, the high expression level of genes in these four pathways in 4 h or 24 h suggests their regulatory role in B. neritina metamorphosis.

Pre-patterning of polypide precursors by axial patterning genes.

In situ hybridizations (ISH) of23 genes in the Wnt, BMP, Notch and Hedgehog signaling pathways were performed. These genes included ligands (BnBMP2, BnBMP3, BnBMP5, BnBMP7, BnDLL, BnShh, BnWntA, BnWnt2, BnWnt5, BnWnt7, BnWnt10b), antagonists (BnChrd, BnGAS, BnTLL), receptors (BnNotch, BnPtc), signal transducers/regulators (BnBAMBI, BnGLi, BnNumb, BnsuH, BnSMO, BnWLS) (nodes highlighted in yellow in Figure 3) of the four signaling pathways. We attempted to perform ISH ofthese genes in SW, preancestrulae at 1 h, 2 h, 4 h, 8 h, 16 h and 24 h, as these time points marks the progression of the lophophore morphogenesis. However, successful results were obtained only for SW, 1 h and 2 h, perhaps due to the development of a lightly calcified chitinous shell11, which inhibits penetration of the RNA probes into the samples from 4 h and subsequent time points. The contigs annotated to these 23 genes, primer pairs used for PCR of each gene are reported in Supplemental Table S2.

All of the genes examined, excluding BnWnt10b and BnNotch, were expressed in the apical blastema (see Figure 4). In SW, BnWntA, BnWnt5, BnWLS (Wntless), BnNumb and BnGAS were

Figure 2 | GFOLD value distribution of the contigs. The top 1% and top 5% of the differentially expressed contigs and the 5 most enriched GO terms in each differentially expressed contig lists in the (A) SW-to-4 h, (B) 4 h-to-24 h and (C) SW-to-24 h comparisons. In the case of duplicated GO terms among the top 5 enriched GO terms, only the term for the highest level of GO hierarchy is listed.

Table 4 | Up-regulated transcription factors

RPKM value GFOLD value

Contig_IDs Gene name Pfam_IDs Pfam_symbol Enterz.Gene.ID SW 4 h 24 h SW-to-4 h 4 h-to-24 h SW-to-24 h

Contig_233733 BnCEBPD PF07716 bZIP 2 100376486 7.89 95.8 1 56 3.63 0.48 4.03

Contig_170052 BnDSM3 PF00751 DM 100329077 0.02 1.4 2.95 3.85 0.44 4.64

Contig_148848 BnDMRTA2 PF00751 DM 100125187 0.10 0.46 3.84 1.55 2.51 4.4

Contig_136956 BnPax6 PF00292 PAX 100892794 0.05 0.12 2.38 0 2.61 3.53

Contig_48229 BnSLP1 PF00250 Fork_head 100650555 0.46 2.57 7.61 2.29 1.22 3.57

Contig_121052 BnFoxD PF00250 Fork head 100329067 0.61 5.14 4.57 2.82 -0.14 2.35

Contig 29550 BnETS PF00178 Ets 445717 0.05 0.04 2.39 0 4.58 4.59

Contig_234791 BnPouP PF00157 Pou 8035365 3.05 42.4 10.1 3.83 -2.16 1.45

Contig_170926 BnESR PF00105 zf-C4 100303418 0 1.26 2.69 3.93 0.47 4.73

Contig_130221 BnISL-1-like PF00046 Homeobox 100906507 0.03 0.66 1.88 2.23 0.71 3.52

Contig_79740 BnHP1 PF00046 Homeobox 8573687 1.88 34 19 4.04 -0.9 2.89

Contig_44749 BnMxp PF00046 Homeobox 410648 0.01 2.56 0.567 5.25 -1.84 2.74

Contig_149788 BnHP2 PF00046 Homeobox 9815194 2.51 29.8 16.5 3.5 -0.91 2.34

Contig_101535 BnHP3 PF00046 Homeobox 100643103 0.02 2.95 0.49 5.17 -2.26 2.23

Contig_13990 BnBarH2 PF00046 Homeobox 100021124 0.38 4.17 1.5 3.13 -1.38 1.37

Contig_148496 BnHLHP1 PF00010 HLH 100115092 0.29 1.99 4.9 2.4 0.89 3.42

expressed in the whole apical blastemal (see Figure 4A, 4C & 4D); BnWnt2, BnWnt7, BnBMP2, BnBMP3, BnBMP5, BnBMP7 and BnGLi were expressed in the mesodermal halfofthe apical blastemal (see Figure 4A, 4B &4D); BnWLS, BnDLL and BnsuH were expressed in both the metasomal sac and the whole apical blastemal (see Figure 4A & 4C); BnShh and BnChrd showed a similar expression pattern to the Wnt signaling antagonist BnsFRP (reported by Wong et al.15) and the transcription factor BnGATA456 (reported by Fuch et al.14), which are expressed in a distinct patch of cells in the posterior region of the apical blastema (see Figure 4B & 4C); and BnTLL, BnBAMBI and BnBMP1 were expressed in the posterior half of the apical blastema (see Figure 4B).

At 2 h preancestrulae (2 h), the disc shaped apical blastema transformed into U-shape. The peripheral and the central part of the disc shaped apical blastemal becomes the upper and lower half of the U-shaped apical blastemal. BnWnt7 was expressed in the upper tip of the U-shaped apical blastema (see Figure 4A), displaying the same expression pattern as BnWnt10 (reported by Wong et al.15). The Hedgehog ligand BnShh, and the BMP protein antagonist BnChrd displayed a similar expression pattern to BnsFRP15, which was expressed in the lower tip of the U-shaped apical blastema (see Figure 4B & 4D). The distinctive expression patterns of these genes suggest their potential role in patterning of the apical blastema. Notch signaling genes BnNotch, BnDLL and Bnsu(H) were expressed

Figure 3 | Graphical representation of the contig expression levels mapped to different genes that participate in the WNT (red nodes), BMP (purple nodes), Hedgehog (blue nodes) and NOTCH (green nodes) signaling pathways. The interaction network between these signaling pathways was built based on the literature (Agilent literature search network). The mean expression patterns of the contigs annotated to each gene is shown in a bar chart beside each node. When multiple contigs were annotated to the same gene, the mean expression (RPKM) value was used. Nodes highlighted in yellow indicate additional analysis using in situ hybridization.

Figure 4 | In situ hybridization of 23 genes involved in the (A) Wnt, (B) BMP, (C) NOTCH and (D) Hedgehog signaling pathways. SW: swimming larval stage; 2 h: 2 h preancestrula stage.

in the basal part of the 2 h preancestrulae, but the negative regulator of Notch signaling pathway, BnNumb, was expressed only in the U-shaped apical blastemal (see Figure 4C). The basal part preancestrula of B. neritina composed of numerous longitudinal fibrous cells15. This finding suggests that Notch signaling may be involved in controlling proliferation of those basal fibrous cells. The expression of BnWnt2, BnWnt10b, BnTLL and BnGAS was not detected in 2 h, while BnSMO and BnPtc expression was detected only in 2 h (see Figure 4D). The expression pattern of BnWnt7 was similar to that of BnWnt10, and both were expressed in the upper tip of the U-shaped apical blastema (see Figure 4A).


Considering that there was neither a reference genome nor a tran-scriptome for any bryozoan, we aimed to make the B. neritina tran-scriptome assembly as accurate as possible. The large numbers of contigs in the pooled assembly compared with other works34,35 is the result of using multiple k-mers and a highly stringent threshold in the removal of redundant contigs. Such an approach enabled us to capture as many transcripts as possible. However, it should be noted that

each non-redundant contig does not necessarily represent a full-length gene because transcriptional variants could be generated by alternative transcription start sites, alternative splicing or partial transcripts. B. neritina genes with low expression level or completely not expressed at the selected developmental stages may also not be detected in the present study. Genome sequencing and appropriate scaffolding will eventually provide information on the complete set of B. neritina genes.

Selective degeneration of larval musculature and neurons is common in metamorphosis of lophotrochrozoans (reviewed in36)12,13,36-39. Although metamorphosis of bryozoans, brachiopods and phoronids are all initiated by the respective larval neuromuscu-lar system, the metamorphic fate of the larval musculature and nervous systems are different. In phoronids, larval musculature and nervous systems associated with the apical organs are selectively lost12,36. Similar selective degeneration of larval tissue occurs in the lecithotrophic larvae of inarticulate brachiopod larvae, where only a few interconnecting neurons are lost during metamorphosis10. In metamorphosis of bryozoans, however, the larval ciliary field and the whole larval nervous system, including neurons innervating the \ >


apical organ and the larval ciliary field, equatorial nerve cord and the central nerve ring, are completely degenerated during meta-morphosis13, leaving no trace of larval musculature and nervous system in the ancestrula. Our results suggest that histolysis of bryozoan larval tissues is an apoptotic process, as in the phoronids actinotroch larvae36, and apoptosis of these larval cells may be mediated by the innate immune system. Based on our transcrip-tomic data on B. neritina, we anticipated that during the early preancestrula stage (4 h), elevation of chitosan and reactive oxygen species such as nitric oxide and hydrogen peroxide may promote macrophage production24,40. Macrophages have been demonstrated to play a central role in tissues remodeling during metamorphosis of ascidians41. In B. neritina, mature macrophages may play a major role in the absorption of internalized larval ciliary field during metamorphosis. Cytokines produced and released by mature macrophage may also trigger apoptosis of internalized larval cells. Macrophage-like cells were also present in the larva of Phoronis pallida and were observed to be released from the blood corpuscle mass at the initial stage of metamorphosis36. Taken together, macrophage mediated absorption and apoptosis of larval tissue during metamorphosis may be a common strategy among Lophophorata, or in taxa that undergo catastrophic metamorphosis. Yet, further experimental validation on the activity of macrophage during metamorphosis of both phoronids and bryozoans taxa is required.

Up-regulation and the distinct expression patterns of the Wnt, BMP, Notch and Hedgehog pathway genes in the apical blastema suggest that the developmental roles of these four conserved signaling pathways in metamorphosis of B. neritina. One of the potential roles of these signaling pathway genes may be on regional specification of the apical blastema. In the Xenopus embryo, the BMP4-chordin/noggin expression gradient dictates dorsal-ventral (DV) axis formation, in which BMP4 serves as a ventral signal while its inhibitors noggin and chordin dorsalizes early embryonic tissues31. Wnt signaling regulates anterior-posterior (AP) axis formation, in which Wnt proteins serve as a posterior signal while the antagonists DKK and sFRP regulates the formation of anterior structures28. Hedgehog signaling regulates left-right axis formation42 and acts upstream of the Wnt and BMP signaling pathways in various developmental systems42,43. Although we obtained ISH data of various BMPs and chordin, we did not observe the formation of a BMP-chordin expression gradient. However, we did observe a Wnt-sFRP axis in the blastemic rudiment of polypide. Based on the ontology of polypide and the ISH results on B. neritina from the previous14,15 and current studies, we speculate that BnWnt7, BnWnt10 and BnsFRP may serve as position determinant of the presumptive polypide, in which expression of BnWnt7 and BnWnt10 define the position of 1) the lophophore or 2) the pharynx (month) and expression of BnsFRP may determine the position of the anus. Bryozoan larvae undergo substantial reorganization at the initial stage such that larval corona, apical organ and gut (in cyphonaute larvae) are degenerated11. While some studies inferred an anterioposterior orientation of bryozoan larva based on the respective embryonic development process44,45, there is so far no study inferring a dorsoventral or anterioposterior axis of the zooid of any bryozoan species. Hence, it is difficult to compare the expression pattern of these axial patterning genes in B. neritina to other bilaterians. Interestingly, from a functional perspective, the lophophore/month Wnt and the anus sFRP expression patterns in the apical blastema ofB. neritina larvae is more similar to the oral Wnt and aboral sFRP expression in cnidarians46 than to the posterior Wnt and anterior sFRP expression in the early embryos of bilaterians28. Although, at this point, we are unable to determine the exact roles of these genes based solely on ISH (at SW and 2 h), our result strongly supported the "pre-patterning" mechanism of bryozoan lophophore development14. ISH experiments on these genes at different time points throughout metamorphosis will definitely

extend our understanding on axial patterning in bryozoans and bilaterians.

BnWnt7, BnWnt10, BnShh, BnChrd, and BnsFRP are ligand genes of the Wnt, BMP, and Hedgehog signaling pathways. Their distinct spatial expression patterns in turn trigger regional expression ofTFs. Up-regulation of these ligand genes at metamorphosis coincide with the up-regulation of Pairbox (PAX) class (BnPax6) and Fork-head (FOX) class (BnFoxD) TFs, which are downstream targets of Wnt and BMP signaling pathways and play major roles in regulating neurogenesis47,48. In the brachiopod Terebratalia transversa, Pairbox class and Fork-head class TFs were expressed in the neurons innervating the apical organ throughout the course oflarval development and were therefore being associated with the development of larval nervous system49. The juvenile nervous system ofbrachiopods is modified from the larval nervous system10. In the case of bryozo-ans, however, larval nervous system is completely independent to the juvenile one. Therefore, PAX and FOX class TFs may regulate the development of juvenile nervous system in both brachiopods and bryozoans.

In conclusion, we reported the transcriptomic changes that occurred during B. neritina metamorphosis. Our results 1) revealed major physiological changes in energy metabolism, 2) provided evidence of the involvement of the innate immune system in the restructuring larval tissues, and 3) provided information about the molecular mechanisms underlying lophophore development. Furthermore, ISH validated the expression of developmental genes during B. neritina metamorphosis and supported the presence of a pre-patterning system in regulating lophophore morphogenesis. Future analyses of the expression patterns of HOX genes and different anterior or posterior transcription factors may provide insight into the establishment of the zooidal body axis of bryozoans during metamorphosis as well as differences in the mechanisms responsible for pattern formation compared with other bilaterians.


Collection of larval samples. Adult B. neritina colonies were collected from the floating rafts of a fish farm in Trio Beach, Hong Kong (22°21' 19N, 114°16' 15E) from February to April 2010. Adult colonies were maintained in a flow-through seawater system at 21°C at the Coastal Marine Laboratory (Hong Kong University of Science and Technology) for no more than 7 days before use.

Free-swimming larvae were collected as described previously15. Approximately 15,000 swimming larvae were divided equally into three groups to sample three developmental stages: SW, 4 h and 24 h. Larvae these time points show the most obvious anatomical and morphological differences. For the SW sample, one group of larvae was collected directly on a mesh sieve and transferred to a 2 mL Eppendorf tube. The larvae were quickly washed with autoclaved filtered seawater (AFSW) and stored at — 80°C. For the 4 h and 24 h larval samples, the remaining two groups of larvae were transferred to Petri dishes containing 10 mL of AFSW and allowed to attach and metamorphose in the dark for 30 min, after which the unattached larvae were discarded. The attached larvae were harvested at 4 h and 24 h and stored in liquid nitrogen for RNA extraction or fixed in 4% paraformaldehyde (PFA) for in situ hybridization.

RNA preparation for Illumina HiSeq2000 sequencing. The frozen larval samples were homogenized using a clean mortar and pestle that had been pre-chilled in liquid nitrogen. Three milliliters of TRIzol reagent (Invitrogen, USA) was added to each homogenized larval sample. Total RNA was extracted according to the manufacturer's protocol and then treated with DNase. The integrity of the resultant total RNA was checked by electrophoresis in a denatured agarose gel and by using a NanoDrop 1000 (Thermo Fisher Scientific, USA). For each stage, three batches (biological replicates) of larval samples were collected. Total RNA extraction, DNase treatment and RNA integrity for these 9 samples (3 stages X 3 batches) were performed independently.

Prior to cDNA synthesis, 10 ug of total RNA from each of the three biological replicates of the same stage were pooled together. The mRNA purification, doublestrand cDNA synthesis, cDNA library construction, adaptor ligation, hydroshearing of cDNA fragments and library validation were performed using the Illumina TruSeqTM RNA Sample Preparation Kit v2 according to the high-throughput protocol. Illumina sequencing was performed using a HiSeq2000 to produce 101 bp paired-end data.

Preprocessing and assembly of Illumina short-read data. Low quality ends of reads (first base with a Q-score <20, corresponding to an error probability of 0.01 and up to the 3' end) and sequencing adapters were trimmed using custom java scripts. Each \ > _

dataset (SW, 4 h and 24 h) was assembled by ABySS v1.3.416 using every odd k-mer from 45 to 81 and then merged using Trans-ABySSv1.4.417 ABySS and Trans-ABySS were also used for consolidation of contigs with various ranges of expression levels. Reads shorter than the k-mer were discarded before ABySS assembly. In addition, we pooled all of the datasets (the pooled data) and followed the same assembly strategy to obtain a comprehensive transcriptome (the pooled assembly). Redundant contigs were removed using CD-HIT-EST v4.6 with a sequence identity threshold of 0.99 (see Table 2)18. The pooled assembly was chosen for the downstream analysis (see Tables 1 and 2). The 454 reads from the previous study19 were assembled using MIRA50.

Functional annotation of the assembled contigs. Coding sequences were obtained using getorf by translating contigs between the start and stop codons of all possible six reading frames51. For functional annotation of the B. neritina contigs, two approaches were used. First, homologous proteins were searched using blastp against the nonredundant (nr) database with a bit score of50 and a match length of20 amino acids as the minimum thresholds52. Second, Interproscan v5.4 was used to search for Pfam domain signatures in the contigs53. Gene ontology terms (as of March 15, 2013) and Reactome pathway ID (as of March 10,2013) information of annotated contigs were retrieved from NCBI and Reactome database54,55. Equivalent information was obtained for the Pfam domain from the Interproscan results. The coding potential of the contigs was calculated using the Coding-Potential Assessment Tool (CPAT) v1.221. Coding sequences for Homo sapiens, Drosophila melanogaster, Caenorhabditis elegans, Strongylocentrotus purpuratus, and Saccharomyces cerevisiae were downloaded from the UCSC genome browser56, and then, 10,000 random sequences were selected. The same number of noncoding sequences (>100 bp) was randomly chosen using NONCODE v3.057. These noncoding sequences were used as the training set for CPAT.

Contig expression analysis. Reads determined from each condition were aligned with the B. neritina contigs using Bowtie2 with the sensitive option58. To examine the correlation between expression level and gene functions, we mapped the short reads of each developmental stage to the contigs, calculated base-level coverage of each contig for each stage and then summed the coverage ofall three stages. The base-level coverage was obtained using genomeCoverageBed of BEDTools59. The number of aligned reads per contig was counted using SAMtools60.

Because we did not have replicates, GFOLD was used to calculate differentially expressed genes (DEGs) between the samples61. Significantly over-represented GO terms and Reactome pathways in all of the DEG lists were calculated using a cumulative hypergeometric distribution.

In situ hybridization. Synthesis of antisense RNA probes for in situ hybridization (ISH) and whole-mount ISH were performed as previously described62. The DNA templates for ISH antisense RNA probe synthesis of each gene were prepared by PCR amplification using carefully designed gene-specific primers (GSP). When multiple contigs were annotated to the same homologous gene, the longest contig was chosen to design the GSP. The lists of GSPs and the contigs for their design are listed in Supplemental Table S2. A T7 promoter sequence was added to the 5' end of the reverse GSPs. Digoxigenin-labeled probes were synthesized from PCR product templates according to the protocol supplied with the DIG RNA Labeling Kit (Roche Diagnostics, USA). Larval samples were dehydrated in 100% methanol at —20°C after fixation in 4% PFA at 4°C overnight. Prior to hybridization, the samples were rehydrated in decreasing concentrations of methanol in PBS-0.1% Triton X100 (PBST). The larval samples were incubated with 20 mg/mL of proteinase K (New England BioLabs, USA) to increase their permeability. Swimming larvae were incubated for 7 min, while the two metamorphosing larval samples were incubated for 15 min. The samples were post-fixed in 4% PFA and washed before pre-hybridization at 56°C for 1 h in a hybridization mixture. RNA probe hybridization was conducted in a 56°C water bath with constant orbital shaking overnight. In each hybridization reaction, no more than 50 ng of antisense RNA probe was used. Subsequent washing steps, anti-DIG immunostaining and stain development procedures have been described previously15. The stained samples were visualized using a microscope (Olympus, Tokyo, Japan) with a dark-field Illumination setting.

1. Hyman, L. H. The Invertebrates. [5] (McGraw-Hill, New York, 1959.).

2. NIELSEN, C. The Phylogenetic Position of Entoprocta, Ectoprocta, Phoronida, and Brachiopoda. Integ. Comp. Biol. 42, 685-691 (2002).

3. Halanych, K. M. et al. Evidence from 18S ribosomal DNA that the Lophophorates are protostome animals. Science. 267,1641-1642 (1995).

4. Dunn, C. W. et al. Broad phylogenomic sampling improves resolution of the animal tree of life. Nature 452, 745-750 (2008).

5. Helmkampf, M., Bruchhaus, I. & Hausdorf, B. Phylogenomic analyses of lophophorates (brachiopods, phoronids and bryozoans) confirm the Lophotrochozoa concept. P. Roy. Soc. Lond. B Bio. 275, 1927-1933 (2008).

6. Hausdorf, B., Helmkampf, M., Nesnidal, M. & Bruchhaus, I. Phylogenetic relationships within the lophophorate lineages (Ectoprocta, Brachiopoda and Phoronida). Mol. Phylogenet. Evol. 55,1121-1127 (2010).

7. Struck, T. H. & Fisse, F. Phylogenetic position of Nemertea derived from phylogenomic data. Mol. Biol. Evol. 25, 728-736 (2008).

8. Nesnidal, M. P. et al. New phylogenomic data support the monophyly of

Lophophorata and an Ectoproct-Phoronid clade and indicate that Polyzoa and

Kryptrochozoa are caused by systematic bias. BMC Evol. Biol. 13, 253 (2013).

9. Herrmann, K. Phoronida. In Harrison, F. W. & Woollacott, R. M. (ed.). Microscopical Anatomy of Invertebrates, Lophophorates, Entoprocta, and Cycliophora. [13] (Wiley-Liss, New York, 1997).

10. Santagata, S. Evaluating neurophylogenetic patterns in the larval nervous systems of brachiopods and their evolutionary significance to other bilaterian phyla.

J. Morphol. 272, 1153-1169 (2011).

11. Woollacott, R. M. & Zimmer, R. L. Attachment and metamorphosis of the cheilo-ctenostome bryozoan Bugula neritina (linne). J. Morphol. 134, 351-382 (1971).

12. Temereva, E. & Wanninger, A. Development of the nervous system in Phoronopsis harmeri (Lophotrochozoa, Phoronida) reveals both deuterostome-and trochozoan-like features. BMC Evol. Biol. 12, 121 (2012).

13. Wanninger, A., Koop, D. & Degnan, B. M. Immunocytochemistry and metamorphic fate of the larval nervous system of Triphyllozoon mucronatum (Ectoprocta: Gymnolaemata: Cheilostomata). Zoomorphol. 124,161-170 (2005).

14. Fuchs, J., Martindale, M. Q. & Hejnol, A. Gene expression in bryozoan larvae suggest a fundamental importance of pre-patterned blastemic cells in the bryozoan life-cycle. EvoDevo. 2, 13 (2011).

15. Wong, Y. H., Wang, H., Ravasi, T. & Qian, P. Y. Involvement of Wnt signaling pathways in the metamorphosis of the bryozoan Bugula neritina. Plos one. 7, e33323 (2012).

16. Simpson, J. T. et al. ABySS: a parallel assembler for short read sequence data. Genome Res. 19, 1117-1123 (2009).

17. Robertson, G. etal. De novo assembly and analysis of RNA-seq data. Nat. Methods. 7, 909-912 (2010).

18. Fu, L., Niu, B., Zhu, Z., Wu, S. & Li, W. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 28, 3150-3152 (2012).

19. Wang, H. et al. Rapid transcriptome and proteome profiling of a non-model marine invertebrate, Bugula neritina. Proteomics. 10, 2972-2981 (2010).

20. Luo, C., Tsementzi, D., Kyrpides, N., Read, T. & Konstantinidis, K. T. Direct comparisons of Illumina vs. Roche 454 sequencing technologies on the same microbial community DNA sample. Plos one. 7, e30087. (2012).

21. Wang, L. etal. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 41, e74-e74 (2013).

22. Bogdanis, G. C., Nevill, M. E., Boobis, L. H. & Lakomy, H. K. Contribution of phosphocreatine and aerobic metabolism to energy supply during repeated sprint exercise. J. Appl. Physiol. 80, 876-884 (1996).

23. Zhang, H. et al. Quantitative proteomics identify molecular targets that are crucial in larval settlement and metamorphosis of Bugula neritina. J. Proteome Res. 10, 349-360 (2010).

24. Peluso, G. et al. Chitosan-mediated stimulation of macrophage function. Biomaterials. 15,1215-1220 (1994).

25. Jovanovic, D. V. et al. IL-17 stimulates the production and expression of proinflammatory cytokines, IL-ß and TNF-a, by human macrophages.

J. Immunol. 160, 3513-3521 (1998).

26. Altaba, A. & Melton, D. A. Involvement of the Xenopus homeobox gene Xhox3 in pattern formation along the anterior-posterior axis. Cell. 57, 317-326 (1989).

27. Schonemann, M. D. et al. Development and survival of the endocrine hypothalamus and posterior pituitary gland requires the neuronal POU domain factor Brn-2. Genes Dev. 9, 3122-3135 (1995).

28. Yamaguchi, T. P. Heads or tails: Wnts and anterior-posterior patterning. Curr. Biol. 11, R713-R724 (2001).

29. Guruharsha, K. G., Kankel, M. W. & Artavanis-Tsakonas, S. The Notch signalling system: recent insights into the complexity of a conserved pathway. Nat. Rev. Genet. 13, 654-666 (2012).

30. Miyazono, K., Kusanagi, K. & Inoue, H. Divergence and convergence of TGF-ß/ BMP signaling. J. Cell Physiol. 187, 265-276 (2001).

31. Piccolo, S., Sasai, Y., Lu, B. & De Robertis, E. M. Dorsoventral patterning in Xenopus: inhibition ofventral signals by direct binding ofchordin to BMP-4. Cell 86, 589-598(1996).

32. Athula, W. H. et al. An ancient role for nuclear ß-catenin in the evolution of axial polarity and germ layer segregation. Nature 426, 446-450 (2003).

33. Cline, M. S. et al. Integration of biological networks and gene expression data using Cytoscape. Nat. Protoc. 2, 2366-2382 (2007).

34. Looso, M. et al. A de novo assembly of the newt transcriptome combined with proteomic validation identifies new protein families expressed during tissue regeneration. Genome Biol. 14, R16 (2013).

35. Adamidi, C. et al. De novo assembly and validation of planaria transcriptome by massive parallel sequencing and shotgun proteomics. Genome Res. 21,1193-1200 (2011).

36. Santagata, S. Structure and metamorphic remodeling ofthe larval nervous system and musculature of Phoronis pallida (Phoronida). Evol. Dev. 4, 28-42 (2002).

37. Barlow, L. A. & Truman, J. W. Patterns of serotonin and scp immunoreactivity during metamorphosis of the nervous system of the red abalone, Haliotis rufescens. J. Neurobiol. 23, 829-844 (1992).

38. Marois, R. & Carew, T. J. The gastropod nervous system in metamorphosis. J. Neurobiol. 21, 1053-1071 (1990).

39. Wanninger, A. Ruthensteiner, B., Dictus, W. J. A. G. & Haszprunar, G. The development ofthe musculature in the limpet Patella with implications on its role in the process of ontogenetic torsion. Invertebr. Reprod. Dev. 36,211-215 (1999).

40. Bogdan, C., Rollinghoff, M. & Diefenbach, A. Reactive oxygen and reactive nitrogen intermediates in innate and specific immunity. Curr. Opin. Immunol. 12, 64-76 (2000).

41. Davidson, B. & Swalla, B. J. A molecular analysis of ascidian metamorphosis reveals activation of an innate immune response. Development 129, 4739-4751 (2002).

42. Tanaka, Y., Okada, Y. & Hirokawa, N. FGF-induced vesicular release of Sonic Hedgehog and retinoic acid in leftward nodal flow is critical for left-right determination. Nature 435, 172-177 (2005).

43. Farzana, L. & Brown, S. J. Hedgehog signaling pathway function conserved in Tribolium segmentation. Dev. Gene Evol. 218, 181-192 (2008).

44. Gruhl, A. Muscular systems in gymnolaemate bryozoan larvae (Bryozoa: Gymnolaemata). Zoomorphol. 127, 143-159 (2008).

45. Nielsen, C. & Worsaae, K. Structure and occurrence of cyphonautes larvae (Bryozoa, Ectoprocta). J. Morphol 271, 1094-1109 (2010).

46. Plickert, G., Jacoby, V., Frank, U., Miiller, W. A. & Mokady, O. Wnt signaling in hydroid development: formation of the primary body axis in embryogenesis and its subsequent patterning. Dev. Biol. 298, 368-378 (2006).

47. Maekawa, M. et al. Pax6 is required for production and maintenance of progenitor cells in postnatal hippocampal neurogenesis. Genes to Cells. 10,1001-1014 (2005).

48. Timmer, J. R., Wang, C. & Niswander, L. BMP signaling patterns the dorsal and intermediate neural tube via regulation of homeobox and helix-loop-helix transcription factors. Development. 129, 2459-2472 (2002).

49. Santagata, S., Resh, C., Hejnol, A., Martindale, M. Q. & Passamaneck, Y. J. Development ofthe larval anterior neurogenic domains ofTerebratalia transversa (Brachiopoda) provides insights into the diversification oflarval apical organs and the spiralian nervous system. EvoDevo 3, 1-21 (2012).

50. Chevreux, B. et al. Using the miraEST assembler for reliable and automated mRNA transcript assembly and SNP detection in sequenced ESTs. Genome Res. 14, 1147-1159 (2004).

51. Rice, P.,Longden, I. &Bleasby, A. EMBOSS: the European molecular biology open software suite. Trends Genet. 16, 276-277 (2000).

52. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403-410 (1990).

53. Mulder, N. & Apweiler, R. InterPro and InterProScan: tools for protein sequence classification and comparison. Methods Mol. Biol. 396, 59-70 (2007).

54. Ashburner, M. etal. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25, 25-29 (2000).

55. Croft, D. et al. Reactome: a database of reactions, pathways and biological processes. Nucleic Acids Res. 39, D691-697 (2011).

56. Kent, W. J. etal. The human genome browser at UCSC. Genome Res. 12,996-1006 (2002).

57. Bu, D. et al. NONCODE v3.0: integrative annotation of long noncoding RNAs. Nucleic Acids Res. 40, D210-215 (2012).

58. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods. 9, 357-359 (2012).

59. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 26, 841-842 (2010).

60. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 25, 2078-2079 (2009).

61. Feng, J. et al. GFOLD: a generalized fold change for ranking differentially expressed genes from RNA-seq data. Bioinformatics. 28, 2782-2788 (2012).

62. Thisse, C. & Thisse, B. High-resolution in situ hybridization to whole-mount zebrafish embryos. Nat. Protoc. 3, 59-69 (2007).


This study was supported by a grant (DY125-15-T-02) from the China Ocean Mineral Resources Research and Development Association, an award from King Abdullah University of Science and Technology [SA-C0040/UK-C0016], and grants from the Research Grants Council of the Hong Kong Special Administrative Region [GRF661611 and GRF662413] to Pei-Yuan Qian.

Author contributions

Y.H.W., T.R.1 (Taewoo Ryu), P.Q. andT.R.2 (Timothy Ravasi) conceived and designed the project, Y.H.W. collected samples, T.R.1, L.S., Y.G., Y.H.W. and S.B. performed the bioinformatics analysis, YHW performed the in situ hybridization experiments and Y.H.W. and T.R.1 wrote the paper. P.Q. and T.R.2 supervised the work.

Additional information

Accession codes: Accessions numbers of raw Illumina paired-end reads of the examined

stages: SAMN02724736 (for SW), SAMN02724737 (for 4 h) and SAMN02724738 (for 24 h).

Supplementary information accompanies this paper at scientificreports

Competing financial interests: The authors declare no competing financial interests.

How to cite this article: Wong, Y.H. et al. Transcriptome analysis elucidates key developmental components of bryozoan lophophore development. Sci. Rep. 4, 6534; DOI:10.1038/srep06534 (2014).

This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder in order to reproduce the material. To view a copy of this license, visit http://