A new transcriptome and transcriptome profiling of adult and larval tissue in the box jellyfish Alatina alata: an emerging model for studying venom, vision and sex

Sample collection

Alatina alata material was collected during a spermcasting aggregation in Bonaire, The Netherlands (April 23-25, 2014, 22:00-01:00) according to the methods in [2]. A single ovulating female medusa (Fig. 1a) was dissected to provide four tissue samples, namely: i. gastric cirri (Fig. 1b), ii. ovaries (Fig. 1c) (within gastrovascular cavity filled with sperm (Fig. 1d)), iii. tentacle (containing nematocysts Fig. 1e) and adjoining pedalium base, which we refer to collectively as the “tentacle sample” below (Fig. 1a), and iv. rhopalium (including rhopaliar stalk) (Fig. 1f). A fifth sample consisted of planulae (Fig. 1h) that developed from blastulae (Fig. 1g) released by females in the lab. Transcripts putatively involved in venom were characterized through their apparent up-regulation in the tentacle and gastric cirri samples; transcripts putatively involved in vision were targeted through analysis of the rhopalium sample; and transcripts putatively involved in sex and embryogenesis were investigated through analysis of the ovaries and planulae samples. We sought to identify the possible onset of expression of candidate genes in the larval planula stage. Given the microscopic size (~150 μm) of the planulae, multiple individuals were pooled to obtain sufficient tissue for RNA isolation and sequencing. Detailed protocol is provided in Methods.

RNA-Seq and bioinformatics

De novo transcriptome assembly

RNA-Seq was performed on five different tissues using the Illumina HiSeq2500 Sequencing System (see Methods). Using Trinity software [50, 51] the ~264 million trimmed raw paired-end sequence reads were assembled de novo into a pooled transcriptome yielding ~126 KTrinity transcripts corresponding to ~84 K Trinity genes with an N50 of 1994 (Table 1). Throughout this paper, we use the term “A. alata gene” to refer to each transcriptome component represented by a unique Trinity gene id, and the term “A. alata transcript” to refer to additional transcriptome components that Trinity assigned as multiple putative “isoforms” of a single unique gene id [50, 51]. Transcriptome completeness was assessed using the subset of 248 widely conserved eukaryote core genes (with low frequency of gene family expansion) using the program CEGMA (Core Eukaryotic Genes Mapping Approach) [52, 53]. We retrieved 242 complete CEGs (98 %) and an additional three partial CEGs (1 %), resulting in 99 % CEG representation. We sought to identify highly expressed transcripts which, by definition, are represented by more reads (than lowly expressed transcripts), and thus have a better chance of being contiguously assembled. Therefore we generated a filtered transcriptome comprising a subset of transcripts expressed above a minimum threshold of 1.5 fragments per kilobase per million fragments mapped (fpkm) [46], based on read quantification and alignment accuracy using RSEM [54]. This filtering step is consistent with our aim of identifying highly expressed transcripts for potential candidate genes across different samples types. The filtered transcriptome, hereafter referred to simply as the A. alata transcriptome, yielded ~32 K transcripts corresponding to ~20 K genes; N50 = 2545 (Table 1). The percentage of sequences above 1000 bp doubled and the percentage of short genes (200–500 bp) was reduced by half (Additional file 1). Figure 2 illustrates the workflow used in this study (modelled after [47]).

	Whole transcriptome	Filtered transcriptome (fpkm = 1.5)
No. of transcripts	126,484	31,776
No. of genes	84,124	20,173
Total assembled bases (bp)	125,647,941	48,556,932
Avg (mean) transcript length (bp)	993	1528
Median transcript length (bp)	456	989
Max transcript length (bp)	9993	9993
N50	1994	2545
GC content (%)	39	40
Percent proper pairs	77	78
Samtools percent mapped and paired	75	77

De novo assembly began with 264,505,922 trimmed raw reads generated from medusa gastric cirri, ovaries, tentacle (with pedalium base), rhopalium, and planulae samples

Gene expression patterns and profiles

We then sought to detect gene expression patterns across the five samples (gastric cirri, ovaries, tentacle with pedalium base, rhopalium, and planulae) with the aim of providing a descriptive analysis of the top expressed gene clusters by sample. In order to estimate transcript abundance we aligned each set of reads back to the A. alata transcriptome and generated an RNA-Seq fragment counts matrix for each sample using RSEM [54]. We subsequently identified differentially expressed genes (see Methods) which we clustered according to their expression profiles using hierarchical clustering analyses within the framework of the EdgeR Bioconductor software package [55], a preferred methodology for studies lacking biological replicates [46, 56]. Of the ~32 K Trinity transcripts (~21 K Trinity genes) identified in the A. alata transcriptome ~10 K transcripts (6676 genes) were found to be differentially expressed, within a broad range, across the five samples (Additional file 4). EdgeR takes the normalized gene counts for all samples (generated from the initial RSEM counts matrix), and then clusters genes with similar mean expression rates across samples [56]. Gene clusters are visualized in the form of a heatmap, permitting pinpointing of genes abundant in certain samples that might be of interest as candidate genes. The results of hierarchical clustering were consistent with our initial RSEM evaluation of abundant genes by sample.

