Sequences of Circadian Clock Proteins in the Nudibranch Molluscs Hermissenda crassicornis, Melibe leonina, and Tritonia diomedea.

While much is known about the genes and proteins that make up the circadian clocks in vertebrates and several arthropod species, much less is known about the clock genes in many other invertebrates, including nudibranchs. The goal of this project was to identify the RNA and protein products of putative clock genes in the central nervous system of three nudibranchs, Hermissenda crassicornis, Melibe leonina, and Tritonia diomedea. Using previously published transcriptomes (Hermissenda and Tritonia) and a new transcriptome (Melibe), we identified nudibranch orthologs for the products of five canonical clock genes: brain and muscle aryl hydrocarbon receptor nuclear translocator like protein 1, circadian locomotor output cycles kaput, non-photoreceptive cryptochrome, period, and timeless. Additionally, orthologous sequences for the products of five related genes-aryl hydrocarbon receptor nuclear translocator like, photoreceptive cryptochrome, cryptochrome DASH, 6-4 photolyase, and timeout-were determined. Phylogenetic analyses confirmed that the nudibranch proteins were most closely related to known orthologs in related invertebrates, such as oysters and annelids. In general, the nudibranch clock proteins shared greater sequence similarity with Mus musculus orthologs than Drosophila melanogaster orthologs, which is consistent with the closer phylogenetic relationships recovered between lophotrochozoan and vertebrate orthologs. The suite of clock-related genes in nudibranchs includes both photoreceptive and non-photoreceptive cryptochromes, as well as timeout and possibly timeless. Therefore, the nudibranch clock may resemble the one exhibited in mammals, or possibly even in non-drosopholid insects and oysters. The latter would be evidence supporting this as the ancestral clock for bilaterians.

Introduction

Most organisms express daily rhythms of activity and physiology that persist in constant conditions and thus are referred to as circadian rhythms. These circadian rhythms are produced by molecular clocks that exhibit cyclical production of various intracellular proteins as the result of transcription-translation feedback loops. The first clock protein to be discovered was PERIOD, initially identified in Drosophila melanogaster by Konopka and Benzer (1971). Since then, homologs of PERIOD have been found in numerous other bilaterian organisms, including humans (Tei et al., 1997). Additional circadian proteins that participate in the molecular clock have been identified and sequenced in a wide array of species, ranging from cyanobacteria to mammals (reviewed in Dunlap, 1999; Zhang and Kay, 2010). Although the individual proteins can vary between widely disparate phylogenetic groups, the presence of negative and positive feedback elements regulating the transcription and translation of each other appears to be a ubiquitous component of circadian pacemakers in all organisms studied to date.

In animals, the molecular mechanisms of the circadian clock are best understood in Drosophila and mammals. While the Drosophila clock is quite derived and different from that of some other insects (Rubin et al., 2006; Zhu et al., 2008; Ingram et el., 2012), we focus on it here because it is one of the best-characterized circadian systems. Briefly, in Drosophila, the canonical core clock proteins consist of circadian locomotor output cycles kaput (CLOCK) and CYCLE, which work together to drive the production of PERIOD and TIMELESS. These latter two proteins feed back on their own transcription, forming a negative feedback loop with a circadian periodicity. A blue-light photoreceptive CRYPTOCHROME (PCRY) is important in providing light input to the clock by triggering the degradation of TIMELESS (Hardin, 2005; Rosato et al., 2006). In mammals, CLOCK and brain and muscle aryl hydrocarbon receptor nuclear translocator-like protein 1 (BMAL1; a homolog of CYCLE) heterodimerize and stimulate transcription and translation of PERIOD and a vertebrate, nonphotoreceptive CRYPTOCHROME (NPCRY). PERIOD and NPCRY then serve as the negative inhibitors of the clock (Reppert and Weaver, 2002; Lowrey and Takahashi, 2004). NPCRY is closely related to a protostome NPCRY, found in many invertebrates, but not present in Drosophila (Lin and Todo, 2005; Özturk, 2016; Michael et al., 2017). A TIMELESS protein was initially identified in mammals (Koike et al., 1998; Sangoram et al., 1998; Takumi et al., 1999), but this is now believed to actually be an ortholog of TIMEOUT (Li et al., 2016), which is a paralog of TIMELESS.

