RNA-Seq of the Caribbean reef-building coral Orbicella faveolata (Scleractinia-Merulinidae) under bleaching and disease stress expands models of coral innate immunity

Sample collection

A concurrent thermal anomaly, coral bleaching event and Caribbean yellow band disease outbreak occurred in 2010 in Puerto Rico. This event provided a unique opportunity to sample colonies of O. faveolata affected by multiple environmental stressors that are known to induce innate immune responses (Mydlarz et al., 2009; Pinzón et al., 2015). Six samples (approximately 25 cm²) from four colonies were collected on a single dive at 10 m depth in October 2010 on Media Luna reef in La Parguera, Puerto Rico (17°56.091 N, 67°02.577 W). Samples were collected under a permit issued by the Department of Natural Resources of Puerto Rico to the Department of Marine Sciences at the University of Puerto Rico at Mayaguez. Reefs in this region experienced ten degree-heating weeks at the time of sample collection. Degree-heating week is a remote sensing metric that estimates accumulated thermal stress in corals during sea surface temperature anomalies (Gleeson & Strong, 1995), and is reported by the National Oceanic and Atmospheric Administration. Five different health conditions were sampled: bleached (sample 1) and asymptomatic tissue (sample 2) of a partially bleached colony; asymptomatic tissue (sample 3) and lesion tissue (samples 4 and 5) from a CYBD-affected colonies; and tissue from a completely asymptomatic colony (sample 6) (Supplemental Information 1). The conditions represented by samples 3 and 4 have not been used for NGS by any previously reported investigation. Photographic examples of each disease condition are presented in Fig. 1. Within one hour of collection and storage at ambient temperature in seawater, tissue samples were transported to the Department of Marine Sciences on Isla Magueyes, flash frozen in liquid nitrogen, photographed while on dry ice, and stored at −80 °C. It was assumed that colonies sampled were non-clonal given their large distances of separation (>10 m) and low clonal levels (3.5%) previously documented for the same species on the same reef (Severance & Karl, 2006).

RNA extraction, sequencing, and de novo transcriptome assembly

The area of each tissue sample was estimated from photographs scaled with millimeter precision using ImageJ software (Schneider, Rasband & Eliceiri, 2012). The ratio of the sample area to volume of Trizol (Life Technologies, CA, USA) was the same for each sample: 2.0 mL of Trizol was added to every 1.00 cm² of sample tissue in 50 ml capped tubes. Tissue was homogenized by vigorous shaking until skeletons were completely denuded. A neutralization reaction occurs between the calcium carbonate skeleton of the coral and the acidic Trizol, so 1 to 5 µl of 6 M hydrochloric acid were added to each sample to minimize DNA contamination of the aqueous phase on the addition of chloroform. For each sample, 700 µl of the aqueous phase was loaded onto and spun through a single RNeasy column, and an on-column DNA digestion step was performed using DNase according to the manufacturer’s protocol (Qiagen, Velno, Netherlands). RNA was eluted from the column with nuclease-free water, and RNA quality and concentration were verified by 2% denaturing agarose gel electrophoresis and a NanoDrop spectrophotometer (Thermo Scientific, Waltham, MA, USA), respectively.

Total RNA was sent for mRNA multiplex sequencing on the Illumina GAIIx platform (Illumina Inc, San Diego, CA, USA) according to standard protocols of the genomic core facility at the Scripps Research Institute in Ft. Pierce, Florida, USA. All samples were sequenced in multiplex on two lanes of an Illumina flow cell. Raw sequence reads from each sample are available on the NCBI Short Read Archive under Bioproject PRJNA236103 and accession numbers for each sample are reported in Supplemental Information 1. From 70,745,177 raw sequence reads of 72 bases in length, adaptor sequences and low-quality bases were trimmed and clipped using cutadapt (Martin, 2011), which resulted in 59,114,519 reads with a mean length of 67 bases (standard deviation of 5 bases). Raw reads were trimmed and clipped to optimize sequence quality based on results of RNA-SeQC analysis (DeLuca et al., 2012). The quality of the data used for sequence assembly before and after processing is reported as a supplementary figure (Supplemental Information 2). The Trinity software suite was chosen for de novo assembly of the O. faveolata transcriptome using 59,114,519 processed sequence reads and default parameters (Grabherr et al., 2011; Haas et al., 2013).

