Developmental Progression in the Coral Acropora digitifera Is Controlled by Differential Expression of Distinct Regulatory Gene Networks.

Corals belong to the most basal class of the Phylum Cnidaria, which is considered the sister group of bilaterian animals, and thus have become an emerging model to study the evolution of developmental mechanisms. Although cell renewal, differentiation, and maintenance of pluripotency are cellular events shared by multicellular animals, the cellular basis of these fundamental biological processes are still poorly understood. To understand how changes in gene expression regulate morphogenetic transitions at the base of the eumetazoa, we performed quantitative RNA-seq analysis duringAcropora digitifera's development. We collected embryonic, larval, and adult samples to characterize stage-specific transcription profiles, as well as broad expression patterns. Transcription profiles reconstructed development revealing two main expression clusters. The first cluster grouped blastula and gastrula and the second grouped subsequent developmental time points. Consistently, we observed clear differences in gene expression between early and late developmental transitions, with higher numbers of differentially expressed genes and fold changes around gastrulation. Furthermore, we identified three coexpression clusters that represented discrete gene expression patterns. During early transitions, transcriptional networks seemed to regulate cellular fate and morphogenesis of the larval body. In late transitions, these networks seemed to play important roles preparing planulae for switch in lifestyle and regulation of adult processes. Although developmental progression inA. digitiferais regulated to some extent by differential coexpression of well-defined gene networks, stage-specific transcription profiles appear to be independent entities. While negative regulation of transcription is predominant in early development, cell differentiation was upregulated in larval and adult stages.

Introduction

Corals belong to the early branching metazoan phylum Cnidaria, which is characterized by diverse life cycles and exceptional regeneration capacity (Steele et al. 2011). This morphogenetic plasticity is a characteristic commonly observed across the lineage; hence, it is likely an ancestral trait present in the last common ancestor of all phylum members (Holstein et al. 2003). Cnidarians are considered the sister group to bilaterian animals (Collins 1998; Medina et al. 2001), and for this reason phylum members have become emerging model organisms to study evolution of developmental mechanisms (Bosch 2007, 2009; Bode 2009). Moreover, corals belong to the most basal cnidarian class, the Anthozoa (Bridge et al. 1992, 1995), making them ideal models to study conservation of developmental processes during animal evolution. There is evidence demonstrating that conserved mechanisms regulate development among all metazoans (Artavanis-Tsakonas et al. 1999; Extavour and Akam 2003; Kusserow et al. 2005; Ingham et al. 2011; Gleason et al. 2014) and it has been suggested that variations from an “ancestral scheme” originated the diversity of body plans observed throughout the animal kingdom (Shubin et al. 2009). Nonetheless, although cell renewal, differentiation, and maintenance of pluripotency are basic biological processes (BPs) shared by multicellular animals, the molecular mechanisms regulating these fundamental processes are still poorly understood.

At the molecular level, clonal reproduction, regeneration, and morphogenesis require that specific cell populations maintain a genetic memory that encodes their pluripotency, while allowing discrete differentiation into specific phenotypes (Chambers and Tomlinson 2009; Thomson et al. 2011). This is achieved via transcriptional networks and interconnected protein–protein and protein–DNA interactions that regulate chromatin structure and gene expression (Meshorer and Misteli 2006; Meshorer et al. 2006). These core networks promote the expression of “pluripotency” genes while repressing the expression of canonical developmental signaling pathways (Orkin and Hochedlinger 2011). Therefore they likely reflect cellular mechanisms underlying morphogenetic plasticity in basal metazoans. Conserved chromatin remodeling factors, such as DNA methyltransferases, histone modifiers, and polycomb proteins, as well as components of canonical developmental pathways, have been identified in Acropora digitifera (Shinzato et al. 2011) and thus likely play important regulatory roles in mediating cellular plasticity and tissue morphogenesis in Acropora.

In order to understand how changes in gene expression regulate morphogenetic transitions at the base of eumetazoa, we performed quantitative RNA-seq analysis during development of the scleractinian coral, A. digitifera. We collected embryonic (blastula, gastrula), larval (sphere, planula) and adult samples to characterize stage specific-transcription profiles, as well as broad expression patterns during developmental progression. Acropora digitifera releases gametes into the water and the first cleavage occurs approximately 2 h (26 °C) postfertilization (HPF), reaching the blastula stage 10–12 HPF (Okubo and Motokawa 2007). Following gastrulation (22–36 HPF) embryos develop to reach a round, motile, larval stage known as a sphere (36–48 HPF), where epithelial cell diversification occurs (Ball et al. 2002; Okubo and Motokawa 2007; Reyes-Bermudez and Miller 2009). At this time, larvae begin to manifest a progressive elongation along the oral/aboral axis until they acquire characteristic planula morphology (48–96 HPF) and abandon the water column to settle on the substrate. Following settlement, planulae metamorphose into primary polyps that originate new colonies (Ball et al. 2002; Okubo and Motokawa 2007; Reyes-Bermudez and Miller 2009) (fig. 1).

F

Coral development. Libraries representing blastula (PC), gastrula (G), postgastrula (S), planula (P), and adult polyps (A) were sequenced in triplicate. Following fertilization developing embryos experience a series of asymmetrical cell divisions that continue until they reach a very distinct blastula stage characteristic of “complexa” claded corals (10–12 h) known as the prawn-chip (PC). Morphogenetic movements during gastrulation (22–36 h) originates ectodermal and endodermal tissues, the blastopore becomes the oral pore (G). Dipoblastic larvae become motile resembling a rotating sphere (36–48 h) at this time cell differentiation of tissue specific lineages begins (S). Progressive elongation of the oral/aboral axis occurs (48–96 h) until larvae acquire the characteristic planula morphology (P). At this stage and under appropriate settlement clues planulae attaches to the substrate, metamorphose into an axial polyp that will originate a new colony (A). O/A within imagines represents the oral/aboral axis. *Presence/absence of Symbiodinium sp.

Although high-throughput quantitative studies have been conducted on corals, most of them have focused on the effect of environmental stressors on coral-specific processes (Meyer et al. 2011; Moya et al. 2012). Gene expression studies of coral development have been restricted to a handful of papers using microarray technology that focused mostly on transcriptional changes occurring during metamorphosis (Grasso et al. 2008; Reyes-Bermudez et al. 2009; Grasso et al. 2011). RNA-seq methods have been used in other cnidarians to study transcript abundance during early stages of development in the related sea anemone Nematostella vectensis (Helm et al. 2013) and differential gene expression between polymorphic polyp types in hydrozoans (Plachetzki et al. 2014), but none have centred on transcription dynamics underlying developmental transitions in corals. In this article, we focused for the first time on changes in regulatory gene networks underlying chromatin structure, cell differentiation, and pluripotency during the A. digitifera life cycle.

Here, we report changes in gene expression between consecutive developmental time points, as well as gene coexpression patterns with emphasis on the molecular mechanisms underlying chromatin regulation, cell differentiation, and morphogenesis. Our results suggest roles for both long noncoding RNA (lncRNA) and coral-specific transcripts during development, and revealed a highly plastic genome able to regulate specific transcriptional circuits at each developmental time point. Likewise, our results identified clear differences in gene expression between early and late developmental transitions and indicate that developmental progression and tissue plasticity in corals are regulated to some extent by differential coexpression of well-defined gene networks. Despite this, stage-specific transcription profiles appear to be independent entities with distinct molecular contexts.

Materials and Methods

Collection of Samples and RNA Extraction

Early coral life history stages from the branching coral, A. digitifera, were raised and collected at the Sesoko island research station, Okinawa, Japan, in 2012, during the annual (June–July) spawning event. Gametes from six colonies were mixed together in six different containers for 2 h until first cleavage was observed. Mixtures were done taking care that: 1) sperm concentration was in the range of 104 and 107 sperms ml−1 (Chui et al. 2014), and 2) that gametes from each colony were represented in each container to ensure biological diversity in the crosses. Developing embryos were maintained in fresh, filtered seawater (1 μm) at approximately 26 °C until they reached the desired developmental stage. Replicates for each stage were collected from different rearing culture vessels. Batches of embryos (∼500) were concentrated in 2 ml cryotubes and seawater was removed using a glass pipet. At this point 1 ml of Trizol was added to each tube and snap frozen in liquid nitrogen for subsequent RNA extraction. Total RNA was extracted from the following key development stages: 1) blastula “prawnchip” (PC) 12 HPF, 2) gastrula (G) 24 HPF, 3) postgastrula “sphere” (S) 48 HPF, 4) planula (P) 96 HPF, and 5) adult colonies (A) (fig. 1). Fragments of branches (∼2 cm) were collected from six different adult colonies and snap frozen to sample adult polyps. Colonies used to collect adult samples were different from the ones used for reproduction.