The molecular circadian clocks in mammals and fruit flies represent two major lineages of animals: deuterostome chordates and protostome ecdysozoans, respectively. As described above, these molecular clocks are quite similar operationally, but they differ in the proteins used for negative feedback and how light input is transmitted to the clock. While circadian genes have been identified in some species in the other major protostome lineage, Lophotrochozoa (Zantke et al., 2013; Bao et al., 2017; Perrigault and Tran, 2017; Schnytzer et al., 2018), we do not yet fully understand the molecular mechanisms underlying clocks in lophotrochozoans as a group. Further elucidating the molecular basis of circadian clocks in this clade may shed light on the evolution of circadian clocks.

Circadian rhythms in gastropod molluscs, which are a lophotrochozoan clade, were first discovered 50 years ago in Aplysia californica (Kupfermann, 1968; Jacklet, 1972; Block and Lickey, 1973; Lickey et al., 1977). Subsequently, circadian rhythms of locomotion, oxygen consumption, and ocular electrical activity have been identified in several other gastropods, including Bulina tropicas (Chaudhry and Morgan, 1983), Bulla gouldiana (Block and Davenport, 1982), Bursatella leachi (Block and Roberts, 1981), Helisoma trivolvis (Kavaliers, 1981), Helix aspersa (Bailey, 1981; Blanc, 1993), Hydrobia ulvae (Barnes, 1986), Limax maximus (Sokolove et al., 1977), Littorina irrorata (Shirley and Findley, 1978), Melanerita atramentosa (Zann, 1973), Melanoides tuberculata (Beeston and Morgan, 1979), and Melibe leonina (Newcomb et al., 2014). However, despite strong interest in gastropod circadian rhythms, and the advantages of this group of animals for investigating the neuronal bases of behaviors, there has been very little progress in identifying circadian genes in gastropods, with the exception of the transcript for period in Bulla gouldiana (Constance et al., 2002), the basic helix-loop-helix (bHLH)-containing proteins BMAL1 and CLOCK in Biomphalaria glabrata, Lottia gigantea, and Patella vulgate (Bao et al., 2017), and some automated annotations on GenBank at the National Center for Biotechnology Information (NCBI, Bethesda, MD). Furthermore, after submission of this paper, Schnytzer et al. (2018) published a study reporting on the transcript sequences for most of the core circadian clock genes in the limpet Cellana rota.

The goal of this study was to identify the core circadian gene products in the central nervous systems (CNS) of three nudibranchs that are common neurophysiological model systems: Hermissenda crassicornis (Eschscholtz, 1831), Melibe leonina (Gould, 1852), and Tritonia diomedea Bergh, 1894 (synonym = Tritonia tetraquetra [Pallas, 1788]), using bioinformatic analyses of recently assembled transcriptomes (Senatore et al., 2015; Tamvacakis et al., 2015). Gene products were identified in all three species for a number of canonical clock genes, as well as some associated sequences, and compared to sequences in clock proteins in other species.

Materials and Methods

Transcriptomes

Previously published CNS transcriptomes were used for Hermissenda and Tritonia (Senatore et al., 2015; Tamvacakis et al., 2015), while the Melibe transcriptome was developed using 15 specimens collected by Monterey Abalone Company (Monterey, CA). Total RNA was extracted from the CNS (the fused cerebropleural and pedal ganglia, plus the buccal ganglia), using an RNeasy Plus Universal Midi Kit (QIAGEN, Hilden, Germany). RNA integrity was confirmed by visualization of electrophoresed RNA on an ethidium bromide-stained agarose gel, as well as with a bioanalyzer, using an Agilent RNA 6000 Pico Kit (Agilent Technologies, Waldbronn, Germany). RNA was then diluted to ~100 ng μl−1 in Trisethylenediaminetetraacetic acid (EDTA) buffer (1 mmol l−1 EDTA, 10 mmol l−1 Tris, pH 7.5) and shipped frozen on dry ice to Beckman Coulter Genomics (Danvers, MA) for cDNA synthesis and Illumina sequencing. At Beckman Genomics, mRNA was isolated via polyA capture, and cDNA was synthesized and amplified using proprietary primers. Paired-end sequencing of cDNA was performed on an Illumina HiSeq2500 platform (2 × 100 bp, Illumina, San Diego, CA), followed by the removal of adaptor sequences and low-quality reads. After delivery of the data, we further trimmed paired reads with the program Sickle (Joshi and Fass, 2011) to remove 3′ and 5′ ends as well as entire reads shorter than 20 bp or with quality scores less than 20 (i.e., default parameters). A de novo transcriptome shotgun assembly (TSA) was then developed using Trinity (release 2013-08-14; Grabherr et al., 2011; Haas et al., 2013; Senatore et al., 2015; Tamvacakis et al., 2015).