Generation of the Orbicella faveolata reference transcriptome

To identify coral host and microbial sequences in the metatranscriptome, a series BLASTn alignments were conducted in parallel. For the first series, transcripts assembled by Trinity were aligned to transcripts previously reported for the photosynthetic endosymbiont, Symbiodinium spp. (Bayer et al., 2012; Shoguchi et al., 2013). Hits with an e-value less than 1E-3 were removed. For the second series, transcripts assembled by Trinity were aligned to O. faveolata expression sequence tags (EST) (E-value cutoff of 1E-6), Acropora digitifera mRNA sequences (E-value cutoff of 1E-3) (Shinzato et al., 2011), Nematostella vectensis mRNA sequences (E-value cutoff of 1E-3) (Putnam et al., 2007), and Hydra magnipapillata mRNA sequences (E-value cutoff of 1E-3) (Chapman et al., 2010). Transcripts with hits lower than the cutoff were assigned as O. faveolata in origin. Sequences downloaded for this analysis from their respective sources are provided as a supplementary file (Supplemental Information 3). Sequences that had significant hits to both Symbiodinium and Cnidarian sequences in both parallel analyses underwent a second-round of filtering using the classifier for metagenomic sequences (CLaMS). The complete genome of N. vectensis and whole transcriptome of Symbiodinium were used as training sets (Pati et al., 2011). Sequences that were binned as Symbiodinium only or both Symbiodinium and Cnidarian in origin were removed, and sequences that were binned as Cnidarian only were retained in the O. faveolata reference transcriptome. The resulting O. faveolata-specific reference transcriptome is provided as a compressed supplementary file (Supplemental Information 4a). A predicted proteome was generated using TransDecoder, a package within the Trinity software suite (Haas et al., 2013). The predicted proteome contains 20,488 sequences and is provided as a supplementary file (Supplemental Information 4b).

Transcriptome annotation and pathway analysis

To assess sequence accuracy of the de novo reference transcriptome, assembled transcripts were aligned to O. faveolata Expressed Sequence Tags (ESTs) in NCBI. O. faveolata ESTs have been submitted to NCBI by various sources. The sequences downloaded and used here are provided as a supplementary file (Supplemental Information 5). Transcriptome annotation was conducted using a combination of methods. Gene ontologies and EMBL/InterProScan protein motifs were assigned using BLASTx alignments on the BLAST2GO platform, which extracts annotations from best hits to non-redundant protein sequences with a maximum e-value of 1E-3 (Boratyn et al., 2013; Conesa et al., 2005; Quevillon et al., 2005). Taxonomic identities of best-hit sequences were extracted for a post-hoc assessment of sequence contamination in the O. faveolata-specific reference transcriptome. Transcripts were also annotated with KEGG Orthologues (KO) using default settings for the KEGG Automatic Annotation Server (KAAS) with a minimum BLAST score of 60 (Moriya et al., 2007). To determine the completeness of the O. faveolata transcriptome, KEGG KAAS annotations were conducted in parallel with mRNA sequences from six other Cnidarians: N. vectensis, H. vulgaris, Porites astreoides, A. millepora, A. digitifera, and P. damicornis (Chapman et al., 2010; Hemmrich et al., 2012; Kenkel, Meyer & Matz, 2013; Moya et al., 2012; Pinzón et al., 2015; Putnam et al., 2007; Shinzato et al., 2011; Traylor-Knowles et al., 2011). The reference transcriptome reported by Pinzón et al. (2015) was also annotated with KEGG KAAS in parallel for direct comparisons to the present data set.

Identification and analysis of putative immunity genes, proteins, and pathways