Total RNA from frozen coral tissues (three samples per stage) was isolated using Trizol lysis reagent (Invitrogen) following product specifications. Modifications from the original protocol, aimed at optimizing RNA extraction from coral tissues were previously described in Reyes-Bermudez et al. (2009). Briefly, two chloroform extractions were performed, followed by isopropanol precipitation, and two washes in 80% ethanol. Pellets were redissolved in 10–15 μl of nuclease-free water. To remove genomic DNA, total RNA samples were DNase I (Invitrogen) treated according to manufacture’s specifications and total RNA was resuspended again in 10 μl of nuclease-free water using 7.5 M LiCl “RNA Precipitation Solution” (Ambiont). RNA quality and integrity were assessed with a Nanodrop ND-1000 spectrophotometer (DNA/RNA ratios) and an Agilent 2100 Bio-analyzer, respectively (ribosomal ratio/RIN number).

Sequencing and Data Analysis

Libraries for sequencing were prepared from total RNA using the protocol described in Aird et al. (2013). This protocol involves the capture of polyadenylated RNA transcripts, followed by template switching at the 3′-end and RACE-PCR amplification of the resulting cDNA (Clontech SMART RACE cDNA Amplification Kit), and then using illumina’s Nextera library preparation kit. This approach preferentially enriches RNAs that undergo posttranscriptional processing, such as polyadenylation and 5′-capping (Harbers et al. 2013), such as mRNAs and lncRNAs. Libraries were sequenced on the illumina GAIIx platform in paired end 50 bp mode. Library construction methods were previously validated using RNA-seq in Aird et al. (2013) and using spike-ins in Aird et al. (2015).

Raw reads were quality trimmed using trimmomatic (v 0.32, Bolger et al. 2014), and mapped to the publically available A. digitifera genome (Shinzato et al. 2011) using Tophat v2 (with the –b2-very-sensitive option, Trapnell et al. 2012). Transcripts were assembled using Cufflinks (default settings; v 2.2.1, Pollier et al. 2013) and annotated using the Shinzato et al. gene models (Shinzato et al. 2011). Assembled transcripts were extracted, and reads were remapped using the RSEM (1.2.1) pipeline (– bowtie-n 3 –paired-end, Li and Dewey 2011) to estimate expected counts transcript fragment counts for each gene. Raw data for all libraries were submitted to the DNA Data Bank of Japan under the bio-project [ID PRJDB3244] and bio-sample IDs [Blastula: SAMD00021035, Gastrula: SAMD00021036, Sphere: SAMD00021038, Planula: SAMD00021037, and Adult: SAMD00021034].

Transcripts were annotated according to the predicted proteome of the coral A. digitifera, (Shinzato et al. 2011) using a local executable copy of BLAST+ v2.2.29 (Camacho et al. 2009). To assess the coding potential of the 6,316 transcripts not mapped to the predicted proteome, we used Coding Potential Calculator (CPC) software with default parameters (Kong et al. 2007). From these, 1,180 transcripts classified as “coding” by CPC were kept for further analysis. Putative peptides provided by CPC for each coding transcript were first matched using BLAST to the A. digitifera proteome. Coding transcripts not represented in the A. digitifera proteome were searched against a custom database containing Nematostella sp., Hydra sp., Danio rerio, Clytia sp., Rattus norvegicus, and Homo sapiens protein sequences consisting of 274,358 entries. Putatively noncoding transcripts were screened using BLASTN against A. digitifera genomic rDNA and significant hits were filtered out from the data set.

Gene Ontology and KEGG Enrichment

The Gene Ontology (GO) and KEGG pathways annotations were constructed using the KEGG orthology-based annotation of A. digitifera (Dunlap et al. 2013) in conjunction with Uniprot database references for KEGG orthologs, release 2014_03. GO enrichment analysis was performed using the GOHyperG function in the GOstats R package (Falcon and Gentleman 2007). The KEGG pathway enrichment was performed using KEGGREST R package (Tenenbaum 2014), and in-house script performing a hypergeometric test (0.05 cutoff value) considering “one versus all” conditions.

EdgeR

Differential gene expression analysis was inferred from the mapped counts using the edgeR R package (Robinson et al. 2010; McCarthy et al. 2012). We filtered out poorly expressed tags using the “filtered_R” function (genefilter package) defining the best quantile value (0.05) for rejection of low abundance tags. To characterize and define transcription profiles for each developmental stage we used the “cpm” function and kept genes with at least 100 counts per million in all replicates (in each one, Robinson et al. 2010) and averaged the remaining transcripts counts. Data were normalized and the common and pairwise dispersion was calculated using generalized linear models (GLM). The GLM model was used to specify probability distributions according to their mean–variance relationship. Top differentially expressed genes (DEGs) were selected with a threshold of P value ≤0.05.

Weighted Gene Coexpression Analysis

We conducted a Weighted Gene Coexpression Analysis (WGCNA) using the R package version WGCNA (v 1.36, Langfelder and Horvath 2008). The analysis was done using fragments per kilobase mapped, which were subjected to a variance stabilizing transformation using the package DESeq2 (Love et al. 2014). We used the “pickSoftThreshold” function to explore soft thresholds from 12 to 46, ultimately choosing a value of 24, which corresponds to an acceptable R2 (>0.8). WGCNA was conducted using signed networks, with a minimal module size set to 30 genes. Module eigengenes were created using default parameters (variation cutoff = 1.0) and merged clusters were formed using a 0.80 similarity. We then conducted module-trait correlations between the module eigengenes, and libraries corresponding to each of the developmental stages. Given the large number of comparisons, we adjusted the P values using FDR correction at a 0.05 family wise significance threshold. GO enrichment analysis was conducted for each module whose expression was significantly correlated with a particular developmental stage.

Results

Sequencing, Mapping and Transcript Abundance

In total, 15 libraries representing blastula (PC), gastrula (G), postgastrula (S), planula (P), and adult polyp (A) stages were sequenced, with three biological replicates each (fig. 1). Illumina results generated 52 Mb (52 million pair reads) of raw data from which 50 Mb (50 million pair reads) passed quality filters and adaptor trimming. High-quality reads were mapped to the publically available A. digitifera proteome and genome data sets (Shinzato et al. 2011) and then assembled into 18,264 unique transcripts. While 69% of all transcripts (12,587) were clearly represented in the predicted Acropora proteome, the remaining 31% (5,677 transcripts) mapped only to the genome. From the latter, 91% (5,187) were identified as putative noncoding transcripts by using CPC software (Kong et al. 2007). Noncoding transcripts were screened using BLASTN against A. digitifera genomic rDNA sequences to validate them as putative lncRNAs. From these, ten transcripts had matches to rDNA and thus removed from the data set. As we captured only poly-adenylated RNA, and the absence of rDNA was confirmed by BLAST, it is reasonable to think that the noncoding transcripts reported in this study are lncRNA. On the other hand, the residual 9% (490) were identified as coding transcripts not represented in the published protein data set.

Stage-specific transcription profiles resulted in 12,223 transcripts expressed in blastula (PC), 13,287 in gastrula (G), 14,148 in postgastrula (S), 11,468 in planula (P), and 11,926 in adults (A). For all stages, coding molecules represented 84–88% of all transcripts and lncRNA the remaining 12–16%. While sequences with known orthologs in other systems represented 71–75% of all transcripts, approximately 13–14% of all transcriptomes are likely to be coral-specific. Transcript abundance in all data sets varied by 3–4 orders of magnitude (fig. 2A) and showed two main tendencies: 1) annotated sequences were slightly more abundant in the upper ranks of the distribution (Q3–4) and 2) 50–55% of all noncoding sequences concentrated in Q1 (fig. 2B). The original count matrix containing number of reads per each stage and replicates can be found in supplementary file S1, Supplementary Material online. Filtered count matrices showing expression levels, coding status, annotation description, and descriptive statistical parameters for each transcriptome can be found in supplementary file S2, Supplementary Material online.

