Various parasite-host interactions that involve adaptive manipulation of host behavior display time-of-day synchronization of certain events. One example is the manipulated biting behavior observed in Carpenter ants infected with Ophiocordyceps unilateralis sensu lato. We hypothesized that biological clocks play an important role in this and other parasite-host interactions. In order to identify candidate molecular clock components, we used two general strategies: bioinformatics and transcriptional profiling. The bioinformatics approach was used to identify putative homologs of known clock genes. For transcriptional profiling, RNA-Seq was performed on 48 h time courses of Ophiocordyceps kimflemingiae (a recently named species of the O. unilateralis complex), whose genome has recently been sequenced. Fungal blastospores were entrained in liquid media under 24 h light-dark (LD) cycles and were harvested at 4 h intervals either under LD or continuous darkness. Of all O. kimflemingiae genes, 5.3% had rhythmic mRNAs under these conditions (JTK Cycle, ≤ 0.057 statistical cutoff). Our data further indicates that a significant number of transcription factors have a peaked activity during the light phase (day time). The expression levels of a significant number of secreted enzymes, proteases, toxins and small bioactive compounds peaked during the dark phase or subjective night. These findings support a model whereby this fungal parasite uses its biological clock for phase-specific activity. We further suggest that this may be a general mechanism involved in parasite-host interactions.

Funding: This work was supported by Marie Curie Actions IOF Project FP7-299501 ( http://ec.europa.eu/research/mariecurieactions ), the Alexander von Humboldt Foundation ( http://www.humboldt-foundation.de ) and funds from the LMUs Friedrich Baur Stiftung ( http://www.klinikum.uni-muenchen.de/Friedrich-Baur-Institut/de/institut/friedrich_baur_stiftung ). All funds were made available to CdB. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Circadian rhythms in ants have been reported [ 20 , 21 , 40 ]. To investigate if parasite circadian clocks could contribute to the observed patterns of manipulated ant behavior, we wished to characterize the temporal program in the parasite, Ophiocordyceps. To this end, we performed RNA-Seq on two-day long time courses of O. kimflemingiae. Rhythmic signals were paired with homology searches against N. crassa to detect candidate core clock genes in the genome. We report rhythmic candidate clock and clock-controlled genes with activity patterns that appear to be synchronized to certain times of day.

Fungal clocks have mainly been studied in the nonpathogenic fungus Neurospora crassa, revealing a negative feedback loop based on transcriptional regulation that is essential for circadian clock function [ 24 ]. In Neurospora, the activation arm of this molecular clock gene mechanism consists of the White Collar Complex (WCC). This complex is made up of the transcription factors White Collar-1 (WC-1) and White Collar-2 (WC-2), which have a second role, namely blue-light perception (important for entraining the clock [ 25 – 27 ]) and regulation of light-induced genes (including many clock-controlled genes, ccgs) [ 28 , 29 ]. As such, the WCC promotes transcription of the frequency (frq) gene. FRQ forms a complex with FRQ-interacting RNA helicase (FRH), which then functions as the repression arm in the transcriptional feedback loop as it inhibits the activity of WCC. After FRQ is synthesized it is post-translationally modified and degraded. As FRQ levels drop, the WCC is reactivated and the cycle can start anew (reviewed in more detail in e.g. [ 30 – 32 ]). Recent genome-wide gene expression studies in N. crassa, performed under constant darkness (free-running) conditions, show that the rhythmic expression of the many ccgs peak in time-of-day-related clusters [ 33 , 34 ]. Catabolic processes appear to be dawn-phased, while anabolic processes are dusk-phased [ 34 ]). Additionally, among the circadian clock-regulated genes, processes such as metabolism, protein synthesis, DNA processing and stress response are highly enriched. and phase-specific, with metabolism primarily peaking in the subjective late night and early morning, and protein synthesis activity rising in the subjective afternoon and early night. Expression levels of DNA processing genes also peaked during the subjective late day/evening [ 33 ]. Unfortunately, little is known about the molecular mechanism and clock-driven processes in other fungi. Saccharomyces cerevisiae shows typical circadian entrainment but no obvious homologs of clock genes are found in its genome [ 35 ]. For Aspergillus spp. the presence of a circadian clock was also shown, though these ascomycetes have different clock components since no homolog for frq was found [ 36 ]. Homologs of the N. crassa “clock genes” have been identified in the basal ascomycete Pyronema confluens and rhythmic expression of the P.c. frq has been observed [ 37 ]. Furthermore, studies have presented phenotypic and molecular evidence for the presence of a functional molecular clock in the phytopathogenic fungus B. cinerea, where it appears to contribute to virulence [ 9 , 10 ]. To a lesser extent, evidence for biological clocks in phytopathogens Cercospora kikuchii and Magnporthe oryzae have also been reported [ 38 , 39 ]. To our knowledge, fungal entomopathogens have not yet been studied from a circadian biology perspective.