Identification of homologous clock sequences

The creation of a bioinformatic pipeline (Pankey, 2018) significantly increased the efficiency with which contiguous sequences from the nudibranch transcriptomes could be identified as potential gene orthologs (Fig. 1). More importantly, this pipeline incorporated phylogenetic analyses to discriminate sequences with orthology to specific clock genes from paralogous proteins that do not serve a circadian clock function. After TSA assembly, candidate coding regions from each nudibranch transcriptome were identified via the in silico translation of the Trinity data sets using TransDecoder (Haas et al., 2013). The translated transcriptome for each nudibranch species could then be imported into the bioinformatic pipeline for the automated discovery of circadian clock gene orthologs.

Figure 1.

Schematic flow diagram of the steps taken to sequence and compare clock proteins. The creation of a bioinformatic pipeline significantly increased the efficiency with which contiguous sequences from the nudibranch transcriptomes could be identified as potential gene orthologs; more importantly, it accurately determined whether they were most closely associated with the products of specific clock genes or more similar to structurally related proteins that do not serve a circadian clock function. To begin, valid reference clock protein sequences (what we termed “BAIT” sequences) were compiled and used to fish out contiguous sequences from the nudibranch transcriptome assemblies. After aligning all sequences, contiguous sequences were then mapped on a phylogenetic tree with the BAIT sequences as well as closely matched protein sequences from automated BLAST searches. This enabled us to identify putative clock protein sequences from each nudibranch transcriptome. ML, maximum likelihood.

Next, various databases (e.g., GenBank, UniProt) were manually queried to locate and compile a valid reference set for each of the following protein sequences: CLOCK, PERIOD, BMAL1, CYCLE, aryl hydrocarbon receptor nuclear translocator-like (ARNTL), PCRY, NPCRY, 6-4 PHOTOLYASE (PHR), CRYPTOCHROME DASH (CRY DASH), TIMELESS, and TIMEOUT. Bait sequences are indicated in blue font in Supplementary Figures 1–5, and accession numbers for these proteins can be found in Supplementary Table 1 (Supplementary Figs. 1–9 and Supplementary Tables 1–3 are available online). These protein sequences were used as queries in BLASTp (NCBIBLAST+) searches of all curated animal genomes on NCBI, as well as the translated nudibranch transcriptomes (Altschul et al., 1990; Fig. 1). Similar to SmartBlast (NCBI), the parallel execution of these tasks efficiently identified protein sequences for targeted circadian clock genes from the nudibranch transcriptomes as well as a diverse array of metazoan taxa that aided with the subsequent rooting of phylogenetic trees and the locating of paralogous gene clades.

Last, the full-length protein sequences for each circadian gene recovered from genome searches, along with the initial query sequences, were aligned using MUSCLE, version 3.8.31 (Fig. 1; Edgar, 2004). Each nudibranch predicted peptide fragment was then iteratively aligned to the full-length alignment, using the “-profile” setting in MUSCLE. With each resulting alignment, RAxML, version 7.3.0, was used to infer an initial gene tree by maximum likelihood (ML) under PROTGAMMAWAG and to assess support for bipartitions via 1000 bootstrap replicates (see Supplementary Figs. 1–5, available online; Yang, 1996; Whelan and Goldman, 2001; Stamatakis, 2014)

Alignment, functional domains, and phylogenetic analyses of homologous clock sequences

The alignments and phylogenetic analyses described above were used as an initial step to identify candidate clock sequences in the three nudibranch transcriptomes. We then used additional alignments and phylogenetic analyses, as described below, for comparing conserved functional domains. To begin, orthologous sequences for each of the 10 circadian clock proteins targeted by this investigation were compiled from different lophotrochozoan, ecdysozoan, and chordate taxa into a reference set (see Supplementary Table 1, available online). Two multiple sequence alignments (MSAs) were sequentially performed on nudibranch protein sequences with orthology to specific clock genes, using Clustal Omega (Sievers et al., 2011). This two-step process prepared the collection of protein sequences associated with each clock gene for phylogenetic analysis by ML.

