A deep transcriptomic resource for the copepod crustacean Labidocera madurae: A potential indicator species for assessing near shore ecosystem health

Coral reef ecosystems of many sub-tropical and tropical marine coastal environments have suffered significant degradation from anthropogenic sources. Research to inform management strategies that mitigate stressors and promote a healthy ecosystem has focused on the ecology and physiology of coral reefs and associated organisms. Few studies focus on the surrounding pelagic communities, which are equally important to ecosystem function. Zooplankton, often dominated by small crustaceans such as copepods, is an important food source for invertebrates and fishes, especially larval fishes. The reef-associated zooplankton includes a sub-neustonic copepod family that could serve as an indicator species for the community. Here, we describe the generation of a de novo transcriptome for one such copepod, Labidocera madurae, a pontellid from an intensively-studied coral reef ecosystem, Kāne‘ohe Bay, Oahu, Hawai‘i. The transcriptome was assembled using high-throughput sequence data obtained from whole organisms. It comprised 211,002 unique transcripts, including 72,391 with coding regions. It was assessed for quality and completeness using multiple workflows. Bench-marking-universal-single-copy-orthologs (BUSCO) analysis identified transcripts for 88% of expected eukaryotic core proteins. Targeted gene-discovery analyses included searches for transcripts coding full-length “giant” proteins (>4,000 amino acids), proteins and splice variants of voltage-gated sodium channels, and proteins involved in the circadian signaling pathway. Four different reference transcriptomes were generated and compared for the detection of differential gene expression between copepodites and adult females; 6,229 genes were consistently identified as differentially expressed between the two regardless of reference. Automated bioinformatics analyses and targeted manual gene curation suggest that the de novo assembled L. madurae transcriptome is of high quality and completeness. This transcriptome provides a new resource for assessing the global physiological status of a planktonic species inhabiting a coral reef ecosystem that is subjected to multiple anthropogenic stressors. The workflows provide a template for generating and assessing transcriptomes in other non-model species.