In an attempt to identify genes that were specifically abundant in each medusa sample, we conducted subsequent hierarchical clustering with the planulae sample excluded. EdgeR analyses of genes from just the four medusa samples (gastric cirri, ovaries, tentacle (with pedalium base) and rhopalium) revealed 2916 differentially expressed genes. To identify patterns of highly expressed gene clusters by medusa sample, we generated a heatmap of the subset of 2916 genes and further partitioned the subset into 10 gene subcluster profiles based on mean expression patterns across samples (Fig. 3a-k; Additional file 5). Furthermore, redoing the hierarchical clustering and subcluster profiling analyses for all five samples using the three subsets of candidate genes (putative venom, vision and sex genes) allowed us to hone in on gene clusters that were relevant to transcriptome functional annotation and profiling of different A. alata samples.

Tissue-specific “core genes”

Our next objective was to identify specific genes that are highly expressed in particular body parts (abundant with respect to other samples), and subsequently identify functional categories informed by the Trinotate report. We constructed a Venn diagram of the 2916 differentially expressed genes from the four medusae samples (Fig. 4a), which revealed that 76 % (2228 genes) were expressed to some degree in all medusa samples (Additional file 1). A subset corresponding to the top 50 most highly expressed genes per sample was plotted in a second Venn diagram (Fig. 4b), and genes unique to each sample’s top 50 were respectively designated as that sample’s “core genes” (Additional file 6). Only one gene (lacking a Trinotate annotation) among each sample’s top 50 most highly expressed genes overlapped in all four samples. An assessment of the core genes revealed that only about 40 % had Trinotate-generated annotations (Fig. 5a-d). Gastric cirri core genes (n = 46) corresponded mostly to putative proteins (26 annotations) implicated in venom and digestion, including metalloproteinases (Fig. 5a) [5, 57, 58]; ovaries core genes (n = 33) corresponded to putative proteins (22 annotations) involved in gametogenesis, including Vitellogenins (Fig. 5b), which, although well studied in bilaterians [59], are newly reported in cubozoans herein; tentacle core genes (n = 34) corresponded to many putative proteins (19 annotations) associated with nematocyst development, including minicollagens (Fig. 5c) [60–62]; and rhopalium core genes (n = 20) corresponded to several putative proteins (12 annotations) identified in vision and the phototransduction pathway, including J-crystallins (Fig. 5d) [28, 33, 63].

Candidate gene profiling

Using the tissue-specific core gene annotations (see above) as a starting point for identifying candidate genes in A. alata medusa body parts, we compiled three candidate-gene lists comprising 148, 109 and 39 terms (related to genes and/or gene functions) broadly associated with “venom”, “vision”, and “sex”, informed primarily by relevant genes and proteins documented in the scientific literature (Additional files 7, 8 and 9: provides lists of terms, references and differential expression matrices for all putative candidate genes). Subsequently, we queried the A. alata Trinotate report separately using each of the lists to identify matching terms among BLASTX, BLASTP and PFAM (protein family) top hits from the Trinotate transcriptome annotation report. This generated three additional targeted Trinotate annotation reports consisting of gene subsets we categorized as “venom implicated genes” (Additional file 7), “vision implicated genes” (Additional file 8), and “sex implicated genes” (Additional file 9). Subsequently, EdgeR hierarchical clustering (see above) was conducted separately on each subset of corresponding candidate genes, generating three new heat maps (Figs. 6, 7 and 8); profiling the respective patterns of expression of each set of candidate genes across all five samples. In all cases, gene cluster patterns were divided into 10 subclusters based on mean expression patterns for genes with potential association with venom (n = 450) (Fig. 6b-k); vision (n = 97) (Fig. 7b-k); and sex (n = 104) (Fig. 8b-k). By comparing planulae samples with samples from body parts of the mature female medusa, an additional aim of this study was to detect the potential onset of expression of these possible candidate genes within developing planulae. Differential expression matrices for all candidate genes and their respective annotations are provided in Additional files 7, 8 and 9. We did not detect expression of most candidate genes in the planulae sample. Instead, we found that transcripts most highly expressed in the planulae comprised mainly: histones (core and early embryonic), nanos transcription factor, and genes implicated in neurogenesis, mitosis, microtubule, and protein processing (Additional file 10).

Putative venom implicated genes

Here we highlight our findings of the 450 transcripts we broadly refer to as “venom implicated genes” based on preliminary candidate gene profiling (above). This effort focused specifically on identifying genes highly expressed in the tentacle (used in prey capture) and gastric cirri (used in digestion). By comparing a body part with tissue abundant in penetrant nematocysts (tentacle and adjoining pedalium base) (Fig. 1e) with one lacking nematocysts (gastric cirri Fig. 1b) our aim was to identify putative site(s) of venom production and nematogenesis in A. alata. Although nematocysts are typically abundant in the gastric cirri of many box jellyfish species [15], only a single individual nematocyst has been documented in the gastric cirri of mature A. alata medusae [2] despite examination of hundreds of mature specimens in several independent studies [2, 15, 38, 40]. Conversely, nematocysts (Fig. 1e) are primarily concentrated in the tentacle (which is contiguous with the pedalium) of A. alata [2], as is the case with all cubozoans [15]. Hierarchical gene cluster profiling (Fig. 6 a-k) revealed that many of the putative venom implicated genes were fittingly highly expressed in the tentacle (Figs. 5c and 6e-i), but were surprisingly also highly abundant in the gastric cirri (Figs. 5a and 6b,c; Additional file 7).

