Differential gene expression during substrate probing in larvae of the Caribbean coral Porites astreoides.

The transition from larva to adult is a critical step in the life history strategy of most marine animals. However, the genetic basis of this life history change remains poorly understood in many taxa, including most coral species. Recent evidence suggests that coral planula larvae undergo significant changes at the physiological and molecular levels throughout the development. To investigate this, we characterized differential gene expression (DGE) during the transition from planula to adult polyp in the abundant Caribbean reef-building coral Porites astreoides, that is from nonprobing to actively substrate-probing larva, a stage required for colony initiation. This period is crucial for the coral, because it demonstrates preparedness to locate appropriate substrata for settlement based on vital environmental cues. Through RNA-Seq, we identified 860 differentially expressed holobiont genes between probing and nonprobing larvae (p ≤ .01), the majority of which were upregulated in probing larvae. Surprisingly, differentially expressed genes of endosymbiotic dinoflagellate origin greatly outnumbered coral genes, compared with a nearly 1:1 ratio of coral-to-dinoflagellate gene representation in the holobiont transcriptome. This unanticipated result suggests that dinoflagellate endosymbionts may play a significant role in the transition from nonprobing to probing behaviour in dinoflagellate-rich larvae. Putative holobiont genes were largely involved in protein and nucleotide binding, metabolism and transport. Genes were also linked to environmental sensing and response and integral signalling pathways. Our results thus provide detailed insight into molecular changes prior to larval settlement and highlight the complex physiological and biochemical changes that occur in early transition stages from pelagic to benthic stages in corals, and perhaps more importantly, in their endosymbionts.

INTRODUCTION

Sessile benthic organisms must find a way to produce offspring that will successfully develop into fully functioning adults. Most marine invertebrates with a benthic adult life have a pelagic life stage, during which they must find the most appropriate substrate to settle and metamorphose. This dispersal step is vital to the species' success, because it decreases competition for resources and inbreeding and increases the chances of species' survival in the event of local extinction events (Mayr, 1970; Pechenik, 1999). Equally important, however, is that larvae select the most suitable habitat to settle and begin colony growth (e.g., Rodríguez, Ojeda, & Inestrosa, 1993), a critical moment in the life cycle of a coral where mortality rates are extremely high (Martinez & Abelson, 2013).

For scleractinian corals, some of the most charismatic ecosystem engineers, two general forms of sexual reproduction are found, brooding—where fertilization and larval development occur internally within the body cavity of maternal polyps—and broadcast spawning—where fertilization occurs externally in the water column after synchronized gamete release. Resulting planulae then disperse through the water column, and the polyp life stage begins after the competent larva settles on the ocean floor and metamorphoses. This metamorphosis from the larval to the polyp stages includes a number of alterations. Larvae have a simple morphology; they are essentially elongated forms with cilia for movement and a mouth structure for passive eating. By contrast, polyps are much more complex; they notably have a gastrovascular cavity, a ring of tentacles, and produce their characteristic rigid calcium carbonate skeleton (Fadlallah, 1983). After the larva settles and becomes the founding polyp of the colony, it secretes a calcium carbonate skeleton and new polyps start budding asexually to form the coral colony, while colony propagation occurs predominantly by sexual reproduction (Baird, Guest, & Willis, 2009; Harrison & Wallace, 1990). Larval settlement and metamorphosis are thus crucial processes that connect the two main life stages of corals: the pelagic planula larva and benthic polyp (e.g., Ritson‐Williams et al., 2009), and it is this step that we investigate here.

Prior to settlement, competent larvae spend time “probing” the benthos (i.e., touching down and crawling) for appropriate substrata. During this time, larvae have been reported to react to a wide range of environmental cues, including light intensity (e.g., Maida, Coll, & Sammarco, 1994; Mundy & Babcock, 1998) and light colour (Strader, Davies, & Matz, 2015), reef sounds mostly generated by fish and crustaceans (Vermeij, Marhaver, Huijbers, Nagelkerken, & Simpson, 2010), chemical signals from crustose coralline algae (CCA) (Heyward & Negri, 1999; Morse, Hooker, Morse, & Jensen, 1988; Ritson‐Williams, Arnold, & Paul, 2016) and microbial biofilms (e.g., Sneed, Sharp, Ritchie, & Paul, 2014; Webster et al., 2004). While the anatomical changes larvae undergo to settle and metamorphose are well understood (e.g., Harrison & Wallace, 1990; Hirose, Yamamoto, & Nonaka, 2008), little is known about the underlying molecular pathways that drive these major changes. For example, in Acropora spp., upregulation in sensory and signal transduction genes was found to accompany the development of larval settlement competency (Strader, Aglyamova, & Matz, 2018) and distinct aboral and oral expression patterns were observed during metamorphosis from settled larvae to polyps (Grasso et al., 2011).

Symbiotic dinoflagellates of the family Symbiodiniaceae, commonly known as zooxanthellae, are essential for survival in many coral species. They live inside host endodermal cells and provide oxygen and other photosynthesis by‐products on which the coral relies. In turn, the coral provides a protected environment and access to compounds required for photosynthesis, such as CO2 and ammonium (Muscatine & Porter, 1977; Yellowlees, Rees, & Leggat, 2008). Zooxanthellae acquisition often correlates with the mode of sexual reproduction. Larvae of brooding species (like Porites astreoides) tend to vertically inherit zooxanthellae from the maternal polyp (Goldberg, 2013). In contrast, larvae of broadcast spawning corals usually acquire zooxanthellae from the surrounding environment later in the larval stage or after metamorphosis (e.g., Adams, Cumbo, & Takabayashi, 2009; Hirose et al., 2008).