In this work, we wish to learn more about how the circadian clock is used in parasitic behavioral manipulation, specifically in stages past the initial host invasion and immune system evasion. Such a behavioral manipulation can be observed in Formica ants infected with the lancet fluke that are manipulated to bite grass. This behavior appears regulated by daily temperature fluctuations [ 14 – 16 ] and its phase-specific activity facilitates parasite transmission from the intermediate ant host to the gut of the ultimate, herbivore host. Another example of daily synchronization of manipulated behaviors is seen in soldier beetles infected with Eryniopsis lampyridarum. Here, infected beetles are thought to grip flower heads with their mandibles in the early morning immediately prior to death. Emergence of conidiophores and conidia (spores; transmission agents) and post-mortem spreading of the wings occurs pre-dawn the following day. This timing is thought to protect the spores via optimal humidity levels [ 17 ]. Another time-of-day specific fungal manipulation of insect behavior is that of fungi of the Ophiocordyceps unilateralis complex as they infect Carpenter ants. Workers foraging for food are infected when they encounter fungal spores [ 18 ]. Infected ants then abandon their tasks as the fungus slowly expands throughout the body [ 19 ]. This behavioral change may represent a disruption of the ant circadian clock, as foraging, and other behaviors, are typically under clock control [ 20 , 21 ]. Ultimately, like in the examples above, the ant displays manipulated biting behavior followed by death at a certain time of day. In field studies in Thailand, Ophiocordyceps camponoti-leonardi induced Carpenter ants to bite vegetation around solar noon [ 22 ]. Laboratory studies with Ophiocordyceps kimflemingiae demonstrated synchronized biting in the early morning [ 23 ]. In both unilateralis species’ interactions, the timing of death after biting is also synchronous. Moreover, mixed transcriptomics data revealed differentially expressed clock gene homologs in O. kimflemingiae-infected ants versus healthy individuals sampled at the same time of day [ 23 ]. This suggests that the fungal infection works at some level through the host circadian clock.

Endogenous temporal programs, which anticipate daily changes in e.g. temperature and light, have evolved in organisms from all phyla [ 1 – 4 ]. These programs function as circadian clocks, organizing physiology and behavior to specific times of day (phases). Circadian clocks share several properties, most notably a free-running rhythm of about 24 h in constant conditions and entrainment or synchronization of these rhythms to exactly 24 h in the presence of highly predictable, regular signals from the environment (zeitgebers, e.g. light and temperature cycles). Various insect behaviors, such as foraging, flight and oviposition show phase-characteristic daily rhythms under entrainment conditions [ 5 , 6 ]. Daily rhythms in immune systems have also been reported in insects, as well as in plants [ 7 , 8 ], an observation that has consequences for timing of inter-species interactions with respect to parasites. For example, it has been proposed that fungi invade plants such that they circumvent rhythms in host immune systems [ 9 ]. Moreover, it has been shown for the plant pathogen Botrytis cinerea that its virulence is regulated by the circadian clock [ 10 ]. Adaptive timing of infection could thus enable a parasite’s life cycle by increasing infection rates, virulence and, eventually, transmission. In vertebrate-infecting parasites, circadian clocks also appear to play a role. For example, circadian rhythms were shown to contribute to fitness in the malaria parasite Plasmodium chabaudi and metabolism in the sleeping sickness parasite Trypanosoma brucei, is under clock control [ 11 – 13 ]. Biological clocks could, therefore, be an important aspect of parasite-host interactions in general.

