Genome-wide profiling of 24 hr diel rhythmicity in the water flea, Daphnia pulex: network analysis reveals rhythmic gene expression and enhances functional gene annotation
- Samuel S. C. Rund†1, 2, 7, 8,
- Boyoung Yoo†3, 9,
- Camille Alam2,
- Taryn Green3,
- Melissa T. Stephens5,
- Erliang Zeng3, 5, 10, 11,
- Gary F. George1, 2,
- Aaron D. Sheppard1, 2,
- Giles E. Duffield1, 2,
- Tijana Milenković1, 3, 4 and
- Michael E. Pfrender1, 2, 6Email author
© The Author(s). 2016
Received: 20 April 2016
Accepted: 5 August 2016
Published: 18 August 2016
Marine and freshwater zooplankton exhibit daily rhythmic patterns of behavior and physiology which may be regulated directly by the light:dark (LD) cycle and/or a molecular circadian clock. One of the best-studied zooplankton taxa, the freshwater crustacean Daphnia, has a 24 h diel vertical migration (DVM) behavior whereby the organism travels up and down through the water column daily. DVM plays a critical role in resource tracking and the behavioral avoidance of predators and damaging ultraviolet radiation. However, there is little information at the transcriptional level linking the expression patterns of genes to the rhythmic physiology/behavior of Daphnia.
Here we analyzed genome-wide temporal transcriptional patterns from Daphnia pulex collected over a 44 h time period under a 12:12 LD cycle (diel) conditions using a cosine-fitting algorithm. We used a comprehensive network modeling and analysis approach to identify novel co-regulated rhythmic genes that have similar network topological properties and functional annotations as rhythmic genes identified by the cosine-fitting analyses. Furthermore, we used the network approach to predict with high accuracy novel gene-function associations, thus enhancing current functional annotations available for genes in this ecologically relevant model species. Our results reveal that genes in many functional groupings exhibit 24 h rhythms in their expression patterns under diel conditions. We highlight the rhythmic expression of immunity, oxidative detoxification, and sensory process genes. We discuss differences in the chronobiology of D. pulex from other well-characterized terrestrial arthropods.
This research adds to a growing body of literature suggesting the genetic mechanisms governing rhythmicity in crustaceans may be divergent from other arthropod lineages including insects. Lastly, these results highlight the power of using a network analysis approach to identify differential gene expression and provide novel functional annotation.
As organisms progress through their 24 h day, they experience a daily cycle of environmental changes, including rhythms in temperature, light, predation risk, and resource availability. To respond to these daily environmental changes, organisms modulate their biology in a rhythmic manner driven by both the coordinated action of an endogenous circadian pacemaker or clock, as well as the direct effect of the environmental light:dark cycle (LD cycle) [1–8]. True ‘circadian rhythms’ are those rhythms that can be observed under constant conditions, as opposed to ‘diel’ (or diurnal) rhythms that are observed under a LD cycle.
Examples of time-of-day specific biology observed in other organisms include increasing olfactory sensitivity prior to times of foraging , rhythmic coordination of detoxification enzymes , and time-of-day specific changes in susceptibility to immune challenge [11, 12]. In the disease vector mosquitoes Anopheles gambiae and Aedes aegypti, gene expression data shows that up to 20 % of the transcriptome is regulated in a circadian and/or diel manner [4, 8, 13]. Similarly, under both circadian and/or diel conditions, as much as 5 % of the transcripts expressed in the fruit fly Drosophila melanogaster head and 4 % in the honeybee Apis mellifera head are rhythmically expressed [14–16].
“Preemptive” or “anticipatory” changes in gene regulation and physiology are possible due to the presence of an endogenous cellular molecular circadian clock that provides crucial timing mechanisms in organisms as diverse as cyanobacteria, plants, fungi, insects, and mammals. The circadian clock is cell autonomous and at the molecular level comprises a series of transcriptional-translational feedback loops (TTFLs) whose completion takes approximately 24 h . The genetic underpinning of the metazoan circadian clock is highly conserved . However, even within insects that typically have a conserved genetic construct with regard to their circadian clock, there is some variation in the arrangement of the TTFL components that are found in at least four different functional configurations [28, 29]. As genomic investigations expand into taxonomically diverse arthropods, there is a greater likelihood to find additional variation as well as gain an increased understanding of the genetic basis of rhythmicity.
The investigation of non-model organisms and their genomes is significantly challenged by a lack of functional annotation for many of their protein coding genes. The specific function of many genes, in even well studied organisms, is incomplete. It is not unusual for upward of a third of the genes in any genome to lack functional annotation [30, 31]. This lack of functional annotation is exacerbated in the diverse genomes of many newly sequenced organisms that often contain lineage-specific genes lacking homology with other assembled genomes [19, 31]. The single genome sequence available for the highly diverse lineage of Branchiopod crustaceans is that of the water flea, Daphnia pulex . In excess of one third of the D. pulex genome comprises lineage-specific genes that lack orthologs in other eukaryote genomes and thus lack any functional annotation . These two factors, diverse clocks and unknown genes, present significant challenges in understanding the rhythmic biology of this species.
Here we describe diel transcriptome analysis performed on female D. pulex. Diel as opposed to constant (circadian/ dark:dark) conditions were chosen, as these conditions are more applicable in the context of existing experimental data and indeed natural, real-world environments. This first genome-wide examination of Daphnia transcriptional rhythms allows us to explore a wide breadth of crustacean biochemistry, physiology, and behavior that may be under rhythmic control. In addition to an analysis of cyclic behavior using available statistical software (JTK_CYCLE ), we employed a comprehensive network modeling and analysis approach to validate and discover both orthologous and novel rhythmic gene regulation in D. pulex, as follows.
Established approaches to the statistical analysis of gene expression patterns can provide valuable insights into cellular functioning. However, these analyses often ignore the functional relationships among genes and miss the opportunity to leverage the patterns of co-expression that result from shared regulatory machinery or complex linkage in networks of biological processes [33–36]. For this reason, we modeled Daphnia gene expression data as a gene co-expression network, in which nodes correspond to genes and edges link genes with similar expressions . Using the resulting network data, we performed comprehensive computational analyses to search for genes that share topological properties (in the sense that they occupy similar network positions or cluster together) with the rhythmic genes identified by independent analysis of the expression data using JTK_CYCLE. Our goals were to illuminate the transcriptional signature of rhythmicity in a model crustacean and develop a network-based analytical approach to validate rhythmic genes and infer functional relationships for genes lacking functional annotations.
Results and discussion
Global transcription analysis and JTK_CYCLE analysis
The medial peak-to-trough fold change of genes determined to be rhythmic by JTK_CYCLE was 2.06. However, many rhythmically expressed genes had much higher daily fold-changes in expression (Fig. 2c). While there were genes with peak expression at all times-of-the day, there were a greater number of genes having a peak expression during mid-day (~ZT 4 – ZT 8) and mid-night (~ZT 16–20) than at other times-of-the day (Fig. 2c). We report peak time-of-day expression in Zeitgeber time (ZT), with ZT 0 defined as the time of lights-on and ZT 12 defined as the time of lights-off.
Network-based prediction of rhythmic genes
For each network, we searched in the given network for the subset of the 1,661 genes that have been identified as rhythmic by JTK_CYCLE analysis of the expression data (Additional file 6). Then, we pursued two directions for identifying new rhythmic genes from network topology.
First, we checked if the subset of the 1,661 JTK_CYCLE-identified rhythmic genes that were in the given network were statistically significantly “central” or “peripheral” in the network [36, 40], and if so, we identified other genes with similar network positions and predicted these genes as additional rhythmic genes (see Methods). Here, we used seven popular node centrality, or peripherality, measures, each of which captures somewhat complementary information on the network position of a node [40, 41]: betweenness, closeness, clustering coefficient, degree, eccentricity, graphlet degree, and k-coreness centrality. In this way, we predicted 3,475 genes by at least one combination of network type and centrality measure to be rhythmic, of which 3,019 were not among the 1,661 JTK_CYCLE-identified rhythmic genes and are thus novel. Each prediction was supported by up to 7 combinations of network type and centrality measure (out of the 10 possible combinations; see Methods).
Second, we checked whether the rhythmic genes statistically significantly grouped (i.e., clustered) together in the network, and if so, we found other genes that clustered with them and predicted these genes as additional rhythmic genes (see Methods) [42–44]. Here, we relied on a popular method called Markov clustering algorithm (MCL) . In this way, we predicted 445 genes from at least one network to be rhythmic by the clustering analysis, of which 133 (Additional file 1) were not among the 1,661 JTK_CYCLE-identified rhythmic genes and were thus novel. Each prediction was supported by up to 3 different networks (out of the 3 possible networks, see Methods).
Validation of the novel network-based rhythmic predictions
We provide as comprehensive as possible bioinformatics validation of our novel predictions. First, we validated the novel predictions via a functional enrichment analysis. Namely, we split genes from our networks into three sets: 1) the 626 genes identified as rhythmic using JTK_CYCLE that were in any of our networks (as positive control or “known” rhythmic genes), 2) our novel predictions (either the 3019, 133, or 116 predictions from the centrality analysis, the clustering analysis, or both analyses, respectively), and 3) the remaining “negative control” genes that have not been predicted as rhythmic by any approach in our study. Then, we measured the enrichment of Gene Ontology (GO) terms in each of the three data sets, with the hypothesis that our novel network-based predictions would be involved in the same processes (i.e., GO terms) as the JTK_CYCLE-determined rhythmic genes, but not as the negative control genes; and this was exactly what we observed (p-value of 4.4×10−5 with respect to the hypergeometric test, Fig. 5b, Additional file 7). Second, we further validated our novel predictions by searching for them in an independent set (i.e. set not considered when making the predictions in the first place) of “known” rhythmic genes. Specifically, we calculated the significance of the overlap between our novel network-based predictions and sequence homologs of rhythmic  mosquito genes as “known” rhythmic genes. Fifty homologs of rhythmic mosquito genes were present in at least one of our networks (anywhere in the given network). Our novel predictions captured a statistically significant portion of these 50 genes (with p-values of 0.0028, 0.0129, and 0.0452 for our novel centrality-based, cluster-based, and both centrality- and cluster-based predictions, respectively), unlike the 626 genes identified as rhythmic using JTK_CYCLE (p-value of 0.2591) or the negative control genes (p-value of 0.999).
These encouraging results (genes in our novel network-based prediction set performing similar functions as, or significantly overlapping with “known” rhythmic genes) imply the effectiveness of the network approach. When combined with the statistical analyses (i.e. JTK_CYCLE) of expression data, this approach was able to uncover additional biological knowledge compared to the statistical analyses alone. Of course, in addition to the above bioinformatics validation of our novel predictions, future experimental validation is certainly of interest.
Network-based prediction of novel functional annotations and their validation
We can use the same network approach as above (i.e. the centrality and clustering analyses) to predict novel functional annotations of currently uncharacterized genes based on how well these genes group with functionally annotated genes. In this context, we first validated the accuracy of the network approach by hiding gene-function associations using leave-one-out cross-validation (see Methods) and measuring how accurately we could predict the hidden associations. We found that the prediction accuracy of the network approach is 95 % (p-value of 2.2×10−16 with respect to the hypergeometric test). Given such a high accuracy, we next applied the network approach to currently uncharacterized genes to predict their functional annotations. In this way, we predicted 477 novel gene-function associations spanning 253 uncharacterized genes and 33 GO terms. Also, we predicted an additional 287 novel associations spanning 121 characterized genes and 22 GO terms (see Additional file 1). Hence, we demonstrated that our network approach complements significantly and with high accuracy the currently limited functional annotation data (See Fig. 3a, Additional file 1, and Additional file 2.). In the sections below, we highlight genes of interest from various rhythmic gene families. Network analysis revealed novel rhythmic genes in only three of these highlighted gene families, peroxidases, C-type lectins, and chitinases. The profiles of these identified genes are visualized in Fig. 3b.
Rhythmic expression of immune genes
A pattern of diel vertical migration (DVM), which serves to reduce exposure to predators and track resources, is a well-documented behavior in many populations of Daphnia living in large bodies of water. This DVM typically manifests as a downward migration during the day and an upward migration at night . While there is a fitness benefit in terms of predator avoidance and resource acquisition, this behavior may also cause increased exposure to pathogens during the daytime since many Daphnia parasites tend to live in water body sediments . This observation suggests that there may be an advantage to increasing the transcriptional activity of genes involved in immune processes during the daytime when Daphnia are lower in the water column to counteract the increased time-of-day risk of infection.
There were three paralogs of Toll-like receptors (TLRs) proteins rhythmically expressed with expression peaking in the morning (JGI_V11_3995, ZT4; JGI_V11_40744, ZT4; JGI_V11_190084, ZT 6; Fig. 6a). These proteins serve as PRRs (but are transmembrane not extracellular proteins) and as components in the Toll signaling pathway. This observation indicates that the Toll pathway may be rhythmically upregulated by means of rhythmic expression of GNBPs and TLRs. Similarly, prophenoloxidase activity/melanization may be rhythmic due to the 24 h rhythmic expression of GNBPs and CTLs. Note the GNBP paralog JGI_V11_329548 had an expression profile that did not peak in the morning, but instead in the middle of the night (ZT 18). TEPs and CTLs also promote phagocytosis, so this immune response may also be regulated in a time-of-day specific manner in D. pulex.
The IMD signal transduction pathway in D. pulex also contained highly rhythmic genes like IMD (JGI_V11_313869, ZT 6), Dredd (JGI_V11_45861, ZT 4), and two Relish homologues (JGI_V11_329057, ZT 6; JGI_V11_62213, ZT 7) (Fig. 6a). This result suggested the IMD immune pathway is more sensitive to activation or will mount a more robust response in the morning. In mice, such a phenomenon exists, with expression of Toll-like receptor 9 (TLR9) and subsequent TLR9 mediated innate and adaptive immunity under control of the circadian clock .
Of the 32 chitinases with rhythmic expression, 27 displayed peak expression in the morning (Fig. 6b). Since chitinases can be extracellular they may contribute to rhythmic immunity by directly hydrolyzing the cell walls of chitin-containing pathogens. We note that network analysis revealed an additional two rhythmically expressed chitinases (Fig. 3b).
The rhythmic nature of many immune genes, which generally peak in expression during the daytime (60 of the 69 rhythmic, as determined by JTK_CYCLE, immune genes discussed in this section peak in expression between ZT2 - ZT8), is consistent with the hypothesis that resistance to infection may be higher during the day than the night. Rhythmic levels of extracellular PRRs, chitinases, and IMD signaling pathway components could drive this time-of-day difference. A time-of-day specific resistance to pathogens has been noted in Drosophila and mice [11, 12, 50], but remains to be investigated in Daphnia.
Metabolic genes, vesicular-type ATPase subunits, and tRNA synthetases are constitutively expressed
It is plausible that D. pulex primarily regulates metabolic gene expression to match available nutrient reserves and/or ambient temperatures, instead of anticipating future nutrient levels. Under the resource-rich, temperature controlled, and oxygen-stable laboratory conditions, it would therefore not be surprising that we did not detect many rhythmic metabolic genes. Rearing Daphnia under cycling temperature conditions or food resource levels in the laboratory will be required to understand if rhythmic metabolic levels indeed exist under cycling environmental condition.
All annotated tRNA synthetases (ligases) displayed constitutive expression (Fig. 7) with the exception of the tRNA synthetase JGI_V11_187913 (ZT 18) (Additional file 1). Rhythms in aminoacyl-tRNA synthetases would suggest an organism has increased protein synthesis activity at particular times-of-day and could indicate there may be rhythmic control at the translational level which produces, enhances, or modifies 24 h rhythms downstream of gene expression. In mosquitoes, 12 of these aminoacyl-tRNA synthetases are rhythmic, and expressed in a similar phase , yet in D. pulex, we find that expression of annotated tRNA synthetases was primarily constitutive, indicating it is unlikely that there are 24 h rhythms in translation. However, there are various other steps involved in protein translation, post-translational modification, and degradation that could still provide time-of-day regulation and modification of protein levels.
Daphnids may also adjust their chemoperception abilities/sensitivities in a time-of-day specific manner, for example to respond to changing needs for predator avoidance. Among the gustatory receptors (Grs) identified by Peñalva-Arana et al. , we found three that were rhythmically expressed (Fig. 8): Gr12 (JGI_V11_327171) and Gr43 (JGI_V11_327170), both peaking at ZT 6; and Gr36 (JGI_V11_329588) which peaks late at night (ZT 22). Gustatory receptors can be very specific in the odorants they detect . It is plausible that D. pulex shows time-of-day specific sensitivity to certain odorants, and not others. A similar olfactory phenomena has been observed in the mosquito, which is more sensitive to certain human host odorants during normal host-seeking times-of-day . D. pulex could, for example, upregulate its ability to sense predator kairomones at the time-of-day when predation risk is highest.
Oxidative stress detoxification
It has been reported that Drosophila has daily rhythms in resistance to oxidative stress . We examined D. pulex genes with known or putative roles in oxidative stress detoxification and found that many were rhythmically expressed (Fig. 8). Catalase had rhythmic expression peaking at mid-day (cat, JGI_V11_308727, ZT 6; Fig. 8), consistent with the peak time of catalase activity described in D. longispina . Four genes encoding putative peroxidases (JGI_V11_304755, ZT 8; JGI_V11_307875, ZT 4; JGI_V11_206041, ZT 10; ZT_V11_314049, ZT 10; Fig. 8) were also rhythmically expressed and peak during the daytime (Fig. 8). Interestingly, sodium oxidase dismutase 1 expression was weakly rhythmic, but with expression peaking much later (SOD1, JGI_V11_329848, ZT20; Additional file 1). We note that network analysis revealed an additional two rhythmically expressed peroxidases (Fig. 3b). These data suggest that response to oxidative stress may be time-of-day specific. During the daytime, when sunlight-induced oxidative stress would be most likely to occur, D. pulex may be more prepared to respond to oxidative stress.
We next looked at D. pulex homologues of arthropod core molecular clock genes primarily as annotated by Tilden et al. . These ‘canonical’ clock genes make up the transcriptional-translation feedback loops (TTFLs) that encompass the molecular circadian clock. TTFLs have been characterized in a variety of organisms including terrestrial arthropods, rodents, and humans [28, 62]. However, the molecular clock has yet to be characterized in Daphnia. Orthologs of the core molecular clock mechanism in Drosophila are found in the D. pulex genome , suggesting that the elements of a conserved arthropod clock are present in Daphnia. Tilden et al.  identified well-characterized circadian clock genes including clock, cycle, period, PAR domain protein 1 ε, vrille, eight timeless paralogs, and both cryptochrome 1 (a Drosophila-like photoreceptor) and cryptochrome 2 (a mouse-like transcriptional repressor) . Finding both forms of cryptochrome suggests that the D. pulex has an endogenous circadian molecular clock that is more similar to that of the monarch butterfly Danaus plexippus and the mosquito An. gambiae than to Drosophila (which lacks the cryptochrome 2 transcriptional repressor) [8, 61, 63, 64].
Based on the observation that a full complement of canonical clock gene orthologues were expressed in D. pulex, and based on current insect clock models, including those of Drosophila and the mosquito, we would predict 24 h sinusoidal rhythmic expression profiles for several of these genes in cells throughout the organism [8, 13–15, 62–66]. This prediction includes at least one, if not both, of the two positive loop components of the TTFL mechanism, clock and cycle, the three negative loop components, period, timeless and cryptochrome 2, and the members of the interlocking loop, pdp1 and vrille.
There are a number of possible explanations for the lack of robust clock gene rhythmicity detected in this experiment. First, there is evidence that specific tissues of crayfish express clock genes in anti-phase of one another; and in mice the phase of clock gene expression can vary between tissues by as much as 4–8 h [67, 68]. As we assayed entire organisms, it may be that differences in peak phases between tissues reduced the amplitude of the oscillation or appeared as multiple peaks of expression; giving the appearance of a lack of rhythmicity. Second, clock gene rhythmicity may be limited only to a subset of tissues and thus would not be detected in a whole-organism assay. For example, rhythmic clock gene expression was not detected in all peripheral tissues of the Zebrafish .
It is also feasible that D. pulex has an alternative, non-canonical, core molecular clock that operates differently from well-characterized insect clocks. In fact, under constant light conditions, D. magna has been reported to have an unusually long 28 h free running period [23, 70]. Alternate clock mechanisms have been suggested in a number of invertebrates including the nematode, Caenorhabditis elegans , the sea squirt, Ciona intestinalis , and in other crustaceans like the speckled sea louse, Eurydice pulchra , and the prawn, Macrobrachium rosenbergii .
In summary, under our laboratory diel conditions and while assaying whole organism transcriptomic data, robust, high fidelity, high amplitude, 24 h sinusoidal expression of Daphnia canonical clock genes was not detected. We suggest several possible explanations, but we cannot distinguish among these alternatives. Future experiments to determine the mechanism of Daphnia circadian clocks should examine tissue-specific expression patterns of clock genes under both diel and constant dark conditions.
Comparison with condition-dependent gene regulation
Comparisons between condition-dependent and temporal regulation of gene expression
Thermal regime a
18 °C vs. 28 °C
Low tolerance genotype C
11 % / 183
Low tolerance genotype E
12 % / 199
High tolerance genotype B
7 % / 116
High tolerance genotype K
7 % / 116
Salinity stress b
5 gL−1 NaCl vs. control
High tolerance genotype
38 % / 631
Low tolerance genotype
33 % / 548
Carbon:Phosphorous ratio c
HiC:LoP (C:P ~800) vs. LoC:HiP (C:P ~100) (3 days)
G1: LoC:HiP tolerance
4 % / 66
G2: HiC:LoP tolerance
8 % / 133
HiC:LoP (C:P ~800) vs. LoC:HiP (C:P ~100) (6 h)
G1: LoC:HiP tolerance
16 % / 266
G2: HiC:LoP tolerance
9 % / 149
We present genome-wide diel transcriptional profiling of Daphnia under controlled laboratory conditions. Our analysis revealed that Daphnia express a great number and diversity of genes in a highly rhythmic manner, including genes involved in sensory processes, response to oxidative stress, and immune-related genes. It is advantageous for an organism to upregulate its ability to combat or avoid environmental stressors at the times-of-day they are most likely to encounter such risks. We highlighted how during the daytime, Daphnia populations may be at greater risk of oxidative damage (from the Sun), predation (from fish, which locate them visually during the day), and exposure to pathogens (which they encounter at greater frequency at the bottom of the water column in Daphnia populations that migrate downwards during the day-time). It therefore should not be surprising that genes that function in detoxifying oxidative stress, chemoperception, and combating pathogens are rhythmically expressed.
We showed that while D. pulex has a full complement of the core clock genes found in other animals, the expected high amplitude rhythmic expression of these genes is not apparent when examined at the whole animal level. Work in other crustaceans suggests that clock rhythmicity may be limited to a subset of tissues or has anti-phasic expression in various tissues that would obfuscate the expression patterns of these genes in a whole-organism assay. It is also possible that D. pulex has an alternative molecular clock. Our work contributes to a growing understanding that there is less conservation than initially reported in the mechanisms of arthropod clocks and that more investigation is prudent in elucidating the core clock of Daphnia.
Finally, we developed a comprehensive network modeling and analysis approach that complements and enhances the standard statistical analyses of differential expression (e.g., JTK_CYCLE). This network analysis, when built on top of the statistical analyses, revealed additional knowledge about rhythmicity that could not be captured by the statistical analyses alone. Importantly, we demonstrated the usefulness of the network approach to identify novel functional annotations for currently uncharacterized genes. As the number of available arthropod genome sequence assemblies increases, and the taxonomic breadth of these genomes widens, methodological approaches that take into account the patterns of co-expression and network configuration will become vital additions to the process of functional annotation.
Future work may include examining rhythms at the protein and functional level and examining rhythms in specific tissues (especially in regards to clock gene expression). Further, because other Zeitgebers (time-givers), such as temperature rhythms, may enhance observed rhythms , their role in governing Daphnia rhythmicity should also be investigated.
Experimental material, RNA collection, extraction, and microarray hybridization
A D. pulex genotype was collected from a vernal pond in southwestern Michigan (42.31979 N 85.35837 W) on May 1, 2011. This pond, referred to as “Roughwood”, has a maximum depth of <2.0 M. Single adult females were isolated in the lab to initiate isoclonal cultures and these lines were maintained through parthenogenetic reproduction. The experimental time course was performed in August of 2011. Stock cultures of this genotype (Roughwood 40-11) were maintained for three generations at 18 °C in 1 quart jars containing COMBO water medium  and a 12 L:12D photoperiod (with abrupt transitions). The media, containing resource in the form of algae Scenedesmus acutus culture at a concentration of 200,000 cells/ml, was replaced every second day to maintain a constant level of resource and individuals were culled to maintain a density of 50 individuals per culture. To establish the experimental generation, female offspring from the second or third clutch of second-generation females were collected within 48 h of birth, pooled, and divided into twelve 1 quart jars containing 750 ml of COMBO media and algae and maintained as above with the same feeding regime and 12 L:12D photoperiod. When these third generation females were 10–12 day old the experiment began and replicate pools of 30 animals per jar were collected, rinsed in distilled water, and snap frozen on dry ice every 4 h over a 44 h period.
To assess patterns of gene expression, we used the D. pulex Expression Array 12x135k GEO Accession GPL11278 (See Colbourne et al. ). The platform is a high-density NimbleGen (Roche-NimbleGen, Madison, WI) microarray that accommodates 12 experiments per glass slide, with each experiment interrogating 137,000 probes. Each predicted and experimentally validated gene is represented by as many as three probes, whereas the remaining probes are designed from transcriptionally active regions whose gene models are not yet known.
Total RNA was purified from a single pool of animals collected at each time point using TRIzol reagent (Invitrogen, Carlsbad, CA) extraction followed by RNA purification using a Qiagen RNeasy Mini Kit with on-column DNAse treatment to isolate total RNA. RNAsecure (Ambion, Austin, TX) was used after RNA purification to inactivate any remaining RNAases. The quality and quantity of resulting RNA were assessed using a Nano Drop (Thermo Scientific, Waltham, MA) and RNA Nano LabChip for the Bioanalyzer (Agilent Technologies, Santa Clara, CA). Beginning with 1.0 μg of total RNA, a single round of amplification using MessageAmpTM II aRNA kit (Ambion) was performed for each RNA sample. The cRNA (10 μg) was converted to double strand cDNA with random primers using the Invitrogen SuperScript Double-Stranded cDNA Synthesis Kit (Invitrogen). From 1 μg double-stranded cDNA, labeled cDNA was generated with NimbleGen’s Dual-colour Labeling Kit (Roche NimbleGen).
RNA from each pooled sample was used to create four technical hybridization replicates. Dual-colour hybridization, post-hybridization washing, and scanning were done according to Roche NimbleGen protocols . Images were acquired using a NimbleGen MS 200 with 2 μm resolution, and GenePix 6.0 software (Molecular Devices, Sunnyvale, CA). The image data from these arrays were processed using NimbleScan v2.5 software (Roche NimbleGen) to extract probe intensity values. These data have been deposited in the NCBI Gene Expression Ominibus and are accessible through GEO accession number GSE67781.
Microarray data analysis
Gene expression values (i.e., gene intensity values) were obtained from a summarization of intensity values of all corresponding probes using the RMA (Robust Multi-array Average) method . The pre-processed microarray data were imported into an in-house analysis pipeline using Bioconductor  for normalization and analysis. All genes were quantile-normalized across arrays, samples, and technical replicates . Unless otherwise noted, data presented is normalized, but not log transformed.
JTK_CYCLE analysis of rhythmic gene expression
To determine the patterns of rhythmic gene regulation we first used the JTK_CYCLE algorithm as previously described [32, 37, 38]. JTK_CYCLE is a nonparametric statistical algorithm designed to identify and characterize cycling variables in large datasets. It applies the Jonckheere-Terpstra-Kendall (JT) test and Kendall’s tau (rank correlation), to find the optimal combination of period and phase that minimize the p-value of Kendall’s tau correction between the experimental time series and each tested cyclical ordering, this being derived from cosine curves. JTK_CYCLE generates period length, phase, and amplitude estimates, as well as corrects for multiple comparisons post hoc. The reported q-value takes into consideration the false discovery rate (FDR) across all genes .
Prior to running JTK_CYCLE, genome features on the microarray not mapping to a gene in the D. pulex v1.1 frozen gene set were excluded, as were all probes that failed to exceed background levels in at least one sample. A background fluorescence cutoff value of 136.5 was used as it was a value that excluded >99.8 % of random probes on the arrays. Genes were scored as being rhythmic if they were found using the JTK_CYCLE algorithm to have a q < 0.1 (a commonly used JTK_CYCLE cutoff) and period length of 22–26 h.
Gene annotations followed the 2013–2014 gene annotations from the Department of Energy – Joint Genome Institute Daphnia Genome Browser  supported by gene annotations from OrthoDb  and euGenes Arthropod 2009.12 annotations . Further gene and gene family identifications were taken from Rivera et al. , Colbourne et al. , McTaggart et al. , Peñalva-Arana et al. , KEGG Pathways Database , and from 2013 to 2014 Drosophila melanogaster homologs listed on FlyBase . For each gene family mentioned in this paper (e.g., chitinases, putative clock genes, glycolysis, etc.), the identities of all genes considered for that family are provided in Additional file 1.
Hierarchical cluster analysis
We performed hierarchical cluster analysis (Fig. 2) using Cluster 3.0 and visualized the results using Java TreeView. Data were mean centered, normalized across the time course for each gene, and centroid linkage clustering performed [87, 88].
Real-time quantitative RT-PCR (qRT-PCR) analysis
Total RNA (as above) was used for cDNA synthesis using a High Capacity cDNA reverse transcriptase kit (Applied Biosystems, Foster City, CA) primed with random hexamers. PCR thermocycling and qRT-PCR were performed using SYBR green reagents per manufacture’s protocol using an Applied Biosystems 7500 Fast Real-Time PCR System. Quantification was based on the generation of standard curves. Dissociation curves to test for primer dimers were generated using dissociation curve software (Applied Biosystems). Normalization of genes was calculated relative to a D. pulex gene deemed constitutively expressed over the 24 h, Alpha tubulin (JGI_V11_301837). See Additional file 3 for primer sequences used.
CircWave v1.4 software, an extension of cosinor analysis, was also used to analyze clock gene expression rhythmicity. This analysis fits a Fourier curve (one sine wave) to the data. The p values reported are the result of F-tests from the software .
We modeled temporal gene expression data with a gene co-expression network as follows: in the network, nodes are genes, and two nodes are connected by an edge if the corresponding genes show similar expression patterns over time. To measure similarity between expression profiles of two genes, we used one of the following three popular edge weight methods: signed Pearson correlation, absolute Pearson correlation, and mutual information, see Rider et al.  for details. Each edge weight method assigned to every pair of genes a score that captures some intuition of how well the genes’ expression levels “correlate” over time; then, we constructed a co-expression network by predicting the top K highest-scoring node pairs as edges in the network. We varied K from N (where N is the number of genes in the expression data) to 10 N in increments of N, and we also tested K = 25 N, 50 N, and 100 N, in order to evaluate the effect of the value of parameter K on the resulting network structure. Intuitively, we wanted to choose this value in a way that keeps only significant edges and provides a meaningful representation as well as interpretation of the data . Namely, for each edge weight method, we aimed to construct a network that ideally linked all genes (i.e., has no isolated nodes), in order to include into the network as much information from the data as possible. At the same time, we aimed to construct a network that is not too dense (where density is defined as the number of edges in the network out of all possible edges), in order to mimic the sparse nature of many real-world networks as well as avoid randomness in network topology [90, 91]. Empirically, by studying the number of non-isolated nodes, edges, connected components, and nodes and edges in the largest connected component, we found that as K increases, the number of non-isolated nodes barely increases, while the network density increases drastically. Therefore, since larger values of K increase network density without introducing many new nodes into the networks, which could unnecessarily increase computation complexity of network analysis methods or include “lower-confidence” interactions, we empirically decided to focus on the following networks (i.e., their largest connected components): signed Pearson correlation with top N interactions, absolute Pearson correlation with top N interactions, and mutual information with top N interactions. Also, since edges that are captured by both Pearson correlation as well as mutual information could be of higher confidence, we also studied (largest connected components of) the intersections of absolute Pearson correlation and mutual information with: 1) top 10 N interactions and 2) top 25 N interactions. Thus, we studied a total of five networks.
Network-based prediction of new rhythmic genes – node centrality-based analysis
For a given network and a given centrality measure, i.e., for each of 5×7 = 35 possible cases (where 5 corresponds to the number of studied networks and 7 corresponds to the number of studied centralities), in 11 cases we observed statistically significant differences (p-values below 0.01) between network positions of the rhythmic genes and the remaining genes. Next, for each of the 11 cases, we ranked all genes in the network from the most central to the least central (or vice versa). We took the top K% most (or least) central genes, where we varied K from 1 to 100 % in increments of 1 %. Then, we computed precision (the portion of the top K% genes that are JTK_CYCLE-identified rhythmic genes), recall (the portion of JTK_CYCLE-identified rhythmic genes from the given network that are among the top K% genes), and F-score (which combines precision and recall into a single value that is easier to interpret than the two individual and typically “contradictory” measures). At the value of K where F-score peaks (meaning that the methodology achieves the highest prediction accuracy at that value of K), we measured the enrichment of the top K% genes in the JTK_CYCLE-identified rhythmic genes (with respect to the hypergeometric test). In 10 out of the 35 possible cases, we observed statistically significant enrichments signal (with p-values below 0.05 after Bonferroni correction for multiple hypothesis testing, which was done using p.adjust package in R). In each of the 10 cases, we took the remaining (non-JTK_CYCLE-identified) genes among the top (or bottom) K% genes and predicted them as network-based rhythmic candidates. We recorded how many of the 10 combinations of network type and centrality measure support each prediction, because the more the combinations, the higher the prediction confidence.
Network-based prediction of new rhythmic genes – network clustering-based analysis
We clustered each network with Markov clustering algorithm (MCL) (with the inflation parameter set to 2, because this value gives empirically optimal cluster size distribution, i.e., not too many very small clusters or too few very large clusters). Then, focusing on meaningful clusters (of size at least two and containing at least two JTK_CYCLE-identified rhythmic genes), for each value of K in 1–100 % range, we found all clusters that have enrichment in the JTK_CYCLE-identified rhythmic genes greater than K% (i.e., clusters in which at least K% of the genes are JTK_CYCLE-identified rhythmic genes), and we measured the statistical significance of the enrichment via the hypergeometric test (at a p-value threshold of 0.05 after Bonferroni correction for multiple hypothesis testing). For each of the resulting significant clusters, we predicted the genes in the cluster as network-based rhythmic genes. We took all predictions from all significant clusters and computed the overall prediction accuracy via precision, recall, and F-score measures. Further, in order to ensure that we could not achieve the same prediction accuracy by chance, we randomly clustered the network 100 times and repeated the above steps to compute “randomized” precision, recall, and F-score. Then, we used these randomized results to determine “optimal” K at which to make predictions from the real clusters. Namely, we selected the value of K at which F-score peaks, but only if that F-score value is statistically significantly high compared to the randomized F-score values. For three of our five networks (absolute Pearson correlation with top N interactions, intersection of absolute Pearson correlation and mutual information with top 10 N interactions, and intersection of absolute Pearson correlation and mutual information with top 25 N interactions), results are statistically significantly better from the actual clusters than from the randomized clusters, and we predicted new rhythmic genes from these three networks (and we left out the other two networks from consideration when making predictions). We recorded how many of the three networks support each prediction, as the higher the number of networks, the higher the prediction confidence.
Network-based prediction of novel functional annotations and their validation
To validate the network approach in the context of predicting novel gene-function associations, we used leave-one-out cross-validation, as follows. When considering all genes in any one network, we hid functional information for one gene of interest at a time. Then, we used functional information about all other genes that group with the gene of interest (with respect to either centrality of clustering analysis) to potentially predict function(s) of the gene of interest. Namely, if the given gene group is statistically significantly enriched in a given function (with p-value below 0.05 after Bonferroni correction for multiple hypothesis testing), we annotated the entire gene group with that function, including our gene of interest. We repeated this procedure for all genes and got the set of predicted gene-function associations for each of centrality and clustering analysis. For increased confidence, we considered association predictions that are in the intersection of the two analyses. For the resulting predictions, we measured their accuracy via precision, which is the percentage of associations that we predicted via the network approach which are in the known functional annotation data. Novel predictions are provided in Additional file 1.
CHT, chitinase; CTL, C-type lectin; DVM, diel vertical migration; FA, fatty acid; GNBP, gram-negative bacteria-binding protein; GO, gene ontology; Gr, gustatory receptor; IMD, immune deficiency; LD, light:dark; MCL, Markov clustering algorithm; PRR, pattern recognition receptor; SIGN N, signed pearson correlation with top N interactions; TCA the citric acid; TEP, thioester-containing protein; TLR, toll-like receptor; TTFL, transcriptional-translational feedback loop; V-ATPase, vesicular-type ATPase
We thank Sayanty Roy for her assistance in the maintenance of experimental stocks; Michael Hughes for provision of and assistance with the JTK_CYCLE algorithm; Enrique Blair and Cameron Houk for their assistance with data processing, Jacqueline Lopez for her assistance with our GEO submission; and we are especially thankful to Tom Little and Seanna McTaggart for sharing their expertise.
S.S.C.R. is funded by a strategic award from the Wellcome Trust for the Centre for Immunity, Infection and Evolution (no 095831) and a Royal Society Newton International Fellowship (NF140517); T.M. by the National Science Foundation (NSF) CAREER CCF-1452795 and CCF-1319469 grants; and G.E.D. by the Eck Institute for Global Health and the NIGMS (R01-GM087508). This project was supported by National Institutes of Health (NIH) grant R24-GM078274 to M.E.P. The funding bodies did not have a role in the design of the study, data collection, analysis, interpretation of data, writing the manuscript, nor the decision to publish.
Availability of data and materials
Microarray expression data has been deposited on NCBI Gene Expression Ominibus and are accessible through GEO accession number GSE67781.
S.S.C.R., C.D.C.A., G.E.D., and M.E.P. collected and prepared the biological material; M.T.S., E.Z., and C.D.C.A. generated and processed the microarray data. C.D.C.A., S.S.C.R., G.E.D., and G.F.G. generated and processed the qRT-PCR data. S.S.C.R performed the JTK_CYCLE analysis. T.G. and B.Y. performed all network-based analyses of the data. T.M. designed and supervised all aspects of the network-based analyses. G.E.D., T.M., and M.E.P. conceived and designed the study. All authors contributed to the analysis of data and writing of the manuscript and approved the final manuscript.
The authors declare they have no competing interests.
Consent for publication
Ethics approval and consent for participation
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Ceriani MF, Hogenesch JB, Yanovsky M, Panda S, Straume M, Kay SA. Genome-wide expression analysis in Drosophila reveals genes controlling circadian behavior. J Neurosci. 2002;22:9305–19.PubMedGoogle Scholar
- Lin Y, Han M, Shimada B, Wang L, Gibler TM, Amarakone A, Awad TA, Stormo GD, Van Gelder RN, Taghert PH. Influence of the period-dependent circadian clock on diurnal, circadian, and aperiodic gene expression in Drosophila melanogaster. Proc Natl Acad Sci U S A. 2002;99:9562–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Albers HE, Liou SY, Ferris CF, Stopa EG, Zoeller RT. Neurochemistry of circadian timing. In: Klein DC, Moore RY, Reppert SM, editors. Suprachiasmatic nucleus: The mind's clock. New York: Oxford University Press; 1991. p. 263–88.Google Scholar
- Rund SSC, Gentile JE, Duffield GE. Extensive circadian and light regulation of the transcriptome in the malaria mosquito Anopheles gambiae. BMC Genomics. 2013;14:218.View ArticlePubMedPubMed CentralGoogle Scholar
- Michael TP, Mockler TC, Breton G, McEntee C, Byer A, Trout JD, Hazen SP, Shen R, Priest HD, Sullivan CM, Givan SA, Yanovsky M, Hong F, Kay SA, Chory J. Network discovery pipeline elucidates conserved time-of-day-specific cis-regulatory modules. PLoS Genet. 2008;4:e14.View ArticlePubMedPubMed CentralGoogle Scholar
- Hastings JW. Biochemical aspects of rhythms - phase shifting by chemicals. Cold Spring Harb Symp Quant Biol. 1960;25:131–43.View ArticlePubMedGoogle Scholar
- Wijnen H, Naef F, Boothroyd C, Claridge-Chang A, Young MW. Control of daily transcript oscillations in Drosophila by light and the circadian clock. PLoS Genet. 2006;2:e39.View ArticlePubMedPubMed CentralGoogle Scholar
- Rund SSC, Hou TY, Ward SM, Collins FH, Duffield GE. Genome-wide profiling of diel and circadian gene expression in the malaria vector Anopheles gambiae. Proc Natl Acad Sci U S A. 2011;108:E421–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Rund SSC, Bonar NA, Champion MC, Ghazi JP, Houk CH, Leming MT, Syed Z, Duffield GE. Daily rhythms in antennal protein and olfactory sensitivity in the malaria mosquito Anopheles gambiae. Sci Rep. 2013;3:2494.View ArticlePubMedPubMed CentralGoogle Scholar
- Balmert NJ, Rund SSC, Ghazi JP, Zhou P, Duffield GE. Time-of-day specific changes in metabolic detoxification and insecticide resistance in the malaria mosquito Anopheles gambiae. J Insect Physiol. 2014;64:30–9.View ArticlePubMedGoogle Scholar
- Lee JE, Edery I. Circadian regulation in the ability of Drosophila to combat pathogenic infections. Curr Biol. 2008;18:195–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Stone EF, Fulton BO, Ayres JS, Pham LN, Ziauddin J, Shirasu-Hiza MM. The circadian clock protein timeless regulates phagocytosis of bacteria in Drosophila. PLoS Pathog. 2012;8:e1002445.View ArticlePubMedPubMed CentralGoogle Scholar
- Leming MT, Rund SSC, Behura SK, Duffield GE, O'Tousa JE. A database of circadian and diel rhythmic gene expression in the yellow fever mosquito Aedes aegypti. BMC Genomics. 2014;15:1128.View ArticlePubMedPubMed CentralGoogle Scholar
- Ueda HR, Matsumoto A, Kawamura M, Iino M, Tanimura T, Hashimoto S. Genome-wide transcriptional orchestration of circadian rhythms in Drosophila. J Biol Chem. 2002;277:14048–52.View ArticlePubMedGoogle Scholar
- Rodriguez-Zas SL, Southey BR, Shemesh Y, Rubin EB, Cohen M, Robinson GE, Bloch G. Microarray analysis of natural socially regulated plasticity in circadian rhythms of honey bees. J Biol Rhythms. 2012;27:12–24.View ArticlePubMedPubMed CentralGoogle Scholar
- Duffield GE. DNA microarray analyses of circadian timing: The genomic basis of biological time. J Neuroendocrinol. 2003;15:991–1002.View ArticlePubMedGoogle Scholar
- Miner BE, De Meester L, Pfrender ME, Lampert W, Hairston Jr NG. Linking genes to communities and ecosystems: Daphnia as an ecogenomic model. Proc Biol Sci. 2012;279:1873–82.View ArticlePubMedPubMed CentralGoogle Scholar
- Ebert D. Genomics. A genome for the environment. Science. 2011;331:539–40.View ArticlePubMedGoogle Scholar
- Colbourne JK, Pfrender ME, Gilbert D, Thomas WK, Tucker A, Oakley TH, Tokishita S, Aerts A, Arnold GJ, Basu MK, Bauer DJ, Caceres CE, Carmel L, Casola C, Choi JH, Detter JC, Dong Q, Dusheyko S, Eads BD, Frohlich T, Geiler-Samerotte KA, Gerlach D, Hatcher P, Jogdeo S, Krijgsveld J, Kriventseva EV, Kultz D, Laforsch C, Lindquist E, Lopez J, Manak JR, Muller J, Pangilinan J, Patwardhan RP, Pitluck S, Pritham EJ, Rechtsteiner A, Rho M, Rogozin IB, Sakarya O, Salamov A, Schaack S, Shapiro H, Shiga Y, Skalitzky C, Smith Z, Souvorov A, Sung W, Tang Z, Tsuchiya D, Tu H, Vos H, Wang M, Wolf YI, Yamagata H, Yamada T, Ye Y, Shaw JR, Andrews J, Crease TJ, Tang H, Lucas SM, Robertson HM, Bork P, Koonin EV, Zdobnov EM, Grigoriev IV, Lynch M, Boore JL. The ecoresponsive genome of Daphnia pulex. Science. 2011;331:555–61.View ArticlePubMedPubMed CentralGoogle Scholar
- Steams SC. Light responses of Daphnia pulex. Limnol Oceanogr. 1975;20:564–70.View ArticleGoogle Scholar
- Starkweather PL. Diel variation in feeding behavior of Daphnia pulex. Influences of food density and nutritional history on mandibular activity. Limnol Oceanogr. 1978;23:307–17.View ArticleGoogle Scholar
- Haney JF, Hall DJ. Diel vertical migration and filter-feeding activities of Daphnia. Arch Hydrobiol. 1975;75:413–41.Google Scholar
- Ringelberg J, Servaas H. A circadian rhythm in Daphnia magna. Oecologia. 1971;6:289–92.View ArticleGoogle Scholar
- Rhode SC, Pawlowski M, Tollrian R. The impact of ultraviolet radiation on the vertical distribution of zooplankton of the genus Daphnia. Nature. 2001;412:69–72.View ArticlePubMedGoogle Scholar
- Loose CJ. Lack of endogenous rhythmicity in Daphnia diel vertical migration. Limnol Oceanogr. 1993;38:1837–41.View ArticleGoogle Scholar
- Dunlap JC, Loros JJ, Decoursey PJ. Chronobiology: Biological timekeeping. Sunderland Mass: Sinauer Associates; 2004.Google Scholar
- Wager-Smith K, Kay SA. Circadian rhythm genetics: from flies to mice to humans. Nat Genet. 2000;26:23–7.View ArticlePubMedGoogle Scholar
- Yuan Q, Metterville D, Briscoe AD, Reppert SM. Insect cryptochromes: Gene duplication and loss define diverse ways to construct insect circadian clocks. Mol Biol Evol. 2007;24:948–55.View ArticlePubMedGoogle Scholar
- Rubin EB, Shemesh Y, Cohen M, Elgavish S, Robertson HM, Bloch G. Molecular and phylogenetic analyses reveal mammalian-like clockwork in the honey bee (Apis mellifera) and shed new light on the molecular evolution of the circadian clock. Genome Res. 2006;16:1352–65.View ArticlePubMedPubMed CentralGoogle Scholar
- Adams M. Open questions: genomics and how far we haven't come. BMC Biol. 2013;11:109.View ArticlePubMedPubMed CentralGoogle Scholar
- Anton BP, Kasif S, Roberts RJ, Steffen M. Objective: biochemical function. Front Genet. 2014;5:210.View ArticlePubMedPubMed CentralGoogle Scholar
- Hughes ME, Hogenesch JB, Kornacker K. JTK_CYCLE: An efficient nonparametric algorithm for detecting rhythmic components in genome-scale data sets. J Biol Rhythms. 2010;25:372–80.View ArticlePubMedPubMed CentralGoogle Scholar
- Sharan R, Ulitsky I, Shamir R. Network‐based prediction of protein function. Mol Syst Biol. 2007;3:88.View ArticlePubMedPubMed CentralGoogle Scholar
- Ideker T, Sharan R. Protein networks in disease. Genome Res. 2008;18:644–52.View ArticlePubMedPubMed CentralGoogle Scholar
- Rider AK, Milenković T, Siwo GH, Pinapati RS, Emrich SJ, Ferdig MT, Chawla NV. Networks’ characteristics are important for systems biology. Network Sci. 2014;2:139–61.View ArticleGoogle Scholar
- Faisal FE, Milenkovic T. Dynamic networks reveal key players in aging. Bioinformatics. 2014;30:1721–9.View ArticlePubMedGoogle Scholar
- Hughes ME, Grant GR, Paquin C, Qian J, Nitabach MN. Deep sequencing the circadian and diurnal transcriptome of Drosophila brain. Genome Res. 2012;22:1266–81.View ArticlePubMedPubMed CentralGoogle Scholar
- Xu K, DiAngelo JR, Hughes ME, Hogenesch JB, Sehgal A. The circadian clock interacts with metabolic physiology to influence reproductive fitness. Cell Metab. 2011;13:639–54.View ArticlePubMedPubMed CentralGoogle Scholar
- Milenković T, Lai J, Pržulj N. GraphCrunch: a tool for large network analyses. BMC Bioinformatics. 2008;9:70.View ArticlePubMedPubMed CentralGoogle Scholar
- Milenković T, Memišević V, Bonato A, Pržulj N. Dominating biological networks. PLoS One. 2011;6:e23016.View ArticlePubMedPubMed CentralGoogle Scholar
- Meng L, Hulovatyy Y, Striegel A, Milenković T. On the interplay between individuals’ evolving interaction patterns and traits in dynamic multiplex social networks. IEEE Trans Netw Sci Eng. 2016;3:32-43.Google Scholar
- Solava RW, Michaels RP, Milenković T. Graphlet-based edge clustering reveals pathogen-interacting proteins. Bioinformatics. 2012;28:i480–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Milenković T, Memišević V, Ganesan AK, Pržulj N. Systems-level cancer gene identification from protein interaction network topology applied to melanogenesis-related functional genomics data. J R Soc Interface. 2010;7:423–37.View ArticlePubMedGoogle Scholar
- Milenković T, Pržulj N. Uncovering biological network function via graphlet degree signatures. Cancer Inform. 2008;6:257–73.PubMedGoogle Scholar
- Van Dongen S. Graph clustering via a discrete uncoupling process. SIAM J Matrix Anal Appl. 2008;30:121–41.View ArticleGoogle Scholar
- Decaestecker E, De Meester L, Ebert D. In deep trouble: habitat selection constrained by multiple enemies in zooplankton. Proc Natl Acad Sci U S A. 2002;99:5481–5.View ArticlePubMedPubMed CentralGoogle Scholar
- McTaggart SJ, Conlon C, Colbourne JK, Blaxter ML, Little TJ. The components of the Daphnia pulex immune system as revealed by complete genome sequencing. BMC Genomics. 2009;10:175.View ArticlePubMedPubMed CentralGoogle Scholar
- Waterhouse RM, Tegenfeldt F, Li J, Zdobnov EM, Kriventseva EV. OrthoDB: a hierarchical catalog of animal, fungal and bacterial orthologs. Nucleic Acids Res. 2013;41:D358–65.View ArticlePubMedGoogle Scholar
- dos Santos G, Schroeder AJ, Goodman JL, Strelets VB, Crosby MA, Thurmond J, Emmert DB, Gelbart WM, FlyBase Consortium. FlyBase: introduction of the Drosophila melanogaster Release 6 reference genome assembly and large-scale migration of genome annotations. Nucleic Acids Res. 2015;43:D690–7.View ArticlePubMedGoogle Scholar
- Silver AC, Arjona A, Walker WE, Fikrig E. The circadian clock controls toll-like receptor 9-mediated innate and adaptive immunity. Immunity. 2012;36:251–61.View ArticlePubMedPubMed CentralGoogle Scholar
- Panda S, Antoch MP, Miller BH, Su AI, Schook AB, Straume M, Schultz PG, Kay SA, Takahashi JS, Hogenesch JB. Coordinated transcription of key pathways in the mouse by the circadian clock. Cell. 2002;109:307–20.View ArticlePubMedGoogle Scholar
- Kanehisa Laboratories. KEGG Pathway Database. [http://www.genome.jp/kegg/pathway.html] Accessed Jan 2014.
- Fleissner G, Fleissner G. Efferent control of visual sensitivity in arthropod eyes : With emphasis on circadian rhythms. Stuttgart, New York: G. Fischer Verlag; 1987.Google Scholar
- Rivera AS, Pankey MS, Plachetzki DC, Villacorta C, Syme AE, Serb JM, Omilian AR, Oakley TH. Gene duplication and the origins of morphological complexity in pancrustacean eyes, a genomic approach. BMC Evol Biol. 2010;10:123.View ArticlePubMedPubMed CentralGoogle Scholar
- Cellier-Michel S, Berthon JL. Rhythmicity of the pigments in the compound eye of Daphnia longispina (Cladocera). J Freshwat Ecol. 2003;18:445–50.View ArticleGoogle Scholar
- Meyer-Rochow VB. The crustacean eye: dark/light adaptation, polarization sensitivity, flicker fusion frequency, and photoreceptor damage. Zool Sci. 2001;18:1175–97.View ArticlePubMedGoogle Scholar
- Peñalva-Arana DC, Lynch M, Robertson HM. The chemoreceptor genes of the waterflea Daphnia pulex: Many Grs but no Ors. BMC Evol Biol. 2009;9:79.View ArticlePubMedPubMed CentralGoogle Scholar
- Hallem EA, Dahanukar A, Carlson JR. Insect odor and taste receptors. Annu Rev Entomol. 2006;51:113–35.View ArticlePubMedGoogle Scholar
- Krishnan N, Davis AJ, Giebultowicz JM. Circadian regulation of response to oxidative stress in Drosophila melanogaster. Biochem Biophys Res Commun. 2008;374:299–303.View ArticlePubMedPubMed CentralGoogle Scholar
- Borgeraas J, Hessen DO. Diurnal patterns of antioxidant activity in alpine and arctic Daphnia under in situ UV-radiation. Archiv für Hydrobiologie. 2002;156:83–95.View ArticleGoogle Scholar
- Tilden AR, McCoole MD, Harmon SM, Baer KN, Christie AE. Genomic identification of a putative circadian system in the cladoceran crustacean Daphnia pulex. Comp Biochem Physiol Part D Genomics Proteomics. 2011;6:282–309.View ArticlePubMedPubMed CentralGoogle Scholar
- Panda S, Hogenesch JB, Kay SA. Circadian rhythms from flies to human. Nature. 2002;417:329–35.View ArticlePubMedGoogle Scholar
- Zhu H, Yuan Q, Froy O, Casselman A, Reppert SM. The two CRYs of the butterfly. Curr Biol. 2005;15:R953–4.View ArticlePubMedGoogle Scholar
- Zhu H, Sauman I, Yuan Q, Casselman A, Emery-Le M, Emery P, Reppert SM. Cryptochromes define a novel circadian clock mechanism in monarch butterflies that may underlie sun compass navigation. PLoS Biol. 2008;6:e4.View ArticlePubMedPubMed CentralGoogle Scholar
- Bell-Pedersen D, Cassone VM, Earnest DJ, Golden SS, Hardin PE, Thomas TL, Zoran MJ. Circadian rhythms from multiple oscillators: lessons from diverse organisms. Nat Rev Genet. 2005;6:544–56.View ArticlePubMedPubMed CentralGoogle Scholar
- Gentile C, Rivas GBS, Meireles-Filho ACA, Lima JBP, Peixoto AA. Circadian expression of clock genes in two mosquito disease vectors: cry2 is different. J Biol Rhythms. 2009;24:444–51.View ArticlePubMedGoogle Scholar
- Escamilla-Chimal EG, Velazquez-Amado RM, Fiordelisio T, Fanjul-Moles ML. Putative pacemakers of crayfish show clock proteins interlocked with circadian oscillations. J Exp Biol. 2010;213:3723–33.View ArticlePubMedGoogle Scholar
- Peirson SN, Butler JN, Duffield GE, Takher S, Sharma P, Foster RG. Comparison of clock gene expression in SCN, retina, heart, and liver of mice. Biochem Biophys Res Commun. 2006;351:800–7.View ArticlePubMedGoogle Scholar
- Whitmore D, Foulkes NS, Strähle U, Sassone-Corsi P. Zebrafish Clock rhythmic expression reveals independent peripheral circadian oscillators. Nat Neurosci. 1998;1:701–7.View ArticlePubMedGoogle Scholar
- Harris JE. The role of endogenous rhythms in vertical migration. J Mar Biol Assoc UK. 1963;43:153–66.View ArticleGoogle Scholar
- der Linden V, Alexander M, Beverly M, Kadener S, Rodriguez J, Wasserman S, Rosbash M, Sengupta P. Genome-wide analysis of light-and temperature-entrained circadian transcripts in Caenorhabditis elegans. PLoS Bio. 2010;8:e1000503.View ArticleGoogle Scholar
- Minamoto T, Hanai S, Kadota K, Oishi K, Matsumae H, Fujie M, Azumi K, Satoh N, Satake M, Ishida N. Circadian clock in Ciona intestinalis revealed by microarray analysis and oxygen consumption. J Biochem. 2010;147:175–84.View ArticlePubMedGoogle Scholar
- Zhang L, Hastings MH, Green EW, Tauber E, Sladek M, Webster SG, Kyriacou CP, Wilcockson DC. Dissociation of circadian and circatidal timekeeping in the marine crustacean Eurydice pulchra. Curr Bio. 2013;23:1863–73.View ArticleGoogle Scholar
- Yang J, Dai Z, Yang F, Yang W. Molecular cloning of Clock cDNA from the prawn, Macrobrachium rosenbergii. Brain Res. 2006;1067:13–24.View ArticlePubMedGoogle Scholar
- Yampolsky LY, Zeng E, Lopez J, Williams PJ, Dick KB, Colbourne JK, Pfrender ME. Functional genomics of acclimation and adaptation in response to thermal stress in Daphnia. BMC Genomics. 2014;15:859.View ArticlePubMedPubMed CentralGoogle Scholar
- Latta LC, Weider LJ, Colbourne JK, Pfrender ME. The evolution of salinity tolerance in Daphnia: a functional genomics approach. Ecol Lett. 2012;15:794–802.View ArticlePubMedGoogle Scholar
- Roy Chowdhury P, Lopez JA, Weider LJ, Colbourne JK, Jeyasingh PD. Functional genomics of intraspecific variation in carbon and phosphorus kinetics in Daphnia. J Exp Zool A Ecol Genet Physiol. 2014;321:387–98.View ArticleGoogle Scholar
- Yoshii T, Vanin S, Costa R, Helfrich-Forster C. Synergic entrainment of Drosophila's circadian clock by light and temperature. J Biol Rhythms. 2009;24:452–64.View ArticlePubMedGoogle Scholar
- Kilham SS, Kreeger DA, Lynn SG, Goulden CE, Herrera L. COMBO: A defined freshwater culture medium for algae and zooplankton. Hydrobiologia. 1998;377:147–59.View ArticleGoogle Scholar
- Lopez J, Colbourne J. Dual-labeled expression microarray protocol for high-throughput genomic investigations. CGB Technical Report. 2011. doi:10.2506/cgbtr-201102.Google Scholar
- Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4:249–64.View ArticlePubMedGoogle Scholar
- Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, Hornik K, Hothorn T, Huber W, Iacus S, Irizarry R, Leisch F, Li C, Maechler M, Rossini AJ, Sawitzki G, Smith C, Smyth G, Tierney L, Yang JY, Zhang J. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5:R80.View ArticlePubMedPubMed CentralGoogle Scholar
- Bolstad BM, Irizarry RA, Astrand M, Speed TP. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003;19:185–93.View ArticlePubMedGoogle Scholar
- Storey JD. A direct approach to false discovery rates. J R Stat Soc B Statistical Methodology. 2002;64:479–98.View ArticleGoogle Scholar
- Department of Energy. Joint Genome Institute Daphnia Genome Browser [http://genome.jgi.doe.gov] Accessed Jan 2014
- Li L, Stoeckert Jr CJ, Roos DS. OrthoMCL: Identification of ortholog groups for eukaryotic genomes. Genome Res. 2003;13:2178–89.View ArticlePubMedPubMed CentralGoogle Scholar
- Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci U S A. 1998;95:14863–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Saldanha AJ. Java Treeview--extensible visualization of microarray data. Bioinformatics. 2004;20:3246–8.View ArticlePubMedGoogle Scholar
- Oster H, Damerow S, Hut RA, Eichele G. Transcriptional profiling in the adrenal gland reveals circadian regulation of hormone biosynthesis genes and nucleosome assembly genes. J Biol Rhythms. 2006;21:350–61.View ArticlePubMedGoogle Scholar
- Hulovatyy Y, DMello S, Calvo RA, Milenkovic T. Network analysis improves interpretation of affective physiological data. J Complex Networks. 2014;2:614–36.View ArticleGoogle Scholar
- Newman M. Networks: an introduction. New York: Oxford University Press; 2010.Google Scholar