The first MSA included a circadian clock protein sequence (e.g., CLOCK) mined from each nudibranch transcriptome—when available—and the relevant orthologs from a small subset of taxa. This facilitated the efficient identification of conserved functional protein domains among lophotrochozoan, ecdysozoan, and/or chordate taxa using Normal SMART (SMART Technologies, Calgary, Alberta, Canada) (Schulz et al., 1998; Letunic et al., 2017). Once the conserved functional domains in each circadian clock protein had been identified and mapped—effectively creating a guide for their rapid identification in orthologous protein sequences from other taxa—sequences flanking these conserved domains were trimmed, and a second MSA was performed.

The refined protein sequences for each clock gene were then subjected to phylogenetic analysis in MEGA7 (Kumar et al., 2016), using ML with 1000 bootstrap replications. The Whelan and Goldman (WAG) model of amino acid substitution, with gamma-distributed and invariant (G+I) sites with five discrete gamma categories, was used for each analysis. The tree inference options were set to nearest-neighbor interchange, with an automatically generated (neighbor-joining [NJ] and BioNJ) initial tree. The tree with the highest log likelihood was considered the best for each phylogenetic analysis.

Results

Sequencing and TSA metrics were previously provided for Hermissenda (Tamvacakis et al., 2015) and Tritonia (Senatore et al., 2015). For Melibe, sequencing resulted in 123,087,512 paired-end reads, 94.5% of which had a Phred score greater than 30. De novo assembly with Trinity yielded 167,841 transcripts with an N50 of 1528 bp; these partitioned into 93,882 nonredundant gene clusters (unigenes). In silico translation of the transcriptome data set with TransDecoder yielded 39,641 predicted protein sequences. Paired-end reads and the TSA were uploaded to NCBI as a BioProject (accession no. PRJNA420367).

Five core clock orthologs were identified in all three nudibranch species (BMAL1, CLOCK, PERIOD, TIMELESS, and NPCRY; Table 1). In addition, we found sequences in the three species for ARNTL, TIMEOUT, and three other CRYPTOCHROMES (PCRY, PHR, and CRY DASH; Table 1). BLASTp analyses provided evidence supporting the similarity of these sequences to other molluscan clock-related proteins and a closer affinity to mammalian orthologs than those of Drosophila (Supplementary Tables 2, 3, available online). These results were further supported by the more rigorous phylogenetic analyses described below and in Supplementary Figures 1–5, available online.

Table 1

Identified core circadian gene transcripts for the three nudibranchs

Designation	Hermissenda crassicornis	Melibe leonina	Tritonia diomedea
clock	MG162587, C, 650	MG189943, C, 651	MG253827, P, 504
bmal1		MG282902, C, 698	MG282903, P, 126
arntl	MG282904, P, 252	MG282905, C, 631	MG282906, C, 633
period	MG282907, C, 1615	MG282908, C, 1491	MG427049, C, 1565
npcry	MG437149, C, 561	MG437150, C, 556	MG437151, C, 559
pcry	MG437152, C, 547	MG470828, C, 543	MG470829, C, 542
phr	MG470830, P, 162	MG470831, P, 236	MG516799, C, 530
cry dash	MG516800, C, 527	MG516801, C, 561	MG516802, C, 550
timeless	MG516803, P, 608	MG516804, P, 696	MG516805, P, 632
timeout	MG549825, C, 673	MG549826, P, 381	MG549827, P, *

In the species columns, each data cell contains the accession number, the designation of a complete (C) or partial (P) sequence, and the number of amino acids in the translated protein. Sequences were determined as complete or partial on the basis of alignment with orthologous sequences and the presence of start and stop codons. Protein translations were done with the standard settings in the ExPASy Translate Tool (Artimo et al., 2012). An asterisk indicates multiple noncontiguous nucleotide sequences, which thus could not be translated.

clock