CaTX/CrTX toxin family genes

We identified eleven different homologs of the CaTX/CrTX toxin family (also annotated individually as CrTX-A or CaTX-A). These were either abundant almost exclusively in the tentacle (n = 4) (Figs. 5c and 6e, f) or in the gastric cirri (n = 7) (Figs. 5a and 6b, c). This gene family consists of pore-forming toxins that cause pain, inflammation and necrosis during human envenomation [64] and prior to this study has exclusively been associated with venom from nematocysts [18, 19, 65]. The taxonomically restricted CaTX/CrTX toxin family [66, 67], previously called the “box jellyfish toxin family” [64, 66, 67], was thought to be restricted to medusozoans [5], but recently a homolog was also identified in the anthozoan coral Acropora [68]. Gene tree reconstruction (Fig. 9) of the eleven CaTX/CrTX gene homologs abundant in either the gastric cirri or the tentacle confirmed homology of the A. alata transcripts with CaTX/CrTX toxin family genes in other cnidarians [5, 19, 68]. The analysis recovered four well-supported groups of A. alata CaTX/CrTX genes, each exclusively containing transcripts with tissue-specific expression patterns, either in tentacle or gastric cirri (Fig. 9). One group includes three A. alata homologs (annotated as CaTX/CrTX, CaTX or CrTX) specific to the tentacle that group with several non-cubozoan medusozoans including the coral Acropora. Three additional groups are all within a well-supported cluster of CaTX/CrTX genes identified from cubozoan taxa. One of these nested groups includes the homolog (CaTX-A) reported more than a decade ago [65] in A. alata (as Carybdea alata). Homologs of this gene are sister to a sub-group comprised of Chironex fleckeri homologs, which have been identified exclusively from tentacle tissue. The two additional nested groups (Fig. 9) include the homolog (CrTX-A) reported more than a decade ago in Carybdea brevipedalia (reported as C. rastonii), consisting of transcripts only identified in our gastric cirri sample, in a well-supported cluster of homologs derived from Carybdea brevipedalia and Malo kingi. Ours is the first report of expression of the CaTX/CrTX toxin family in a medusozoan body part that lacks nematocysts.

Nematocyst structural genes

Genes encoding putative nematocyst structural proteins were also characterized in this study. This is in line with our aim to characterize molecular components of cubozoan nematocysts and pinpoint putative regions of nematogenesis, given the current view that venom deployment in medusozoans is exclusively controlled by nematocysts. Our findings revealed that three minicollagens, key components in nematocyst capsule development [62], were abundant almost exclusively in the A. alata tentacle (with adjoining pedalium base) (Figs. 5c and 6e, h), with slight expression signal in all other medusa samples. All three of the minicollagen genes identified from A. alata possess the characteristic collagen-like domain of short repeated tripeptides of the form Gly-X-Y flanked on both sides by proline repeats by N- and C-terminal cysteine rich domains (CRD). The CRDs for two of the three genes are of the regular form (CXXXCXXXCXXXCXXXCC) and thus would be classified as Group 1 minicollagens [69]. The third has a regular CRD at the C-terminus but a variant form at the N-terminus that we refer to as Group 2 variant. Gene tree reconstruction of minicollagen genes for cnidarian taxa (Fig. 10) show that the minicollagens identified for A. alata cluster primarily with other Group 1 minicollagens. Of the additional nematogenesis related genes [17, 61, 62, 70–74] expressed in this study, most were almost exclusively abundant in the tentacle (with adjoining pedalium base). They include “nematocyst outer wall antigen” (NOWA) (Fig. 6e), “chondroitin proteoglycan 2” (Fig. 6f), “nematoblast-specific protein nb035-sv2/nb035-sv3/nb012a” (Fig. 6e, f), “nematogalectin-related protein” (Fig. 6e), and “Dickkopf-related protein 3” (Fig. 6e, f, h). Five different “Dickkopf-related protein 3” homologs were abundant in A. alata tentacle sample; with only two having notable expression in other body parts (gastric cirri, ovaries or rhopalium).

Putative vision implicated genes

Here we highlight our findings of the 97 transcripts we broadly refer to as “vision implicated genes” based on preliminary candidate gene profiling (above). This effort focused specifically on genes expressed in the rhopalium of A. alata, which bears a pair of lens eyes with cornea and retina, two pairs of simple ocelli-comprising photoreceptors, and a statocyst (Fig. 1f). By comparing the rhopalium with its visual capabilities and planulae with its known eye-spot photoreceptors, against the medusa samples that lack known photoreceptors (gastric cirri, ovaries and tentacle), our aim was to identify the expression of opsins and other vision implicated genes in the rhopalium of A. alata, as well as in putative extraocular photoreceptors in A. alata. Hierarchical gene cluster profiling (Fig. 7a–k) revealed that most of the 97 putative vision implicated genes (see “gene-profiling” above) were abundant in the rhopalium (Fig. 5d), but in many cases they were more highly expressed in other samples, in particular in medusa samples (Fig. 7c, e-j; Additional file 8).