Correct open reading frames and domain architectures for predicted protein sequences were verified by hmmscan in the HMMER web server (Finn, Clements & Eddy, 2011). For phylogenetic analysis of predicted Wnt protein sequences, assignment to a specific Wnt family member was based on phylogenetic tree construction and clustering with previously described Wnt proteins for N. vectensis (Kusserow et al., 2005). Whole-length sequences were aligned by MUSCLE (Edgar, 2004), conserved regions were curated using Gblocks (Talavera & Castresana, 2007), maximum likelihood phylogenies were estimated using PhyML with 100 bootstraps and the Dayhoff substitution model (Dayhoff, Schwartz & Orcutt, 1978), and trees were constructed using TreeDyn (Chevenet et al., 2006). This pipeline was executed using the Phylogeny.fr platform (Dereeper et al., 2008). Immune signaling pathways were constructed by mapping assigned KO terms to KEGG pathway maps. Pathways were modified to illustrate the presence or absence of essential signaling components in the O. faveolata transcriptomes reported here and by Pinzón et al. (2015). Components of the miRNA and siRNA pathway gene list were selected based components reported to be conserved across invertebrates, including Cnidarians (Ding & Voinnet, 2007; Moran et al., 2013).

Transcriptome assembly and quality

The metatranscriptome sequences assembled by Trinity included 35,967 and 47,760 transcripts that were identified by BLASTn alignments as O. faveolata and Symbiodinium in origin. However, 3504 transcripts were identified as both coral and Symbiodinium by CLaMS analyses. Those sequences were removed thus producing an O. faveolata-specific transcriptome that contains 32,463 sequences. BLASTx alignment of these sequences to non-redundant protein sequences in SwissProt by BLAST2GO revealed most frequent hits to N. vectensis followed by other metazoans (Fig. 2). This provided a post hoc confirmation that contaminating sequences were successfully removed from the O. faveolata reference transcriptome. The size of the reference transcriptome is similar to previous studies that have used the Illumina GAII platform, which report between 33,000 and 48,000 unique coral transcripts (Barshis et al., 2013; Libro, Kaluziak & Vollmer, 2013). The GC content of the reference transcriptome is 44%, which is comparable to previous reports for corals (Sabourault et al., 2009; Soza-Ried et al., 2010; Vidal-Dupiol et al., 2013). Sequence accuracy was high with reference transcriptome sequences sharing 96% identity with corresponding O. faveolata ESTs in NCBI (Supplemental Information 6). The N50 was 1736 bp, which is comparable to recent studies that produced de novo coral transcriptomes (Barshis et al., 2013; Burge et al., 2013; Lehnert, Burriesci & Pringle, 2012; Libro, Kaluziak & Vollmer, 2013; Moya et al., 2012; Pinzón et al., 2015; Pooyaei Mehr et al., 2013; Shinzato, Inoue & Kusakabe, 2014; Sun et al., 2012). This value is low compared to the most recent sequencing efforts for corals, but this is likely an indication of the low number of total sequence reads rather than sequence quality. High quality of short reads is demonstrated by RNA-SeQC results, which compare the raw sequence reads to trimmed and clipped sequences used as the input for assembly by Trinity (Supplemental Informations 1 and 2). Collectively, these results demonstrate that the methods used were sufficient to assemble an accurate de novo reference transcriptome for O. faveolata.

Transcriptome annotation and novel sequences

BLAST2GO annotation assigned Gene Ontology (GO) terms and protein domain identities to 20,913 (64%) sequences with a maximum E-value of 1E-3 (Supplemental Information 7). A summary of the most abundant gene ontology terms for biological process, molecular function and cellular component is presented in Fig. 3. As expected, there was an abundance of transcripts associated with essential processes, such as metabolism, transcription, translation and protein complexes. In addition, there are transcripts associated with processes related to the physiological state of stress and disease, such as response to oxidative stress, death, immune system process, symbiosis, and response to wounding. Complete KO term and KEGG pathway annotations are provided as supplementary files (Supplemental Informations 8 and 9, respectively). To assess transcriptome completeness, a number of other Cnidarian transcriptomes were annotated in parallel for direct comparisons to the present data set (Table 1). Representation of metabolic and protein complex pathways for the present reference transcriptome was similar or slightly lower than data sets for other Cnidarian species. The results indicate that limited sequencing depth may underrepresent the full repertoire of possible O. faveolata transcripts. However, the present data set has a number of sequences associated with immune system-related pathways that are on par with or exceed the other data sets for Cnidarians (Table 1).