F

Stage-specific trancriptomes. To characterize stage-specific transcription profiles we kept transcripts with at least 100 counts per million in all replicates and averaged the remaining transcripts counts. Expression levels displayed differences of 3–4 orders of magnitude between minimum and maximum values (A). In all cases, while annotated sequences were slightly more abundant in the upper ranks of the distribution, noncoding transcripts concentrated in the low expression quartile (Q1) (B). Stage-specific transcription profiles reconstructed coral development, revealing two main expression clusters. The first cluster grouped PC and G and the second grouped subsequent developmental stages. S was the most distinct stage among all developmental time points (C).

Stage-Specific Transcriptomes Show Overrepresentation of Distinct Molecular Functions

To identify changes in key BPs and/or molecular functions (MFs) underlying developmental progression, we performed GO enrichment analysis of complete stage-specific transcriptomes. We identified stage-specific enriched MFs that reflected cellular complexity and molecular context at each developmental time point. For example, glutathione binding (GO:0043295) was overrepresented only in PC, GDP-binding (GO:0019003) did so in G, phosphogluconate dehydrogenase activity (GO:0004616) was the only MF category uniquely enriched in P and histone deacetylase (GO:0004407) and calmodulin (CaM)-dependent protein kinase (GO:0004683) were overrepresented only in A.

Interestingly, while repressing transcription factor (TF) binding (GO:0070491) was over represented in all stages but S, TF binding (GO:0008134) was enriched only in this stage indicating the usage of very distinct gene networks at this developmental time point. The fact that a number of MF categories with developmental regulatory roles such as beta-catenin binding (GO:0008013, PC and S) and GTPase regulator activity (GO:0005083) were coenriched in more than one stage, suggests coexpression of similar transcriptional networks during development. A similarity matrix based on stage-specific transcription profiles reconstructed development and revealed two main clusters. The first one grouped early developmental stages (PC and G) and the second grouped the stages following gastrulation (S, P, and A). This matrix postulates S as the most distinct transcriptome, as P and A grouped together within this cluster (fig. 2C). A summary of GO terms mentioned in this section can be found in table 1.

Table 1

GO Enrichment Summary in Complete Data Sets and Consecutive Stages Developmental Progressions

	Stagea	GO ID	Node Size	Sample Match	P Adj	Term	Ontology
General GO enrichment	PC	GO:0043295	36	29	3.82.E−02	Glutathione binding	MF
	G	GO:0019003	97	67	2.81.E−02	GDP binding	MF
	S	GO:0008134	1105	641	3.21.E−05	Transcription factor binding	MF
	P	GO:0004616	25	22	7.30.E−03	Phosphogluconate dehydrogenase (decarboxylating) activity	MF
	A	GO:0004407	27	23	4.82.E−02	Histone deacetylase activity	MF
		GO:0004683	406	223	3.64.E−02	Calmodulin-dependent protein kinase activity
	PC_A	GO:0005083	586	317	7.20.E−03	Small gtpase regulator activity	MF
	PC_S	GO:0008013	155	102	1.58.E−02	Beta-catenin binding	MF
	PC_G_P_A	GO:0070491	108	301	2.51.E−03	Repressing transcription factor binding	MF
	PC_S_P_A	GO:0007030	404	263	3.33.E−12	Golgi organization	BP
	All	GO:0030258	677	416	9.48.E−13	Lipid modification	BP
		GO:0006907	2,415	1,469	3.90.E−51	Pinocytosis
		GO:0006396	28,24	1,669	1.91.E−40	RNA processing
aEnriched in stages shown in column. Node size = Total number of Go terms in node. Sample match = number of transcripts with GO terms associated to specific nodes.

Enriched in stages shown in column. Node size= Total number of GO terms in node. Sample match = number of transcripts with GO terms associated to specific nodes. GO categories shown in this table were selected based on: 1) significant P values (0.05 cutoff value) and 2) evidence in the literature of their involvement in the regulation of developmental processes. A complete list of enriched GO categories can be found in supplementary file S4, Supplementary Material online.

Distinct Gene Expression Patterns Were Identified between Early and Late Developmental Transitions

To understand transcription dynamics throughout development we identified DEGs (P value ≤0.05) between consecutive stages. We performed pairwise comparisons based on developmental progression and searched for transcripts that were differentially expressed between data sets. Overall common dispersion was 0.0381 and the biological coefficient of variation was 0.1919. To classify DEGs based on their expression levels, we selected all significant DEGs in each developmental transition and then identified the subset of transcripts that showed log-fold changes (FCs) ≥1.5. We selected this value, which differs from the popular and arbitrary cutoff of ≥ 2 as we consider that FC values of ≥2 would overestimate the number of low expressed DEGs (Dalman et al. 2012).

We identified a peak of DEGs and FC’s during early development, especially in the G to S progression. Although the transition P to A spans settlement and metamorphosis, this progression reported the lowest number of DEGs as well as lowest FC values (fig. 3). Interestingly, lncRNAs were slightly more abundant during transitions involving gastrulation, indicating roles for lncRNA during early development (fig. 3A). A summary of all DEGs and GO terms reported in this study can be found in supplementary files S3 and S4, Supplementary Material online, respectively. GO terms and DEGs of developmental relevance mentioned in the discussion section can be found in tables 1 and 2, respectively.

F

Differential gene expression during developmental progressions. A peak of DEGs was identified during early developmental transitions especially in the G to S progression. This transition also showed the highest FC’s across comparisons. Annotated transcripts were more abundant during late developmental transitions and noncoding DEGs more abundant in the comparisons involving gastrulation (A and B). Nonparametric regressions (LOWESS) identified a tendency for medium and high abundant DEGs in PC and S to have lower FCs relative to G and P, respectively. Yet, medium and high abundant DEGs in G and P tended to increase FCs relative to S and A, respectively (B). Although the transition P to A spans settlement and metamorphosis, this progression reported the lowest number of DEGs as well as the lowest FC values (A and B). Overall, we observed clear differences in gene expression between early and late developmental transitions, with higher numbers of DEGs and FCs around gastrulation.

Table 2

Summary of DEGs of Developmental Relevance Identified During Consecutive Developmental Progressions

Comparison	Up ina	ID	FC	P value	Definition
PC vs. G	PC	adi_v1.09698	5.17	1.25. E − 08	Cadherin EGF LAG seven-pass G-type receptor 2
		adi_v1.01477	1.16	3.91. E − 05	WNT inhibitory factor 1
		adi_v1.09632	2.46	2.06. E − 28	Polycomb protein SCMH1
		adi_v1.08662	2.45	4.99. E − 06	Cadherin EGF LAG seven-pass G-type receptor 1
		adi_v1.04586	2.13	1.24. E − 18	Groucho
		adi_v1.00777	6.06	2.69. E − 06	Homeobox protein MOX
		adi_v1.05704	4.31	1.32. E − 08	Homeobox protein OTX
		adi_v1.06807	4.1	1.69. E − 03	SWI/SNF-related matrix-associated
		adi_v1.24238	3.15	2.87. E − 38	SWI/SNF-related matrix-associated
		adi_v1.04163	4.93	1.27. E − 04	BTB/POZ domain-containing protein 3/6
		adi_v1.03255	2.29	1.52. E − 08	BTB/POZ domain-containing protein 9
	G	adi_v1.02467	−0.77	1.64. E − 02	POU domain transcription factor, class 3
		adi_v1.01257	−2.42	2.19. E − 20	POU domain transcription factor, class 3
		adi_v1.07955	−0.79	2.33. E − 03	Polycomb protein EED
		adi_v1.22403	−1	7.30. E − 09	Transcription factor YY
		adi_v1.05095	−2.02	2.03. E − 30	Krueppel-like factor 8/12
		adi_v1.05096	−4	6.58. E − 50	Kueppel-like factor 5
		adi_v1.16694	−3.31	4.29. E − 08	Forkhead box protein C
		adi_v1.24494	−3.42	1.26. E − 27	Forkhead box protein J1
		adi_v1.22785	−3.77	1.26. E − 18	N-Myc proto-oncogene protein
		adi_v1.22791	−4.83	7.54. E − 45	Myc proto-oncogene protein
		adi_v1.00241	−3.87	1.25. E − 30	Protein sprouty homolog 1
		adi_v1.00237	−7.2	6.22. E − 55	Protein sprouty homolog 1
		adi_v1.10004	−4.2	4.70. E − 04	Collagen, type I/II/III/V/XI/XXIV/XXVII, alpha
		adi_v1.09766	−4.96	2.99. E − 05	Collagen, type I/II/III/V/XI/XXIV/XXVII, alpha
		adi_v1.07457	−5.64	3.38. E − 87	Homeobox protein DLX, invertebrate
		adi_v1.10929	−8.3	9.30. E − 34	Homeobox protein GSH
		adi_v1.04989	−1.34	2.33. E − 05	BTB/POZ domain
		adi_v1.02105	−2.42	4.91. E − 23	BTB/POZ domain (germ cell-less protein-like 1)
		adi_v1.12091	−7.52	2.80. E − 105	Hairy and enhancer of split, invertebrate bHLH
		adi_v1.14589	−6.95	4.61. E − 184	Hairy and enhancer of split, invertebrate bHLH
		adi_v1.00267	−2.36	9.53. E − 14	SOX group C
		adi_v1.03401	−3.17	6.23. E − 47	SOX group E/F
		adi_v1.11949	−3.93	5.70. E − 17	Jagged-1
		adi_v1.20013	−3.87	4.41. E − 41	Jagged-2
		adi_v1.01978	−3.65	2.31. E − 09	Notch-1
		adi_v1.14878	−5.57	1.08. E − 05	Notch-like
		adi_v1.08519	−3.42	1.06. E − 48	Snail, invertebrate
		adi_v1.11963	−2.7	1.21. E − 36	Snail, invertebrate