Opsins

In this study the Trinotate report corresponding to the A. alata transcriptome contained a total of 41 transcripts with PFAM annotations corresponding to homologs of the “7 transmembrane receptor (rhodopsin family)”. Of the rhodopsin family, opsins are considered universal light sensitive proteins associated with photoreceptor cells of animal retinas. We found eleven opsin homologs to be variably expressed across A. alata medusa samples, with only six homologs most abundant in the rhopalium (Figs. 6g, h and 7i, h). A gene tree reconstruction of all rhodopsin family cnidarian genes (Fig. 11), which we rooted on the group that includes all previously known medusozoan opsins, recovered two of the three previously identified cnidarian opsin groups [30]. Group A includes only anthozoan taxa, while Group B includes all opsins previously known from medusozoans, eleven transcripts from Alatina and a few from anthozoans. However, cnidarian opsins previously identified as Group C fell into two groups—one of which includes thirty A. alata transcripts. This appears to be the first example of medusozoan opsins outside group B cnidarian opsins.

Among the top 10 most abundant genes in the rhopalium of A. alata was a homolog of the cubozoan lens-eye opsin for Carybdea rastonii (=C. brevipedalia) (Figs. 5d and 7i). The Carybdea lens eye opsin was also highly expressed in all A. alata samples, including planulae which have eye spots (Fig. 1h). Normalized counts revealed three additional opsin genes that were not differentially expressed across all samples, were almost exclusively found in the planulae sample (Fig. 11). Only two A. alata putative rhodopsin family homologs were expressed almost exclusively in the rhopalium (Trinotate top BLAST hits: “dopamine receptor 2” and “visual pigment-like receptor peropsin”) (Fig. 11). Among the putative “rhodopsin family” genes expressed in medusa samples, including those annotated in the Trinotate report as non-opsin based photoreceptors, were: “compound eye opsin bcrh1/d(1b) dopamine receptor” (Fig. 7e) and “opsin rh1/mu-type opioid receptor” (Fig. 7h), “blue-sensitive opsin” (Fig. 7k), “visual pigment-like receptor peropsin” (Fig. 7d), “melanopsin-b” (Figs. 7b and 8b; Additional files 8 and 9).

Oogenesis and embryogenesis

We found that homologs for the putative large lipid transfer protein Vitellogenin-2 were most abundant in A. alata ovaries (Figs. 6i and 8e), though highly expressed in all medusa samples. Conversely, genes annotated as Vitellogenin-1 or simply Vitellogenin were expressed in all medusa samples, but most abundant in the tentacle (Fig. 6i, g). Likewise, genes annotated as Apolipophorin or Apolipoprotein B-100, the other major animal protein group involved in lipoprotein processing [86], were expressed across all medusa samples but most abundant in the tentacle (Fig. 6g-i). Genes of the Vitellogenin family are responsible for lipid transfer from ovarian follicle cells to oocytes, providing nutrition during embryogenesis in bilaterians and some cnidarians [86, 87], and have a documented role as egg yolk protein precursors in the ovaries of animals, including anthozoans [88]. Gene tree reconstruction of A. alata transcripts with known cnidarian Vitellogenin and Apolipophorin-like proteins (Fig. 12) confirmed their homology with other cnidarian large lipid transfer proteins.

We report the abundance of multiple creatine kinase isozymes in the ovaries (Fig. 5b), all of which were expressed to some degree in all samples including planulae (Figs. 6d and 8h). Creatine kinase activity has a documented role in oogenesis and early embryogenesis in mammals [89]. More broadly though, creatine kinase is important in cells with variable rates of energy turnover, such as muscle, neurons, photoreceptors, and primitive spermatozoa [90], which is consistent with its expression in all samples in this study.

Prey capture and defense

In cubozoans, and more broadly in all cnidarians, prey capture and defense are based on nematocyst (stinging organelles) and associated venom. By comparing a body part abundant in penetrant nematocysts (tentacle and adjoining pedalium base) with one lacking nematocysts (gastric cirri) our aim was to identify putative site(s) of nematocyst development (nematogenesis) and venom production in A. alata.

We found that transcripts corresponding to a number of putative nematocyst structural proteins (minicollagens, nematogalectin, NOWA, chondroitin, and Dickkopf homologs) were abundant in the tentacle (and adjoining pedalium base). Although, putative nematogenic transcripts were detected primarily in the tentacle, some were also detected in non-tentacle medusa samples. We expect that this signal stems from the abundant adherent nematocysts covering the medusa bell. Together these findings are consistent with nematogenesis in A. alata occurring primarily, but not solely, in the region comprising the tentacle and adjacent pedalium base. Future in situ hybridization studies employing genes identified in nematogenesis in this study can help pinpoint more precisely nematogenic regions in A. alata.