To normalize samples within a time course, candidate reference genes were selected from genes displaying non-cycling expression patterns in the same RNA-Seq dataset, as well as low variation expression patterns in the RNA-Seq dataset from [ 23 ]. The pool of non-cycling genes (JTK-CYCLE p >0.057) was first restricted by selecting for genes with a low fold-change, as done in [ 47 ] and by setting a cutoff for the coefficient of variation, as described in [ 48 ]. Reference candidates displayed a <2-fold-change between the maximum and minimum RPKM expression values under both LD and DD conditions. In addition, log2-transformed expression values displayed a coefficient of variation ≤5%. Reference candidates were further restricted to genes that lacked differential expression in a previously published transcriptomics dataset that compared stages of the fungal parasite in culture, inside the manipulated host and inside the dead host [ 23 ]. From this restricted pool, five genes were amplified using RT-QPCR after which the two most stable ones were selected with RefFinder [ 49 ] ( Table 1 ).

RNA-Seq results were verified through RT-QPCR by determining the relative expression of five clock gene homologs. Two biological replicates from independently performed time courses and technical triplicates for each time point and amplified gene were used. Primers were designed to span exon-exon boundaries when possible ( Table 1 ). Efficiencies were calculated from a standard curve of serial dilutions. RNA was converted to cDNA with the QuantiTect Reverse Transcription Kit (Qiagen). RT-QPCR reactions were performed on a ViiA 7 (Applied Biosystems) running in Fast Cycling Mode followed by a melt curve analysis, with 10 ng cDNA per 10 μL reaction and PowerUp SYBR Green Master Mix (Applied Biosystems). Amplification data was analyzed using the Relative Quantification and the Design and Analysis Application Apps of the ThermoFisher Cloud Software.

Functional annotation enrichment was determined with a Fisher Exact test with the Benjamini-Hochberg procedure and a corrected p-value cutoff of 0.05 to correct for multiple testing. Hierarchical clustering was performed with Cluster 3.0 and resulting heat maps were visualized with Java TreeView 1.16r4 [ 46 ]. For hierarchical clustering, normalized gene expression values were increased by 1 prior to log2-transformation to avoid negative values, and scaled per gene by centering around the mean. Next, genes were hierarchically clustered using the average linkage method.

RNA reads were mapped to the latest version of the O. kimflemingiae genome [ 44 ] with CLC Genomics Workbench 9. Reads were also mapped to inter-genic regions and other parameters were kept at default. Expression levels were calculated for all genes and normalized to Reads Per Kilobase of exon model per Million mapped reads (RPKM).

Total RNA was extracted as previously reported [ 43 ]. Sample quality and quantity were checked with a Bioanalyzer and RNA 6000 Nano Kit (Agilent Technologies), and a Qubit with the RNA BR Assay Kit (Invitrogen). Sequence libraries were constructed using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England BioLabs) and manufacturer’s instructions. Library quality and quantity were checked with a Bioanalyzer and the High Sensitivity DNA Analysis Kit (Agilent Technologies) and a Qubit with the ds DNA HS Assay Kit (Invitrogen). Samples were sequenced on a HiSeq1500 (Illumina) using 50 bp single-end sequencing in rapid-run mode. The dataset is available at the NCBI’s Gene Expression Omnibus under the accession code GSE101312.