In all three nudibranchs, we identified a single RNA transcript for the gene clock (Table 1). The length of the predicted amino acid sequence was similar for Hermissenda (650 amino acids) and Melibe (651), whereas in Tritonia, CLOCK was >100 amino acids shorter (504), indicating that it was likely a partial sequence. A BLAST search indicated that Melibe CLOCK had the highest BLAST score with a predicted CLOCK sequence in the sea hare Aplysia californica (63% identity; Supplementary Table 2, available online). The nudibranch CLOCK proteins contained the conserved regions present in orthologs of other species, including a helix-loop-helix (HLH) domain, two Per-Arnt-Sim (PAS) domains, and aPAS-associated C-terminal (PAC) motif (see Fig. 2 for example of the alignment; alignments for subsequent circadian transcripts are provided in Supplementary Figs. 6–9, available online). Similar to BLASTp comparisons (Supplementary Tables 2, 3), phylogenetic analysis indicated that the nudibranch CLOCK proteins evolved from an ancestral lophotrochozoan ortholog with a closer relationship to chordate CLOCK proteins than arthropod (ecdysozoan) orthologs, suggesting a dynamic evolutionary history of duplication and losses of this homolog (Fig. 3; Supplementary Fig. 1, available online).

Figure 2.

Alignment of conserved regions for CLOCK. The nudibranch CLOCK proteins contained conserved domains present in orthologs for Drosophila melanogaster and Mus musculus: blue indicates the helix-loop-helix domains, light green the Per-Arnt-Sim (PAS) domains, and red the PAS-associated C-terminal domains. Symbols immediately below the aligned sequences indicate functional similarity of residues among the CLOCK protein orthologs from the five different taxa. An asterisk indicates identical residue, a colon indicates residues with strongly similar properties, and a period indicates residues with weakly similar properties.

Figure 3.

Maximum likelihood analysis indicating the evolutionary relationships between CLOCK orthologs. The nudibranch CLOCK proteins grouped together with orthologs in other molluscs and shared a more recent common ancestor with chordate CLOCK proteins than arthropod orthologs. The tree was rooted on CLOCK from the cnidarian Nematostella vectensis. Numbers at branch points indicate bootstrap values for maximum likelihood analysis with 1000 replicates.

bmal1

A single RNA transcript for the bmal1 gene was identified in Melibe, which encoded a protein of 698 amino acids (Table 1). Only a partial bmal1 transcript could be determined in Tritonia (encoding 126 amino acids), and the bioinformatics pipeline did not identify a bmal1 sequence in Hermissenda. The predicted Melibe BMAL1 protein contained the conserved regions present in orthologs of other species, including the HLH, PAS, and PAC domains (Supplementary Fig. 6, available online). A BLAST search indicated that Melibe BMAL1 was most similar to a predicted ARNTL1 protein in the freshwater snail Biomphalaria glabrata (69% identity; Supplementary Table 2, available online).

In addition to bmal1, complete RNA transcripts of the related arntl gene were identified in all three species, and the Melibe ARNTL protein was used to root the phylogenetic analysis of BMAL1 (Fig. 4). Results of this analysis suggest that Melibe BMAL1 shares a sister grouping with the other molluscan BMAL1 (Crassostrea gigas) and that these arose from a homolog in a lophotrochozoan ancestor. As with CLOCK and BLASTp analyses (Supplementary Table 3, available online), the Melibe BMAL1 shared greater sequence similarity with chordate BMAL1 proteins than crustacean BMAL1 proteins or the related CYCLE found in insects.

Figure 4.

Evolutionary relationships of BMAL1/CYCLE orthologs, based on maximum likelihood analysis. In this study, a complete sequence for BMAL1 was identified only in Melibe leonina. The putative Melibe BMAL1 was closely related to orthologs in other lophotrochozoans (e.g., molluscs and annelids) and, similar to CLOCK, shared a more recent common ancestor with chordate BMAL1 proteins than arthropod orthologs. The tree was rooted on the newly identified ARNTL protein in Melibe.

period

A single transcript for the period gene was identified for all three nudibranchs (Table 1). Melibe PERIOD had the highest BLAST score to a predicted, but uncharacterized, protein in Aplysia, based on a BLASTp search (53% identity; Supplementary Table 2, available online). All three nudibranch PERIOD proteins contained the conserved HLH, PAS, and PAC domains present in orthologs of other species (Supplementary Fig. 7, available online). Phylogenetic analysis of these conserved regions suggested that, as with CLOCK and BMAL1, the nudibranch PERIOD proteins evolved from a lophotrochozoan ancestor (Fig. 5). In contrast to BLASTp comparisons that suggested that Melibe PERIOD was more similar to Mus than Drosophila (Supplementary Table 3, available online), phylogenetic analysis suggested that molluscan PERIOD proteins shared a common ancestor more recently with protostomes other than chordates (Fig. 5). However, low bootstrap support at deep nodes, together with the absence of non-bilaterian PERIOD orthologs (Supplementary Fig. 3, available online) to serve as an outgroup, hinder efforts to resolve the evolutionary history of this key player.