The scleractinian genus Porites includes some of the most recognizable and dominant reef builders and has a circumtropical distribution (Veron, 2000). Porites astreoides Lamarck, 1816 has recently increased in abundance to become a locally dominant shallow‐water species on many Caribbean reefs (Green, Edmunds, & Carpenter, 2008) and serves as a coral model system (e.g., Kenkel, Goodbody‐Gringley, et al., 2013; Kenkel, Meyer, & Matz, 2013; Serrano et al., 2016). Due to their relative abundance and predictable larval brooding and release schedule, P. astreoides colonies are good candidates to study larval settlement and metamorphosis. In addition, P. astreoides planulae are well developed upon release, including formed mesenteries, stomodaeum (mouth and pharynx), gastrodermis and nematocysts in the ectoderm and endoderm (Chornesky & Peters, 1987) and endosymbiont‐rich, which allows for the symbiont's role throughout the coral life cycle to be studied. Previous research on this species focused on larval settlement induction via microbial biofilms isolated from crustose coralline algae (CCA) (Sneed et al., 2014), larval settlement response in the presence of macroalgae (Campbell, Sneed, Johnston, & Paul, 2017; Kuffner et al., 2006; Olsen, Sneed, & Paul, 2016; Ritson‐Williams et al., 2016), response to reduced pH (Campbell et al., 2017), larval phenotypic variability implications for settlement success (Putron et al., 2017) and effects of thermal stress and elevated nitrate levels on larval settlement and mortality (Serrano et al., 2018). However, little has been done to understand the genetic mechanisms of settlement.

Studying differential gene expression of probing and nonprobing larvae in controlled environments has the potential to illuminate the mechanisms that drive and guide larval settlement. Differential physiological and molecular responses are exhibited during this transition (e.g., Grasso et al., 2011), so studying planulae without distinguishing between probing and nonprobing behaviours could therefore impact differential gene expression results. In addition, previous studies were conducted on aposymbiotic larvae (e.g., Grasso et al., 2011; Hayward et al., 2011; Meyer, Aglyamova, & Matz, 2011; Meyer, Aglyamova, et al., 2009; Meyer, Davies, et al., 2009; Polato et al., 2010; Reyes‐Bermudez, Lin, Hayward, Miller, & Ball, 2009; Reyes‐Bermudez, Villar‐Briones, Ramirez‐Portilla, Hidaka, & Mikheyev, 2016; Strader et al., 2018). However, no larval gene expression study to date has been conducted on corals that have endosymbionts throughout larval development. Thus, this study builds upon previous research as the first transcriptomic analysis of coral larvae and their endosymbionts to focus on probing versus nonprobing behaviour cued by the presence of a preferred settlement substratum. For this, we explored general gene expression patterns, putative genes and corresponding protein functions, and major biological pathways associated with the transition from nonprobing to probing larvae. We expected to observe significant differences between probing and nonprobing larvae, with mostly significant upregulation in probing samples, since this period marks the onset of a distinct transition and previous work has demonstrated expression shifts during transitional stages in the development (e.g., Reyes‐Bermudez et al., 2016). We also expected to find differentially expressed genes corresponding to cell growth and differentiation, signalling and transport, and metabolism and that these genes would tend to be upregulated in probing larvae, to signify increasing complexity throughout development (e.g., Reyes‐Bermudez et al., 2016; Strader et al., 2018). Additionally, we expected to find that most differentially expressed genes within the coral holobiont could be attributed to the coral host rather than the endosymbiont, as had been previously demonstrated in adult corals (e.g., Kaniewska et al., 2015), and endosymbiont genes might not be significantly affected by substrate probing. Overall, this study examines underlying molecular components of the shift from nonprobing to substrate probing, a discrete and potentially major transitional period during larval development.

MATERIALS AND METHODS

Sample collection and storage

In the Florida Keys, peak larval release for Porites astreoides occurs during the night on days near the new moon in April and May (McGuire, 1998), and larvae are competent to settle upon release and remain competent for a week or more (see Ritson‐Williams et al., 2016). Porites astreoides adult coral heads were collected from Wonderland Reef in the lower Florida Keys (24.56069°, −81.50135°) and maintained in flow‐through seawater aquaria at Mote's Tropical Research Laboratory in Summerland Key, Florida. Following their regular larval release cycle (McGuire, 1998), 19 colonies released between 2 and 2,500 coral larvae on the nights of 9–10 May 2013. Larvae released on the night of 9 May and 10 May were kept separate and maintained in flow‐through seawater in buckets with 180 µm mesh at the bottom, that is without any benthic substrate. Buckets were cleaned daily to avoid build‐up of microbial biofilms that might induce settlement in these larvae. On 11 May 2013, 30 larvae were selected randomly and transferred to individual wells within 6‐well plates containing 10 ml of 0.2 µm filter‐sterilized seawater (FSW). Approximately 1 cm2 of the crustose coralline alga Titanoderma prototypum (CCA) was added to some wells to initiate probing behaviour (Ritson‐Williams et al., 2016), while others contained only FSW with no settlement cue. After 3–4 hr of exposure, 20 larvae were collected that were actively probing the CCA in preparation for settlement and metamorphosis. Twenty larvae were also collected from the wells without CCA, that is they were not exposed to a settlement cue and were not probing. Collected larvae were kept together in their experimental groups, and then, individuals were separated and immediately fixed in RNAlater and stored at −80°C until further processing.

RNA‐Seq library preparation

mRNA extraction

To extract mRNA, either individual larvae or pools of ten larvae (Table S1) were transferred to tubes containing 500 μl of TRIzol (Invitrogen) to break down the tissue and preserve RNA and DNA. After grinding the larvae in TRIzol solution with an RNase‐free plastic pestle, another 500 μl of TRIzol was added. To capture as much RNA as possible from such small tissue samples, 10 μl of glycogen was added. The solution was then incubated at room temperature (RT) for 10 min. To isolate RNA, 100 μl of bromochloropropane (BCP) was added, followed by another 10 min of incubation at RT. Afterwards, the solution was centrifuged at 16,000 rcf for 15 min at 4°C. The upper aqueous layer, which contains the total RNA content, was recovered and transferred to a new tube. A mixture of 250 μl isopropanol and 250 μl high salt plant solution (Clontech Laboratories, Inc.) was added and stored overnight at −20°C for RNA precipitation.

The next day, samples were centrifuged at 16,000 rcf for 15 min at 4°C, and the supernatant was removed. Total RNA precipitation was performed by washing the remaining pellet twice in 1 ml of 75% ethanol and centrifuging it at 7,600 rcf for 5 min at 4°C. The dried pellet was eluted in 100 μl of RNA Storage solution (Invitrogen). The purification of mRNA was performed with the Dynabeads mRNA Purification Kit (Invitrogen) following the manufacturer's instructions, and mRNA was eluted in 13 μl of 10 mM Tris‐HCl. The supernatant was immediately transferred to a 0.5‐ml cryovial and preserved at −80°C. Remaining rRNA contamination was assessed with pico‐RNA assays in an Agilent 2100 Bioanalyzer (Agilent Technologies), and the mRNA extraction protocol was repeated if necessary.