G vs. S	G	adi_v1.09989	7.06	7.56. E − 11	Chromodomain helicase DNA binding protein 8 [EC:3.6.4.12]
		adi_v1.15318	1.09	2.53. E − 07	Chromodomain-helicase-DNA-binding protein 3 [EC:3.6.4.12]
		adi_v1.07285	7.88	1.41. E − 18	Chromobox protein 1_polycomb family
		adi_v1.16661	4.08	5.19. E − 34	Chromatin modification-related protein EAF7
		adi_v1.16299	3.14	1.56. E − 47	Histone acetyltransferase HTATIP [EC:2.3.1.48]
		adi_v1.12821	5.38	1.92. E − 91	BTB/POZ domain-containing protein 3/6
		adi_v1.12822	5.93	1.64. E − 97	BTB/POZ domain-containing protein 3/6
		adi_v1.13155	4.03	5.75. E − 29	WNT inhibitory factor 1
		adi_v1.21181	1.16	1.29. E − 06	WNT inhibitory factor 1
		adi_v1.01479	5.62	6.26. E − 65	Dickkopf
		adi_v1.08519	2.54	5.04. E − 29	Snail, invertebrate
		adi_v1.11963	2.2	2.37. E − 22	Snail, invertebrate
		adi_v1.20687	1.89	1.47. E − 06	Brachyury protein-like
		adi_v1.21035	2.16	2.14. E − 02	Forkhead box protein L
		adi_v1.17538	2.06	3.59. E − 28	Forkhead box protein N
		adi_v1.14488	3.05	6.75. E − 08	Forkhead box protein P
		adi_v1.12738	−2.65	6.42. E − 28	Chromatin modification-related protein YNG2
		adi_v1.12357	−4.46	5.16. E − 07	Histone deacetylase 6/10 [EC:3.5.1.98]
		adi_v1.06801	−6.33	1.11. E − 34	Chromobox protein 6-polycomb family
		adi_v1.06066	−5.81	5.62 E − 39	Fibroblast growth factor receptor
		adi_v1.09510	−4.75	8.98 E − 05	Fibroblast growth factor receptor
		adi_v1.07835	−4.35	5.00 E − 50	Fibroblast growth factor
		adi_v1.00239	−3.4	6.55 E − 10	Fibroblast growth factor
		adi_v1.20515	−2.81	1.46 E − 05	Fibroblast growth factor receptor
		adi_v1.09253	−6.42	2.80. E − 71	Forkhead box protein L
		adi_v1.00195	−2.99	1.12. E − 33	Forkhead box P3
		adi_v1.16190	−4.68	2.93. E − 18	Forkhead box protein A2
		adi_v1.10031	−1.23	2.04. E − 05	Forkhead box protein O3
		adi_v1.05516	−4.99	8.77. E − 17	Forkhead box protein P
		adi_v1.16084	−7.83	4.81. E − 26	SOX group B
		adi_v1.23073	−8.8	8.27. E − 54	SOX group B
		adi_v1.13187	−6.19	2.18. E − 75	SOX group E/F
		adi_v1.23415	−7.46	6.02. E − 46	SOX group E/F
		adi_v1.05035	−4.57	3.38. E − 27	Transcription factor Sp, invertebrate
		adi_v1.07373	−2.62	2.02. E − 39	Transcription factor Sp2
		adi_v1.05034	−1.13	3.32. E − 07	Transcription factor Sp5
		adi_v1.02807	−3.61	9.70. E − 03	Transcription factor Spi-B
		adi_v1.04124	−0.61	1.16. E − 02	Brachyury protein
		adi_v1.14801	−2.33	9.55. E − 16	Dishevelled associated activator of morphogenesis
		adi_v1.07373	−2.62	2.02. E − 39	Transcription factor Sp2
		adi_v1.05035	−4.57	3.38. E − 27	Transcription factor Sp, invertebrate
		adi_v1.06348	−6.71	5.51. E − 12	Homeobox protein HoxA/B2
		adi_v1.11050	−6.87	1.10. E − 30	Matrix metalloproteinase-23 (CA-MMP) [EC:3.4.24.-]
		XLOC_003365	−7	7.23. E − 16	Dishevelled associated activator of morphogenesis
		adi_v1.16084	−7.83	4.81. E − 26	Transcription factor SOX1/3/14/21 (SOX group B)
		adi_v1.14274	−7.97	4.37. E − 67	Matrix metalloproteinase-23 (CA-MMP) [EC:3.4.24.-]
		adi_v1.23073	−8.8	8.27. E − 54	Transcription factor SOX1/3/14/21 (SOX group B)
		adi_v1.06125	−9.77	4.03. E − 59	Homeobox protein HoxA/B/C6
S vs. P	S	adi_v1.19555	6.54	5.96. E − 06	Homeobox protein aristaless-related
		adi_v1.06348	6.71	1.31. E − 07	Homeobox protein HoxA/B2
		adi_v1.05443	6.06	6.55. E − 45	Homeobox protein GSH
		adi_v1.16084	7.83	9.78. E − 16	Transcription factor SOX1/3/14/21 (SOX group B)
		adi_v1.23415	8.91	3.87. E − 29	Transcription factor SOX7/8/9/10/18 (SOX group E/F)
		adi_v1.05035	6.03	2.00. E − 23	Transcription factor Sp, invertebrate
		adi_v1.07373	3.99	5.80. E − 54	Transcription factor Sp2
		adi_v1.23437	3.68	1.07. E − 09	Alkaline phosphatase [EC:3.1.3.1]
	P	adi_v1.01571	−1.77	6.69. E − 17	Transcription factor SOX1/3/14/21 (SOX group B)
		adi_v1.00157	−1.75	2.16. E − 17	Transcription factor Sp, invertebrate
		adi_v1.05150	−3.46	7.32. E − 38	Transcription factor Sp, invertebrate
		XLOC_010624	−6.25	1.17. E − 06	Homeobox protein MOX
		adi_v1.20693	−4.12	2.09. E − 27	Homeobox protein aristaless-like 4
		adi_v1.00574	−3.94	1.44. E − 05	PAX 3/7/D
		adi_v1.06991	−5.37	1.02. E − 25	PAX 3/7
		adi_v1.22792	−1.4	1.30. E − 06	Myc proto-oncogene protein
		adi_v1.22785	−1.67	1.85. E − 09	Myc proto-oncogene protein
		adi_v1.20978	−1.79	1.47. E − 03	Transcription factor AP-2, invertebrate
		adi_v1.15539	−2.77	1.38. E − 35	Transcription factor AP-1
		adi_v1.04124	−1.18	2.42. E − 07	Brachyury protein
		adi_v1.13155	−2.8	5.28. E − 10	WNT inhibitory factor 1
		adi_v1.21181	−0.58	1.92. E − 02	WNT inhibitory factor 1
		adi_v1.16974	−1.05	1.35. E − 05	Bone morphogenetic protein receptor type-2 [EC:2.7.11.30]
		adi_v1.15796	−5.77	1.63. E − 05	Bone morphogenetic protein 2/4_isoform
		adi_v1.12454	−1.11	1.96. E − 02	Frizzled 10-like
		adi_v1.12455	−4.25	1.73. E − 104	Frizzled 8-like
