Identification and temporal expression of putative circadian clock transcripts in the amphipod crustacean Talitrus saltator

Animal husbandry, tissue sampling and RNA extraction

Animals were caught by hand from Ynyslas beach, Wales, UK at night during their active period and returned to glass tanks containing damp sand at ambient temperature (17 °C) and a 12:12 LD regime. Lighting at dawn and dusk was ramped/dimmed over 30 min to approximate natural conditions. Animals were fed fish food flakes ad libitum. After seven days acclimation, the rhythmic emergence and locomotor activity of a sub-sample of 60 animals was established. T. saltator has been shown to express more robust and less variable locomotor rhythms in small groups (Bregazzi & Naylor, 1972). Therefore, animals were housed in groups of five in a glass tank containing 10 cm-deep damp sand and compartmentalized with Plexiglas dividers. Across each compartment, infrared beams were passed via bespoke recording apparatus fabricated by Trikinetics (Waltham, MA, USA). All activity, registered as interruptions to infrared beams was recorded via proprietary software on a stand-alone PC. Activity data were recorded in one minute bins and was analysed and plotted using ClockLab software (Actimetrics, Wilmette, IL, USA) run via Matlab^® v6.2. T. saltator brains were rapidly dissected from animals drawn from the population shown to be behaviourally rhythmic (see above). Dissections were done in ice-chilled DEPC treated physiological saline. Ten brains were pooled to form each replicate at 3-hour intervals and snap frozen in liquid nitrogen before storing at −80 °C until use.

Total RNA was extracted using Trizol^® (Invitrogen ™, Thermo Fisher, UK) according to the manufacturer’s instructions except that an additional wash step in 75% ethanol of the RNA pellet was introduced prior to drying and rehydration in 30μl DEPC-treated water. All RNA samples were treated with Turbo DNA-free ™(Ambion^®, Thermo Fisher, Hampshire, UK) to remove contaminating DNA. RNA was then pooled for degenerate PCR, RACE PCR and Illumina HiSeq2500 sequencing, whilst RNA time-course samples were prepared separately for temporal expression analysis.

Degenerate PCR

Initially we used a strategy of degenerate PCR and 5^′ and 3^′ Rapid amplification of cDNA ends (RACE) to identify full-length core canonical clock genes. Full details of the PCR conditions used can be found in File S1 and primers used are tabulated in Table S1.

Illumina RNAseq protocol

The TruSeq cDNA library preparation protocol (Illumina, Cambridge, UK) was carried out on eight time-point samples according to manufacturer’s instructions. Amplified cDNA was run on a 1% agarose gel to validate correct library size range. The cDNA libraries were sequenced on an Illumina RNA TruSeq sequencer (Illumina Inc.) and quality checked at Aberystwyth University’s Core Genomics Facility. Read library quality was investigated using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), thereafter all samples were pooled and de novo transcriptome assembly was performed with Trinity software, version 2012-10-05 (Grabherr et al., 2011). Transcripts smaller than 300 bases were removed from the pooled assembly. This consisted of 186,495 contigs with a mean length of 1,215 bp pairs and an assembly N50 statistic of 2,371 bp.

Transcriptome mining and sequence analysis