cDNA library construction and next‐generation sequencing

cDNA libraries were constructed in the Apollo 324 automated system using the PrepX mRNA 8 Protocol Kit (IntegenX) with ~1–2 pg/μl mRNA per larvae. cDNA libraries were amplified using the KAPA Library Amplification Kit (Kapa Biosystems) according to the manufacturer's instructions, with 16–20 PCR cycles. During the PCR process, samples were marked with different 6 bp indices for multiplexing. Final library quality and size distribution were measured with a DNA high‐sensitivity assay in an Agilent 2100 Bioanalyzer (Agilent Technologies). Library concentrations were measured using a qPCR with the KAPA Library Quantification Kit (Kapa Biosystems). Samples were paired‐end sequenced (2x 150 bp) on an Illumina HiSeq 2500 platform at the Faculty of Arts & Sciences (FAS) Center for Systems Biology at Harvard University.

For this study, seventeen different RNA‐Seq libraries were sequenced. Five individual larvae were sequenced for probing day 1 (PD1a‐e), control day 1 (CD1a‐e) and control day 2 (CD2a‐e), and two samples of ten pooled larvae were sequenced for probing day 2 (PD2A&B).

Data processing, de novo transcriptome assembly and gene expression analyses

Sequences were processed on the Odyssey cluster at Harvard University Faculty for Arts & Sciences. Raw reads were trimmed using Trim Galore (Babraham Bioinformatics, 2015), removing adapters, indices, low‐quality positions (Q < 30) and reads shorter than 25 bp. Mitochondrial DNA (mtDNA) and rRNA were filtered out by using Bowtie 2.0 (Langmead, Trapnell, Pop, & Salzberg, 2009) against a database of known cnidarian mtDNA and rRNA sequences from GenBank.