P vs. A	P	adi_v1.12103	1.92	1.97. E − 07	Chromobox protein 2
		adi_v1.24534	1.79	1.04. E − 16	Chromodomain-helicase-DNA-binding protein 4 [EC:3.6.4.12]
		adi_v1.18668	6.35	1.47. E − 07	Chromobox protein 6
		adi_v1.06504	2.15	4.74. E − 09	Calmodulin-like5
		adi_v1.01102	1.4	2.93. E − 08	Calmodulin CaM
		adi_v1.05564	1.13	1.66. E − 02	Calmodulin-like4
		XLOC_013827	0.93	1.60. E − 05	Calmodulin-7
		adi_v1.03833	0.73	9.25. E − 03	Calmodulin-8
	A	adi_v1.20187	−2.32	8.44. E − 06	Histone-lysine N-methyltransferase MLL1 [EC:2.1.1.43]
		adi_v1.12357	−2.4	2.49. E − 03	Histone deacetylase 6/10 [EC:3.5.1.98]
		adi_v1.01196	−1.07	2.30. E − 04	Polycomb group RING finger protein 3
		adi_v1.19438	−1.46	8.88. E − 03	Polycomb group RING finger protein 5
		adi_v1.09632	−1.93	2.77. E − 11	Polycomb protein SCMH1
		adi_v1.10053	−1.62	1.01. E − 03	Polycomb protein SCMH1
		adi_v1.14017	−1.22	5.58. E − 06	Polycomb protein SCMH1
		adi_v1.00878	−1.45	4.39. E − 06	Polycomb protein SUZ12
		adi_v1.16355	−0.97	3.58. E − 03	Carbonic anhydrase [EC:4.2.1.1]
		adi_v1.11313	−2.61	9.17. E − 10	Carbonic anhydrase [EC:4.2.1.1]

UP regulation in each stage shown. Transcripts shown in this table were selected based on: 1) significant P values (0.05 cutoff value) and 2) evidence in the literature of their involvement in the regulation of developmental processes. A complete list of DEGs can be found in supplementary file S3, Supplementary Material online.

Coexpression of Distinct Gene Networks Regulates Developmental Progression in Acropora

To identify associations between gene expression patterns during development, we performed a weighted correlation network analysis (WGCNA). Gene modules of coexpressed genes usually reflect common functional and regulatory relationships between differentially expressed transcripts. In this analysis, our 18,264 unique sequences were assigned to 37 different gene modules that ranged from 30 to 3,400 transcripts. Modules could be partitioned into three major coexpression clusters, with eigengenes showing similar expression patterns (C1–C3). In most cases module eigengene expression was significantly correlated with specific stages (fig. 4). A summary of transcripts in each coexpression cluster can be found in S5. GO terms mentioned for coexpression clusters are found in table 3.

F

Coexpression gene networks. Transcripts (18,264) were assigned to 37 different gene modules that ranged from 30 to 3,400 transcripts and grouped in three main coexpression clusters. In some cases modules were differentially expressed between two stages indicating negative/positive regulation of specific processes at different developmental time points. Eigengenes were calculated for each module and although we were able to identify discrete gene expression patterns, in most cases significant module-trait correlation were observed in a stage specific fashion. * P value ≤0.05, ** P value ≤0.01, *** P value ≤0.001.

Table 3

GO Enrichment Summary in WGCNA

	Stagea	GO ID	Node size	Sample match	P adj	Term
WGCN clusters’ GO enrichment	PC	GO:0034660	1571	341	0.0030	ncRNA metabolic process
		GO:0007049	6345	1270	0.0082	Cell cycle
		GO:2000036	26	10	0.0165	Regulation of stem cell maintenance
		GO:0048866	8	4	0.0477	Stem cell fate specification
		GO:0006349	77	28	0.0002	Regulation of genetic imprinting
		GO:0043697	12	6	0.0152	Cell dedifferentiation
		GO:0040034	125	35	0.0093	Regulation of development, heterochronic
		GO:0001510	321	93	0.0000	RNA methylation
		GO:0043045	49	18	0.0028	DNA methylation involved in embryo development
		GO:0097198	14	6	0.0351	Histone H3-K36 trimethylation
		GO:0061085	16	7	0.0205	Regulation of histone H3-K27 methylation
		GO:0034401	84	24	0.0219	Regulation of transcription by chromatin organization
		GO:0070870	4	3	0.0236	Heterochromatin maintenance chromatin silencing
		GO:0006338	476	113	0.0050	Chromatin remodeling
		GO:0036093	23	9	0.0201	Germ cell proliferation
		GO:0016332	41	15	0.0064	Maintenance of polarity of embryonic epithelium
		GO:0002070	38	12	0.0447	Epithelial cell maturation
		GO:0044334	6	4	0.0141	Wnt signaling pathway regulation of epithelial transition
		GO:0008543	439	101	0.0196	Fibroblast growth factor receptor signaling pathway
	G	GO:0001703	95	16	0.0051	Gastrulation with mouth forming first
		GO:0046529	77	14	0.0043	Imaginal disc fusion, thorax closure
		GO:0003384	18	6	0.0026	Apical constriction involved in gastrulation
		GO:0048615	10	3	0.0446	Embryonic anterior midgut (ectodermal) morphogenesis
		GO:0007370	44	11	0.0008	Ventral furrow formation
		GO:0045919	29	8	0.0020	Positive regulation of cytolysis
		GO:1902337	11	4	0.0099	Apoptotic process involved in morphogenesis
		GO:0045167	85	15	0.0042	Protein localization involved in cell fate determination
		GO:0097193	491	57	0.0065	Intrinsic apoptotic signaling pathway
		GO:0007191	48	13	0.0001	Dopamine receptor signaling pathway
		GO:0030521	128	20	0.0046	Androgen receptor signaling pathway
		GO:0042921	16	5	0.0080	Glucocorticoid receptor signaling pathway
		GO:0048384	40	7	0.0455	Retinoic acid receptor signaling pathway
		GO:0007210	82	12	0.0384	Serotonin receptor signaling pathway
		GO:0006836	627	65	0.0399	Neurotransmitter transport
		GO:0043409	421	45	0.0490	Negative regulation of MAPK cascade
		GO:0016360	36	7	0.0271	Sensory organ precursor cell fate determination
		GO:0030218	336	38	0.0339	Erythrocyte differentiation
		GO:0070849	81	12	0.0358	Response to epidermal growth factor
		GO:0006122	17	4	0.0477	Electron transport, ubiquinol to cytochrome c
		GO:0015980	1031	107	0.0105	Energy derivation by oxidation of organic compounds
		GO:0051901	32	9	0.0009	Positive regulation of mitochondrial depolarization
		GO:0016042	765	87	0.0018	Lipid catabolic process
		GO:0005996	688	87	4.21 E − 05	Monosaccharide metabolic process
	S	GO:0001510	298	76	0.0003	RNA methylation
		GO:0034660	914	184	0.0157	ncRNA metabolic process
		GO:0007530	188	43	0.0364	Sex determination
		GO:0009994	523	107	0.0473	Oocyte differentiation
		GO:0048232	1123	232	0.0035	Male gamete generation
		GO:0045815	78	24	0.0031	Positive regulation of gene expression, epigenetic
		GO:0016441	566	115	0.0487	Post-transcriptional gene silencing
		GO:0007601	494	108	0.0079	Visual perception
		GO:0009785	21	8	0.0211	Blue light signaling pathway
		GO:0009639	61	17	0.0297	Response to red or far red light
		GO:0048665	88	25	0.0081	Neuron fate specification
		GO:0021587	126	32	0.0172	Cerebellum morphogenesis
		GO:0048854	181	43	0.0210	Brain morphogenesis
		GO:0021692	34	11	0.0269	Cerebellar Purkinje cell layer morphogenesis
		GO:0021527	24	8	0.0469	Spinal cord association neuron differentiation
		GO:0031987	128	31	0.0351	Locomotion involved in locomotory behavior
		GO:0021692	124	31	0.0233	Mechanosensory behavior
		GO:0009950	191	47	0.0087	Dorsal/ventral axis specification
		GO:0048263	40	14	0.0063	Determination of dorsal identity
		GO:0060811	97	26	0.0152	imRNA localization anterior/posterior axis specification