Introduction
Copepods are ubiquitous in aquatic and semi-aquatic habitats, living in marine, estuarine, freshwater and interstitial environments from the deepest ocean trenches to the top of mountain peaks [1].Labidocera madurae is in the family Pontellidae, which are free-living surface dwelling planktonic copepods that are particularly abundant in coastal marine environments [2].The genus Labidocera is a key member of oligotrophic waters surrounding coral reefs in the Pacific and Indian Oceans including Kāne'ohe Bay, Oahu, Hawai'i [3,4].Kāne'ohe Bay has a thriving coral reef community, which has shown significant resilience and the ability to recover from major environmental perturbations, including pollution, eutrophication, high temperatures, and low salinities [5][6][7].It is also one of the best-studied coral reef ecosystems, and serves as a natural laboratory for experimental research on coral reef habitats [8][9].Equally important are the pelagic regions that surround coral reefs, which serve both as a source of food and habitat for reef dwellers.Fishes, corals and other invertebrates have bi-phasic lifestyles: their larvae spend days to months in the plankton before settling nearshore, often within 100 m of their parents [10].Furthermore, planktivorous reef-dwelling fishes and invertebrates depend on the abundant supply of zooplankton brought to them by currents [11][12].Thus, the coral reef ecosystem includes both coral reef areas and the surrounding open water.
The zooplankton community in Kāne'ohe Bay is dominated by copepods including two cyclopoid species in the genus Oithona, two paracalanid species (Bestiolina similis and Parvocalanus crassirostris), and L. madurae [3].Genetic barcoding indicates that while the L. madurae present in Kāne'ohe Bay is genetically unique, it is clearly a member of the L. madurae species complex [3,4].Because L madurae occurs throughout Kāne'ohe Bay and its surrounding inshore waters, and it is moderately abundant year-round, it has the potential to be an indicator species for the pelagic regions of this estuarine system [13,14].As one of the larger copepod species in Kāne'ohe Bay (Fig 1 ), its physiology and behavior has been investigated [15][16][17][18][19].However, L. madurae has been inaccessible to the genetic and genomic research tools that, applied to model organisms, have yielded so much insight into basic biology.As a group, copepods and other crustaceans are under-represented in the number of sequenced genomes and genomic resources.Thus, much of the basic understanding of the taxon potentially available from such resources is lacking.While there are several crustacean genome projects in progress (e.g., the water flea Daphnia magna, the copepods Tigriopus californicus, Tigriopus kingsejongensis, Eurytemora affinis and the amphipod Hyalella azteca), the genome of Daphnia pulex stands out as the only crustacean genome thus far completed, curated, fully annotated, and accessible through a searchable web portal (wFlea-Base; http://wfleabase.org/)[20].
Transcriptomes can be reconstructed with high-throughput sequencing technologies.However, the quality of de novo assemblies is variable [21,22], and poor quality limits their usefulness in physiological and cellular studies that use gene expression profiles.Thus, the goal of this study was to generate a deep and high-quality de novo transcriptome for L. madurae.Furthermore, multiple workflows were used to provide complementary indicators for assessing its quality and depth.RNA was obtained from multiple developmental stages to increase the representation of transcripts, since a significant percentage of genes are silent in any particular stage [23].Bioinformatics tools described in the Methods section were used to assemble and provide an initial evaluation of quality based on assembly, mapping and annotation statistics.This analysis was followed by targeted searches for transcripts encoding proteins of interest: green fluorescent proteins (GFP; Fig 1), the voltage-gated sodium channel (Na V ), and the proteins involved in circadian signalling.All of the proteins are highly conserved across eukaryotes and possess stereotypical structural domains that were used to vet the completeness of the identified sequences.

Sample preparation and RNA sequencing
Total RNA was obtained from two developmental stage groups of L. madurae: mix of copepodites (CIII to CV) and adult females (CVI) (S1 Table ).All animals used here were collected in summer 2015 from central Kāne'ohe Bay (Hawai'i) (Lat: 21˚4'N; Long: 157˚7'W) using surface net tows with a 0.25 m diameter, 125-μm mesh plankton net.The field collection did not require any permits or approval and was performed by PHL and DKH using a personal watercraft.Zooplankton collections were immediately diluted into a bucket containing 5-10 L of seawater and returned to the laboratory.Adult female and copepodite L. madurae were sorted from samples under the microscope, rinsed in filtered seawater, transferred onto a sieve to remove excess seawater and either preserved in RNAlater (Ambion) (adult females) or prepared for immediate RNA extraction (copepodites).The copepodites were inspected for stage distribution prior to total RNA extraction.Three biological samples were obtained for each group with 5 to 6 pooled individual females and approximately 15 to 26 pooled copepodites for each replicate sample (S1 Table ).
Total RNA was extracted using the QIAGEN RNeasy Plus Mini Kit (catalog # 74134) with Qiashredder (catalog # 79654) following the instructions of the manufacturer and stored in a -80˚C freezer.For each sample, RNA concentration and quality were checked using an Agilent model 2100 Bioanalyzer (Agilent Technologies).Total RNA samples were shipped on dry ice to the Georgia Genomics Facility (University of Georgia, Athens, GA; dna.uga.edu) for library preparation and sequencing.Double-stranded cDNA libraries were prepared using the Kapa Stranded mRNA-seq kit (KK8420) following manufacturer's instructions with a mean library insert size of 201-300 bp.Briefly, RNA samples were first purified with two oligo-dT selection (poly(A) enrichment using oligodT beds), and then fragmented and reverse transcribed into double-stranded complementary cDNA.Each sample was tagged with an indexed adapter and paired-end sequenced (151 bp, 300 cycles) using a High Output Flow Cell in a single lane using an Illumina NextSeq instrument (NextSeq 500) (S1 Table ).

De novo assembly and functional annotation
Prior to assembly, raw sequencing reads were assessed for quality using FASTQC (v1.0.0;Illumina Basespace Labs).The six RNA-Seq libraries were quality filtered using FASTQ Toolkit (v.2.0.0;Illumina Basespace Labs) by trimming the first nine bp, removing Illumina adapters (TruSeqLT universal primer) and low quality reads ("Phred" cutoff score !30), and setting the minimum read length to 50 bp.This led to the removal of an average of 11% of reads, leaving from 79 to 85 million reads per sample for the de novo assembly.The resulting reads from the six libraries were combined and assembled using Trinity (v.2.0.6)[24] on the National Center for Genome Analysis Support's (NCGAS; Indiana University, Bloomington, IN, USA) Mason Linux cluster.The initial parameters of Trinity were set to:-seqType fq-CPU 32-max_memory 200G -min_contig_length 300 -normalize_max_read_cov50. The minimum sequence length in the assembly was set to 300 bp.A summary of the assembly statistics was obtained using the script TrinityStat.pl(v2.0.6).Quality-filtered reads were mapped back to the reference using Bowtie2 software (v2.1.0)[25].
Functional annotation was performed in different steps.First, we predicted transcripts with coding regions (CDS) using TransDecoder (v3.0.0) with default settings (minimum open reading frame [ORF] length 100 amino acid and multiple ORFs per transcript) [24].Then, all predicted transcripts with coding region were automatically annotated using a local BLAST webserver on a Beowulf cluster running the NCBI BLAST algorithm [26].The BLASTx algorithm was used to search against the SwissProt protein database [27] (downloaded on 18 th September, 2015 from NCBI) employing a maximum E-value for annotation of 10 −3 .As a third step, the resulting BLAST annotations were mapped against the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database using UniProt [28].The transcripts with GO terms were classified under three categories: biological process, molecular function and cellular component, which are hierarchically organized into levels.Lastly, "Bench-marking universal single-copy orthologs" (BUSCO) software (v1.22) was used to identify core genes: a set of singlecopy genes highly conserved among eukaryotes and thus expected to be present in a complete assembly [29].BUSCO analysis was performed using the Arthropoda dataset consisting of 2,675 single-copy orthologs.

Transcriptome mining and confirmation of protein identification
In addition to the automated annotation step, a targeted approach was used to identify and vet transcripts encoding Na V s and GFPs and circadian signalling system proteins.The complete assembled transcriptome was downloaded to a local Beowulf cluster running the NCBI BLAST algorithm [26], and queried using known protein sequences for transcripts encoding putative homologs of the target groups (GFP: the copepod Pontella mimocerami; Na V and circadian system: fruit fly Drosophila melanogaster, monarch butterfly Danaus plexippus or the copepod Calanus finmarchicus).
Nucleotide sequences with low E-value hits were translated (TranSeq or ExPASy) and then aligned (MAFFT, (v7) [30]) with and checked for homology to the query protein (typically better than 50% identity).Each deduced protein was used to query the NCBI non-redundant proteins (nr) to confirm the annotation.For Na V channels, conserved regions were located in the MAFFT alignments with the C. finmarchicus predicted proteins as a check on the identification.Protein identity was confirmed by the presence of the characteristic four amino acid (DEKA) selectivity filter [31].For GFP proteins, the online program Pfam (v 29.0) [32] was used to check for the presence of a GFP domain.BLAST searches for transcripts encoding putative circadian signaling system proteins including those for core clock, clock-associated, clock input pathway and clock output pathway proteins [33][34][35].The circadian proteins were identified as "full-length" if they exhibit a functional signal sequence (including a "start" methionine) and were flanked on their C-terminus end by a stop codon, while "partial" proteins either lacked a start methionine (referred to as C-terminal partial proteins), or a stop codon (referred to as N-terminal partial proteins), or both of these features (referred to as internal fragment proteins).Next, each predicted L. madurae protein was used as the input query in a BLAST search of the annotated Drosophila protein dataset present in FlyBase (v FB2016_05) [36], except for CRY1 and CRY2.For these two proteins, the extant D. plexippus protein dataset present in GenBank was used for the reciprocal BLAST.The arthropod protein most similar to each L. madurae sequence was subsequently determined by conducting a BLAST search of the non-redundant arthropod protein dataset (taxid:6656) curated at NCBI.Finally, protein structural motifs were analyzed for each of the L. madurae proteins using the online program Pfam (v 29.0) [32].This manual annotation was compared with the KEGG pathway annotation (map0471).
A key member of the circadian system is pigment dispersing hormone (PDH), which undergoes post-translational modification.Thus, the mature structures of L. madurae PDH and several other peptides derived from the PDH preprohormone were deduced using a workflow employed previously for peptide structural prediction in crustaceans, including copepods [37,38].Specifically, the precursor protein in question was assessed for the presence of a signal peptide using the online program SignalP 4.1 [39]; the D-cutoff values of SignalP 4.1 were set to "Sensitive".Prohormone cleavage sites were identified based on homology to known arthropod PDH preprohormone processing schemes.Carboxyl (C)-terminal amidation at glycine residues were predicted by homology to known peptide isoforms, while the sulfation state of tyrosine residues was predicted using the online program "Sulfinator" [40].

Reference transcriptomes and differential gene expression
Four different transcriptomes were constructed and assessed for differential gene expression between copepodites and adult females.In addition to the full transcriptome ("Full") consisting of 211,002 transcripts, three "reference" transcriptomes were generated and searched: 1) "Trinity predicted genes", consisting of unique TR#_c#_g# and the longest "i"; 2) "Full-CDS", which included only transcripts with predicted coding regions using TransDecoder [24] on the full transcriptome; 3) "Pred.genes-CDS", which was derived from the Trinity predicted gene transcriptome and included only transcripts with predicted coding regions using Trans-Decoder [24].
Mapping and statistical analysis were performed using the pipeline described for "Differential expression using a Trinity assembly" [24] employing kallisto for mapping and edgeR for the statistical analysis.We compared these analyses to a second approach using Bowtie as the mapping program, followed by edgeR.Briefly, the quality filtered reads from the six RNA-Seq libraries were mapped against each reference transcriptome using either Bowtie (default settings; v. 2.0.6)[25] or kallisto (default settings; v.0.43.1) [41].Each dataset generated by the mapping program was then tested for statistical significance using the BioConductor package edgeR [42].As implemented by edgeR, prior to statistical testing, RNA-Seq libraries were normalized using the TMM methods (trimmed means of M values), followed by the removal of transcripts with expression levels below 1 count per million (1 cpm).Transcripts with a Benjamini-Hochberg corrected p-value <0.05 were considered differentially expressed (DEGs).Venny (v.2.1) and BioVenn were used to generate Venn diagrams of the DEGs identified using kallisto and Bowtie [43,44].Differential expression of the target genes was analyzed and compared across transcriptomes.