A reference transcriptome was assembled using Trinity (Haas et al., 2013) with default parameters and read normalization. The seventeen RNA‐Seq libraries from this study were combined with sixteen additional RNA‐Seq libraries from different pools of P. astreoides larvae and spat to generate the reference transcriptome (Table S2). Holobiont transcriptome assembly completeness was assessed through BUSCO v3, using a set of 303 single‐copy orthologs across Eukaryota and 987 single‐copy Metazoan orthologs (Waterhouse et al., 2018). BUSCO v3 was also used to separately assess the completeness of the coral transcriptome and symbiont transcriptome that were separated from the holobiont transcriptome (see below). Further, holobiont transcriptome and separated coral and symbiont transcripts quality were evaluated by mapping back the raw reads of each library. Holobiont transcriptome isoform and gene abundance were quantified through eXpress (Roberts & Pachter, 2012). Isoform and gene‐level abundance estimates were used to generate normalized expression values (TPM, transcripts per million transcripts) based on transcript length, mapped reads per transcript and total reads in eXpress. Differentially expressed genes between probing and nonprobing samples were identified using cross‐sample normalized expression values (TMM, trimmed mean of M values, from TPM values) in DESeq2 (Love, Huber, & Anders, 2014). The most significantly differentially expressed genes were extracted and arranged based on expression pattern relatedness (Love et al., 2014). A p‐value (p) of ≤ .01 was chosen for differential expression analyses to employ a stringent threshold and represent greater biological significance (e.g., Bay & Palumbi, 2015; Mayfield, Wang, Chen, Lin, & Chen, 2014). Differential expression analyses were conducted for both genes and isoforms but are presented here only for the gene level. Here, we are using the Trinity definition of “isoform” and “gene”, that is “isoforms” are unique assemblies of transcripts while a “gene” is the sum of all isoforms that share transcripts (https://github.com/trinityrnaseq/trinityrnaseq).

Taxonomic identification

Holobiont genes were taxonomically identified as coral host or endosymbiont genes in BlastN with an E‐value cut‐off of ≤ 10–5 as in De Wit et al. (2012). Since there is no sequenced P. astreoides genome, we used the symbiont‐free transcriptome of the closely related coral species Porites lobata (Bhattacharya et al., 2016) as a taxonomic reference for the coral host. Since Porites astreoides is most commonly associated with the genus Symbiodinium (formerly called Clade A) (Hauff, Haslun, Strychar, Ostrom, & Cervino, 2016), we used the genome of Symbiodinium microadriaticum (Aranda et al., 2016) as a taxonomic reference for its symbionts. Genes were labelled as coral if they matched solely or with a more significant E‐value with P. lobata than S. microadriaticum; genes were labelled as symbiont in the same manner if they matched solely or with a more significant E‐value with S. microadriaticum. However, genes were only included in downstream annotation, gene ontology and KEGG analyses if they matched solely to either the P. lobata or S. microadriaticum references. Genes identified by Trinity do not necessarily represent unique genes in the Porites astreoides or endosymbiont dinoflagellate genome. To quantify this inconsistency between actual genes and the putative genes identified here by Trinity, we assessed how many Trinity genes matched genes in the Porites lobata and Symbiodinium microadriaticum references.

Annotation of transcriptome and differentially expressed genes, gene ontology and KEGG analysis

Genes were annotated by BlastX (E‐value cut‐off of ≤ 10–5) homology searches using public databases in the UniProt Knowledgebase (UniProt Consortium, 2011): Swiss‐Prot and TrEMBL. Since TrEMBL is not curated, it was only used if there was no available Swiss‐Prot annotation. Annotations were carried out at the isoform level; the isoform with the most significant E‐value was chosen to represent a particular gene.

Gene Ontology (GO) annotations for genes were obtained from the Gene Ontology Annotation database (Camon et al., 2004; Huntley et al., 2015; http://www.geneontology.org). The GO Annotation database was used to identify putative functions of transcriptome genes and the subset of differentially expressed genes. A background of Swiss‐Prot annotated transcriptome genes was used for functional enrichment annotation analysis of annotated differentially expressed genes using DAVID (Database for Annotation, Visualization and Integrated Discovery; https://david.ncifcrf.gov) Bioinformatics Resources 6.8 (Huang, Sherman, & Lempicki, 2009). Through DAVID's functional annotation platform, the false discovery rate (FDR) correction was calculated by using the Benjamini–Hochberg algorithm, and enrichment was determined with an adjusted FDR p ≤ .01. Additionally, KEGG orthology‐based annotations were identified for coral genes using default parameters from the KEGG Automatic Annotation Server (KAAS) (Moriya, Itoh, Okuda, Yoshizawa, & Kanehisa, 2007).

RESULTS

Transcriptome assembly

Almost 650 million 150 bp paired‐end sequencing reads were generated for this study from 17 different samples (15 individual specimens and 2 sample pools), and on average, ~38 million raw reads were generated per specimen/pool (full range: 16–81 M; Table S2a). After trimming poor‐quality reads (~1.5% on average) and discarding rRNA reads (~12% on average), 575 million raw reads were normalized and then used to generate the holobiont reference transcriptome; 34 million mRNA reads (range: 15–75 M) were available on average per sample for differential gene expression analyses (Table S2a).

The assembled Porites astreoides holobiont reference transcriptome consisted of 918,990 isoforms and 446,849 putative genes (Table S2b), which is comparable with other coral holobiont transcriptomes (e.g., Mansour, Rosenthal, Brown, & Roberson, 2016; Pinzón et al., 2015; Shinzato, Inoue, & Kusakabe, 2014). Holobiont transcriptome completeness assessments using BUSCO v3 identified 244 out of 303 complete (80.5%), 38 fragmented (12.5%) and 21 missing (7.0%) Eukaryote BUSCO orthologs and 746 out of 978 complete (76.3%) as well as 113 fragmented (11.6%) and 119 missing (12.1%) Metazoan BUSCO orthologs (Table S2c). Transcriptome quality checks matched 312 million raw reads to the reference holobiont transcriptome (54%; 6.5–45 million per sample), and 71% of matched reads aligned in the orientation expected from paired‐end sequencing.

Characterization of the holobiont transcriptome

Out of the 446,849 putative genes in the holobiont transcriptome assembly, 159,722 (36%) matched either one or both references: 83,195 putative genes matched the Porites lobata reference transcriptome and 86,007 matched the Symbiodinium microadriaticum reference genome. Among transcriptome genes that returned matches, there was approximately a 6% overlap with 9,480 genes matching both references. 73,715 genes matched only to P. lobata and 76,527 genes matched only to S. microadriaticum (Figure 1a). Of the genes that matched both genomes, 3,201 matched better to P. lobata and 6,279 matched better to S. microadriaticum so these genes were considered coral and dinoflagellate, respectively. In summary, 76,916 genes were considered to be derived from the coral host and 82,806 genes were considered to be from endosymbionts (i.e., 17% and 19% of the 446,849 putative de novo genes; Figure 1b).

Figure 1

Taxonomic identification of genes in the holobiont transcriptome (a and b) and differentially expressed genes (c and d). In a and c, genes that matched both references, Porites lobata and Symbiodinium microadriaticum, are identified as dual matches. In contrast, in b and d, genes that matched both references were identified as either coral or dinoflagellates genes, depending on the better match

The 83,195 putative de novo coral genes aligned to 18,500 P. lobata genes, so we covered 88% of the 21,062 genes currently identified in the P. lobata transcriptome (Bhattacharya et al., 2016). The 86,007 putative de novo dinoflagellate genes aligned to 29,930 S. microadriaticum genes, so we recovered 61% of the 49,109 genes currently identified in the S. microadriaticum genome (Aranda et al., 2016). Separate BUSCO analyses for coral and symbiont transcriptome subsets identified 193 out of 303 complete (63.7%) as well as 61 fragmented (20.1%) and 49 missing (16.2%) Eukaryote BUSCO orthologs in the coral subset compared with 187 complete (61.7%), 40 fragmented (13.2%) and 76 missing BUSCO orthologs (25.1%) in the Porites lobata reference transcriptome (Bhattacharya et al., 2016) (Table S2c). In the symbiont transcriptome subset, we found 182 complete (60%), 43 fragmented (14.2%) and 78 missing BUSCO Eukaryote orthologs (25.8%) compared with 60 complete (19.1%), 19 fragmented (6.3%) and 224 missing (73.8%) BUSCO Eukaryote orthologs in the S. microadriaticum genome (Aranda et al., 2016). Comparisons with the BUSCO eukaryotes, protists and alveolates data sets confirm (Table S2c) that our holobiont, coral and symbiont transcriptomes are qualitatively similar or better than the reference data sets. The results for the BUSCO analyses with the Metazoan subset are similar (Table S2c). The assembled holobiont transcriptome thus putatively contains a large proportion of genes from both the host and the symbiont.

The coral subset of our de novo transcriptome matched to 112.7 million sequence reads in probing and control samples (i.e., 47.7% of all reads)—approximately 46.4 million reads aligned to probing and 66.3 million reads to control samples—on average 6.6 million reads aligned per sample (Table S2a). The symbiont subset matched to 58.5 million reads (24.8% of all reads)—approximately 24.1 million reads aligned to probing and 34.4 million reads to control samples— with an average of 3.4 million reads per sample (Table S2a). This indicates that even though less reads were attributed to symbionts, we achieved sufficient sequencing depth for both coral and symbiont differential gene expression analyses (e.g., Liu, Zhou, & White, 2014).

When evaluating read depth across isoforms, we found that there were 5,940 symbiont isoforms out of 274,863 (2.2%) that had more than zero TPM (transcripts per million transcripts), and there were 28,720 coral isoforms out of 443,025 (6.5%) that were above the zero TPM cut‐off. The difference between the proportion of coral versus symbiont expressed isoforms is highly significant (z test statistic, p‐value ≤ .01). This finding could indicate that the majority of symbiont isoforms are relatively lowly expressed in coral larvae. However, this result might also imply that a significant proportion of symbiont isoforms were not sequenced deeply enough for detection through the method developed in this study, although this seems unlikely. Further, more highly expressed isoforms range from 5 × 10−27 to 123,155 TPM in symbionts (on average, 170 TPM) and from 1 × 10−96 to 6,160 TPM in corals (on average, 0.262 TPM). This demonstrates that although expressed coral isoforms outnumber symbiont isoforms, highly expressed symbiont isoforms were expressed significantly more than coral isoforms. Few gene expression studies report read count statistics for individual isoforms or genes; however, one recent study similarly found thousands of contigs (contiguous gene sequences) that surpassed a minimum read count threshold to compare the transcriptomic response to recurrent bleaching events between two acroporid coral species (Thomas, López, Morikawa, & Palumbi, 2019).

Differential gene expression analysis

Out of the 446,849 putative genes and 918,990 isoforms in the reference transcriptome, 321,170 putative genes (71.9%) and 615,301 isoforms (67.0%) were identified among probing and nonprobing larval samples, respectively. Out of those identified in the probing subset transcriptome, 860 putative genes (0.3%) and 236 isoforms (0.04%) were identified as significantly differentially expressed (SDE) between probing and nonprobing samples (p ≤ .01; Figure S1). Only one gene was found to be differentially expressed among probing samples between Day 1 (D1) and Day 2 (D2) (it was upregulated on D2), and nine genes were found to be differentially expressed among nonprobing samples (two were upregulated on D1 and seven were upregulated on D2). These results strongly suggest that the observed differences in gene expression between probing and nonprobing larvae are not due to experimental artefacts but due to biological and behavioural differences (Figure S2).

In general, gene expression patterns were significantly different between probing and nonprobing samples but highly similar among samples within these two groups—at the isoform (Figure 2a) as well as the gene level (Figure 2b). Gene expression profiles were more similar among nonprobing samples than among probing samples, due to increased gene expression heterogeneity among overexpressed genes in probing compared with nonprobing samples (Figure 3). Most significantly differentially expressed (SDE) genes were upregulated in probing compared with nonprobing samples, indicating increased and additional activity; out of the 860 SDE genes, 837 (97.3%) were significantly upregulated in probing samples compared with 23 genes (2.7%) that were significantly downregulated (Figure 3). On average, 609 (72.8%) significantly upregulated genes were found in each probing sample and 19 (82.6%) upregulated genes were found in each nonprobing sample. One probing larva (sample PD1e) had only 165 significantly upregulated genes and fell in between nonprobing and the other probing samples, which might be due to biological variability or indicate that the larva was not actively probing at the time of collection (Figure 3). The most upregulated gene in probing samples was c267535_g1 (log2 fold change of 11.4), which matched solely to the P. lobata reference but returned no BlastX hit. The most upregulated gene in nonprobing samples was c260000_g2 (log2 fold change of 8.6), which did not match to either reference, and no BlastX hits were found.

Figure 2

Sample correlation matrix heat map of relatedness between samples, based on differential expression patterns of Trinity isoforms (a) and Trinity genes (b). Differential expression for sample comparison was determined based on p ≤ .01. The colour key represents a spectrum of most different (green) to most similar (red) in terms of relatedness. Sample relatedness trees are featured on the heat maps' left y‐axes and top x‐axes; the red and blue lines, respectively, correspond to probing and nonprobing samples

Figure 3

Samples versus gene heat map representation of differential expression patterns. Differential expression was determined at a significance level of p ≤ .01. The colour key represents a spectrum of lowest (purple) to highest (yellow) expression. Sample relatedness trees are featured on the heat maps' top x‐axes (red and blue lines, respectively, correspond to probing and nonprobing samples). Gene relatedness trees are on the heat maps' left y‐axes. Right y‐axes are not included but are available upon request

Characterization of differentially expressed genes

Out of the 860 significantly differentially expressed (SDE) genes, 674 (78%) could be matched to either one or both references. Interestingly, out of the 674 identified genes, only 41 SDE genes (6.1%) matched the P. lobata transcriptome exclusively (36 genes) or better (5 genes) than the S. microadriaticum genome and were thus considered coral genes. In contrast, 633 SDE genes (93.9%) matched the S. microadriaticum genome exclusively (613 genes) or better (20 genes) and were thus considered dinoflagellate genes. That means more than fifteen times more symbiont genes than coral genes were found to be significantly differentially expressed between probing and nonprobing samples (73.6% vs. 4.8% of all SDE genes; Figure 1d). This imbalance between coral and dinoflagellate genes remains if the 25 ambiguous genes that matched both references are excluded (71.3% vs. 4.2%; Figure 1c), and only genes that matched solely to one reference (36 coral and 613 symbiont genes) were included in downstream analyses.

Importantly, nearly half of all genes in our reference transcriptome were conservatively identified as coral (73,715 vs. 76,527), and there was higher coverage of predicted coral genes that matched to P. lobata compared with predicted symbiont genes that matched to S. microadriaticum (88% vs. 61%). Further, the majority of coral genes matched solely to one P. lobata reference gene (30 out of 36), and all but one P. lobata reference gene matched solely to one putative P. astreoides coral gene. Symbiont genes matched to approximately two S. microadriaticum reference genes, and S. microadriaticum genes, on average, matched to one putative symbiont gene. The imbalance between coral and symbiont SDE genes is therefore unlikely to be due to analytical artefacts but is likely a realistic representation of differential gene expression profiles during substrate probing.

The majority of differentially expressed coral genes (n = 36) were upregulated in probing samples (30 genes, 83.3% of SDE coral genes) compared with six upregulated genes in nonprobing samples (Table S3). Out of 36 SDE coral genes, 17 were successfully annotated (Table S4). Putative protein functions included macromolecule binding (9 genes, 53%), metabolism (4 genes, 24%), and signalling and transport (3 genes, 18%). A subset of coral genes was also linked to development (3 genes) and transcription regulation (2 genes). The most notable GO term for a coral contig is for xenobiotic metabolism (GO:0,006,805)—involved in the environmental stress response and proposed as a proxy for coral health (Rougee, Richmond, & Collier, 2014). All SDE coral genes but one (c240963, linked to nucleotide binding and signalling) were upregulated in probing larvae.

Similarly, all 613 SDE Symbiodinium genes were upregulated in probing samples (Table S3) and 523 were successfully annotated (85.3%) (Table S4). In addition to macromolecule binding (185 genes, 35%), signalling and transport (58 genes, 11%) and metabolism (20 genes, 4%), putative protein functions of note include responses to environmental stimuli (21 genes), transcription (10 genes), translation (8 genes), cell division and the cell cycle (6 genes), photosynthesis (5 genes) and immune response (5 genes). Notable GO terms for symbiont genes include heat shock protein binding (GO:0031072), glutathione transferase activity (GO:0004364), glutathione peroxidase activity (GO:0004602), glutathione metabolic process (GO:0006749)—glutathione is linked to regulation of oxidative stress (e.g., Asada, 1999)—and signal transduction by p53 class mediator (GO:0043516)—linked to cell cycle regulation (e.g., Evan & Littlewood, 1998).

Further, when compared in the DAVID enrichment analysis, significant functional enrichment was detected among holobiont SDE genes. Functional enrichment was found for proteins linked to intra‐ and extracellular transport, signalling, macromolecule recycling and subcellular localization to the plasma membrane and other membrane‐bound organelles (Table S5).

Overall, the transition from nonprobing to substrate probing is accompanied by significant expression‐level changes for a variety of coral and symbiont genes. Many of these differentially expressed genes are involved in nucleotide and protein binding, suggesting an increase in the activation of certain cellular processes, such as signalling pathways, transcription and translation, and cell division.

KEGG analysis of differentially expressed genes

The 17 and 523 SDE and annotated coral and symbiont genes, respectively, were subjected to KEGG pathway analyses. Of the 540 annotated genes, 130 (9 coral and 121 symbiont) had KEGG entries and a subset of 40 (4 coral and 36 symbiont) were found to be involved in 80 unique KEGG pathways.

The putative coral gene Magi3, which encodes for a membrane‐associated guanylate kinase that functions as a scaffolding protein at cell–cell junctions and is linked to cell signalling, was identified in the Rap1 signalling pathway (KEGG:rno04015)—linked to signal transduction and environmental information processing. Putative coral gene CYP3A4, belonging to the cytochrome P450 protein family of monooxygenases and linked to biosynthesis and metabolic functions, was identified in several KEGG metabolic pathways. It occurs for example in the linoleic acid metabolism (KEGG:hsa00591), as identified in previous studies of coral lipid composition (e.g., Latyshev, Naumenko, Svetashev, & Latypov, 1991; Lõhelaid & Samel, 2018), in the retinol metabolic pathway (KEGG:hsa00830; e.g., Tarrant et al., 2009) and in the metabolism of xenobiotics by cytochrome P450 (KEGG:hsa00980), which was also represented by putative symbiont gene (GSTT1). Two putative coral genes (FOLH1 and ddo‐1) and one symbiont gene (glnA3) were found in the KEGG pathway for alanine, aspartate and glutamate metabolism (KEGG:ssc00250, KEGG:cel00250 and KEGG:ddi00250), respectively. Putative coral oxidase gene ddo‐1 and symbiont transporter gene (SLC25A17) were also identified in the peroxisome KEGG pathway (KEGG:cel04146).

Three putative symbiont genes (nat, AK812_SmicGene4920 and ACS) were found in KEGG pathways for biosynthesis of secondary metabolites (KEGG:mbo01110, KEGG:syn01110 and KEGG:ath01110). Further, there were two putative symbiont genes (Ank1 and Ank2) found in the KEGG pathway for proteoglycans in cancer (KEGG:hsa05205 and KEGG:mmu05205). Proteoglycans have been found in plants and algae and are involved in a variety of cell surface processes, including cell adhesion and signalling (Knox, 2016). Two putative genes (TNXB and MYLK) were found in the KEGG pathway for focal adhesion (KEGG:hsa04510 and KEGG:ocu04510) (Table S4), which is involved in cell surface interactions, including cell anchoring and cell surface‐extracellular matrix communication. These cell surface proteins could play a vital role in host–symbiont communication during substrate probing in coral larvae.

DISCUSSION

The number of raw reads, raw reads per sample and clean reads per sample in the current study is exceptionally high compared with other published studies. It is therefore not surprising that the overall number of holobiont genes and isoforms (446,849 and 918,990, respectively; Table S2) is on the higher end of the spectrum commonly reported in coral holobiont studies (e.g., Mansour et al., 2016; Pinzón et al., 2015; Shinzato et al., 2014).

The number of putative coral and dinoflagellate genes in the present study (73,715 and 76,527, respectively) is higher than predicted in transcriptomic and genomic studies for corals (e.g., ~21,000; Bhattacharya et al., 2016) and for species in the Symbiodiniaceae family (37,000–49,000; reviewed in Aranda et al., 2016). However, the number of putative coral and symbiont genes is well within the range of genes assembled in other coral RNA‐Seq studies. For example, Pinzón et al. (2015) reported 178,943 coral contigs and 130,217 Symbiodiniaceae spp. contigs in their Orbicella faveolata holobiont transcriptome. The two published transcriptomes for P. astreoides range from 31,663 holobiont isogroups identified by Kenkel, Meyer, et al. (2013) to 95,294 coral genes and 145,570 dinoflagellate genes by Mansour et al. (2016).

The ratio of coral to dinoflagellate genes in coral holobiont transcriptomes is usually roughly 1:1, which is very similar to the ratio observed here (1:1.04). This observed balance between identified putative coral and symbiont genes is remarkable given the particularly strong imbalance between significantly differentially expressed (SDE) coral and symbiont genes in probing versus nonprobing larvae. Out of the 860 SDE genes, only 36 were conservatively identified as coral genes, while 613 were identified as symbiont genes, that is a ~1:17 ratio. This was completely unexpected, but it indicates that endosymbiotic dinoflagellates may play a major role in substrate probing and/or respond directly to larval probing behaviour.

Analytical artefacts may have partly influenced this dramatic shift in dinoflagellate gene expression. For example, the ratio of putative de novo genes to “reference” genes was 3.9:1 for coral genes but only 1.8:1 for symbiont genes, which means that on average every “coral transcriptome gene” was split into 3.9 de novo genes, while every “symbiont genome gene” was split into only 1.8 putative de novo genes. This difference could be due to regular molecular processes such as alternative splicing, which leads to multiple mRNA transcripts per gene (Berget, Moore, & Sharp, 1977). However, the difference between these two ratios could also be due to a more fragmented coral transcriptome, which may lead to biases in differential expression analyses if individual coral genes were split across more transcriptome contigs than symbiont genes. Our BUSCO analyses did indicate slightly more fragmented and duplicated genes in the coral (20.1% and 29.4%) compared with the symbiont data set (14.2% and 25.7%), but the difference was small (18 and 11 out of 303 genes, 5.9% and 3.6%, respectively) and indicates that gene fragmentation in the de novo transcriptome is unlikely to contribute much to the observed difference in differential gene expression. In addition, our ability to identify coral and symbiont genes might have been affected by having to rely on references of closely related species. Moreover, although most of the transcriptome was generated from individual larvae, the symbiont data still represent large pools of specimens since every coral larval hosts thousands of single‐cell endosymbionts, which can lead to pooling biases and associated inflation of differentially expressed symbiont genes. Variation in read depth per gene in the host versus the symbiont could be another bias of the DEG analysis. While our analyses indicate that both symbiont and coral reads were sufficient for differential gene expression analyses, almost twice the number of coral reads compared with symbiont reads were recovered (on average, 6.6 million vs. 3.4 million). We also found that coral isoforms above a set minimum threshold of zero transcripts per million transcripts (TPM) greatly outnumbered symbiont isoforms (28,720 = 6.5% of all coral isoforms vs. 5,940 = 2.2% of symbiont isoforms), but the average expression was higher for symbionts than corals (170 vs. 0.262, respectively). Nonetheless, these potential analytical biases should not affect SDE more than other genes and are thus unlikely to fully account for the observed 17‐fold difference in SDE dinoflagellate versus coral genes.

Our results indicate a significant increase in gene expression in endosymbiotic dinoflagellates during planula substrate probing, which suggests an influence of larval symbionts on substrate choice. Zooxanthellae are unequally distributed in a coral larva and are most abundant in the oral region (e.g., Huang, Wang, Chen, Fang, & Chen, 2008), while the larva's aboral part is most actively probing the substrate (Grasso et al., 2011; Tran & Hadfield, 2013). This reinforces our idea that zooxanthellae play a role in environmental sensing. For example, the photosynthetic dinoflagellates might assess the light environment at potential settlement sites, which is a vital characteristic for the holobiont's growth and survival, and relay this information to their larval host. For example, the putative symbiont SDE gene ELD1 is involved in “response to light stimulus” (GO:0009416) and “cellular response to light stimulus” (GO:0071482, Table S4). The role of the symbiont in substrate and settlement site choice might have been overlooked in the past and should be considered in future studies.

Holobiont genes and cellular pathways of interest

Differentially expressed genes between nonprobing and substrate‐probing larvae were linked to general cellular processes, such as macromolecule binding and transport, intra‐ and extracellular signalling, and increased metabolic function and regulation. Putative functions were also involved in the development and environmental sensing.

Previous studies on larval settlement have led to significant advances in our understanding of the molecular mechanisms that drive metamorphic changes in response to environmental cues in corals. Studies have demonstrated differential gene expression between early and late larval developmental stages (Grasso et al., 2008; Reyes‐Bermudez et al., 2016; Strader et al., 2018), including upregulation in genes linked to sensing and signalling (Strader et al., 2018) and calcification (Grasso et al., 2008). During metamorphosis, aposymbiotic larvae are found to be enriched in genes involved in transcription regulation, oxidative stress response and microtubule‐based movement (Meyer et al., 2011), and express different calcification‐linked genes before and after settlement (Grasso et al., 2011; Hayward et al., 2011). Putative genes involved in apoptosis, immune response (Grasso et al., 2008, 2011; Meyer, Aglyamova, et al., 2009; Reyes‐Bermudez et al., 2016) and integral signalling pathways (e.g., Notch, Wnt and G‐protein coupled signalling pathways) (Meyer et al., 2011; Meyer, Aglyamova, et al., 2009; Reyes‐Bermudez et al., 2016; Strader et al., 2018) have also been identified in sequenced planula larvae.

In addition, we identified putative holobiont genes of interest based on high differential expression levels and/or known links to signalling pathways, biological functions or environmental responses relevant to coral or more broadly marine larval settlement. For example, galaxin was significantly upregulated in probing samples. This gene was not included in downstream expression or annotation analyses since it matched both references, but it matched orders of magnitude better to P. lobata (E‐value ≤ 1 × 10−180) than to S. microadriaticum (E‐value = 3.21 × 10−05). Galaxin is a well‐known coral skeleton and skeletal extracellular matrix protein (Fukuda et al., 2003). Differential expression of Galaxin‐related genes has been identified in Acropora digitifera adults and postsettlement polyps (Fukuda et al., 2003), in A. millepora adults (Ramos‐Silva et al., 2013) and in larvae during settlement and metamorphosis (Hayward et al., 2011; Reyes‐Bermudez et al., 2009). Due to its function as a skeleton‐associated protein, the upregulation of galaxin in probing samples is likely in the preparation for postsettlement metamorphosis and skeleton formation. Adhesive secretions have also been observed in the early stages of skeleton formation, which may serve to attach the larva to the substratum (Goreau & Hayes, 1977), so galaxin could be a component of these secretions and play an important role in the attachment process.

Coral genes PAX3‐b and AdiPAXD2 were also significantly upregulated in probing samples. PAX (paired box) genes play a critical role in embryonic development (Blake & Ziman, 2014) and are known to associate with Wnt signalling pathways (e.g., Gan et al., 2014). PAXD genes are considered cnidarian equivalents to bilaterian PAX3 genes (Breitling & Gerber, 2000; Hadrys, DeSalle, Sagasser, Fischer, & Schierwater, 2005), which are linked to segmentation and neurogenesis in diverse Metazoan groups (Davis, Jaramillo, & Patel, 2001). PAX3 was found to be upregulated in coral larvae when responding to settlement inducers (Meyer et al., 2011) and in postsettlement polyps when compared to planula larvae (Hayward et al., 2011). Further, the coral gene glutaredoxin was among the most upregulated in probing samples. Glutaredoxin proteins regulate cellular redox state and redox‐dependent signalling pathways (Lillig, Berndt, & Holmgren, 2008; Meister, 1994), and have been shown to directly interact with Wnt signalling pathway components (Hanschmann, Godoy, Berndt, Hudemann, & Lillig, 2013; Hirota et al., 2000; Meyer, Belin, Delorme‐Hinoux, Reichheld, & Riondet, 2012; Murdoch et al., 2014; Sandieson, Hwang, & Kelly, 2014). The significant upregulation of glutaredoxin and PAX genes in probing samples suggests a complex regulation of the Wnt signalling pathway in probing larvae and could be the result of late‐stage larval maturation and development. The Wnt signalling pathway has also been found to potentially play a role in other marine larval groups. For example, Wnt proteins are involved in larval development in the sea anemone Nematostella vectensis (Marlow, Matus, & Martindale, 2013) and the marine annelid Platynereis dumerilii (Marlow et al., 2014), and during metamorphosis in the bryozoan Bugula neritina (Wong, Wang, Ravasi, & Qian, 2012).

The holobiont gene arylsulfatase G (ARSG) exhibited the most significantly downregulated gene in probing samples. It matched to both references but orders of magnitude better to S. microadriaticum (E‐value ≤ 1 × 10−180) compared with the P. lobata reference (E‐value = 1.43 × 10−16). The protein encoded by ARSG is a lysosomal enzyme (Kowalewski et al., 2014) and member of the sulfatase enzyme family (ARS), which are involved in hydrolysing sulphate esters of macromolecules for degradation (Ferrante, Messali, Meroni, & Ballabio, 2002). Arylsulfatase C enzyme activity has previously been monitored in adult corals to better understand environmental stress impacts (Rougee et al., 2014). The downregulation of ARSG suggests a decrease in symbiont lysosomal activity during substrate probing, potentially to account for increased integral compound retention and demand for other cellular functions. Arylsulfatases have also been identified in many genera of bacteria (e.g., Gao, Jin, Yi, & Zeng, 2015; Kim et al., 2004) and during larval development in molluscs and echinoderms (e.g., Degnan & Morse, 1995; Yang, Angerer, & Angerer, 1989). Additionally, gene DNAJ was significantly upregulated in probing samples and matched to the S. microadriaticum genome. DNAJ chaperone proteins (also known as Hsp40s) are well known for both preventing stress‐denatured protein aggregation and breaking down existing aggregates during the heat shock response (Qiu, Shao, & Wang, 2006). DNAJ and other related proteins have been identified in Symbiodiniaceae species (e.g., Aranda et al., 2016), and expression of DNAJ family and associated genes has been reported in coral‐associated endosymbiont species in response to heat stress (Baumgarten et al., 2013; Gierz, Forêt, & Leggat, 2017; Leggat, Hoegh‐Guldberg, Dove, & Yellowlees, 2007). In coral‐associated Symbiodinium, a DNAJ‐related gene was found to be upregulated under thermal and chemical stress (Yuyama et al., 2012) and DNAJ was found to be upregulated in planula larvae, compared with postsettlement polyps (Hayward et al., 2011). In other marine invertebrates, heat shock proteins have been implicated in morphologic changes (e.g., Bishop, Bates, & Brandhorst, 1998; Bishop, Bates, & Brandhorst, 2001; Gunter & Degnan, 2007; Kroiher, Walther, & Berking, 1992), which might indicate that the observed increase in DNAJ expression is not directly involved in substrate probing but rather in simultaneously occurring morphogenetic processes. For examples, HSP90 has been frequently related to cell growth and differentiation (e.g., Deane & Woo, 2003; Krone, Evans, & Blechinger, 2003; Krone, Sass, & Lele, 1997; Mayer & Bukau, 1999; Queitsch, Sangster, & Lindquist, 2002; Rutherford & Lindquist, 1998). However, in Acropora digitifera larvae, relative expression of Hsp90 decreased after the second day postfertilization, while Hsp70 did not change significantly (Nakamura, Morita, Kurihara, & Mitarai, 2012).

Individual larval transcriptomes

In this study, no significant differences in gene expression patterns between individual coral larvae and two pooled samples were detected (Figure 3). For example, only one gene was found to be significantly differentially expressed between probing day 1 (5 individual larvae) and probing day 2 samples (2 pools of 10 larvae each). Our results and conclusions are therefore unaffected by including two pooled and 15 individual larval samples.

It is however worth noting that this is the first study to date that generated transcriptome data sets for individual coral larvae, adapting a protocol designed and optimized for microarthropods (Fernández, Edgecombe, & Giribet, 2018). Previous transcriptomic studies pooled coral larvae (e.g., Mansour et al., 2016; Meyer, Aglyamova, et al., 2009) to extract enough mRNA for RNA‐Seq library generation. While pooling samples may reduce biological variability and thus increase the analytical power to detect SDE genes, it can also generate pooling bias (i.e., differences between pool measures and means of individual measures), which is particularly sensitive to increasing the number of false positives (Rajkumar et al., 2015). In addition, individual outlier samples (e.g., sample PD1e) cannot be identified and may introduce distortions and misleading information. The ability to generate RNA‐Seq libraries from individual coral larvae therefore offers a valuable new option for future differential gene expression studies of coral larvae and other small samples.

CONCLUDING REMARKS

This study characterizes differential gene expression between nonprobing and substrate‐probing Porites astreoides larvae. It finds significant and complex processes and shifts in holobiont gene expression during this important transitional life history stage, highlights the importance of distinguishing between nonprobing and substrate‐probing larvae in future larval development and settlement research, and connects observed gene expression patterns with broader trends seen in other marine organisms with a pelagic larval stage. It was predicted that significantly more genes would be upregulated in probing versus nonprobing larvae, due to the presumed increased complexity of the substrate‐probing process, compared with the free‐swimming stage. We found a sizeable subset of genes that exhibited significant differential expression, and the majority of these genes were upregulated in substrate‐probing larvae. Differential gene expression was also predicted to predominantly occur in cellular processes involved in signalling and metabolic pathways and cell growth and division, which was observed among both coral and symbiont genes. It was also found that differentially expressed holobiont genes were linked to environmental sensing and response, which would likely aid in identifying a suitable substrate for settlement. Additionally, the discovery that the overwhelming majority of identified differentially expressed genes were attributed to endosymbionts was unexpected and presents interesting and particularly significant implications for the potential roles of the endosymbiont in larval probing and settlement. Future analyses conducted on the coral holobiont should put more emphasis on symbiont–coral interactions and dynamics. An important future direction would be to conduct experiments that seek to definitively identify and characterize genes and their protein products between nonprobing and substrate‐probing larvae, with genes identified in this study serving as good candidates for further exploration, and to test whether different environments (e.g., depth and light variation) and environmental cues may trigger expression of alternative genes.

AUTHOR CONTRIBUTIONS

Larval settlement cue experiments were designed and conducted by V.J.P. and J.M.S. The field samples were preserved by V.J.P., J.M.S. and G.G. The study was designed by N.S.W., G.G. and D.C. mRNA extraction protocols were optimized by R.F. and N.S.W. The laboratory research was conducted by N.S.W, and the data were analysed by N.S.W. and D.C. The manuscript was written by N.S.W. and D.C., along with contributions from G.G. All authors read and approved the final manuscript.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.