Figure 5.

Phylogenetic tree indicating evolutionary relationships of PERIOD orthologs across taxa. The nudibranch PERIOD proteins were most closely related to other molluscs but, unlike CLOCK and BMAL1, shared more recent common ancestry with arthropods than chordates, though with low bootstrap support for this relationship. The tree was rooted on NEURONAL NPAS DOMAIN PROTEIN (NPAS) from the oyster Crassostrea gigas (Perrigault and Tran, 2017).

Cryptochrome

We identified four distinct cryptochrome transcripts in all of the nudibranchs: npcry, pcry, phr, and cry dash (Table 1; Supplementary Fig. 4, available online). All of the translated protein sequences were complete, with the exception of the Hermissenda and Tritonia CRY DASH proteins. BLASTp searches (Supplementary Table 2, available online) indicated that Melibe NPCRY was most similar to a predicted cryptochrome-1-like protein in Aplysia californica (85% identity). Melibe PCRY was also most similar to a predicted cryptochrome-1-like protein, but in Biomphalaria glabrata (69% identity). NPCRY is a core clock protein in mammals, and possibly many invertebrates as well, so we compared conserved regions of the nudibranch NPCRY proteins with those of the bivalve Crassostrea (Supplementary Fig. 8, available online). The nudibranch NPCRY proteins all contained the conserved DNA photolyase and flavin adenine dinucleotide (FAD) binding 7 domains, as indicated by analysis in SMART software (SMART Technologies). PCRY is a core clock protein in Drosophila and also contains the same conserved DNA photolyase and FAD binding 7 domains as NPCRY, so our ML analysis included protein sequences from both of these cryptochrome families (Fig. 6). The nudibranch NPCRY and PCRY sequences were clustered together in their respective clades and grouped with other molluscan proteins. Similar to findings in other studies, the lophotrochozoan NPCRY proteins shared a common ancestor with vertebrate NPCRY proteins. PCRY formed a clade separate from NPCRY, and molluscan PCRY proteins were homologous to insect PCRY proteins.

Figure 6.

Phylogenetic tree of NPCRY and PCRY orthologs, based on maximum likelihood analysis of conserved regions. The nudibranch NPCRY proteins were similar to other molluscs and shared more recent ancestry with orthologs in chordates than arthropods. Nudibranch PCRY proteins were also grouped with molluscs. Chordates do not have orthologs for PCRY. The tree was rooted on a CRY protein from the sponge Amphimedon queenslandica (Rivera et al., 2012).

timeless and timeout

We identified a complete RNA transcript for the timeout gene in Hermissenda, but only partial sequences for timeout in the other two species, as well as only partial transcripts for the timeless gene in all three nudibranchs (Table 1). Hermissenda TIMEOUT contained the conserved TIMELESS domain present in orthologs in other species (Supplementary Fig. 9, available online). Phylogenetic analysis with the Hermissenda TIMEOUT suggests that it shares homology with nematode and chordate orthologs (Fig. 7). Although there are no TIMEOUT sequences yet reported in other lophotrochozoans, we have recovered a clade of putative TIMEOUT orthologs that includes annelids and molluscs (Supplementary Fig. 5, available online). TIMELESS and TIMEOUT proteins fell out into respective clades, with the exception of Drosophila TIMEOUT, which was in the midst of the TIMELESS clade. However, this and numerous other branches had low bootstrap support (<50%).

Figure 7.

Evolutionary relationships between TIMELESS and TIMEOUT orthologs. A complete TIMEOUT sequence was identified in Hermissenda crassicornis, while only partial sequences of TIMEOUT were identified in the other two nudibranchs. Similarly, only partial sequences for TIMELESS were identified in all three nudibranchs. Therefore, only the H. crassicornis TIMEOUT was included in this maximum likelihood analysis. The putative H. crassicornis TIMEOUT shared most recent ancestry with orthologs in Caenorhabditis elegans and chordates. TIMEOUT and TIMELESS largely grouped together, with the exception of arthropod TIMEOUT proteins (see Drosophila melanogaster, Danaus plexippus, and Tribolium castaneum). However, bootstrap values for more ancestral branch points were below the typical 50%–75% threshold for high confidence. The tree was rooted on TIMEOUT in the cnidarian Nematostella vectensis.