Results and discussion
To date, the majority of publications describing de novo transcriptomes of calanoid copepods have targeted a single genus, Calanus [23,[45][46][47][48].The individuals used in the current study are from the coastal region of Oahu, Hawai'i: they belong to the L. madurae species complex [3,4].Illumina sequencing of six libraries yielded 528 million paired-end reads ranging in length from 50 to 151 base pairs (bp) and over 92% of these reads were of high quality (Phred score !30).These reads were de novo assembled using the Trinity software package (see Methods)(Table 1).The first step in quality assessment was to generate the battery of standard statistical measures characterizing the results.The assembly produced 211,002 transcripts with an average length of 872 bp, a maximum of 23,836 bp and an N50 value of 1,184 bp (Table 1).It contained 153,604 "Trinity predicted genes" that is transcripts with unique "TR# | c#_g#" identifiers (Table 1).Of the "Trinity predicted genes", the majority (127,025) were singletons (83%), with the remaining genes (26,579) possessing from two to 71 "Trinity predicted isoforms" (TR#|c#_g#_i#).This is similar to the percentage reported for C. finmarchicus [23].
For the L. madurae assembly, the mapping rate was high, ranging from 88 to 92% for the six individual samples (Table 1, S1 Table), which is above the suggested cut-off at 80% mapping rate for a successful assembly.Ambiguous mapping, which was ~30% (31-37% of reads that aligned >1 time; S1 Table ), is likely due to the large number of multiple isoforms assembled by Trinity.
The complete de novo transcriptome containing 211,002 transcripts was used in three separate workflows to further assess the quality of the assembly (Fig 2, see methods).First, bioinformatic tools were used to identify transcripts with coding regions (CDS), which were then annotated against SwissProt, Gene Ontology and KEGG databases, followed by BUSCO analysis.Next, targeted protein discovery was focused on large conserved and complex proteins ("giant proteins" and Na V s), proteins of interest of this copepod group (GFPs and crystallins), and proteins members of a complex pathway (circadian signalling system).Finally, several approaches were tested for generating a representative reference transcriptome for gene expression studies.

Functional annotation of the transcriptome
TransDecoder (see Methods) identified 72,391 transcripts with coding regions (CDS; length ! 100 amino acids) in the de novo assembly.Nearly 87% of the CDS retrieved Table 1.De novo assembly and annotation statistics.Labidocera madurae RNA-Seq data from six samples were combined, quality filtered and trimmed and assembled using Trinity software [24].Missing genes (%) 12 * Trinity's hierarchical nomenclature ("TR# | c#_g#_i#") classifies assembled sequences by similarity."TR#" corresponds to gene "families"; unique "TR# | c#_g#" corresponds to predicted "genes".** Minimum sequence length of > 300 bp was set as one of the assembly parameters ***"Complete" is defined as a gene with a predicted length that is within two standard deviations of the BUSCO group mean length that get annotation against the "Eukaryotes databases"."Complete duplicate" indicates that multiple transcripts annotated to the same core gene such as transcripts with predicted isoforms."Fragmented genes" refers to transcripts that encode partial proteins.

Sequencing and Quality Filtering
https://doi.org/10.1371/journal.pone.0186794.t001 significant hits with E-values of 10 −3 or lower when blasted against the SwissProt database, and over 95% of these were further annotated with gene ontology terms (Table 1).Within the "biological process" category, L. madurae transcripts covered broadly conserved eukaryotic processes with "cellular process", "metabolic process" and "single-organism process" representing more than 60% of the annotated transcripts (Fig 3).Eighty percent of transcripts with GO terms were annotated within the KEGG database (Table 1), indicating good coverage of transcripts encoding proteins/enzymes involved in lipid, amino acid and energy metabolism pathways (S1 Fig) .BUSCO analysis identified 76% (2,036) complete orthologs of 2,679 core eukaryotic genes in the CDS with <1% of these genes present in more than one copy (duplicated).An additional 11% of fragmented core genes were found among the CDS, with only 12% of core genes missing completely (Table 1).
The assembly and annotation statistics of the L. madurae de novo transcriptome were compared with those of other non-model arthropods: three insect species and five other copepods (Table 2) [23,47,[49][50][51][52][53].The number of assembled transcripts is quite variable across de novo transcriptomes with the number in the L. madurae transcriptome (~200K) being among the highest (Tables 1 and 2).The number of transcripts with coding regions is higher in copepods, including L. madurae, than that reported for the insect, Lygus hesperus (Western tarnished plant bug) [49].Interestingly, the L. madurae annotation rates (87% of transcripts with coding regions) were higher than those reported in the other copepods which can in part be attributed to limiting annotation to protein encoding transcripts (Table 2).The number of predicted core proteins was similar across the transcriptomes with an approximate coverage of 80 to 90% based on the BUSCO analysis (Table 2).Overall, the annotation statistics suggests that the L. madurae transcriptome is at least as good in quality and depth as the others with which it was compared.
The large number of putative lncRNA transcripts in L.madurae suggests that there may be more lncRNA loci in this crustacean than in D. melanogaster [54][55].However, a shotgun assembly only produces predicted transcripts, and further analyses are needed to confirm which transcripts are indeed lncRNAs, as opposed to genes coding for very small proteins (<100 amino acids long), incomplete transcripts, or assembly artifacts (e.g.fragmented UTRs which have been found in this transcriptome).

Searches of target genes based on automated annotation
"Giant" proteins.The presence of transcripts encoding "giant" proteins (those >4,000 amino acids) was used as an indicator of quality of the assembly.The L. madurae assembly included 23 transcripts that exceeded 15,000 bp in length.The lengths of these transcripts are comparable to those reported for six of the transcriptomes listed in Table 2.The majority of the long transcripts encoded "giant" proteins belonging to titin/connectin family, such as "twitchin", and proteins involved in cellular architecture/cytoskeleton such as "nesprin".Examples of long transcripts, all of which are predicted to be full-length and annotations are given in Table 3.
Crystallins.An unusual feature of Labidocera and other pontellids is a sophisticated frontal eye structure [56,57].Unlike most copepods, the pontellid eye includes a clear lens, which requires structural proteins that are both stable and transparent.However not much is known about the structure of invertebrate lenses [58].In vertebrates, the structural proteins of lenses include crystallins, which have been well characterized.A search of the L. madurae list of automated annotated transcripts identified 20 putative crystallins.Fifteen of these encode putative α-crystallins, with others encoding putative members of the β-crystallin (2), the γ-crystallin (1) and λ-crystallin (1) families (S2 Table ).The β-and γ-crystallins, which form a partnership with α-crystallins, are the primary structural proteins of the vertebrate lens [59,60].Thus, one or more of these transcripts might be involved in lens formation in L. madurae.* BUSCO analysis was performed in 2017 using publicly accessible NCBI "transcriptome shotgun assembly".TSA data were first processed using transdecoder, followed by BUSCO (v.1.22)specifying the "Arthropoda" dataset, which included 2,675 core genes-analysis.TSA accession numbers: GAXK00000000 (C.finmarchicus), GCJT01000000 (P.nana), GCHA01000000 (T.japonicus), GDFW00000000 (T.kingsejongensis) ** # of transcripts given is the number of isotigs, N50 value is the isotig N50.L. madurae de novo assembly included a significant number of contigs (>100K), which lacked an open reading frame.Many of these non-coding sequences could belong to a class of transcripts called "long (>200 nucleotides) non-coding RNAs" (lncRNAs).While these sequences are often omitted from de novo transcriptomes, they are unlikely to be "assembly artifacts". https://doi.org/10.1371/journal.pone.0186794.t002

Manual sequence annotation using targeted gene discovery
Green fluorescent proteins (GFPs).Pontellids are well known for the presence of GFPs, which include some of the brightest GFPs currently known [61].In L. madurae, GFPs are concentrated at the base of the appendages as seen in the side view of an adult female (Fig 1C and  1D).Three transcripts were found that putatively encode GFPs (S2 Table ).Two of the predicted proteins, both full lengths, shared 90% amino acid identity with a pair of GFPs identified in a closely related species, Pontella mimocerami [61].The third L. madurae GFP is most similar to a jellyfish (Aequorea victoria) GFP with which it shares 90% amino acid identity (S2 Table ); this protein appears to represent a new class of copepod GFP.These putative transcripts encoding crystallins could serve as a starting point for any study investigating lens formation in copepods, specifically the pontellids, which possess modified naupliar eyes.
Large proteins with splice variants: voltage-gated sodium channels (Na V ).Large proteins that belong to families with closely-related members and which possess multiple splice sites or other regions of variation can be challenging to assemble and group dependably.One such protein family comprises the Na V s.In arthropods and in particular copepods de novo transcriptomes, incomplete or fragmented genes are common within this family (e.g.see publicly accessible transcriptomes in the following references: [23,45,48,52] and NCBI Bioprojects PRJEB20069, PRJNA231234).Thus, as a stringent test of transcriptome quality, we assessed the assembly of the L. madurae Na V s proteins (Labma Na V s), comparing it with that from our previously published well-vetted transcriptome for C. finmarchicus [23,38,62,63].We examined whether expectations were met in: 1) the number and completeness of predicted Na V genes, identified by their expected characteristics (match statistics, conserved motifs, length); 2) the occurrence and nature of predicted splice variants; 3) how well Na V s were grouped into the Trinity hierarchy; 4) the occurrence and nature of irregularities (incorrect or incomplete sequences).
Characteristics of Na V s expected to be present in an invertebrate transcriptome include occurrence of contigs from two families of orthologous genes, designated Na V 1 and Na V 2 [64].However, in L. madurae three predicted gene families (TR#) were identified as Na V s by the automated annotation.This is one more than expected (Table 4).These had low E-values (<8e-156) and were identified either as para or 60E, the D. melanogaster designations for Na V 1 and Na V 2 respectively.Querying the full transcriptome with a well-vetted arthropod sodium-channel sequence from D. melanogaster (SwissProt SP3500) retrieved 13 sequences from the same three gene families with E-values < 1e-88.Sequences with the next higher Evalues had features of voltage-gated calcium channels.The retrieved sequences are shown diagrammatically in Fig 4 .ReBLASTing each of the Na V contigs into Flybase returned either para Table 3. Giant proteins.Four transcripts encoding "giant" proteins assembled using Trinity software in Labidocera madurae transcriptome.For each transcript, transcript length, predicted protein length, annotation name (NCBI), Accession No. of top blast hit (NCBI), E-value annotation (NCBI), protein family and protein function are listed.

TR75346|c7_g2_i1
TR27483|c2_g1_i1 TR79107|c1_g1_i1 TR75290|c0_g1_i1 Transcript length (bp) or 60E (Table 4).To further resolve the identity of the contigs, they were used to query the C. finmarchicus transcriptome [23], retrieving top hits for 7 isoforms corresponding to Na Full-length proteins of the Na V family are expected to be around 200 kD in size.Completeness of predicted proteins was verified for one or more contigs from each Labma Na V 1 gene as well as from the single reconstructed Labma Na V 2 gene (see below).Start and stop codons as well as 5' and 3' UTRs are present in all three.When all optional sequence segments (putative exons) are included, predicted proteins 2072 and 2069 amino acids long result for Labma Na V 1.1 and Labma Na V 1.2, respectively.These match the lengths predicted for corresponding genes of C. finmarchicus (2094 and 2079 respectively) [23], for D. melanogaster Na V 1 (2131; UniProtKB P35500), and for human Na V 1.1 (2009; UniProtKB P35498).Similarly, the 2533 residue length of the reconstructed Labma Na V 2 was within 2% of that for C. finmarchicus (2485aa) and 10% of that for D. melanogaster (2821aa).Thus, three Na V genes, with appropriate characteristics, are well assembled in the L. madurae transcriptome.Two or more sites of splice variation separated by more than a cDNA-insert-length of identical bridging sequence 2 Top BLASTp result from Flybase annotated proteins; "para" = Na V 1; "NaCP60E" = Na V 2 3 Original identification based on automated annotation 4 Sodium channel not fully characterized The Drosophila melanogaster Na V 1 sequence (sp|P{35500) para was used as a query in a tBLASTn probe of the Labidocara madurae 2015 transcriptome (column Drome e-value) The top hits (Trinity ID number column), with e-values < e-84, were translated into protein sequnces and reblasted using the tBLASTn tool against the Calanus finmarchicus Gulf of Maine transcriptome [23].The top hits from that BLAST are indicated in the column "C.finmarchicus top hit," with e-values given in the column "Calfi e."These are used to identify the protein (column "Labma name") using the correspondence of comp222993 and comp299307 with Na V 1.1, comp44060 and comp233807 with Na V 1.  value), so the associations implied by the contigs assembled that include those two regions are unreliable.This does not imply a poorer quality of assembly compared with other paired-end assemblies of cDNA inserts of the same length: it is intrinsic to the shotgun approach.This caveat applies to four of the seven contigs of Labma Na V 1.1 (Fig 4), but as well to the long contigs (18 in all) of Calfi Na V 1.1 (see Fig 10 of Lenz et al [23]).Despite this ambiguity, the Labma Na V 1.1 contigs gave solid evidence for the presence of four optional segments at Site I and one alternative segment at Site II, which is qualitatively similar to the pattern found in C. finmarchicus.No clear evidence for splice variants was found for Labma Na V 1.2 (i2 is an anomalous fragment, possibly artifactual), the same being the case for Calfi Na V 1.2.For Labma Na V 2, the two members of each pair of fragments (TR65477 and TR68660) differ in the presence of "optional" segments in each, a feature not found in Calfi Na V 2 (arrows in the Na V 2 diagram of Fig 4).Thus, aside from this last case, the L. madurae transcriptome showed splicing features expected from the C. finmarchicus assembly.Most differences in the details (see below) are likely species differences.Hierarchical transcript grouping by Trinity, as outlined in Methods, enables classifying assembled sequences into likely gene families, genes and isoforms.It performed well on the Labma Na V 1 genes, separating them correctly into two genes nested within a single family.In contrast, transcripts for the same Calfi Na V 1 genes are more broadly assigned, spanning four "Chrysalis components" (comps = gene proxies; Table 4) [24].Reassembly of the C. finmarchicus transcriptome using Trinity 2.0.6only reduced this number from four to three and failed to include them in the same gene family.Thus the L. madurae transcriptome is of higher quality in this respect than that of C. finmarchicus.On the other hand, a single transcript coded for Calfi Na V 2, while Labma Na V 2 was present as two fragments assigned to different Trinity (2.0.6) predicted gene families (Table 4).Still, these fragments had overlapping ends and could be amalgamated to form a full-length predicted protein with all of the expected properties.Thus the overall structure of the three Na V genes was successfully assembled in the L. madurae transcriptome with about the same quality as for that of the C. finmarchicus.
Irregularities in the L. madurae assembly were of several types, described in more detail in S5 Fig.To summarize, the number of Labma Na V s assembled was smaller (three vs. six) than for C. finmarchicus.This is likely in part a species difference.Anomalous sequences of various origins were also noted.These include a short contig (TR25803|c0_g1_i1) that may represent an additional Labma Na V 1 (Table 4) and a sequence with a frame-shift that is probably an error.In addition, several issues appear to have arisen from the ambiguity in assembling regions of variation bridged by segments with identical sequences that are longer than one cDNA-insert-length: 1) isoforms, especially within the Labma Na V 1.1 gene, code for partial rather than full-length proteins (Table 4); 2) Calfi Na V 1.1s have many more full-length contigs (18 vs 2) perhaps reflecting a greater leniency of Trinity 1.0 for matching variable regions; 3) genetic variability within the population may have increased the number of variable regions, possibly contributing to premature truncation of sequences.
Overall, the L. madurae transcriptome assembled Na V s as well as or better than that of C. finmarchicus [23].However, it highlighted the limitations inherent in matching variant segments separated by stretches of identical sequence longer than a cDNA-insert-length.
Translation of the identified transcripts revealed that the vast majority encoded full-length proteins (Table 5, S3 Table ), with just two encoding partial sequences (Table 5).For many protein groups, multiple variants, all likely derived from a common gene, were predicted.These variants were most likely derived from alternative splicing, as well as single nucleotide polymorphisms (e.g., the five CYC variants shown in S2 Fig) .In addition, for a number of groups, proteins derived from multiple genes were identified (e.g., the four distinct PP1s shown in S3 Fig) .PDP1 was represented with four predicted genes, one with a splice variant, as shown in Fig 5 .While parts of the molecule were very conserved, there were significant differences between the predicted proteins, which may reflect diversity in function.In the case of the CRY2 protein, 12 distinct transcripts were identified, and while they differed in length (Table 5), the predicted proteins were all identical.These transcripts differed in the two untranslated regions (5'UTR and 3'UTR), which may be related to differential processing and/ or tissue-specific expression.
In addition to vetting the completeness/quality of the L. madurae transcriptome, the mining of this resource for circadian protein-encoding transcripts has shed light on the clock system of this species, and for that matter, those of crustaceans in general.The large suite of proteins predicted from the Labidocera transcriptome (Table 5), include, among others, the canonical core clock proteins CLK, CYC, PER and TIM, all showing significant homology to those of D. melanogaster (S4 Table ).They possess structural domains consistent with their fruit fly homologs, domains required for normal function (S1 Fig) .Moreover, putative L. madurae homologs of both CRY1 and CRY2 were identified (Table 5), a finding that suggests that the Labidocera circadian system is organized more similarly to the "ancestral-type" clock proposed for lepidopteran/mosquito species than to that of D. melanogaster [66].Specifically, CRY2, which is missing in Drosophila, but participates in the core clock itself, is likely to be a repressor of CLK-CYC-mediated transcription, while CRY1 functions as a photoreceptor, putatively providing photic input to the core clock.This result is consistent with the "ancestral-type" circadian systems described in other crustaceans that have been examined via genome/ transcriptome analyses [33][34][35]69], suggesting that this type of clock organization is broadly conserved within members of this arthropod subphylum.
The mining of the Labidocera transcriptome resulted in the discovery of the first PDP1s from a member of the Copepoda.The results suggest the presence of multiple genes from several protein families: DBT (three genes), PDP1 (four genes), PP1 (four genes), MTS (two genes), TWS (two genes) and SGG (two genes).No members of PDP1 had been identified previously from either C. finmarchicus or T. californicus [33,34].The identification of the L. madurae PDP1 genes allowed for the revisitation of the C. finmarchicus and T. californicus transcriptomes for putative homologs.Using the Labidocera PDP1 predicted proteins as queries, related proteins have now been discovered in these two copepod species (A.E. Christie, unpublished).Moreover, mining of the assembly led to the prediction of a novel isoform of PDH, NSEMLHILRSMPKDMGKIIRNamide, which is just the second member of this peptide family identified from a copepod [37], a peptide that may serve as an output signal from the Labidocera clock for controlling its physiology and behavior.Comparison between the results from the targeted gene discovery workflow with the results from the automated annotation is shown in Fig 6 .The circadian system pathway retrieved from the KEGG database (map0471) resulted in the identification of five of the eight expected genes (Fig 6).The automated annotation programs failed to identify VRI, PDP1 and PER among the L. madurae transcripts with coding regions (CDS).These results underscore the value of targeted gene discovery in combination with the automated bioinformatics tools to obtain a complete annotation for a de novo transcriptome.