The fungal parasite used in this study, O. kimflemingiae [ 41 ], was isolated from an infected and manipulated Camponotus castaneus ant [ 42 ]. This species resides within the O. unilateralis s.l. complex and was recently named [ 41 ]. Previous data for this species’ interaction with C. castaneus suggests that biological clocks could be involved in the observed manipulation [ 23 , 42 ]. Cultures were maintained in Grace’s insect medium (Sigma) supplemented with 10% fetal bovine serum (Invitrogen). For time course studies, cultures were grown as fungal blastospores ( Fig 1A ) in 25m 2 tissue culture flasks in 20 mL media and shaken at 60 rpm at 28 o C. Fungal hyphae ( Fig 1B ) were grown in a 250 mL Erlenmeyer flask holding 100 mL media shaken at the same speed. To synchronize fungal cells, cultures were kept in constant light (light intensity 300 Lux) for 48 h and then transferred to 24 h LD (12 h:12 h) cycles for 5 days at 28 o C. To avoid starvation conditions during sampling, on the fourth day at the end of the light phase, cultures were briefly spun down in a swing-out rotor (3000 rpm for 5 min), and old media was removed and replaced with fresh media. Each culture was subsequently split in two 20 mL cultures: one to keep in LD during sampling, and one to be transferred to DD. After an additional 24 h under LD conditions, half of the cultures were transferred to constant darkness (DD samples). The other half of the flasks were maintained under LD conditions (LD samples) ( Fig 1C ). During this time, fungal cultures had reached an OD600 of ~ 1 and material was collected every 4 h under red light in a climate controlled room set to 28 o C. Cultures were kept in incubators running zeitgeber cycles in antiphase, which allowed harvesting 48 h time courses over 12 h. Fungal cells were pelleted by spinning down 2 mL of the culture for 3 min at 10,000 rpm. The supernatant was removed and the samples were snap frozen in liquid nitrogen and stored at -80 o C until processing for RNA extraction.

Results and discussion

O. kimflemingiae blastospores: Growth and sampling conditions Fungal entomopathogens have a yeast-like (blastospores) and a vegetative (hyphae) growth state. Daily foraging behavior is disrupted and phase-specific biting behavior is induced when the infecting fungus is in a blastospore state (Fig 1A). When secured in the biting position, the host dies and the fungus quickly switches to a vegetative, hyphal growth (Fig 1B). For this study, we identified conditions for O. kimflemingiae that would induce morphological growth resembling that in situ (blastopores in an infected, live C. castaneus ant). Culturing O. kimflemingiae on solid media results in vegetative (hyphal) growth. Hyphae are also formed in liquid shaking cultures grown in Erlenmeyer flasks. In order to trigger blastospore growth in liquid cultures, we tried several media compositions (S1 Table). Vegetative growth persisted in all conditions. We next cultured fungi in tissue culture flasks with rectangular edges shaken gently at 60 rpm. In the presence of Grace’s medium with 10% FBS, the fungi switched to a yeast-like growth and started forming blastospores. Blastospores reverted back to hyphal growth when they were placed into Erlenmeyer flasks. We do not quite understand the reason for this but speculate that this could be due to different aeration rates. The greater turbulence in tissue culture flasks might also decrease flocculation and, therefore, lead to favored blastospore growth. For the experiments in this study, fungal cells were grown in tissue culture flasks throughout time course experiments. Experiments were performed as described in the Materials and methods section and in Fig 1C.

Mining the O. kimflemingiae genome for clock gene homologs Clock and clock-controlled genes of fungi, cyanobacteria, plants, insects and mammals have been reported [1]. We surveyed the O. kimflemingiae genome for homologs by searching for protein family (Pfam) domains [50] that are encoded in those genes. Per Arnt Sim (PAS) domains (PF00989) are found in all kingdoms of life [51] and have been shown to be involved in eukaryotic circadian clocks [52]. A frequency (FRQ) domain (PF09421) has been identified in the clock gene frequency of the fungal model organism N. crassa [53] as well as other fungi [10, 37]. PHY domains (PF00360) are found in phytochrome photoreceptors that can sense (far-)red light. Phytochromes have been mainly studied in plants where they regulate development via photoperiodism. To a lesser extent, they have been reported for (cyano)bacteria and fungi as well. In Aspergillus nidulans this receptor was shown to repress sexual spore formation [54] as it forms a light-regulator complex with receptors similar to the white collar proteins in N. crassa [55]. In N. crassa, however, phy transcripts are not regulated by light. There, phytochrome-1 mRNAs levels appear under the control of the circadian clock [56]. We mined the O. kimflemingiae genome for the presence of these functional domains and found eight candidate homolog genes (Table 2). The Pfam search was followed by BlastP alignments against the N. crassa OR74A genome [57]. These alignments confirmed O. kimflemingiae homologs of the N. crassa clock genes frequency (frq), white collar 1 and 2 (wc-1 and wc-2), and vivid (vvd). In Neurospora, VVD is only expressed in the light and regulates gating of light input to the circadian clock [58]. Four additional homologs of PAS domain-containing genes were found through alignment, one of which was a putative phytochrome 1 (phy-1) (Table 2). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Table 2. Identified clock(-controlled) gene homologs. https://doi.org/10.1371/journal.pone.0187170.t002