GO:2000738		14	6	0.0245	Positive regulation of stem cell differentiation
GO:0045165		1858	355	0.0387	Cell fate commitment
GO:0003263		8	5	0.0059	Cardioblast proliferation
GO:0061325		37	13	0.0080	Cell proliferation involved tract morphogenesis
GO:0016055		1123	224	0.0190	Wnt signaling pathway
GO:0000188		54	16	0.0205	Inactivation of MAPK activity
GO:0046426		56	16	0.0286	Negative regulation of JAK-STAT cascade
GO:0007173		643	131	0.0350	Epidermal growth factor receptor signaling pathway
GO:0032502		12791	2291	0.0327	Developmental process
GO:0048611		10	5	0.0193	Embryonic ectodermal digestive tract development
GO:0002119		798	167	0.0073	Nematode larval development
GO:0007394		15	6	0.0348	Dorsal closure, elongation of leading edge cells
	GO:0019954	831	83	0.0000	Asexual reproduction
	GO:0007276	3318	228	0.0441	Gamete generation
P	GO:0019953	778	67	0.0023	Sexual reproduction
	GO:0036166	80	16	0.0000	Phenotypic switching
	GO:0010172	179	18	0.0294	Embryonic body morphogenesis
	GO:0009949	36	6	0.0221	Polarity specification of anterior/posterior axis
	GO:0010085	13	3	0.0426	Polarity specification of proximal/distal axis
	GO:0008258	164	20	0.0029	Head involution
	GO:0048580	238	23	0.0205	Regulation of postembryonic development
	GO:0042488	40	10	0.0001	Regulation of odontogenesis of dentin-containing tooth
	GO:0030860	116	15	0.0054	Regulation of polarized epithelial cell differentiation
	GO:0042481	77	10	0.0202	Regulation of odontogenesis
	GO:0045880	104	12	0.0273	Positive regulation of smoothened signaling pathway
	GO:0035481	4	2	0.0212	Positive regulation of Notch signaling pathway
	GO:0031930	20	5	0.0064	Mitochondria-nucleus signaling pathway
	GO:0008543	374	40	0.0005	Fibroblast growth factor receptor signaling pathway
	GO:0007173	385	40	0.0008	Epidermal growth factor receptor signaling pathway
	GO:0060853	6	2	0.0487	Notch signaling pathway involved in cell commitment
	GO:0038032	152	19	0.0028	t G-protein coupled receptor signaling pathway
	GO:0008277	498	45	0.0073	G-protein coupled receptor protein signaling pathway
	GO:1901297	3	2	0.0110	Wnt signaling pathway involved in cell fate commitment
	GO:0016441	566	52	0.0030	Post-transcriptional gene silencing
	GO:0030466	155	17	0.0158	Chromatin silencing at silent mating-type cassette
	GO:0022008	5808	405	0.0022	Neurogenesis
	GO:0045165	1858	135	0.0274	Cell fate commitment
	GO:0007626	1158	87	0.0350	Locomotory behavior
A	GO:0009059	10086	632	0.0414	Macromolecule biosynthetic process
		GO:0006006	881	76	0.0008	Glucose metabolic process
		GO:0009060	256	25	0.0112	Aerobic respiration
		GO:2001171	6	3	0.0037	Positive regulation of ATP biosynthetic process
		GO:0051280	9	3	0.0137	Negative release of sequestered calcium into cytosol
		GO:0051562	9	3	0.0137	Negative regulation of mitochondrial calcium ion
		GO:0051926	28	5	0.0237	Negative regulation of calcium ion transport
		GO:0046534	32	5	0.0400	Positive regulation of photoreceptor cell differentiation
		GO:0042304	59	8	0.0237	Regulation of fatty acid biosynthetic process
		GO:0016042	765	60	0.0196	Lipid catabolic process
		GO:0060290	10	3	0.0187	Transdifferentiation
		GO:0043046	9	3	0.0137	DNA methylation involved in gamete generation
		GO:0035093	8	3	0.0096	Spermatogenesis, exchange of chromosomal proteins
		GO:0002478	603	47	0.0389	Presentation of exogenous peptide antigen
		GO:0051023	55	8	0.0160	Regulation of immunoglobulin secretion
		GO:0044029	2	2	0.0036	Hypomethylation of CpG island
		GO:0071514	88	11	0.0156	Genetic imprinting
		GO:0000183	48	7	0.0230	Chromatin silencing at rDNA
		GO:0044337	8	3	0.0096	Wnt signaling pathway in regulation of apoptotic
		GO:0032483	90	13	0.0027	Regulation of Rab protein signal transduction
		GO:0038096	282	27	0.0110	Fc-gamma receptor signaling pathway
		GO:0033209	123	13	0.0328	Tumor necrosis factor-mediated signaling pathway
		GO:0042753	42	6	0.0377	Positive regulation of circadian rhythm
		GO:0070257	8	3	0.0096	Positive regulation of mucus secretion
		GO:0003331	18	5	0.0034	Regulation of extracellular matrix constituent secretion
		GO:0070278	23	4	0.0451	Extracellular matrix constituent secretion
		GO:0009101	962	80	0.0018	Glycoprotein biosynthetic process
		GO:0018279	185	20	0.0077	Protein N-linked glycosylation via asparagine
		GO:0018242	16	4	0.0130	Protein O-linked glycosylation via serine
		GO:0006891	134	14	0.0300	Intra-Golgi vesicle-mediated transport
		GO:0032956	1041	76	0.0416	Regulation of actin cytoskeleton organization

Enriched in Cluster shown in column. Node size= Total number GO terms in node. Sample match= number of transcripts with GO terms associated to specific nodes. GO categories shown in this table were selected based on: 1) significant P values (0.05 cutoff value) and 2) evidence in the literature of their involvement in the regulation of developmental processes. A complete list of enriched GO categories can be found in supplementary file S4, Supplementary Material online.

Cluster 1: Downregulated in Blastula and Adult but Upregulated in Planula

C1 consisted of 3,050 transcripts and showed two main trends in coexpression. The first one (203 genes) contained modules upregulated in P and mostly downregulated in all other stages, including two modules that showed a significant module-trait correlation in P (76 and 97 genes, respectively). The second one (2,847 genes), contained modules mainly downregulated in PC and A with a tendency to be upregulated in all other stages. This coexpression pattern includes two modules (469 and 575 genes, respectively) that showed a significant module-trait correlation in PC as well as two modules (313 and 374 genes, respectively) that showed a significant module-trait correlation in A (fig. 4).

Cluster 2: Downregulated in Gastrula/Planula and Upregulated in Sphere/Adult

C2 (8,474 genes) grouped molecules that were downregulated in G and P but upregulated in S and A. This cluster showed three main trends in coexpression. The first one (5,129 genes) consisted of molecules mostly downregulated in PC and G, but upregulated in S. This pattern included one module that that showed a significant module-trait correlation in PC, (1,458 genes), two modules that did so in G (66 genes each), and two that showed a significant correlation in S, (3,400 and 68 genes, respectively). Likewise, the second coexpression pattern (1,325 genes) grouped molecules that were mostly upregulated in PC, G and S but downregulated in P and A. This pattern included two modules that showed significant correlation in P, (371 and 199 genes, respectively).

The third coexpression pattern consisted of 2,020 transcripts that were upregulated in A, but with a strong tendency to be downregulated in all other stages. However, this pattern contained a module (676 genes) that was upregulated in PC and showed significant module-trait correlation in both PC and P stages. In a similar way, this coexpression pattern included two modules that showed significant module-trait correlation in G (405 and 263, respectively) from which, one (405 genes) also showed a significant module-trait correlation in S and the other (263 genes) did so in A. This pattern also grouped two modules that showed only a significant module-trait correlation in A (194 and 482 genes, respectively) (fig. 4).

Cluster 3: Down Regulated in Sphere but Upregulated in Gastrula