Assembled contigs were initially interrogated for canonical clock gene transcripts using the programme BioEdit (Hall, 1999), Basic local alignment (BLAST) searches (http://blast.ncbi.nlm.nih.gov/Blast.cgi) were done using Drosophila melanogaster protein sequences as search terms and tblastn search functions set with default parameters (e-value threshold at 1e⁻¹⁰⁰ except for PDH, where stringency was relaxed to 1000). Where no sequence information was available for Drosophila (e.g., cryptochrome 2) alternative arthropod sequences were used, details of which are given in Table 1. Contigs were translated using EXPASY Translate tool (http://web.expasy.org/translate/) and the coding sequences checked manually using blastP for homology to known circadian proteins.

EMBL SMART (Letunic, Doerks & Bork, 2009) servers were used to detect and analyse the conserved functional domains and motifs of clock proteins. Deduced T. saltator protein sequences were used to query the FlyBase (http://flybase.org) and NCBI non-redundant protein databases using the blastp algorithm. The identity and similarity between protein sequences were calculated using EMBOSS Pairwise Alignment Algorithms (http://www.ebi.ac.uk). The server SignalP v.4 (Petersen et al., 2011) was used to predict the presence and location of signal peptide cleavage sites.

Quantitation of transcript abundance

Using the Trinity downstream analysis tools, reads from each time-point were mapped to the pooled assembly and transcript abundance estimation was carried out using RSEM (Li & Dewey, 2011). In this way the FPKM (Fragments Per Kilobase of exon per Million fragments mapped) values were calculated for each time point. Following this TPM values (transcripts per million) were calculated following Wagner, Kin & Lynch (2012). Quantitative PCR methods for TPM validation are given in File S1 and primers in Table S1. All data files for time-course RNAseq can be found on the NCBI SRA data base, bioproject Accession No. 297565.

Blast2GO analysis- transcriptome annotation

Before running Blast2GO, the number of transcripts was further reduced by clustering them using CD-HIT-EST (version 4.5.4) with default options except for the sequence identity threshold, which was set to 95% (Li & Godzik, 2006; Li, Jaroszewski & Godzik, 2001). The clustering step removes very similar splice variants, which may include some misassemblies, and it significantly reduces the number of transcripts to be included in the time-consuming BLAST analysis. In total 156,766 transcripts with a mean length of 1,534 bp and an assembly N50 statistic of 968 bp were run through BLAST. Using the BLAST output the transcripts were annotated using BLAST2GO software (Conesa et al., 2005; Gotz et al., 2008) including reports for Gene Ontology (GO) terms and EC numbers for the KEGG pathway. The BLAST2GO cut-off parameters used to filter out poor quality BLAST hits for the annotation were as follows: Annotation rule cut-off = 55; E-value = 1e–6; Hit-HSP overlap = 0; and the GO weight = 5.

Determining cycling transcripts

Rhythmicity in clock gene expression was determined using the JTK_CYCLE software (http://openwetware.org/wiki/HughesLab:JTK_Cycle) developed by Professors Michael E. Hughes, Karl Kornacker and John Hogenesch (Hughes, Hogenesch & Kornacker, 2010; Miyazaki et al., 2011) following the JTK_CYCLE Users guide. Changes in gene expression values (TPM) were also analysed by one-way ANOVA.

Data availability

Raw data files have been deposited in public sequencing databases as indicated in the text or at https://figshare.com/s/a2513243c63bf557b720.

Identification of putative Talitrus saltator circadian proteins

Our principle objective was to describe the cDNAs encoding the core elements of the circadian clock system in T. saltator, initially adopting homology cloning and standard sequencing approach but, subsequently superseded by RNAseq strategies. The Illumina RNAseq generated 141,769,456 reads, 128,386,193 of these were assembled into 156,766 clustered contigs (minimum length 300 bp). This T. saltator Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the Accession No. GDUJ00000000. The version described in this paper is the first version, GDUJ01000000. From this transcriptome we identified contigs encoding putative circadian clock genes (Table 1). The coding regions and conserved functional domains of identified ‘core’ clock genes period (Talper), timeless (Taltim), cryptochrome2 (Talcry2) and bmal1 (Talbmal1) and clock (Talclk) are diagrammatically represented in Fig. 2. In the interests of space, alignments of all genes (core and clock-related) are shown in Figs. S1–S26 and all BLAST results are tabulated in Tables S2–S4.

Cryptochrome 2

Degenerate PCR coupled with RACE PCR yielded a 1,843 bp cDNA sequence coding for a 565 amino acid protein; a putative cryptochrome 2 assigned TalCRY2. An identical sequence was revealed in the transcriptome data but assembled as two separate contigs (Table 1). High fidelity PCR amplification and sequencing confirmed that Talcry2 is expressed as one contiguous transcript. The deduced protein sequence contains two SMART identified domains, a DNA photolyase domain and a FAD binding 7 domain (Figs. 2 and S1) and shared an identity of 46.4% and a similarity of 56.0% with Danaus plexippus CRY2 (Accession No. ABA62409; Zhu et al., 2005). The highest sequence identity protein in the NCBI non-redundant database is that for the closely related isopod crustacean, Eurydice pulchra CRY2 sequence (Accession No. AGV28717; Zhang et al., 2013).

Period and timeless

For Talper we mined an 8,001 bp contig encoding a putative full-length protein of 1557 amino acids aligning to dPER with an identity of 29% and a similarity of 50%. The deduced TalPER contains characteristic features of PER proteins including two PAS domains and a PAC domain as well as nuclear export signal, cytoplasmic localization and nuclear localization signal motifs. In common with the isopod E. pulchra, TalPER also has a CK1ε binding site sequence and Period C domain (Figs. 2 and S2). The putative TalPER sequence aligned most closely to dPER-PB in FlyBase (FlyBase No. FBpp0304590) and in the NCBI database to Eurydice pulchra PER (Accession No. AGV28714; Zhang et al., 2013). We also used the Drosophila melanogaster TIM (dTIM) protein (Accession No. AAC46920; Myers et al., 1995) as a search term to identify a 1,209 bp contig with a coding region for a 402 amino acid protein sharing sequence similarity to dTIM (identity of 24.3% and similarity of 41.0%, Fig. S3). SMART analysis of this sequence identified a 228 amino acid Timeless domain (Figs. 2 and S3). The partial protein sequence, assigned TalTIM aligned most closely in the NCBI database to the partial Timeless-like protein of the termite Zootermopsis nevadensis (Accession No. KDR17447; Terrapon et al., 2014) and in FlyBase to Timeout-PA (FlyBase No. FBpp0082180).

Clock

Using the D. melanogaster CLK (dCLK) sequence (Accession No. AAC62234; Bae et al., 1998) as search term, a single 5,723 bp contig was identified as a candidate TalCLK encoding sequence. However, this sequence lacked a stop codon and we were unable to deduce the C-terminal end of the protein, which is missing a C-terminal transactivation domain (Figs. 2 and S4). Nevertheless, our partial sequence contains SMART identified domains including one PAS domain, one PAC domain and coiled coil domains, whilst missing a bHLH domain and a second PAS domain present in dCLK. The TalCLK PAS domain is likely homologous to the PAS-2 domain of dCLK as it aligned better (identity of 64.2% and a similarity of 80.6%) than it did to the PAS-1 domain (identity of 17.3% and a similarity of 36.0%). Moreover, the PAC domains align closely with an identity of 75.0% and a similarity of 90.9% to that of dCLK. The most significant hit for TalCLK in the NCBI database was CLK of the prawn Macrobrachium rosenbergii (Accession No. AAX44045; Yang et al., 2006). In FlyBase the most similar protein is the D. melanogaster CLK orthologue (FlyBase No. FBpp0306710).

Cycle/Bmal1

By searching the transcriptome using the Drosophila melanogaster CYC (dCYC) sequence (Accession No. AAF49107; Adams et al., 2000), we identified two contigs as potential TalCYC/BMAL encoding sequences each of which aligned closely to separate regions of dCYC (Figs. 2 and S5). In fact these two contigs, when correctly oriented, overlap to form one contiguous sequence encoding a characteristic bHLH domain, two PAS domains a PAC and a C terminal transactivation domain. This sequence was confirmed by standard, Hi-fidelity PCR, cloning and sequencing. TalBMAL1 aligned with identity of 43% and similarity of 57% to BMAL1 of signal crayfish Pacifasticus leniusculus (Accession No. AFV39705). TalBMAL1 aligned with 44% identity 58% similarity to E. pulchra BMAL1.

In addition to the ‘core’ clock elements, we also isolated a comprehensive suite of ‘clock-associated’ transcripts from contigs that conceptually translate to proteins known to contribute to circadian clock functioning (Table 1). In the interests of space, we refer the reader to the Supplemental Information; herein the majority of detail is omitted from the text. Sequence alignments are displayed in Figs. S6–S26 whilst BLAST search outputs of identified contig sequences to nucleotide, protein and FlyBase repositories are given in Tables S2–S4.

Casein Kinase II (CKII) α and β subunits

We identified single contigs coding for putative CK11α and CKIIβ. For CK11α we found a 5,147 bp contig enveloping a coding sequence for a 353 amino acid protein including a highly conserved serine/threonine protein kinase catalytic domain (Table 1; Fig S6). This putative TalCKIIα sequence aligned in the NCBI database most closely to human CK2 Chain A (Accession No. 1NA7_A; Niefind et al., 1998) and CKIIα from the copepod crustacean, Paracyclopina nana (Accession No. AII16523). For CKIIβ we identified a 1,567 bp coding sequence for a 220 amino acid protein, including a casein kinase regulatory subunit domain (Table 1; Fig. S6). This sequence aligned most closely in the NCBI database to parasitic wood wasp Orussus abietinus CKIIβ (Accession No. XP_012287730).

Clockwork orange (CWO)

A candidate cwo contig was mined that encodes a partial protein sequence of 167 amino acids (Table 1; Fig. S7) and included a SMART-identified bHLH domain and the Orange of the Hairy/E (SPL) family domain. TalCWO matched most closely in the NCBI database to an uncharacterized protein from the mite Metaseiulus occidentalis (Accession No. XP_003744690); however, the same sequence aligned with CLOCKWORK ORANGE from the firebug Pyrrhocoris apterus (Accession No. AGI17571; Bajgar, Jindra & Dolezel, 2013) with an e value of 8e⁻⁴².

Doubletime (DBT)

A putative 1,092 bp dbt contig coding for a deduced a full-length 310 amino acid protein was identified containing a serine/threonine protein kinase catalytic domain (Table 1; Fig. S8). The deduced sequence aligned most closely in the NCBI database to DBT of the isopod crustacean Eurydice pulchra (Accession No. AGV28719; Zhang et al., 2013) and in FlyBase to the Drosophila melanogaster orthologue DCO (FlyBase No. FBpp0306615).

Par domain protein 1ε (PDP1ε)

A single putative 3,423 bp pdp1ε contig with a coding sequence for a 508 amino acid protein was identified containing a basic leucine zipper domain in the C-terminal region (Table 1; Fig. S9). The deduced TalPDP1ε protein sequence aligned most closely in the NCBI database to the hepatic leukaemia factor of the clonal raider ant Cerapachys biroi (Accession No. EZA50108; Oxley et al., 2014); however PDP1ε from the mosquito Culex quinquefasciatus (Accession No. XP_001865130) was also a top ranking BLAST hit aligned with an e value of 1e⁻³⁵. In FlyBase the sequence aligned most closely to Drosophila melanogaster PDP1 (FlyBase No. FBpp0289727).

Protein phosphatase 1 (PP1)

A 1,725 bp contig was identified coding for a 357 putative PP1 containing a protein phosphatase 2Ac catalytic domain of the serine/threonine phosphatase family (Table 1; Fig. S10). The translated sequence aligned most closely to PP1 of the jumping ant Harpegnathos saltator (Accession No. EFN86649; Bonasio et al., 2010) in the NCBI database. In FlyBase the sequence aligned most closely to Drosophila melanogaster PP1 (FlyBase No. FBpp0306442).

Protein phosphatase 2A (PP2A)

PP2A catalytic subunit TWINS (TWS)

A unique 1,633 bp contig coding for a putative 445 amino acid protein with high identity to Drosophila TWS, (Table 1; Fig. S13) was identified. The deduced sequence contained seven WD40 domains identified using SMART. The protein aligned most closely in the NCBI database to PP2A subunit B of the mud crab Scylla paramamosain (Accession No. AFK24473) and to Drosophila melanogaster TWINS (FlyBase No. FBpp0081671).

Shaggy (SGG)

A 4,413 bp transcript incorporating a coding region for a 418 amino acid protein that contained a SMART identified serine/threonine protein kinase catalytic domain (Table 1; Fig. S14). In the NCBI database, the sequence aligned most closely to glycogen synthase kinase-3 of the turnip sawfly Athalia rosae (Accession No. XP_012256017) and in FlyBase to Drosophila melanogaster SGG (FlyBase No. FBpp0070450).

Supernumerary limbs (SLIMB)

A candidate 588 amino acid putative SLIMB protein encoded within a 2,121 bp contig was mined. This sequence included one D domain of beta-TrCP, one F box and seven WD40 domains (Table 1; Fig. S15). The deduced sequence aligned in the NCBI database most closely to F-box/WD repeat-containing protein 1A of the termite Zootermopsis nevadensis (Accession No. KDR19729; Terrapon et al., 2014) and in FlyBase to Drosophila melanogaster SLIMB (FlyBase No. FBpp0303082).

Vrille (VRI)

A single 3,949 bp contig incorporating coding sequence for a 509 amino acid protein was identified as a putative vrille transcript. The putative TalVRI protein contains a SMART identified basic region leucine zipper domain (Table 1; Fig. S16) and aligned most closely in the NCBI database to nuclear factor interleukin-3-regulated protein of the termite Zootermopsis nevadensis (Accession No. KDR19729; Terrapon et al., 2014) and in FlyBase to Drosophila melanogaster VRI (FBpp0312171).

Ebony

A single 4,380 bp contig coding for a putative 974 amino acid protein with one AMP binding domain, one AMP binding C domain and one PP binding domain was mined (Table 1; Fig. S17). The deduced protein most closely aligned in the NCBI database to the β-alanyl conjugating enzyme of the cockroach Periplaneta americana. P. americana EBONY (Accession No. CAI26307; Blenau & Baumann, 2005) has been shown to have β-alanyl-dopamine (DA) synthase (BAS) enzymatic activity and the Drosophila melanogaster EBONY sequence (Accession No. ABO27280) aligned with the TalEBONY candidate sequence with an e value of 8e⁻⁹⁰, supported by EBONY (FlyBase No. FBpp0083505) being the most closely aligned protein to the TalEBONY query in FlyBase.

Pigment dispersing hormone 1 and II (PDHI and PDHII)

In the Trinity assembly, a 2,430 bp contig was identified by searching for the conserved NSELINS domain. The identified contig included a coding region for a putative 129 amino acid TalPDH prepropeptide (TalPDH-I, Table 1, Fig. S18). This was subsequently extended to 2,471 bp by the addition of a 5^′ UTR by RACE PCR. The TalPDH-I contains a signal peptide between residues 20–21, a 77 amino acid PDH-precursor-related peptide (PPRP) sequence ending in a K-R dibasic cleavage site. Unusually however, the deduced mature peptide sequence is 32 amino acids long and lacks an amidation signal. SMART identified PDH domains are present in the TalPDH. The deduced 32 amino acid Tal-PDH-I mature peptide has an identity of 40.6% and a similarity of 53.1% with the Uca pugilator β-PDH consensus sequence.

A second putative TalPDH encoding transcript was found in the transcriptome of 3,392 bp in length that included a coding region for an 89 amino acid prepropeptide that we assigned TalPDH-II (Table 1; Fig. S18). This prepropeptide sequence was predicted to include a 23 residue signal peptide. The 43 amino acid PPRP thus extends from residue 24 to a K-R dibasic cleavage site at residue 66. The 23 amino acid mature peptide terminates in an amidation signal and has an identity of 47.8% and a similarity of 69.6% with Uca pugilator β-PDH consensus sequence. No domains were identified by SMART but a corresponding PDH domain region was identified.

RORA

One 2,217 bp partial contig, coding a 599 amino acid containing a C4 zinc finger domain was identified (Table 1; Fig. S19) as a predicted candidate for a T. saltator homologue of the hormone receptor RORA. The 3^′ HOLI ligand binding domain present in other RORA proteins is not identified in TalRORA. It is possible that the HOLI ligand binding domain is present in the unsequenced 3^′ section. The most closely aligned protein in the NCBI database is the house fly Musca domestica nuclear hormone receptor HR3 (Accession No. XP_011290218). The Drosophila hormone receptor-like in 46 protein (FlyBase No. FBpp0297438) was the closest protein to TalRORA in FlyBase.

Reverb

A 5,385 bp full length contig was identified, coding for a putative 1110 amino acid protein containing both a C4 zinc finger domain and a HOLI ligand binding domain (Table 1; Fig. S20). The sequence aligned most closely in the NCBI database to nuclear hormone receptor E75, the non-mammalian REVERB homologue in the carpenter ant Camponotus floridanus (Accession No. XP_011259848). The closest aligning protein in FlyBase was the Drosophila melanogaster E75 protein (FBpp0297726).

Sirt 1, 2, 4, 6 and 7

Five contig sequences of lengths 4,033 bp, 2,275 bp, 2,977 bp, 1,209 bp and 5,180 bp all coding for partial proteins of length 955, 376, 354, 402 and 948 amino acids, respectively were identified in the transcriptome as homologues for SIRTUIN proteins 1, 2, 4, 6, and 7. Each sequence contained one SMART identified SIR2 domain (Table 1 and Figs. S21–S25). TalSIRT1 aligned most closely to the trematode Schistosoma mansoni SIRT1 (Accession No. ABG78545) in the NCBI database and Drosophila melanogaster SIRT1 (FBpp0080015) in FlyBase. TalSIRT2 was most closely aligned to the red flour beetle Tribolium castaneum hypothetical protein (Accession No. EFA06770) in the NCBI database and D. melanogaster SIRT2 (FBpp0310647) in FlyBase. The putative TalSIRT4 protein is most closely aligned in the NCBI database to SIRT4 of the Asian citrus psyllid Diaphorina citri (Accession No. XP_008480918) and in FlyBase to the D. melanogaster SIRT4 (FBpp0070817). The TalSIRT6 sequence most closely aligns to the water flea Daphnia pulex hypothetical protein (Accession No. EFX74386) in the NCBI database and in the FlyBase database to D. melanogaster SIRT6 (FBpp0293897). TalSIRT7 most closely aligns in the NCBI database to the SIRT7 protein of the leafcutter bee Megachile rotundata (Accession No. XP_012143211) and aligns most closely in the FlyBase database to the D. melanogaster SIRT7 protein (FBpp0084733).

Jetlag

One 2,454 bp contig coding a partial putative 458 amino acid protein was identified. This candidate putative protein sequence contains one F-box domain and multiple leucine-rich repeat domains (Table 1; Fig. S26). The sequence was most closely aligned in the NCBI database to the red flour beetle Tribolium castaneum F-box/LRR-repeat protein (Accession No. XP_008193983). In the FlyBase database the TalJET sequence aligned most closely to a D. melanogaster protein from the F-box and leucine rich repeat region group (FBpp0111980).

Blast2Go analysis

Gene ontology (GO) annotations were categorised into the three groupings: ‘Biological Process,’ ‘Molecular Function’ and ‘Cellular Component’. These categories of annotation were most easily displayed graphically at ontology level 2 (Fig. 3) where the most frequent Molecular Functions were “Binding” and “Catalytic activity”; the most common Biological Processes were “Cellular and Metabolic”; and the most common Cellular Components were “Cell Components” and “Organelle components.” These GO term outputs were subsequently compared with results from transcriptomes previously generated and analysed including the larva of the moth Plutella xylostella (Xie et al., 2012), the silkworm Bombyx mori, the parasitic nematode Teladorsagia circumcincta (Menon et al., 2012) and the North Atlantic copepod Calanus finmarchicus (Lenz et al., 2014). Proportions of GO terms were all found to be at broadly similar levels in all organisms transcriptomes and for each category analysed; for example, in the category Molecular Function ontology level 2 ‘binding’ terms for T. saltator, P. xylostella, B. mori and T. circumcincta were 51%, 48%, 41% and 40% respectively whilst level 2 ’catalytic activity’ terms were 31%, 36%, 38% and 40% respectively.

Temporal expression of clock genes

JTK_CYCLE analysis applied to TPM values indicated oscillations (limited to a period of between 21–27h) in the core clock genes Talper (q = 1.69⁻⁰⁷) and Talbmal1 (q = 0.001) (Fig. 4). Alternative analyses by ANOVA showed significant differences in mean values at each interval for Talcry2 (F_7,30 = 2.745, P = 0.032), Talper (F_7,30 = 4.86, P = 0.002) and Talbmal1 (F_7,30 = 2.445, P = 0.049). Of the clock-related genes Talpdh-II (q = 5.43⁻⁰⁴), TalckII β (q = 0.0017), Talpp1 (q = 0.023), Talshaggy (q = 0.021, Talebony (q = 0.028), Talsirt1 (q = 0.001), Talsirt7 (q = 0.037) and Taljetlag (q = 0.027) exhibited oscillatory expression profiles by JTK_CYCLE analysis. Interrogation of these data by ANOVA showed TalCK11β (F_7,30 = 3.284, P = 0.014) Talpp1 (F_7,30 = 3.99, P = 0.005), Talsirt1 (F_7,30 = 4.91, P = 0.002) and Talslimb (F_7,30 = 3.19, P = 0.016) to be differentially expressed over time. Talper mRNA accumulation increased approximately 6-fold during early night and peaked towards the middle of the subjective night (CT18-21) before diminishing at CT24 reaching a nadir at CT6. This profile is approximately antiphasic to TalckIIβ, Talpp1, Talsirt1 and Talslimb that all showed peaks in expected daytime. Talcry2 and Talbmal1 respective contigs also showed antiphasic relationships to Talper. However, JTK_CYCLE analysis did not reveal significant rhythmicity in the Talcry2 transcript.

Given the robust rhythm of Talper expression we chose this transcript with which to validate our RNAseq expression data by quantitative PCR. RNA levels of Talper in samples taken from an independent time-course experiment show peak expression (∼3-fold change from CT3) at CT19 (F_7,32 = 4.54, P = 0.01, see Fig. S27). These data support those revealed by RNAseq with peak expression in the early to mid night albeit with a slightly lower amplitude.

Reads for Talcwo were below levels acceptable for analysis and excluded from further investigation.