Detection of candidate genes with cycling transcripts The identification of clock gene homologs suggests that O. kimflemingiae may have clock-controlled daily oscillations in gene expression. To detect genes that oscillate under LD and DD conditions, we generated genome wide gene expression profiles of samples collected over a 48 h time course using RNA-Seq. Reads were mapped to the latest version of the O. kimflemingiae genome, which is 23.92 Mb in size and has been annotated to encode 8,629 genes [44]. Sequencing yielded an average of 20.5 million single-end reads per sample with an average mean quality score of 38. This read depth is around the maximum read depth that was tested for D. melanogaster [59]. It should thus be sufficient to detect cycling transcripts in our fungal model since its genome is about five times smaller. An average of 95.9% of the reads mapped uniquely to this genome. Resulting expression profiles were normalized to Reads Per Kilobase of exon model per Million mapped reads (RPKM). Normalized expression profiles were analyzed with the JTK_CYCLE algorithm [45]. Our sampling method (single samples at 4 h intervals) did not allow us to assign false discovery rates. We, therefore, ranked genes based on the Bonferroni-adjusted p-values. In our ranking approach, we relaxed the p = 0.05 threshold, which has been the topic of debate for many years now [60, 61]. We set the threshold for rhythmicity at the p-value of an Ophiocordyceps homolog of a rhythmically expressed Neurospora clock gene. In N. crassa, ca. 24 h rhythmic expression levels have been observed for frq under both LD and DD conditions. Under LD conditions, vvd is rhythmic, but when transferred to DD, vvd oscillation is only seen in the first 24 h. For wc-1, rhythmicity has been observed in LD and though its promoter shows rhythms in DD, its mRNA levels are not predicted to be highly rhythmic [26, 33, 58, 62, 63]. In addition, we used the expression patterns from fungal tissue cultured in LD conditions. JTK_CYCLE returned Bonferroni-adjusted p-values of 0.002, 0.057, and 0.009, respectively, for the above-mentioned homologs under LD conditions. As such, we set our threshold for rhythmicity at the adjusted p-value for wc-1 under these conditions (i.e. rhythmic candidates p ≤ 0.057). This led to the detection of 333 candidate genes with oscillating transcripts under LD conditions with a median period of 24 h (S1 Fig and S1 Data File). Under DD conditions, adjusted p-values for these genes were markedly higher (0.30, 0.46, and 1 for frq, wc-1, vvd, respectively). This suggests either that endogenous rhythmicity at the mRNA level damps over 48 h in DD and constant temperature conditions or that rhythms detected in LD are simply driven by the zeitgeber cycle. To minimize false discovery, we set the threshold for rhythmicity in DD to the same adjusted p-value as for the LD samples. This led to the detection of 154 genes with suggested circadian rhythms with a median period of 26 h (S1 Fig and S2 Data File). Studies using other model organisms also report fewer rhythmic genes under DD than under LD conditions [64–68]. In addition, lower numbers of cycling transcripts for cells taken out of their natural context and grown under laboratory cell culture conditions have previously been shown [69]. Moreover, a decreased amplitude under DD conditions for the significantly rhythmic gene frq is observed in B. cinerea [10]. Similarly, a decreased amplitude in O. kimflemingiae under DD conditions, combined with single 4h-interval data points could have resulted in poor detection of endogenously rhythmic genes in this study. Among the rhythmic candidate genes in LD, we identified 14 transcription factors (TFs). One of these TFs, as identified above, is the O. kimflemingiae homolog of wc-1. Two of these TFs appear to be involved in regulation of Neurospora conidiation (homologs of N. crassa znf-21; NCU02671 and csp-2; NCU06095) and one regulates ascus development (N. crassa homolog asd-4; NCU20921). ChIP sequencing revealed that N. crassa TFs znf-21 and csp-2 are among the direct targets of the White Collar Complex (WCC, formed by WC1 and WC2) [70], a blue-light photoreceptor that is the key transcriptional activator of the circadian oscillator in this fungus. Additionally, we identified a putative blue-light photoreceptor cryptochrome DASH homolog (cry; NCU00582) (S1 Data File). In Neurospora, CRY was shown to be an essential photoreceptor for the correct functioning of a novel frq-less-oscillator (FLO) named the CRY-dependent oscillator [71]. We performed enrichment analyses on the rhythmic genes. For both LD and DD conditions, secreted proteins (SPs) and small secreted proteins (SSPs) were significantly over-represented (i.e. LD: 41 SPs and 25 SSPs, Fisher-exact corrected p-value 3.88E-07 and 0.0029, respectively; and DD: 15 SPs and 15 SSPs, Fisher-exact corrected p-value 0.0008 and 0.00160, respectively). One of the SSPs was homologous to clock-controlled gene-6 of N. crassa (ccg-6; NCU01418). Rhythmic expression has been demonstrated for this conidiation-related gene in N. crassa, which also appears to be photoinducible and developmentally regulated [72]. Among the rhythmic SPs under LD conditions were five genes encoding enterotoxins and six encoding predicted proteases. These, as well as the potentially bioactive SSPs, could be involved in O. kimflemingiae interactions with its ant host. Moreover, we found rhythmicity in genes encoding secreted enzymes that have previously been identified as candidate manipulation genes: tyrosinase and protein tyrosine phosphatase, [23]. Twenty-six genes had rhythmic transcripts under both cycling and free-running conditions (S1A Fig). A relatively low overlap between LD and DD rhythmic transcripts has previously been observed [65–68, 73]. To obtain insight into the biological processes that were rhythmic regardless of LD or DD conditions, we again performed an enrichment analysis. Despite the small number of endogenously cycling genes (26), we again found an over-representation of SPs (six genes) and SSPs (five genes; Fisher-exact corrected p-value 1.95E-05 and 0.0030, respectively). Nine out of the eleven genes with encoded secretion signals appeared to peak during the dark phase in LD and in the subjective night in DD (Fig 2A). These comprised four unknown SSPs that appear unique to O. kimflemingiae since they only aligned with deposited sequences of this species in the NCBI Database. Night-active SPs comprised a chloroperoxidase, an exo-beta-D-glucosaminidase, a metallocarboxypeptidase, an enterotoxin, and a tyrosinase. The other two genes peaked during the (anticipated) light phase and comprised another unique SSP with unknown function and a homolog of the Aspergillus fumigatus allergen Asp f9, which is a cell wall glucanase (Fig 2A). Two fungal transcription factors (TFs) were also among the endogenously rhythmic genes. We additionally identified a rhythmically expressed histidine phosphotransferase, a gene encoding a predicted two-component regulatory system protein, which would facilitate responses to changing stimuli in the environment [74]. This histidine phosphotransferase, as well as one of the fungal TFs appeared to peak during the (subjective) day, while the other TF showed peaks in expression during the (subjective) night (Fig 2B). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 2. Expression patterns of endogenously cycling genes. (A) Rhythmic mRNA levels of secreted proteins (SPs) and small secreted proteins (SSPs) under cycling (LD) and free-run (DD) conditions. (B) Rhythmic mRNA levels of genes involved in regulatory systems under cycling (LD) and free-run (DD) conditions. Normalized expression levels are plotted as their standard scores (Z-scores) and represent the amplitude of the gene expression observed at that time point versus the mean expression measured over the entire time course of 48 h. https://doi.org/10.1371/journal.pone.0187170.g002