Venom is a complex cocktail of bioactive compounds (e.g., protein and/or peptides called toxins, salts and neurotransmitters) secreted by one animal that is delivered to another animal by an infliction [57, 95]. Venom disrupts physiological and biochemical molecules of prey and predators, thus facilitating feeding and defense [57]. Nematocysts have long been considered the sole secretory structure for venom deployment in cnidarians [96]. However, we found preliminary evidence for venom production in the gastric cirri, where nematocsyts are lacking in mature A. alata. Futhermore, we found that the gastric cirri and tentacle express distinct groups of homologs of a major family of cnidarian venom proteins, the CaTX/CrTX toxin family [66, 67]. This suggests that venom plays an important, and possibly different role in the gastric cirri and tentacle. Venom components likley differ between the nematocyst-bearing tentacle, with a primary role in immobilizing prey and warding off predators, and the gastri cirri, with a primary role in killing and digesting prey [13].

Based on our findings, we hypothesize that A. alata has gland cells that secrete toxins associated with the gastric cirri. Evidence was recently presented for toxin-secreting gland cells in the ectoderm of the sea anemone Nematostella in regions containing nematocysts as well as areas that may lack nematocysts [97, 98], but our findings represent the first putative case in a cubozoan. Future morphological studies examining the ultrastructure of the stomach and gastric cirri, and venom gene candidate localization studies, will permit testing of the hypothesis of toxin-secreting gland cells associated with the gastric cirri of A. alata.

Although differences exist in the exact complement of putative bioactive toxins between the gastric cirri and tentacle sample, the venom cocktail in each body part includes transcripts from similar digestive enzyme families. A recent review of jellyfish toxins lists a number of toxin-like digestive enzymes that are deployed as components of nematocyst venom to disable homeostatic processes in prey or predators [58], as has been noted in animals possessing venom glands [57, 99, 100]. These bioactive proteins function in cytolytic, paralytic and hemolytic roles, thereby facilitating prey digestion [58, 64, 101]. Specifically we note the abundance of several enzyme groups primarily in either the tentacle or gastric cirri in A. alata that have been well studied in venomous animals [57, 99, 102, 103], namely astacin-like metalloproteinase and serine proteinase (and inhibitors), and more broadly cysteine-rich secretory proteins (CRISPs). Metalloproteinase and serine proteinase (and inhibitors) are a common component of the venom of animals with venom glands either activating toxins or acting as toxins themselves [18, 99, 100]. In particular, cysteine-rich secretory proteins identified in snake venoms are thought to inhibit smooth muscle contraction in bite victims [104]. Both metalloproteinase and CRISPs have previously been characterized in the tentacles of cubozoan [5, 105] and other cnidarians [19, 20, 61, 68, 106]. The abundance of multiple isozymes of astacin-like metalloproteinase and serine proteinase (and inhibitors) and CRISPs in the gastric cirri and tentacle of A. alata suggest a dual role in venom and digestion. Further studies are required to test this hypothesis given the broad involvement of these bioactive proteins in other biological processes [5, 103].

Vision and the phototransduction pathway

Cubozoans are the earliest diverging animal clade to have image-forming lens eyes, which are part of specialized sensory organs called rhopalia. By comparing a medusa body part bearing conspicuous eyes (the rhopalium) and planulae with eye spots (rhabdomeric photoreceptors) against the medusa samples that lack documented photoreceptors (gastric cirri, ovaries and tentacle), our aim was to profile the molecular components of the opsin-regulated phototransduction pathway and identify additional regions of putative extraocular sensory perception in A. alata.

We found that many genes with conserved roles in vision (opsins and crystallins) were abundant in the rhopalium. Although transcripts with putative roles in light-mediated phototransduction pathway were detected primarily in the rhopalium where eyes are present, their expression was broadly detected across the medusa samples, and in some cases in planulae. We expect this signal stems from the presence of additional photoreceptors (yet undescribed) throughout the body of this cubozoan. Together these findings are consistent with a vision-related role for opsins and crystallins in the lens-eye of the rhopalium, as well as a role in putative photoreceptors within non-rhopalium tissues and in planulae eye spots.

The animal phototransduction pathway is mediated by photopigments in photoreceptors consisting of two parts: a membrane protein (apoprotein) “opsin” and a chromophore “retinal” (vitamin A derivative) [107]. Opsins mediate light as phototypical G protein-coupled receptors in both visual and non-visual systems [108]. Currently more than 1000 types of opsin are known across Metazoa, with three subfamilies recognized in bilaterians: rhabdomeric (r-opsins), Go-coupled plus retinochrome retinal G protein-coupled receptor (Go/RGR) and ciliary (c-opsins) [107]. Studies characterizing opsins in cnidarians have raised the possibility that cnidarian opsins form a monophyletic clade referred to as “cnidops” that is sister to the c-opsins [26–28, 109, 110]. Other metazoan-wide analyses of opsins have categorized cnidarian opsins into three groups, A, B and C, in which each of these cnidarians opsin groups has been found, albeit with limited support, to be sister to each of the respective bilaterian opsin groups [30]. One study has revealed support for r-opsins in cnidarians [111], which is consistent with the identification of planulae eye spots as rhabdomeric photoreceptors [24, 25]. Although a consensus is lacking on the relationships between cnidarian opsins and other metazoan opsins, our study identified a number of transcripts with molecular characters corresponding to diverse metazoan opsins (rhodopsin family) in A. alata, adding to the known diversity of this gene family within cnidarians.