Reference transcriptome analysis
Identification of differentially expressed genes between L. madurae developmental stages.The generation of a transcriptome that provides robust results for gene expression profiling is key for application to physiological ecology.While sequenced and annotated genomes are used as reference in model species, de novo assembled transcriptomes, in combination with bioinformatic tools for annotation and statistical testing, provide a powerful alternative.However, for a transcriptome of a non-model species to be used as an alternative for a genome, it needs to be of high quality and complete.Here, we compare four strategies for obtaining a reference for read mapping and identification of differentially expressed genes (DEGs).While the full transcriptome (211,002 transcripts) is optimal for targeted gene discovery, including the identification of genetic variants (i.e., splice variants, indels, SNPs), it also generates a large percentage of ambiguous mapping that could affect statistical testing.In addition to the full transcriptome ("Full"), we generated three alternative "reference" transcriptomes from the "Full" assembly by: 1) selecting the longest transcript for Trinity predicted genes (unique TR#_c#_g#; "Pred.genes"); 2) selecting only transcripts with coding regions (CDS) ("Full-CDS"); and 3) selecting only transcripts with coding regions (CDS) from the "Trinity predicted genes" transcriptome ("Pred.genes-CDS").
Table 6 shows the effects of applying these filters.The number of transcripts decreased from 211K to 45K in the smallest "reference".Nevertheless, the four transcriptomes were comparable with respect to the number of core eukaryotic proteins, which declined only by 3% between the full and the Trinity-predicted "unique" gene transcriptomes ("Pred.gene", "Pred.gene-CDS").With the exception of the full transcriptome, the number of duplicated genes (genes with more then one copy) was low (< 0.5%).The percentage of mapped reads using Bowtie decreased from 91% to 68% between the Full and Pred.genes-CDS references Furthermore, the three derived reference transcriptomes had fewer ambiguous reads than the full transcriptome, and the "unique gene" approach led to the lowest number of reads mapped more than once (14% and 6% for "Pred.genes" and "Pred.genes-CDS", respectively).Differences among these potential "reference transcriptomes" were further evaluated by testing for differential gene expression between copepodites and adult females.Thus, we mapped reads to the four "reference transcriptomes" using two different bioinformatics tools (Bowtie and kallisto) followed by statistical testing to DEGs (edgeR).While the number of counts (= mapped reads) associated with each transcript is higher in Bowtie than in kallisto, this did not affect the number of transcripts tested for relative gene expression after applying the 1 cpm filter (Table 6).The number of DEGs identified by edgeR using counts generated by Bowtie varied by more than a factor of two among the references used.Nevertheless, 8,970 DEGs were shared among the four references (S4 Fig) .In contrast, the number of DEGs identified with kallisto was similar for all four transcriptomes (Table 6), with 6,229 shared among all references ( Fig 7).A comparison between Bowtie and kallisto of the shared DEGs identified 5,438 common DEGs (S4 Fig) .The smallest reference transcriptome ("Pred.genes-CDS") had best agreement between Bowtie and kallisto with 9,827 shared DEGs, which represented approximately 89% (kallisto) and 77% (Bowtie) of identified DEGs, which is not surprising given that this transcriptome had the smallest number of ambiguous reads (S4 Fig) .In general, mapping by kallisto is more conservative, making it the preferred mapping program for the identification of DEGs, in particular in association with an assembly program like Trinity, which is designed to preserve isoform variants [24].
While the large number of shared DEGs regardless of mapping program or reference transcriptome (5,438 DEGs) was reassuring, there were still many of DEGs that were identified in one or two references but not the others as shown in the Venn diagram for DEGs generated from kallisto mapped reads ( Fig 7).There was good agreement between the Full and Full-CDS (9,637 DEGs) and the Pred.genes and Pred.genes-CDS (9,028 DEGs; Fig 7).
Differential expression of transcripts identified through targeted gene discovery.To gain further insight into differences in expression, we examined expression results for the targeted genes identified in the previous sections (Tables 3, 4 and 5).For all investigated Table 6.Comparison across four possible reference transcriptomes generated from the de novo assembly for gene expression studies.Reference transcriptomes-"Full": complete de novo Trinity assembly; "Pred.genes": retained a single (longest) isoform each Trinity-defined unique genes; "Full-CDS": de novo Trinity assembly filtered using TransDecoder with only transcripts with predicted coding regions retained; "Pred.genes-CDS": "Pred.genes" transcriptome filtered using TransDecoder with only transcripts with predicted coding regions retained.Number of transcripts, Bowtie mapping statistics and BUSCO analysis is given for each reference.Differential gene expression results include the number of transcripts that were included in the statistical analysis (expression level: > 1 cpm) and number of identified differentially expressed genes (DEGS) using either Bowtie or kallisto software as the mapping program.transcripts, expression rate was higher such that the transcripts did pass the 1 cpm filter and were considered for the statistical test.The transcripts encoding "giant" proteins were represented in all reference transcriptomes, and two transcripts, fibrillin-1 and nesprin-1, were consistently identified as differentially expressed (Table 7).Other target genes that contributed to the shared DEGs (6,229) included a Na V (Labma1.2) and one transcript each of PER-v1, CWO-v1 and VRI (Table 7; Fig 7).The transcriptomes differed in the number of Na V transcripts given the presence of isoforms.Thus, Na V 1.1 and 1.2 had seven and two isoforms, respectively, in the Full and Full-CDS transcriptomes, while the two unique gene transcriptomes (Pred.genes and Pred.genes-CDS) had single transcripts representing each of these two genes (Table 7).Two isoforms (i3, i4) of the Na V 1.1 transcripts were differentially expressed in the Full and Full-CDS transcriptomes, however, the single Na V 1.1 in the other two transcriptomes was not among the DEGs (Table 7).
Several CRY2 isoforms were identified as differentially expressed in the Full and Full-CDS transcriptomes, but not in the Pred.genes and Pred.genes-CDS references (Table 7).In these two transcriptomes the "i6" isoform was among the DEGs (Table 7).A similar pattern was observed with JET-two out of three isoforms were among the DEGs, while the third isoform The references transcriptomes are defined as: "Full" with 211K transcripts (purple), "Pred.genes" consisting of longest transcript for Trinity predicted genes (yellow), "Pred.genes-CDS"consisting of transcripts with coding regions (CDS) from the "Pred.genes"(green) and "Full-CDS" consisting of transcripts with coding regions (CDS) from "Full" (pink).Relative transcript abundance as determined using kallisto, and DEGs were identified by statistical analysis using edgeR with P<0.05 and false discovery rate (FDR) cutoff at 5%. https://doi.org/10.1371/journal.pone.0186794.g007Table 7.Comparison among reference transcriptomes in the identification of differentially expressed genes (DEGs) between L. madurae copepodites and adult females among transcripts encoding for "giant" proteins, voltage-gated sodium channels and circadian system proteins.Transcripts were identified as DEGs using a Benjamini-Hochberg corrected p-value <0.05.

Transcript
Reference transcriptomes Protein name Trinity identification # "Full" "Pred.genes" "Full-CDS" "Pred.genes-CDS" "Giants"  was identified as differentially expressed in the references with single isoforms (Table 7).While these are examples of disagreements between the references (Fig 7 ), the results are consistent in identifying at least one isoform of CRY2 and JET as differentially expressed.

Voltatge-gated sodium channel
Another pattern that occurred was the inclusion of multiple isoforms among the DEGs in both Full and Full-CDS transcriptomes, but not in the "unique gene" ones (e.g., CYC, Tim, PDP1-III and PDP1-IV).The reverse, differentially expressed according to the "unique gene" transcriptomes, but not the other two, occurred for transcripts of one doubletime (Labma-DBT-II), one PAR-domain protein 1 (Labma-PDP1-II) and pigment-dispersing hormone (Labma-PDH).Four DEGs were identified in a single reference (3 in Full-CDS and 1 in Pred.genes-CDS), while one DEG was shared between the two CDS-based reference transcriptomes (Table 7).In summary, comparing DEGs identified with four reference transcriptomes for the target genes indicated good concordance between the Full and the CDS-based transcriptomes (29/33) and the two "unique gene" references (Pred.genes and Pred.genes-CDS: 11/13).Agreement between all four transcriptomes regardless of isoform was observed in eight out of 13 genes.Inconsistent results across reference transcriptomes are typically associated with transcripts belonging to genes with multiple isoforms, such as those with predicted splice variants.Thus, independent of the method used for generating a reference transcriptome, it is important to assess the number of isoforms for each differentially expressed gene.

Conclusions
High-throughput sequencing in combination with bioinformatics tools has made transcriptomic approaches accessible to non-model species, including those of ecological interest.Thus, transcriptomics can now be used to investigate the eco-physiology of key species within the context of life history strategies, population cycles and ecosystem dynamics.However, these types of studies, which involve gene expression profiling, depend on good reference transcriptomes.Application of multiple workflows to evaluate the quality and completeness of a transcriptome generated for the copepod L. madurae demonstrates that no single criterion is sufficient to assess a de novo assembly.High-throughput bioinformatics tools were used to identify transcripts with protein coding regions and provide annotations.Targeted gene discovery provided information on completeness of individual genes, identified possible sources of fragmentation, established predicted gene variants, and provided additional annotations.The analysis of four different strategies for generating a reference for gene expression studies suggest good agreement among references when a predicted gene assembled into a single isoform.However, many predicted genes include a multiplicity of isoforms, and when these are included in the reference they contribute to ambiguous mapping.Thus, one source of disagreement among transcriptomes in the identification of DEGs is related to which genes are regulated, and weather they are represented by multiple isoforms.The workflows developed in this study if used in a routine assessment of de novo transcriptomes would enhance the reliability of gene expression studies.

Fig 1 .
Fig 1.Light micrographs of Labidocera madurae copepodite (A, B) and adult female (C,D).(A) Copepodite stage CIII, dorsal view (magnification: 4x).(B) Same copepodite as in A under fluorescent light showing expression of green fluorescent protein (GFP) (magnification 10x).(C) lateral view of the anterior portion of an adult female showing one dorsal and the ventral ocelli, feeding appendages and GFP expression (magnification 10x).(D) Lateral view of the same individual as in C under fluorescent light showing GFP expression at the base of the swimming legs (magnification 10x).Scale bar: 0.5 mm.https://doi.org/10.1371/journal.pone.0186794.g001

Fig 2 .
Fig 2. Diagram of the workflow used to generate the de novo transcriptome for Labidocera madurae and the three approaches used to test for completeness and quality of the assembly.https://doi.org/10.1371/journal.pone.0186794.g002

V 1. 1 (
TR7852|c0_g1), 2 corresponding to Na V 1.2 (TR7852|c0_g2) and 4 in two TR groupings corresponding to the Na V 2 gene (TR65477_c0_g1 and TR68660_c0_g1).The motifs expected of Na V s, shown diagrammatically at the top of Fig 4 (see caption for details), were validated through sequence alignment (MAFFT) and hence the various groups have been designated Labma Na V 1.1, 1.2 and 2.

Fig 4 .
Fig 4. Labidocera madurae voltage-gated sodium channel sequences assembled by Trinity.Diagram at top shows the four wellconserved domains (DI-DIV) bridged by less-well-conserved loops.Conserved domains are depicted vertically expanded to show approximate locations of six trans-membrane α-helical segments (colored bands labeled S1, S2-S6).Sodium-selectivity of the Na V 1 transcripts (but not Na V 2) is confirmed by the occurrence of four characteristic amino acids (aspartic acid, glutamic acid, lysine and alanine [DEKA]) in specific locations termed the "P-loops"[31].Coverage by variants of three putative genes, Labma NaV 1.1 Labma Na V 1.2 and Labma Na V 2 indicated by bars labeled with the i number assigned by Trinity.For Labma Na V 1.1, no one sequence possessed all of the pieces (putative exons), so the overall span across the diagram represents a manual reconstruction generated by including all of the pieces from the different i's.Gaps in sequences are indicated by fine dotted lines.Identical 5' (504 nucleotide) UTRs for i1-i7 have been omitted, as have the identical 3' UTRs (1518 nucleotides) of i1 and i2.Within each gene, corresponding residues across different i's were identical (reflected in the same coloration of the bars) in almost all cases, except for the splice variant indicated in red for Na V 1.1 i3.Sequences representing partial predicted proteins not initiated by an M at the N-terminal or terminated by a stop codon ("X" above the bar) at the C-terminal are indicated with a short diagonal bar.Positions of the domains for Na V 2 differ somewhat from those of Na V 1 shown in the top diagram and are indicated by thickening of the bars.Two sites of putative splice variation (Site I and II) are indicated below the Na V 1.1 diagram, and one non-optional segment within Site I is designated "1" (96aa).Arrows in the Na V 2 diagram indicate short optional pieces (gaps in the horizontal bars), and the overlap region between the two pairs of isoforms of 44 identical amino acids (aa) is indicated.https://doi.org/10.1371/journal.pone.0186794.g004 Protein type abbreviations: F, full-length protein; N, amino (N)-terminal partial protein; C, carboxyl (C)-terminal partial protein.Proteins used as queries in tblastn searches: CLK, Drosophila melanogaster CLK (Accession No. AAC62234); CRY2, Danaus plexippus CRY2 (Accession No. ABA62409); CYC, D. melanogaster CYC (Accession No. AAF49107); PER, D. melanogaster PER, isoform A (Accession No. AAF45804); TIM, D. melanogaster TIM (Accession No. AAC46920); CKII α, D. melanogaster CKIIα, isoform A (Accession No. AAN11415); CKIIß, D. melanogaster CKIIß, isoform B (Accession No. AAF48093); CWO, D. melanogaster CWO, isoform A (Accession No. AAF54527); DBT, D. melanogaster discs overgrown, isoform A (Accession No. AAF57110); JET, D. melanogaster JET, isoform A (Accession No. AAF52178); PDP1, D. melanogaster PDP1, isoform B (Accession No. AAN12022); PP1, D. melanogaster PP1 (Accession No. CAA39821); MTS, D. melanogaster MTS, isoform A (Accession No. AAF52567); TWS, D. melanogaster TWS, isoform A (Accession No. AAF54498); WDB, D. melanogaster WDB, isoform A (Accession No. AAF56720); SGG, D. melanogaster SGG, isoform A (Accession No. AAN09082); SLIMB, D. melanogaster SLIMB, isoform A (Accession No. AAF55853); VRI, D. melanogaster VRI, isoform A (Accession No. AAF52237); CRY1, D. plexippus CRY (Accession No. AAX58599); PDH, Eucyclops serrulatus Prepro-PDH I (deduced from Accession No. GARW01021210); PDHR, D. melanogaster pigment dispersing factor receptor, isoform A (Accession No. AAF45788).https://doi.org/10.1371/journal.pone.0186794.t005

Fig 6 .
Fig 6.Predicted gene mapping to the circadian rhythm pathway obtained through KEGG annotation.Circadian rhythm pathway shown represents a map for Drosophila melanogaster (map04711).Highlighted boxes (green) represent L. madurae transcripts with coding regions (CDS) automatically annotated against the Kyoto Encyclopedia of Genes and Genomes (KEGG).PER, VRI, PDP1 were not identified by the automated annotation (white boxes).https://doi.org/10.1371/journal.pone.0186794.g006

Fig 7 .
Fig 7. Non-proportional Venn diagram for the number of differentially expressed genes (DEGs) identified using four different transcriptomes as a reference for mapping of reads.The references transcriptomes are defined as:"Full" with 211K transcripts (purple), "Pred.genes" consisting of longest transcript for Trinity predicted genes (yellow), "Pred.genes-CDS"consisting of transcripts with coding regions (CDS) from the "Pred.genes"(green) and "Full-CDS" consisting of transcripts with coding regions (CDS) from "Full" (pink).Relative transcript abundance as determined using kallisto, and DEGs were identified by statistical analysis using edgeR with P<0.05 and false discovery rate (FDR) cutoff at 5%.
in the reference transcriptome but not differentially expressed YES Transcript present in the reference transcriptome and differentially expressed X Transcript not present in the reference transcriptome In bold Transcripts resulting differentially expressed in all 4 reference transcriptomes https://doi.org/10.1371/journal.pone.0186794.t007
https://doi.org/10.1371/journal.pone.0186794.t004cannotbeassembledreliablywithout additional information.Labma Na V 1.1 has two sites with variant segments at opposite ends of the molecule.SiteI is an N-terminal region of optional segments (putative exons;Fig 4); site II is an alternatively spliced segment nearer the C-terminal end.Both sites correspond to ones in Calfi Na V 1.1 (Table4).The two sites are separated by a minimum of 630 residues in L. madurae, well over a cDNA-insert-length (200-300 bp mean

Table ) ,
and vetted via reciprocal BLAST searches (S4Table and S5 Table) and protein structural motif analysis (S6 Table

Circadian signaling system protein Transcript/protein identifications Transcript Deduced protein Clock component Family Trinity identification number Length* Name Length + Type
*Length in nucleotides.+ Length in amino acids.