To highlight the unique contributions of the present investigation to the bioinformatic data available for O. faveolata, an in-depth comparison to the (Pinzón et al., 2015) O. faveolata transcriptome is presented. The 70,745,177 raw sequence reads from which the present reference transcriptome is derived is small in comparison to the 387,512,512 raw sequence reads reported by Pinzón et al. (2015). When annotated in parallel, we estimate that the two projects share 5,618 unique KO terms (79% of total). The present study contributes 524 unique KO terms (7% of total) and the latter contributes 995 unique KO terms (14% of total). Greater sequencing depth for the (Pinzón et al., 2015) reveals better coverage of KEGG pathways used conventionally to assess transcriptome completeness (e.g. metabolic and protein complex pathways) (Table 1). However, the present reference transcriptome had a greater number of transcripts mapped to the following stress and immunity-related pathways: MAPK, Ras, Notch, phagosome, Rig-i-like, bacterial invasion, nuclear factor kappa beta, complement, and leukocyte migration Pathways (Fig. 4). Therefore, the transcriptome presented here contributes a considerable number of sequences not reported previously that are associated with evolutionarily conserved pathways of the innate immune system.

The coral innate immune system

Models of coral immune response to disease are generally characterized by inflammation that involves the production of antimicrobial peptides and reactive oxygen species, the production of antioxidants to reduce self-harm, the migration of phagocytic cells to sites of infection, and the accumulation of melanin to prevent the spread of infection (Mydlarz et al., 2008; Palmer, Modi & Mydlarz, 2009; Tchernov et al., 2011; Vidal-Dupiol et al., 2011; Weis, 2008). In the arms race between invading pathogens and the coral host, a breakdown of host homeostasis leads to the activation of apoptosis and ultimately tissue mortality (Weis, 2008). Signaling pathways that control the coral immune system can be organized into 4 levels: (1) pathogen sensing by pattern recognition receptors, (2) downstream signaling cascades, (3) activation of inflammatory cytokine expression, and (4) effector mechanisms to lead to survival or death (Palmer & Traylor-Knowles, 2012). The present investigation provides new data for the study of these pathways in O. faveolata by revealing components of evolutionarily conserved immune signaling pathways that regulate immune cell development, migration and host-microbe interactions.

Wnt and Notch pathways in immune cell development, migration and communication

Gene families, such as Wnt and Notch, have been studied with respect to genome evolution and larval development in Cnidarians (Käsbauer et al., 2007; Kusserow et al., 2005; Marlow et al., 2012; Miller, Ball & Technau, 2005; Radtke, Fasnacht & MacDonald, 2010). It is also known that these genes are important regulators of the immune system through the control of immune cell differentiation, migration and communication (Duncan et al., 2005; Radtke, Fasnacht & MacDonald, 2010; Staal, Luis & Tiemessen, 2008). However, few studies have explored hypotheses about the role of Wnt and Notch in coral immune responses to disease. If one considers that pluripotent, phagocytic immune cells (i.e. coral amoebocytes) are at the front lines of wound healing and pathogen removal, it is likely that these genes regulate coral immune responses disease and stress. Therefore, pathways that have dual roles in the regulation of development and innate immunity should be researched further. To that end, we report novel members of the Wnt-like protein family and a complete gene set for the Notch pathway in the O. faveolata transcriptome.