Our opsin gene tree only included the known cnidarian opsins and thus does not address the question of cnidarian opsin monophyly. However, our analysis recovered homologs within two of the three previously identified cnidarian opsin groups, namely group A and B [30]. Our analysis also recovered many opsin homologs within a large group that also contains cnidarian opsin group C, previously thought to only be present in anthozoans. The presence of this opsin group in both anthozoans and medusozoans suggests that it was present in the cnidarian ancestor.

We found several opsin genes in A. alata to be highly expressed in samples other than the rhopalium, and similar results have been reported for opsins in another cubozoan, T. cystophora [27, 28]. Based on these findings we hypothesize that cubozoans have opsin-mediated extraocular photoreception activity possibly related to phototaxis, circadian rhythm or light-mediated spawning, such as has been demonstrated in other animals, including anthozoans [29, 35, 76, 112]. These findings are also suggestive of extraocular photosensitivity [25, 108, 110, 113] that has a documented role in rhythmic behaviors and physiological processes in vertebrates and invertebrates, including nematocyst firing in cnidarians [28, 30, 32, 110, 111, 113]. Such suggested extraocular photoreceptor cells may also comprise anatomically dispersed light sensitive neurons, in addition to ciliary or rhabdomeric morphotypes, possibly functioning in dispersed photoreception, also called the “dermal light sense” (for a review see [114]). Future characterization of the absorbance spectra for different opsin types in cubozoans, and visualization of the precise locality of expression using in situ hybridization, will help elucidate their potential functions in different medusa body parts and planulae.

In this study the expression of some of the components of the retinal photoisomerization pathway in all samples including planulae suggests that A. alata metabolizes the universal chromophore retinal [35, 107]. However, transcripts for a putative blue-sensitive photoreceptor protein and circadian clock regulator cryptochrome homolog suggest an additional putative chromophore in A. alata that might function in non-rhopalium related blue-light mediated processes (e.g., phototaxis); such a function has previously been documented in other metazoans [35, 76]. Determining the precise chromophore utilized by A. alata must await future functional studies.

Crystallins are multifunctional proteins often related to stress or metabolic enzymes that serve as important lens components controlling optical properties [115]. The dual role crystallins play in eye lens as well as non-eye related tissues is known as “gene sharing” [80]. We found that all three types of J-crystallins previously reported in T. cystophora [115] were present in A. alata and that these were typically most highly expressed in the rhopalium, with J2 crystallins showing more variable expression across samples.

We also identified transcripts corresponding to the developmental transcription factors Six and eyes absent (Eya), representing the first homologs of these genes identified from cubozoans. Genes in the Six-Eya homolog complex have known functions in eye development, including during embryogenesis and regeneration, in both non-bilaterians and bilaterians [79, 84, 115]. Six-Eya complex genes have been shown to act downstream of Pax genes [84], and PaxB expression has been reported in both adult and larval eyes of T. cystophora, where it is inferred to promote J-crystallin expression [28]. We also identified transcripts corresponding to the developmental transcription factors Six and eyes absent (Eya), representing the first homologs of these genes identified from cubozoans. Genes in the Six-Eya homolog complex have known functions in eye development, including during embryogenesis and regeneration, in both non-bilaterians and bilaterians [79, 84, 115]. Six-Eya complex genes have been shown to act downstream of Pax genes [84], and PaxB expression has been reported in both adult and larval eyes of T. cystophora, where it is inferred to promote J-crystallin expression [28]. Conversely, in the scyphozoan Aurelia, development of simple eyes is mediated by Six-Eya complex genes independent of PaxB expression [116]. Although the Trinotate report for the filtered A. alata transcriptome did not contain any transcripts annotated as PaxB, we identified two transcripts (comp95018 and comp20156) annotated as other homeobox genes that appear to be putative PaxB homologs based on sequence identity (tBLASTx) with Nematostella vectensis PaxB mRNA. Whether eye development and Six-Eya expression in A. alata are dependent or independent of PaxB expression remain open questions. Future studies determining the spatial localization of gene expression during eye development in A. alata may be useful for further elucidating the gene regulatory networks functioning in eye development in cubozoan rhopalia and planulae eyes spots.

Sexual reproduction and embryogenesis

Cubozoan lifecycles alternate between an asexually reproducing sessile polyp stage and a sexually reproducing motile medusa stage. By profiling the transcripts from an adult body part (ovaries) abundant in developing oocytes, our aim was to characterize the molecular components of oogenesis and early embryogenesis in A. alata. Additionally, because we found that sperm are internalized and interact with newly ovulated eggs within the gastrovascular cavity of A. alata females, our ovaries tissue sample also provided the opportunity to identify genes that might be involved in fertilization.