Discussion

We report the sequences of putative core clock gene products in three nudibranchs: Hermissenda crassicornis, Melibe leonina, and Tritonia diomedea. The majority of the predicted proteins have higher percent similarity with orthologs in Mus than those in Drosophila (Supplementary Table 3, available online). This similarity to chordate sequences was supported by most phylogenetic analyses, with the inclusion of numerous orthologs in other species (Figs. 3–7). The group of putative clock proteins identified in the nudibranchs suggests that the gastropod clocks are likely to function more like the mammalian clock than the Drosophila clock.

Bioinformatics pipeline

The nudibranch transcripts were identified using an automated bioinformatics pipeline that used manually determined bait sequences (usually from canonical species such as Drosophila and Mus) to fish out contiguous sequences from the transcriptome assemblies. These were then mapped on a phylogenetic tree with the bait sequences, as well as additional closely matched protein sequences from automated BLAST searches. Compared to simple BLAST-style searches of a transcriptome, this pipeline significantly facilitated the ease with which contiguous sequences could be identified as potential orthologs, and, more importantly, it accurately determined whether they were most closely associated with the products of specific clock genes or whether they were more similar to paralogs that likely do not serve a circadian clock function. For example, in the early stages of our transcriptome searches, we identified a gene product in Melibe using BLAST-style queries with BMAL1 and CYCLE bait sequences. However, it was not until we used the automated pipeline that we realized that the Melibe sequence that we had identified as a putative BMAL1 or CYCLE ortholog actually fell out with non-circadian ARNTL proteins. We were then further able to identify a transcript that more closely matched BMAL1 from related species (Fig. 4). The pipeline is a LINUX-based sequence of bioinformatics modules, using freely available online databases and the transcriptome of a selected organism. We have made this pipeline available for download (Pankey, 2018).

Nudibranch clock sequences

In mammals, CLOCK and BMAL1 work together to drive the transcription of other clock genes, including period and npcry (Gekakis et al., 1998), whereas in flies, CLOCK works with CYCLE, a homolog of mammalian BMAL1, to perform a similar function (Rutila et al., 1998). Even though a complete sequence for BMAL1 could only be found in Melibe, a partial sequence was identified in Tritonia. Bioinformatics analysis indicated that these proteins shared common ancestry with other BMAL1 proteins, as opposed to the homologous insect CYCLE (Fig. 4). This mirrors the situation in other lophotrochozoans investigated to date (Arendt et al., 2004; Zantke et al., 2013; Bao et al., 2017; Perrigault and Tran, 2017), as well as in crustaceans (Zhang et al., 2013) and some insects (Rubin et al., 2006).

Another key difference between the mammalian clock and the dipteran clock is the molecular partner for PERIOD. In Drosophila, PERIOD heterodimerizes with TIMELESS, and the pair act as a brake on their own transcription, as well as that of other clock-associated genes (Gekakis et al., 1995). Although it is hypothesized that TIMELESS and its paralog TIMEOUT were both present early in animals, TIMELESS has been lost in lineages leading to Caenorhabditis elegans and deuterostomes (Li et al., 2016). Thus, chordate proteins initially named TIMELESS have now proven to actually be orthologs of TIMEOUT. In regard to the mammalian circadian clock, PERIOD works with NPCRY, instead of TIMELESS, to perform a similar negative feedback function (Kume et al., 1999). Drosophila has lost the gene npcry, so TIMELESS may have been co-opted to serve a similar function. All three nudibranchs in this study had both TIMELESS and TIMEOUT, as well as NPCRY. While there is no clear evidence for a role for TIMEOUT in circadian clocks, it remains to be determined whether TIMELESS and/or NPCRY are involved in the negative feedback loop of the nudibranch circadian clock.