C3 (2,588 genes) grouped molecules that were consistently downregulated in S with a tendency to be upregulated in G, and either up or down regulated in PC, P and A. Similar to C1, this cluster showed two main coexpression trends. The first one (4,152 genes), consisted of molecules upregulated in PC and downregulated in S, including three modules that showed a significant module-trait correlation in PC (599, 1,761, and 306 genes, respectively) as well as other three that showed a significant module-trait correlation in S (916, 336, and 234, genes, respectively).

The second coexpression pattern (2,588) grouped transcripts that were downregulated in PC and S, upregulated in G and either up- or downregulated in P and A. This trend included one module (247 genes) that showed a significant module-trait correlation in P, one that did so in PC (283 genes) as well as four that that showed significant module-trait correlations in G (101, 345, 1,451, and 161, respectively). From these, one (161) also showed a significant module-trait correlation in S (fig. 4).

GO Enrichment Analyses of Modules Significantly Correlated with Developmental Stages Revealed Stage-Specific Overrepresentation of Distinct BPs

To identify changes in the BPs underlying developmental progression, we performed GO enrichment analysis of modules, which had expression patterns that were significantly correlated with any of the developmental stages. This analysis identified stage-specific enriched BPs that reflected cellular complexity and molecular context. In some cases modules were differentially expressed between two stages indicating negative/positive regulation of specific processes at different developmental time points. For example, PC showed eight significant modules from which only one (676 genes in C2) was shared between PC and other stage (P). This module was upregulated in PC and downregulated in P (fig. 4).

In a similar way, G showed eight significant modules. From these, three modules also displayed significant correlation with other stages (G, S, and A). While two modules (C2) were downregulated in G and either upregulated in S (405 genes) or A (263 genes); the remaining module (C3) was upregulated in G and downregulated in S (161 genes). Likewise, S showed seven significant modules from which, only two displayed significant correlation with other stage (G) (described above). P showed six significant modules from which, only one displayed significant correlation with other stage (PC) (described above). Finally, A showed five significant modules from which only one was shared with other stage (G) (also described above) (fig. 4).

Among BPs overrepresented in PC’s significant modules, we found ncRNA metabolic process (GO:0034660), cell cycle (GO:0007049), stem cell maintenance (GO:2000036), stem cell fate specification (GO:0006349), regulation of gene expression by genetic imprinting (GO:0006349), and establishment of embryonic epithelium (GO:0016332, GO:0002070). Overrepresented categories in G included: gastrulation and morphogenesis (GO:0001703, GO:0046529), apoptosis (GO:0045919), and ncRNA metabolic processes (GO:0034660). Interestingly, we also observed enrichment of signaling pathways that seem to be of vital importance during gastrulation, such as dopamine and serotonin receptor (GO:0007191, GO:0007210), retinoic acid receptor (GO:0048384), androgen receptor (GO:0030521), and glucocorticoid receptor (GO:0042921) signaling pathways (table 3).

In a similar way, overrepresented BPs in S’s significant modules included: nervous system development (GO:0021692), positive regulation of stem cell differentiation (GO:2000738), larval development (GO:0032502), axial specification processes (GO:0009950), sex determination (GO:0007530), posttranscriptional gene silencing (GO:0016441), light perception (GO:0007601), locomotion (GO:0031987), mechanosensory behavior (GO:0007638) as well as RNA methylation (GO:0001510) and ncRNA metabolic process (GO:0034660). Enriched categories in P included: larval development (GO:0002119), axial specification (GO:0009949), epithelial cell differentiation (GO:0030860), neurogenesis (GO:0022008), cell fate commitment (GO:0045165), locomotory behavior (GO:0007626), and chromatin silencing (GO:0030466, GO:0016441) (table 3).

Finally, overrepresented BPs in A’s significant modules included: glucose (GO:0006006) and lipid (GO:0042304), metabolic processes, aerobic respiration (GO:0009060), photoreceptor differentiation (GO:0046534), gamete generation (GO:0043046, GO:0035093), negative regulation of development (GO:0045992), transdifferentiation (GO:0060290), positive regulation of circadian rhythm (GO:0042753), chromatin silencing at rDNA (GO:0000183), genetic imprinting (GO:0071514), hypomethylation (GO:0044029), immune response (GO:0002478, GO:0051023), and positive regulation of mucus secretion (GO:0070257). Interestingly, we found overrepresentation of BPs that might reflect skeleton deposition such as: calcium ion homeostasis (GO:0051280), extracellular matrix secretion (GO:0003331), intra-Golgi vesicle-mediated transport (GO:0006891), protein glycosylation (GO:0018279), and regulation of actin cytoskeleton (GO:0032956) (table 3).

Discussion

Overall, our results showed clear differences in gene expression between early and late developmental transitions that likely reflect changes in the regulatory gene networks underlying the shift between embryonic to larval/adult life stages. While S was identified as the most distinct transcriptome P and A clustered together, which was consistent with previous coral developmental gene expression studies (Reyes-Bermudez et al. 2009). Moreover, the finding that in all transcriptomes approximately 12–16% of all transcripts were identified as putative lncRNAs—slightly more abundant around gastrulation—, suggests roles for these molecules during coral development. This idea is supported by the fact that ncRNA metabolic processes (GO:0034660) were enriched in gene modules that showed significant correlations in PC, G, and S stages. lncRNA molecules represent a poorly understood level of genome regulation able to control chromatin architecture, epigenetic imprinting, and gene expression (Mattick and Makunin 2006; Mercer et al. 2009), strongly implying conserved analogous roles in embryonic coral cell populations. Comparative analyses with lncRNA from other organisms are necessary to understand the role of these molecules during coral development and to further characterize the type of lcnRNA molecules present in our data set.

Furthermore, while the concentration of annotated transcripts in the high ranks of the distribution probably reflects conserved, fundamental cellular processes occurring globally in developing embryos, the low abundance of noncoding transcripts suggests functions for these molecules in specific cell populations (fig. 2B). It has been shown that lncRNA molecules are cell type-specific, with distinct cellular localizations and functions (Mattick and Makunin 2006). In situ hybridizations will be necessary to test this idea. On the other hand, the fact that 13–14% of all transcripts had no orthologs in other systems indicates taxa-specific modifications of fundamental developmental processes. Newly evolved and/or highly divergent taxon-restricted genes with roles in axial patterning and endoderm formation have been recently identified in the hydrozoan cnidarian Clytia sp. (Lapebie et al. 2014), supporting the idea of coral-specific modifications of ancestral developmental signaling cascades.

Blastula to Gastrula: Negative Regulation of Transcription and Early Imprinting of Cell Lineages

Whereas transcripts encoding brick-a-brack, tramtrack, broad-complex/poxvirus zinc finger (BTB/POZ) proteins were abundantly upregulated in PC, molecules encoding basic helix-loop-helix (bHLH) TFs were exclusively upregulated in G. BTB/POZ proteins are known repressors of transcription (Collins et al. 2001) and some bHLH proteins are transcription repressors that influence cell proliferation and differentiation during embryogenesis (Atchley and Fitch 1997), suggesting distinct levels of transcription repression between blastula and gastrula stages. Upregulation in PC of a transcript encoding the transcriptional corepressor, groucho (Jennings and Ish-Horowicz 2008), and inhibitors of wnt signaling indicate to some extent negative regulation of cell fate at blastula prior gastrulation. However, enrichment of GO categories: stem cell maintenance, and specification (GO:2000036, GO:0006349), as well as regulation of gene expression by genetic imprinting and chromatin organization (GO:0006349, GO:0034401) in PC’s significant modules, suggests early imprinting of cell lineages in blastula. Upregulation in PC of an AP-2 ortholg, which is first expressed in the primitive ectoderm of eumetazoans (Eckert et al. 2005) supports this idea (tables 2 and 3).

On the other hand, upregulation in both PC and G of distinct hox TFs implies complex regulatory networks underlying axial specification and morphogenetic gradients early in development (Botas 1993; Deschamps and van Nes 2005; DuBuc et al. 2012). Upregulation in G of TFs known to regulate stem cell differentiation such as hes1 (Kobayashi et al. 2009), myc (Ambrosone et al. 2012; Tansey 2014; Zinin et al. 2014), pou-domain (Millane et al. 2011) as well as a diverse array of sox (Jager et al. 2011), krupple (McConnell and Yang 2010), and fox (Zaret and Carroll 2011) are likely to reflect germ layer lineage differentiation during gastrulation. Moreover, enrichment in G’s significantly correlated modules of components of dopamine and serotonin receptor (GO:0007191, GO:0007210), retinoic acid receptor (GO:0048384), androgen receptor (GO:0030521, GO:0033574), and glucocorticoid receptor (GO:0042921) reveals key transcriptional regulatory roles for these signaling pathways in the gastrula stage (tables 2 and 3).