Wnt proteins are extracellular ligands of transmembrane receptors, Frizzled (FZD) and low-density Lipoprotein Receptor (LRP), that transduce signals to control the expression of target genes required for development in a β-catenin-dependent and independent manner (Gordon & Nusse, 2006; Logan & Nusse, 2004). Components of the Wnt pathway were first reported for Cnidarians in the species Hydra vulgaris (Hobmayer et al., 2000), from which time their role in Cnidarian larval development has been studied extensively (de Jong et al., 2006; Guder et al., 2006; Kortschak et al., 2003; Kusserow et al., 2005; Miller, Ball & Technau, 2005; Randall & Szmant, 2009; Technau et al., 2005a). Therefore, it is not surprising that a set of genes putatively involved in Wnt signaling is present in the O. faveolata transcriptomes (Supplemental Information 10). However, the role of Wnt-signaling in coral amoebocyte development during migration and differentiation has not been investigated. A number of Wnt-like predicted protein sequences were identified by annotation of the reference transcriptomes reported here and by Pinzón et al. (2015). The parallel annotation of other Cnidarian species identified Wnt-like sequences in the transcriptomes of the corals Acropora millepora (Kortschak et al., 2003; Moya et al., 2012) and Acropora digitifera (Shinzato et al., 2011), and the sea anemone, Aiptasia pallida (Lehnert, Burriesci & Pringle, 2012). Twelve Wnt proteins were originally identified in the N. vectensis genome (Kusserow et al., 2005). These sequences were included in the phylogenetic analysis of Wnt sequences to guide homology prediction for Wnt sequences in the O. faveolata transcriptome (Supplemental Information 11). Seven unique O. faveolata transcripts are reported here as putative homologues of the following N. vectensis Wnt proteins: Wnt1, Wnt5, Wnt6, Wnt7, Wnt8, Wnt10, and Wnt11 (Fig. 5). The extent to which Wnt protein expression is restricted to specific developmental life stages or environmental stimuli is not well understood in corals. Therefore, the proteins described here may only represent a subset of the complete repertoire encoded in the genome of this species. A complete genome sequence or a reference transcriptome derived from a more diverse set of tissue types and environmental conditions is required for an exhaustive assessment of O. faveolata Wnt-like genes. The present study is the first to characterize O. faveolata Wnt-like protein sequences. This new data can be used in future mechanistic studies on the role that Wnt protein family plays in the coral innate immune system.

Notch signaling is involved in the regulation of cell identity, proliferation, differentiation and apoptosis (Gazave et al., 2009). Notch proteins are transmembrane receptors that detect delta-like ligands expressed on the surface of adjacent cells. On engagement of its ligand, the intracellular domain of Notch is cleaved and translocates to the nucleus where it serves as a transcription factor. Notch receptors have been described previously for various Cnidarians, including Nematostella, Hydra, and Acropora (Dunlap et al., 2013; Käsbauer et al., 2007; Marlow et al., 2012; Münder et al., 2010; Technau et al., 2005a). Inhibition of the Notch pathway was recently shown to disrupt wound healing in N. vectensis thus establishing a role for Notch in the regulation of innate immunity (DuBuc, Traylor-Knowles & Martindale, 2014). A role of Notch signaling in coral immune responses to disease is also supported by an RNA-Seq investigation that reported the differential expression of a Notch-like gene in response to acute exposure to bacterial Pathogen Associated Molecular Patterns (PAMPs) (Weiss et al., 2013). A putative Notch signaling pathway based on the presence or absence of pathway components in the O. faveolata transcriptome is presented in Fig. 6. This data can be used to investigate the role of Notch signaling in wound healing and immune responses to disease in O. faveolata.

Intracellular host-microbe interactions

Signaling through innate immune pathways is activated on detection of extracellular PAMPs by Pattern Recognition Receptors (PRR). PRR-induced signals are transmitted by cytosolic kinases, such Mitogen Activated Protein Kinase (MAPK). MAPKs make up a highly conserved family of protein kinases that regulate immune signaling pathways (Chen et al., 2001; Schwarz et al., 2008). Signals transmitted by kinases activate pro-inflammatory gene expression by transcription factors that regulate effector responses, such as Nuclear Factor Kappa Beta (NF-κB) and Jun (Palmer & Traylor-Knowles, 2012). A complete prototypical inflammatory signaling cascade has been documented in Hydra in the context of host-associated microbiota. In this model, PAMP recognition by PRRs, such as Toll-Like Receptors (TLRs), activates MAPK signaling cascades through MyD88-dependent phosphorylation of Jun-Kinase (JNK). JNK in turn activates Jun and leads to the expression of pro-survival factors, such as Bcl-2 (Franzenburg et al., 2012). In O. faveolata, the expression of genes associated with this pathway has been detected under conditions of environmental stress (Schwarz et al., 2008; Voolstra et al., 2009). Mechanistic studies to confirm a functional role for these genes in O. faveolata immune responses to disease can now be conducted with the use of full-length gene sequences reported here and by Pinzón et al. (2015).