Animal cryptochromes are thought to have evolved from ancestral photolyases associated with light-dependent DNA repair (Mei and Dvornyk, 2015) and now include regulators of circadian rhythms (Özturk et al., 2007; Michael et al., 2017). Animal cryptochromes can be photoreceptive (PCRY), providing direct light input into the circadian clock, or non-photoreceptive (NPCRY), possibly acting as transcriptional repressors of clock genes. In this study, in addition to all three nudibranchs expressing related DNA repair genes phr and cry dash, they also exhibited gene products for both pcry and npcry. The presence of both pcry and npcry is similar to many invertebrates, although Drosophila has only pcry (Emery et al., 1998; Stanewsky et al., 1998); and vertebrates (Hsu et al., 1996; Todo et al., 1996) and some insects (Rubin et al., 2006; Ingram et al., 2012) have only npcry. It has been hypothesized that in animals that have both pcry and npcry, both of these genes may be important in circadian clock function, with PCRY acting as a blue-light photoreceptor to provide light input to the clock and NPCRY serving a transcriptional repressor function for the clock (Yuan et al., 2007). Further research will be necessary to determine whether both of these cryptochromes are involved with the circadian clock of nudibranchs.

While many animals appear to have single orthologs of the canonical core clock genes seen in Drosophila and Mus, there have also been instances during animal evolution where these clock genes have duplicated one or more times. For example, mammals have three period genes (period1, period2, and period3) (Shearman et al., 1997; Zylka et al., 1998), although only period1 and period2 play a significant role in the central clock mechanism (Bae et al., 2001); and they have two paralogs for npcry (referred to as cryptochrome1 and cryptochrome2) (Hsu et al., 1996). In invertebrates, the horseshoe crab Limulus polyphemus has two paralogs for cycle and three gene copies for period (Chesmore et al., 2016); and the silk moth Antheraea pernyi has two paralogs for period, located on different sex chromosomes (Gotter et al., 1999). In our study, we found only single copies of transcripts for all of the canonical clock genes. Therefore, we did not find any evidence of clock gene duplication, such as occurred in ecdysozoans (arthropods) and deuterostomes (mammals).

Evolution of circadian clocks

The two best-studied clock gene networks in the animal kingdom are those in Drosophila and mammals. Therefore, we compared the Melibe proteins with orthologs in Drosophila and Mus (Supplementary Table 3, available online). In almost all cases, the Melibe proteins had a higher similarity with those found in the mouse. The lone exception was TIMELESS, but that is likely due to the fact that TIMELESS has been lost in mammals, and the timeless gene in Mus is actually timeout (Li et al., 2016). Based on evolutionary relationships, this similarity of Melibe clock sequences to those of mammals may be due to recent divergence and/or rapid gene evolution of Drosophila. While the overall feedback loop of circadian clocks in Drosophila and mammals is similar, and there are numerous homologous genes in these clocks, there are distinct differences as well. For example, TIMELESS is the transcriptional repressor partner for PERIOD in Drosophila (Gekakis et al., 1995), whereas NPCRY serves that function in mammals (Kume et al., 1999). In addition, PCRY provides light input to the clock in Drosophila (Emery et al., 1998; Stanewsky et al., 1998), whereas MELANOPSIN plays this role in mammals (Hattar et al., 2002; Panda et al., 2002; Ruby et al., 2002). These differences raise the question of what the ancestral clock was like before the divergence of protostomes and deuterostomes. Lophotrochozoans, such as molluscs, represent an animal clade that can potentially shed some light on this question. Recent research on the oyster Crassostrea gigas found that the bivalve circadian clock may be intermediate to that seen in mammals and Drosophila, with both TIMELESS and NPCRY possibly acting as transcriptional repressors and PCRY providing light input to the clock (Perrigault and Tran, 2017). This is actually similar to the proposed clock mechanism in some non-drosopholid insects, such as butterflies (Zhu et al., 2008), and the chelicerate Limulus (horseshoe crab) (Chesmore et al., 2016). Considering that the nudibranchs in this study expressed a suite of clock genes most similar to butterflies, horseshoe crabs, and oysters, this potentially increases the number of clades with this type of circadian clock.

Going further back in evolutionary time, cnidarians also express orthologs of clock, bmal1 and cycle, and both pcry and npcry (reviewed in Reitzel et al., 2013). However, to date there is no evidence of period or timeless homologs existing in cnidarians. Therefore, the metazoan ancestral clock may resemble that present in non-drosopholid insects, horseshoe crabs, oysters, and nudibranchs, with the exception that NPCRY may have been the sole transcriptional repressor of clock and bmal1 and cycle. Additional studies in other lophotrochozoans and cnidarians, as well as other phylogenetically informative clades near the base of Metazoa (e.g., placozoans), should shed additional light on the evolution of these animal circadian clocks.