We identified several apparent homologs of Vitellogenin and Apolipoprotein, which have documented roles in oogenesis and embryogenesis, [86–89] and found these to be most abundant in the ovaries of A. alata. Vitellogenin is an animal egg yolk protein that is synthesized in somatic cell lineages and subsequently incorporated into developing oocytes (by receptor mediated endocytosis), eventually serving as a nutrition source during embryogenesis [86–88]. In medusozoans little is known about the characteristics of Vitellogenins as they have only been documented as egg yolk proteins in two coral species [88, 116] and the model sea anemone Nematostella vectensis [87]. Vitellogenin proteins are expressed in both ovarian (or putative ovaries in anthozoans, e.g., [87]) and extra-ovarian somatic cells, consistent with their important roles in processing large lipoproteins in a broad range of complex biological processes among metazoans [86], including their distinct role as honey bee venom allergens [117]. Consistent with this, we found that in A. alata, apparent homologs of Vitellogenin-2 were expressed most highly in the ovaries, yet they and other Vitellogenins and Apolipoprotein-like homologs were detected in all medusa samples. In this study we also found a number of creatine kinase genes to be most abundant in the ovaries, but many were also detected (though at much lower expression levels) in all of our samples. Creatine kinases play an important role in oogenesis and early embryogenesis in mammals [89], having a broad enzymatic function in yielding ATP by catalyzing the reversible transfer of phosphate from creatine phosphate to ADP in cells with high activity (e.g., photoreceptors, primitive-type spermatozoa) [90, 91].

We did not recover any genes characteristic of meiosis in the ovaries sample of A. alata, suggesting that the tissue was composed exclusively of mature ova at the time of sampling. It is also possible that the expression levels of putative meiosis transcripts were too low to be detected by our analyses, given our conservative transcriptome analysis protocol (see Methods). However, few studies exist that characterize the molecular aspects of sexual reproduction in cnidarians [88, 116], limiting the number of potential gametogenic candidate genes targeted in this study. Future transcriptome and proteome profiling studies of the gonads of A. alata and other cubozoans during medusa maturation are needed to shed light on the molecular underpinnings of the processes controlling gametogenesis in cubozoans.

In this study we also detected the expression of genes with putative roles in sperm flagella activation, pro-acrosomal vesicles and sperm capacitation, with many of these being most abundant in the ovaries of A. alata. These morphological and biochemical changes to the sperm are necessary for the sperm to reach and fertilize an oocyte, and their occurrence has been documented within the female reproductive tract in many animals [92]. Sperm capacitation was previously thought to occur exclusively in mammals, but more recently it has been documented in several invertebrates [59]. Our study is the first to suggest that sperm capacitation might occur within the gastrovascular cavity (putative female reproductive tract) of a cnidarian. We note that although sperm storage structures have been reported in a single family of cubozoans (Tripedaliidae) [1, 118], we do not expect sperm storage to occur in A. alata. Morphological observations during the course of this study as well as previous studies in A. alata have identified no structure(s) with a putative role in sperm storage in either male or female medusae [2, 38, 119]. Future histological studies of A. alata medusae undergoing internal fertilization should elucidate the ultrastructure of the female reproductive tract and provide further insight into fertilization dynamics in this species.

Based on our observations in this and a prior study [2], monthly spermcasting aggregations of A. alata medusae consist entirely of males and females with mature gonad morphology. We therefore hypothesize that gonad development occurs offshore in response to environmental and molecular cues related to the lunar cycle that may instigate inshore migrations. During these monthly nearshore aggregations, which span three to four consecutive days, both sexes exhaust their entire gamete reserves in a process known as “controlled gonad rupture” [2, 120]. Male gonads completely disintegrate over the course of several hours, and females simultaneously ingest massive quantities of sperm for internal fertilization. The interaction of sperm and eggs witnessed in the gastrovascular cavity, followed by release of blastulae into the surrounding water by females within hours, along with the abundance of sperm and fertilization-related transcripts detected in the ovaries sample, corroborate previous observations [2] that fertilization occurs immediately following sperm ingestion and ovulation, adjacent to the ovaries within the gastrovascular cavity. Future molecular studies characterizing expression in sperm and eggs prior to and during fertilization will provide further insight into the dynamics of fertilization in cubozoans.

Specimen vouchers

Male and female A. alata medusae were collected in Bonaire, The Netherlands. All proper collection and export permits were obtained. Medusae were kept within a glass aquarium in filtered seawater for several hours. Males shed sperm into the water that was taken up by the female manubrium. Using light microscopy (1000x) to observe the ovaries, which are located within the gastrovascular cavity, confirmation was made of ovulation and sperm and egg interaction (i.e., putative fertilization) (Fig. 1c, d). No prey items were present within the stomach and associated gastric cirri or attached to the tentacle. A live female A. alata medusa undergoing internal fertilization was placed on ice, and using a sterile RNase-free disposable scalpel tissue samples were quickly excised from the gastric cirri, ovaries, tentacle, rhopalium (Fig. 1). A fifth sample consisting of thousands of swimming planulae that had developed from blastulae released from different females in the lab was also collected. All samples were placed in 2 ml cryovials and flash frozen with liquid nitrogen. Frozen samples were shipped via Cryoport to the Smithsonian Biorepository. Additionally, a single spawning A. alata female medusa was collected at Oil Slick Leap, Kralendijk on April 22, 2014, relaxed in 7.5 % Magnesium chloride, fixed and preserved in 8 % formalin, and deposited into the collection of the National Museum of Natural History, Washington, D.C. as a morphological voucher (USNM 1248604). No specific permissions were required from an ethics committee to conduct the research described herein as no humans or protected species were used.