Host cells can also detect PAMPs associated with intracellular pathogens by Nod-Like Receptors (NLRs) and Rig-Like Receptors (RLRs) (MacKay, Wang & Kurt-Jones, 2014; Ting et al., 2008; Yoneyama & Fujita, 2007). NLRs detect components of bacterial cell walls and RLRs detect cytosolic viral nucleic acids. Various bacterial pathogens of corals, including the causative agents of CYBD, invade the host cytoplasm (Cervino et al., 2008; Cervino et al., 2004; Kushmaro et al., 2001). Therefore, NLRs should be investigated in coral immune responses to infectious disease. Recent investigations have also demonstrated the important roles that viruses play in pathogenesis and symbiosis (Atad et al., 2012; Barr et al., 2013; Davy et al., 2006; Soffer et al., 2013; Weynberg et al., 2015; Wilson et al., 2005). Therefore, RLRs, and the related dicer-like protein family, should be investigated to elucidate their role in the control of host interactions with symbiotic, commensal and pathogenic viruses. To support future studies along these lines, the present investigation identifies components of NLR, RLR and Dicer signaling pathways present in the O. faveolata transcriptome.

Nod-like receptors and intracellular bacterial detection

The NLR protein family has conserved roles in regulating innate immunity through the recognition of intracellular bacteria by Leucine-Rich Repeat (LRR) domains. LRR domains are also hallmarks of membrane-bound TLRs, however, NLRs are located exclusively in the cytoplasm (Ting et al., 2008). After pathogen recognition, NLRs can promote caspase-mediated cell death through the macromolecular assembly of the inflammasome, which converts inactive procaspase to active caspase by proteolysis. Caspase-mediated cell death has been reported in corals and Cnidarians and has been described in depth by multiple reviews (Cikala et al., 1999; Lasi et al., 2010; Tchernov et al., 2011; Weis, 2008). Alternatively, NLRs can mediate inflammation through the aforementioned MAPK, JNK, or NFκB pathways. Genes homologous to NLRs have been identified in basal metazoans including Hydra, Nematostella, Acropora, Amphimedon (Augustin, Fraune & Bosch, 2010; Bosch, 2012; Hamada et al., 2013; Yuen, Bayes & Degnan, 2014). In the Hydra model system, there is evidence to support a role for NLR in bacterial detection and caspase activation (Bosch et al., 2011; Lange et al., 2011). Therefore, it is not surprising that we also identified NLR homologues in the O. faveolata transcriptomes. A putative O. faveolata NLR signaling pathway is constructed to guide future investigations on coral immune responses to intracellular bacterial pathogens (Fig. 7).

Detection of viral nucleic acids

To date, interactions between the coral host and associated microbes have largely focused on coral-fungal and coral-bacterial interactions (Kim et al., 2000; Kvennefors et al., 2008; Sutherland, Porter & Torres, 2004). However, viruses have emerged as important members of the coral holobiont and can have beneficial or harmful affects on the coral host (Van Oppen, Leong & Gates, 2009, Weynberg et al., 2015). While bacteriophages have been demonstrated to eradicate coral pathogens and prevent infection (Barr et al., 2013; Efrony et al., 2007), herpes-like viruses have been associated with virulent coral diseases (Soffer et al., 2013; Thurber et al., 2008). Therefore, recognition of viruses and regulation of their activities by the coral immune system is hypothesized to be essential in maintaining homeostasis.

Dicer-like proteins, which share protein domains with RLRs, are essential members of the Micro RNA (miRNA) and Small Interfering RNA (siRNA) pathways. They have only recently gained attention in corals, and they may have conserved roles in antiviral immunity (Liew et al., 2014; MacKay, Wang & Kurt-Jones, 2014; Moran et al., 2013). The present study identifies full-length Dicer-like proteins and the essential components for siRNA and miRNA signaling in the O. faveolata transcriptomes (Fig. 8). The dicer-like protein sequences used in this study are provided in a supplementary file (Supplemental Information 12). RLRs have evolutionarily-conserved roles in recognition of viral nucleic acids (Loo & Gale, 2011; Mukherjee, Korithoski & Kolaczkowski, 2014; Zou et al., 2009). Homologues of RLR-like genes have been reported for Nematostella in studies of RLR evolution from basal metazoans to chordates (Zou et al., 2009), but their function has not been investigated in Cnidarians. The O. faveolata reference transcriptome has many of the evolutionarily conserved components required for RLR signaling, and a putative pathway is presented in Fig. 9 to guide future research on coral immune responses to viral infection.