Gastrula to Sphere: Cell Diversification and the Initiation of Larval Life

Upregulation in G of several repressors of wnt signaling, such as dickkopf (Niehrs 2006) together with over representation in S of DEGs with diverse roles in cell diversification (GO:0030855, GO:0021575, GO:0048732, GO:0048644, GO:2000738, GO:0048665, GO:0045165), larval development (GO:0002119), and locomotory behavior (GO:0007626), suggests terminal differentiation of cellular phenotypes in S. This is consistent with the idea that developmental progression is characterized by a decline in undifferentiated cell populations, followed by increased committed cell types (Chambers and Studer 2011). Interestingly, both G and S upregulated distinct and diverse “pioneer” TFs from the fox family (Hannenhalli and Kaestner 2009; Zaret and Carroll 2011). Fox TFs are transcription regulators able to bind condensed chromatin during cell differentiation and thus primed loci for gene expression (Hannenhalli and Kaestner 2009; Zaret and Carroll 2011); indicating the initiation of cell-specific transcriptional circuits during larval body morphogenesis (tables 1–3).

Likewise, G and S upregulated distinct brachyury isoforms suggesting the usage of complex transcriptional networks during blastopore and endodermal specification (Kispert and Hermann 1993). Brachyury expression was reported in the blastopore and developing mesenteries of developing embryos from the related anthozoan, N. vectensis (Scholz and Technau 2003) indicating similar expression domains and conserved function in A. digitifera (table 2). Finally, cellular diversification and emergence of a larval body plan in S was also reflected in the upregulation of diverse components of FGF signaling, several matrix metalloproteinases as well as two disheveled coding transcripts. FGF signaling controls pluripotency and lineage segregation during development (Lanner and Rossant 2010), metalloproteinases have been associated with tissue remodeling in Drosophila (Page-McCaw et al. 2003) and disheveled proteins are regulators of the actin cytoskeleton during morphogenetic processes (Li et al. 2011) (table 2).

Sphere to Planula: Motile Life and Tissue Imprinting for Settlement and Metamorphosis

Enrichment of GO category larval development (GO:0002119) in significant modules for S and P stages suggests that establishment of mature larval morphology is achieved by the use of intricate and complementary transcriptional circuits. This is consistent with upregulation of distinct sox, hox, and sp TFs in both developmental time points. Our results are consistent with the idea that Acropora sp. planulae is transcriptionally primed for habitat switch (Grasso et al. 2011) as FCs and number of DEGs in the transition spaning settlement and metamorphosis (P to A) were the lowest in the data set (fig. 3). Exclusive upregulation of pax, myc, and AP-1/2/4 TFs in P, postulate these molecules as key regulators of cellular imprinting prior to the switch from a pelagic to a sessile existence (tables 2 and 3). Although Myc proteins regulate chromatin structure, proliferation, and terminal cell differentiation (Ambrosone et al. 2012; Tansey 2014; Zinin et al. 2014), Pax and AP molecules are regulators of cell-fate specification and tissue regionalization (Chi and Epstein 2002).

Moreover, upregulation in P of the oral/aboral axial determinant, brachyury, which also has roles in specification of endodermal structures (Kispert and Hermann 1993), suggests that axial polarity and imprinting of cell populations is indeed actively occurring at the onset of metamorphosis. Interestingly, upregulation in P of two frizzled receptors and two inhibitors of wnt signaling suggests that similar to the hydrozoan Clytia hemisphaerica, axial polarity prior metamorphosis in A. digitifera is determined by asymmetric activation of the wnt pathway (Momose and Houliston 2007). Furthermore, upregulation in P of a BMP-2/4 ligand and its receptor, which determines tissue boundaries (Hayward et al. 2002),—and are considered calcifying epithelium markers (Zoccola et al. 2009)—, indicates roles for BMP signaling in tissue reorganization during metamorphosis. Likewise, overrepresentation of molecules with roles in nucleosome assembly (GO:0006334), RNA methylation (GO:0001510), chromatin silencing (GO:0030466, GO:0016441), and mitochondrial ATP synthesis (GO:0042775) in the subset of DEGs and modules differentially expressed by P, suggest changes in chromatin structure and energy metabolism at the onset of settlement an metamorphosis (tables 1 and 3).

Planula to Adult: Habitat Switch and Responses to the Environment

DEGs with roles in cell-cell junctions (GO:0045216), immune responses, (GO:0002433), oocyte generation (GO:0007292), and asexual reproduction (GO:0019954) were over represented in the subset of DEGs upregulated by A. Likewise, BP categories macromolecule biosynthetic process (GO:0009059), positive regulation of circadian rhythm (GO:0042753), and positive regulation of mucus secretion (GO:0070257) were enriched in A’s significant modules. These results together with the fact that two carbonic anhydrases orthologs,—which are enzymes known to regulate pH and skeleton deposition in corals (Moya et al. 2008)—were also upregulated in A, indicate the initiation of transcriptional circuits underlying adult specific processes such as skeleton deposition. Overrepresentation in A of GO categories: extracellular matrix secretion (GO:0003331, GO:0070278), intra-Golgi vesicle-mediated transport (GO:0006891), protein glycosylation (GO:0009101, GO:0018279, GO:0018242), and regulation of actin cytoskeleton (GO:0032956) support this idea (tables 1–3).

The finding that five DEGs encoding CaM-like molecules were exclusively upregulated in P in this comparison, indicates diversification of calcium signaling pathways at the onset of settlement and metamorphosis. CaM-like molecules have been previously reported as key regulators of settlement and metamorphosis in corals (Reyes-Bermudez et al. 2009, 2012) (table 2). Furthermore, overrepresentation of GO categories respiration (GO:0022904, GO:0006120), RNA methylation (GO:0001510), and chromatin silencing (GO:0030466, GO:0016441) in the subset of transcripts differentially expressed by P, indicates high metabolic rates prior to the lifestyle switch, and suggests that epigenetic regulation prior metamorphosis in Acropora might be happening. To test this idea, more research is necessary (tables 1–3).

Conclusions

This study revealed clear differences in gene expression between early and late developmental transitions, with higher numbers of DEG and FCs in the progression involving gastrulation. These differences might reflect transcriptional changes underlying the transition between embryonic to larval and adult life stages, revealing a highly active and plastic A. digitifera genome. During early transitions, transcriptional networks seemed to regulate cellular fate and morphogenesis of the larval body. In late transitions, these networks are likely to play important roles preparing planulae for the lifestyle switch and in colonial polyps to regulate adult processes. Although, development and tissue plasticity in corals are likely to be regulated to some extent by differential coexpression of well-defined gene networks, the fact that some modules were restricted to specific developmental time points indicates that stage-specific transcription profiles are independent entities with very distinct molecular contexts.

Similar to vertebrates, developmental networks in corals appeared to be linked to changes in chromatin architecture (Meshorer and Misteli 2006; Meshorer et al. 2006) and to control cell differentiation during early development by repressing the expression of canonical developmental signaling pathways (Orkin and Hochedlinger 2011). Furthermore, imprinting of embryonic cell populations by widely conserved TFs is likely to reflect ancestral regulatory pathways underlying cell differentiation in eumetazoa. Despite this, the finding that approximately 13% of all transcripts in our data set are coral-specific suggests taxa-specific modifications of fundamental developmental processes. In situ hybridization studies and functional experiments are necessary to fully characterize the role of both coral-restricted and widely distributed molecules during coral development.

Comparison of gene expression between sequential developmental stages and coexpression gene network analysis provided different resolution of transcriptional dynamics underlying coral development, yet the results obtained by the two approaches were consistent and complementary. This work provides a quantitative perspective on global transcriptional dynamics during A. digitifera development, but lacks spatial resolution. It constitutes a framework for future studies; thus we encourage researchers to use the data set and to examine in detail genes and gene expression patterns that were out of the scope of this study.

Supplementary Material

Supplementary files S1–S5 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).