Sequencing

The five frozen tissue samples (gastric cirri, ovaries, tentacle (and adjoining pedalium base), rhopalium, and planulae) were sequenced at the University of Kansas Medical Center—Genomics Core (KUMC), where total RNA (0.5 ug) was used for library preparation for each sample. Illumina HiSeq 2500 Sequencing System was used to generate FASTQ files, which were de-multiplexed into individual sequences for further downstream analysis.

Transcriptome assembly and post-assembly analyses

The 278 M paired end (100 bp) raw reads from five samples were analyzed on the Smithsonian Institution High Performance Cluster, SI/HPC, and filtered using the program TrimGalore! [121] with the adaptor trimming tool Cutadapt [122] and FastQC [123] (--quality 30 --phred33 --length 25) to remove Illumina lane and multiplex adaptors (overlapping by 1 bp). ALLPATHSLG error correction software [124] was used on the 265 M trimmed paired end reads (PAIRED_SEP option was set to 100), and unpaired reads following trimming to predict and correct sequencing errors (see [125]) and mitigate potential errors in transcriptome assemblies. All five samples (i.e., gastric cirri, ovaries, tentacle, rhopalium, and planulae) were pooled and assembled de novo into a reference transcriptome (FASTA format) for A. alata using Trinity (version trinityrnaseq_r20131110) [50, 51], with the following additional flags: --no_bowtie --normalize_reads--path_reinforcement_distance 75.

Differential expression estimates and analyses

RNA-Seq by Expectation Maximization (RSEM) was run on each of the five samples separately to estimate transcript abundance (read counts). A single matrix was generated corresponding to expression values for all samples as normalized Trimmed Mean of M-values (TMM) [126]. EdgeR was then used to identify differentially expressed genes in the counts matrix (--dispersion 0.1) [55]; followed by differential expression analysis to extract all genes most significantly expressed, i.e., with p-values < =0.005 and with at least a fourfold change of differential expression (--matrix iso_r123456.TMM.fpkm.matrix -P 1e-3 -C 2). This EdgeR step generated a single expression matrix of the results of all pairwise comparisons between the five samples. Further, hierarchical clustering generated a heatmap indicating clustering of similarly expressed genes (vertical axis) plotted by sample type (horizontal axis), while maintaining column order by sample (--order_columns_by_samples). This was done for all differentially expressed transcripts for all five samples (Additional file 4); just for the medusa samples (Fig. 3a); and for each of the three subsets of differentially expressed candidate genes (Figs. 6a, 7a and 8a). Color-coding on the vertical access of each heatmap indicates gene clusters with similar mean expression levels. Gene cluster patterns were further subdivided into 10 K-mean subclusters, which were visualized as subcluster profile plots (Figs. 3b-k, 6, 7 and 8b-k). In the absence of biological replicates in this study, the specific significance of fold-change expression levels of each of the differentially expressed genes was of limited value, and we therefore chose to not further filter transcripts based on additional statistical analyses. Instead, all differentially expressed genes were targeted as candidates for narrowing our search for genes of interest by sample type. Furthermore, redoing the hierarchical clustering analysis on just the three subsets of candidate genes (putative venom, vision and sex genes) allowed us to hone in on gene clusters that were relevant to transcriptome functional annotation and profiling of A. alata samples types.

Additional analyses

Venn diagrams were constructed using Venny [127]. Trinotate reports for each of the three sets of candidate genes investigated in this study (venom, vision and sex) were generated by filtering the original A. alata Trinotate report using this custom Python script: https://github.com/pbfrandsen/SI_scripts/blob/master/cheryl_trinotate.py.

Gene tree reconstruction

Amino acid sequences corresponding to predicted ORFs (TransDecoder), or translated nucleotide sequences, from the A. alata transcriptome were aligned using MUSCLE (default parameters with 5 iterations) against other cnidarian homologs from NCBI Genbank for the respective candidate genes of interest. ProtTest v. 3.2 was used to determine the most appropriate model of amino acid evolution (i.e., LG + G or WAG + G + F, or BLOSS62 I-G-F) for each alignment. Shimodaira-Hasegawa-like branch support indices [128] and are shown at each node of the ML topology. All ORF alignments (.nex files) predicted from TransDecoder, and the corresponding gene tree reconstructions (.tre files) are available at: https://figshare.com/articles/Supplemental_Information_for_A_new_transcriptome_and_transcriptome_profiling_of_adult_and_larval_tissue_in_the_box_jellyfish_Alatina_alata_an_emerging_modelfor_studying_venom_vision_and_sex/3471425.