Non-Pleiotropic Coupling of Daily and Seasonal Temporal Isolation in the European Corn Borer

Speciation often involves the coupling of multiple isolating barriers to produce reproductive isolation, but how coupling is generated among different premating barriers is unknown. We measure the degree of coupling between the daily mating time and seasonal mating time between strains of European corn borer (Ostrinia nubilalis) and evaluate the hypothesis that the coupling of different forms of allochrony is due to a shared genetic architecture, involving genes with pleiotropic effects on both timing phenotypes. We measure differences in gene expression at peak mating times and compare these genes to previously identified candidates that are associated with changes in seasonal mating time between the corn borer strains. We find that the E strain, which mates earlier in the season, also mates 2.7 h earlier in the night than the Z strain. Earlier daily mating is correlated with the differences in expression of the circadian clock genes cycle, slimb, and vrille. However, different circadian clock genes associate with daily and seasonal timing, suggesting that the coupling of timing traits is maintained by natural selection rather than pleiotropy. Juvenile hormone gene expression was associated with both types of timing, suggesting that circadian genes activate common downstream modules that may impose constraint on future evolution of these traits.


Introduction
Speciation typically occurs through the accumulation of multiple reproductive barriers [1,2].The evolution of reproductive isolation is tied to the coupling of different barrier traits that produce increased isolation when acting together when compared to each acting separately [3,4].Although we understand how coupling between prezygotic and postzygotic isolation occurs through indirect selection during reinforcement [5], the mechanisms that produce coupling between different prezygotic barriers are largely unexplored [3].It is possible that some prezygotic traits could be coupled through a shared genetic architecture rather than by direct or indirect selection.For example, coupling could occur due to pleiotropy, with a change in a single locus producing multiple types of isolation [6].Coupling may also emerge from physical linkage or reduced recombination on a chromosome.Such genetic architecture would lead to the rapid evolution of prezygotic reproductive isolation, but it may increase evolutionary constraint of barrier traits and potentially slow the speed of future evolution [7][8][9][10].In contrast, different and unlinked genes underlying coupled barrier traits could allow for faster, independent evolution of traits, but require selection to generate and maintain coupling of reproductive isolation.Information on pleiotropy's role in coupling is also needed to estimate the effect on reproductive isolation for any single "speciation" gene [11].Thus, genetic mechanisms contributing to coupling among barrier traits will have important consequences for the speed, stability, and direction of the evolution of reproductive isolation.
Variation in the timing of mating often contributes to prezygotic reproductive isolation between sympatric populations by reducing the encounter rates between potential mates (eels [12]; salmon [13]; corals [14,15]; mosquitos [16]; moths [17][18][19][20]).Daily (24 h) and seasonal (time of year) forms of allochrony can be coincident and can both limit gene flow between incipient species [21,22], but it is unknown how they become coupled.Given that both forms of allochrony involve endogenous time-keeping, it is possible they rely on the same circadian clock pathway genes [23].However, there have been few direct comparisons between the genes that are associated with different forms of allochrony in the same species.Studying the mechanisms that produce daily timing variation in a species with seasonal timing variation would allow us to determine if genetic mechanisms underlie coupling, improving our understanding of how these allochronic isolating barriers evolve and lead to speciation [22].
Variation in daily mating time is the result of variation in the timing of release of sexual signals to attract mates, receptivity to these signals, and the release of gametes [24].The timing of these traits may often be controlled by the circadian clock pathway, which regulates sexual traits, such as courtship song [25] and pheromone release [24].The circadian clock synchronizes biological rhythms with 24-h cycles of light and dark using light-sensitive proteins and negative feedback loops (Figure 1) [26].The heterodimer formed by CLOCK (CLK) and CYCLE (CYC) activates the transcription of a number of circadian genes by binding to E-box sequences in their promoters; these include period (per), timeless (tim), vrille (vri) and PAR domain protein 1 (pdp1) [27][28][29][30].These activated genes feedback to regulate clk and cyc transcription (Figure 1).Mutations in these core circadian genes have been associated with changes in daily timing of mating, particularly in insects (Drosophila [31,32], Bactrocera [33,34], Spodoptera [35], and Clunio [36]).Although evidence suggests that the differences in daily timing are often associated with the changes in regulation or expression of circadian clock genes, whether daily and seasonal biological rhythms both rely on the circadian clock pathway remains a subject of debate.Seasonal cues, such as increases in temperature and day length (photoperiod), strongly influence the time of year at which insect mating occurs by triggering dormancy break and reproductive development [37].The light-sensitive circadian clock pathway could potentially sense changes in photoperiod and pleiotropically cause both circadian and seasonal variation in mating time [38].Diapause response often shows circadian rhythmicity [39] and circadian clock genes do alter the photoperiodic responses that are associated with seasonal timing in species, such as bean bugs (Riptortus pedestris), where cycle RNAi impedes diapause break and period RNAi affects diapause entry [40].However, a functioning endogenous clock does not seem to be necessary for diapause in Drosophila [6].Photoperiodic response can evolve separately from circadian rhythm in pitcher plant mosquitos Wyeomyia smithii [41,42], providing further evidence that seasonal timing may be independent of the circadian clock.
Natural variation in daily and seasonal mating time occurs in the European corn borer moth (ECB; Ostrinia nubilalis) and causes incomplete reproductive isolation, providing an ideal system to explore coupling and the genetic basis of these traits [4].Two ECB strains are defined by the use of different pheromones that females use to attract males.In the Z strain, females produce and males preferentially respond to a 3:97 ratio of E and Z isomers of 11-tetradecenyl acetate, while the E strain produces and responds to a 99:1 pheromone ratio [43].In some locations in North America, sympatric E and Z strains differ in seasonal mating time, with the E strain having a mating flight in May and the Z strain mating flight occurring in July [4,44].Seasonal timing alone provides incomplete reproductive isolation, limiting 66% of possible matings between strains [4].Daily temporal isolation has also been reported, but it is unknown if this is coupled with seasonal timing.Using populations with unknown seasonal mating times, Liebherr & Roelofs (1975) observed a two-hour difference in the mean mating time during the night between strains, with Z strain from Ontario mating two hours earlier than E strain from New York [45].A reversed pattern has been documented in European populations; E-strain females release pheromone and mate earlier than Z strain females [46].
In ECB, seasonal timing is at least partially determined by the time that is needed for overwintering larvae to terminate winter diapause (dormancy) in the spring and summer months [47].Diapause termination is triggered by a combination of photoperiodic and temperature cues [48] and can be fast (~14 days after receiving diapause breaking cues) or slow (~44 days).These time windows ultimately determine the time of mating flight and the number of generations produced per summer [44,47,49,50].In New York, E strain moths exhibit fast diapause termination while sympatric Z strain individuals have slow termination [4,44,47].Per is physically linked to the major quantitative trait locus (QTL) controlling the termination time [51][52][53].Wadsworth & Dopman (2015) quantified gene expression prior to and during diapause termination [54].Genes that were differentially expressed between strains during diapause termination and/or differing in amino acid sequence included circadian clock genes (pdp1 and per) and several insulin-signaling genes [54].
In this study, we investigate differences in the daily mating time between E and Z ECB strains that have been previously characterized for differences in seasonal timing.Strains had been maintained in a common lab environment for several generations to ensure that observed differences among strains are genetic rather than environmentally induced.We test if behavioral nighttime differences in mating are phenotypically coupled to differences in seasonal timing and estimate how much reproductive isolation these differences cause.If cis-acting genetic changes in gene regulation underlie differences in mating time, we would expect causal loci to differ in the expression at the peak mating time for each strain.We therefore examine patterns of gene expression associated with changes in times of peak mating activity to identify genetic and regulatory factors that are associated with and potentially contribute to divergence in circadian mating time in ECB.Furthermore, we compare the overlap between genes differentially expressed during mating to candidate genes for seasonal timing (those differentially expressed during diapause termination [54]) to explore the hypothesis that genetic mechanisms produce coupling and that common molecular changes explain the variation in both types of reproductive timing.

Mating Trials
The stocks of bivoltine E (BE) and univoltine Z (UZ) strain ECB were donated by Charles Linn at the New York State Agricultural Experiment Station in Geneva, NY.The BE stock was originally derived from bivoltine E strain individuals collected from Geneva, NY (42.8680 • N, 76.9856 • W) and the UZ stock from univoltine Z strain individuals that were collected from Bouckville, NY (42.8892 • N, 75.5513 • W).The stocks were maintained via mass rearing (>200 individuals per generation) at Tufts University for several generations and were the same lines that were previously characterized for seasonal timing [54].ECB larvae were reared on corn borer diet (Southland Products, Lake Village, Arkansas, USA) under 16:8 L:D, 26 • C conditions until pupation.300 pupae per strain were removed from rearing containers and placed in individual 44.4 mL plastic cups with lids, each with a 3.8 cm dental wick soaked in water, and returned to the incubator.We replicated methods that were used in a historical study of ECB mating time [45].One day after eclosion, two males and one female of the same strain were randomly chosen and transferred from individual 44.4 mL cups to one 355 mL transparent plastic container for mating observations.Transfers were completed within the last hour of photophase.During scotophase (16:00-23:00), all of the mating groups were checked for mating activity every 20 min.Each mating group was used only once, with initiation of mating recorded over a 7-h period, for a total of 76 E-strain mating trials and 79 Z-strain trials.
We calculated the absolute strength of the reproductive isolation stemming from mating time following Dopman et al. [4].We generated the expected number of hybrid and pure-strain offspring after one generation of random mating using Hardy-Weinberg proportions.For each one-hour window (i, from 1 to w), the expected frequency of hybrid (2p i q i ) and pure-strain offspring (p i 2 + q i 2 ) offspring was calculated based on the number of mating moths that were Z strain (p i ) or E strain (q i ).The expected number of hybrid and pure-strain moths each hour was then calculated by multiplying the total number of moths (n i ) by their expected frequencies.Summing the expected numbers of intra-and inter-strain matings across hours (w) provided an estimate of the total number of hybrid and parental offspring produced during scotophase.Thus, the strength of reproductive isolation ranges from 0 (random mating) to 1 (complete reproductive isolation), and it is given by Equation (1) below:

Sampling Time Points
In order to identify changes in expression related to circadian control of female mating receptivity, independent of male presence or the physical act of mating, a second experiment was conducted to examine RNA expression levels for the two strains.Additional individuals were reared, as described above, but because mating time is predicted to be female-controlled, individuals were sexed as pupae and only female pupae were kept in isolated 44.4 mL cups until eclosion.
One-day-old virgin females were randomly selected for sacrifice at one of three time points.The first time point was one hour before the end of photophase and provided a daylight baseline for expression change during scotophase (photophase time point).The two other time points corresponded to the median observed mating time for each strain, as determined during the mating trial experiment: 1.3 h into scotophase (median mating time E strain) and 4 h into scotophase (median mating time Z strain).At designated sacrifice time points, the containers were placed in triple-layered black plastic bags and were kept at −20 • C for 10 min to sedate the moths.Under dark conditions, female moths were transferred to 1.5 mL microcentrifuge tubes containing 1.0 mL of RNAlater (Qiagen, Hilden, Germany).Moths in RNAlater were stored at −20 • C for the duration of the experiment, before long-term storage at −80 • C. A total of 72 females of each strain were preserved, and 60 of these were ultimately used for sequencing.

Library Sequencing
Adult moths were decapitated and heads were pooled according to time point and strain.Heads were used due to known circadian gene expression in antennal and neural tissue [55][56][57][58].Tissue from five combined heads was used, with four replicate pools for each strain and time point, for a total of 24 pools used for RNA extraction.Pooled heads were used to get sufficient RNA for RNA-seq.RNA was extracted from pooled samples using RNeasy kits (Qiagen, Hilden, Germany).RNA samples were quantified with a Nanodrop (Thermo Scientific, Wilmington, DE, USA) and Qubit Broad Range RNA assays (Life Technologies, Carlsbad, CA, USA).
cDNA libraries were prepared from mRNA using the TruSeq Sample Prep Kit v2 Set A (Illumina Inc., San Diego, CA, USA) using 1 mg total RNA, and prepared libraries were quantified using the Qubit High Sensitivity DNA assay.Libraries were quantified a second time on an Agilent Bioanalyzer (Santa Clara, CA, USA).Two replicate libraries for each strain and time point were run on each of the two lanes of an Illumina HiSeq 2500, located at the Tufts University Core Facility for Genomics (Boston, MA, USA) to generate 100 bp single-end reads from 23 libraries.One UZ 1.3 h library failed.Sequences are available in the GenBank Sequence Read Archive (SRA, accession: SRP135924).

Transcriptome Assembly
Single-end Illumina sequencing reads were assessed for quality using the FastQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Sequences were then trimmed using Trimmomatic version 0.32 to remove adapter sequences, bases with low sequence quality, and any reads that were shorter than 36 base pairs [59].FastQC reports were generated for each file again to confirm post-trimming quality.Mitochondrial DNA and ribosomal RNA sequences were removed using Bowtie2 version 2.1.0[60] by aligning against known mtDNA sequences and identical reads were collapsed prior to assembly (but counts retained) using the FastX Toolkit version 0.013 (http://hannonlab.cshl.edu/fastx_toolkit).The transcriptome was assembled de novo using Trinity and a k-mer length of 25 from all reads (both strains and all time points) [61].The longest transcript for each component were retained using custom scripts and the R package SeqinR [62].This transcriptome is available in Dryad (doi:10.5061/dryad.k7t20).

Differential Expression Analysis
The set of longest transcripts was indexed as a reference transcriptome, and the individual reads were mapped back to this transcriptome using Bowtie2 version 2.1.0[60].The resulting SAM file was parsed by individual libraries to generate count data, which were merged and filtered for weakly expressed genes, defined as genes that did not have at least one count per million reads in at least seven libraries (the minimum number of libraries representing the two strains at a given time point).Libraries were normalized for GC content and between-library sequencing depth using the R package EDASeq [63] and input as offsets along with the raw counts into edgeR [64].
Differential expression analysis was done by fitting negative binomial generalized log-linear models in edgeR [64].Specific contrasts were done between strains at the photophase, 1.3 h, and 4 h time points, as well as within-strain between each pair of time points (Figure 2).P-values were corrected for multiple testing using a false-discovery rate (FDR) cutoff of 0.05.Additionally, we emphasized transcripts that were specifically upregulated during at the peak mating time for each strain.We compared genes that were significantly differentially expressed between strains at any mating time point (photophase, hour 1.3, hour 4; FDR q < 0.05) to genes that were differentially expressed between strains during diapause break (day 1 and day 7) from the data set in Wadsworth & Dopman [54].Orthologous transcripts were identified via blastn and retaining only transcripts with blast scores that were greater than 100.

Annotation
Transcripts were annotated by BLASTing sequences against known Ostrinia nubilalis-Drosophila melanogaster gene pairs [65].Supplemental annotations were derived from blasts against the Bombyx mori [66] and Danaus plexippus proteomes [67].For circadian transcripts, we list the location in the homologous B. mori protein by amino acid number (BMAA).Transcripts were also blasted against known B. mori gene pairs to estimate chromosomal locations.Enriched gene ontology (GO) terms were identified for differentially expressed transcripts.This was done using the FlyBase D. melanogaster gene IDs [68] derived from the annotated gene pairs and GOrilla [69].GO enrichment was performed only for transcripts that were also differentially expressed between strains during diapause break.We also searched for previous identified olfactory receptor transcripts from Ostrinia using BLASTn searches against previous published mRNA sequences from Ostrinia furnacalis [70,71].
Correlations in expression among annotated circadian genes across time points were explored using data from both strains together (six time points total, three E and three Z).For each transcript, normalized read counts were averaged within strain and time point.Correlations in mean expression were calculated across time points and strains using custom scripts and the corrplot package in R [72].

Mating Trials
For E strain, a total of 21 of the 76 females (27.6%) mated at any point during the scotophase period.For the Z strain, 25 of the 79 females (31.6%) were observed to mate at any time during the night.
The proportions of females mating from each strain were not significantly different (Fisher's exact test, p = 0.602).The median mating times for each strain were 1.33 h into scotophase for the E (mean = 2.41 h) and 4 h for the Z (mean = 3.96).An independent two-group Mann-Whitney-Wilcoxon test indicated that the distributions of the two strains were also significantly different (W = 395, p = 0.0035, Figure 3).When mating times were compared with Liebherr & Roelofs [45], Mann-Whitney-Wilcoxon tests showed that the current distribution for the E strain differed from the historical distribution (W = 371, p < 0.0001), as did the Z strain distribution (W = 420.5,p = 0.014).The absolute strength of circadian mating time as a reproductive isolating barrier between the current E and Z strain populations was calculated to be 0.393, whereas the historical absolute strength of circadian mating time was 0.525.

Transcriptome Assembly
After trimming of raw Illumina reads and removal of mtDNA and rRNA sequences, 61.5 million unique reads remained, representing a total of 78,236 transcripts and 40,811 components.The mean read length was 1136 base pairs, with a minimum of 201, a maximum of 17,430, a median of 691, and an N50 length of 1918 (Table S1).

Differential Expression Analysis
Of transcripts for which differential expression was calculated, a total of 4449 unique transcripts (24.25%) showed patterns of significant differential expression in any contrast.The number of differentially expressed (DE) transcripts varied greatly between the different contrasts tested (Figure 2).Generally, far more DE transcripts were found in between-strain comparisons than any within-strain comparisons.
Sex peptide receptor is related to the termination of pheromone calling and receptivity to mating [74].Five transcripts of sex peptide receptor were differentially expressed between strains; four of which were upregulated in the Z strain in at least one contrast.Two transcripts of ecdysis-triggering hormone receptor (comp36183) were upregulated in the E strain during scotophase (transcript a: log FC 4 h = 1.363; transcript b: log FC 1.3 h = 1.023, 4 h = 1.370), while an ecdysone receptor (comp16498, log FC 1.3 h= 1.829, 4 h = 1.583) and ecdysteroid-regulated protein (comp28452, log FC 1.3 h = 1.167; 4 h = 1.265) were upregulated Z strain throughout scotophase (Figure S1).
Thirteen olfactory receptors (ORs) previously identified in Ostrinia were detected in our female head/antennal transcriptome (ORs numbered following [71]).Of these, five were differentially expressed between strains at one or more time points.OR15 was more highly expressed in E strain at photophase (comp28365, log FC = 1.43) and hour 4 (log FC = 1.27).OR37 (comp15171) was more highly expressed in Z strain at photophase (log FC = −8.59)and hour 4 (log FC = −2.91).OR2 (the olfactory co-receptor, comp34293) was more highly expressed in E strain at hour 4 (log FC = 1.17).OR38 was more highly expressed in Z strain at photophase (log FC = −1.80).OR19 was more highly expressed in the Z strain at hour 1.3 (comp12551, log FC = −2.06).
Finally, we identified ten transcripts that fit a predicted pattern of significant differential expression between the strains at both hour 1.3 and hour 4. Of these, six were successfully annotated, including four with immunological functions: two transcripts of cecropin A, and two of attacin.The remaining annotated transcripts were identified as an IQ motif protein and a prickle-like protein.
A total of 6,430 transcripts (15.75%) were annotated with D. melanogaster gene IDs.In the between-strain contrasts, enriched GO terms found amongst all of the assembled transcripts that were significantly differentially expressed included catalytic and hydrolytic activity during photophase, structural constituents at hour 1.3, and thiolester hydrolytic activity and the regulation of triglyceride metabolism at hour 4. When separated by strain, additional enriched GO terms were identified.In the E strain, light absorption was enriched at hour 1.3 and pigmentation-related terms and melatonin defense response regulation were enriched at hour 4 (Table S3).In the Z strain, the regulation of wound healing was enriched during photophase.
Fewer enriched GO terms were identified in within-strain contrasts.Within the E strain, enrichment of voltage-gated ion activity was found in transcripts significantly differentially expressed between hour 1.3 and 4. Within the Z strain, the only GO terms enriched between any of the contrasts were microvillus and actin-based cell projection in the photophase vs. hour 4 contrast.We also performed GO enrichment on 157 annotated joint outliers from the diapause time course.This analysis found that only macromolecule depalmitoylation, pyrimidine nucleobase metabolic process, and palmitoyl-(protein) hydrolase activity were significantly enriched.

Discussion
Our results indicate that differences in daily mating time are coupled with differences in the seasonal timing in ECB.E strain females that break diapause early in the season mated 2.6 h earlier during scotophase than Z strain females that break diapause later in the season.This coupling will enhance the total isolation among these sympatric ECB strains.Dopman et al. (2010) estimated that seasonal timing differences in New York can limit up to 66% of inter-strain mating [4], based on 10 years of trapping data.We estimate daily mating time further reduces inter-strain mating by 39%.Daily timing is a weaker reproductive isolating barrier on its own when compared to seasonal timing.However, cumulative reproductive isolation is enhanced when later acting barriers limit the portion of gene flow remaining after earlier barriers [4].When daily mating time is coupled with seasonal timing differences, it could limit an additional 39% of the 34% of gene flow not limited by seasonal timing.Together, these two barriers can limit 79% of inter-strain mating and lead to considerable reproductive isolation in ECB populations.
We found differences in the daily mating time measured among ECB populations in this study and those that were measured in the historical Liebherr & Roelofs study [45].Our E strain population originated from a site 50 km distant from the historical E site and our Z strain site was 501 km distant from the historical Z site.In our E and Z strain populations, we found a weaker contribution of mating time to reproductive isolation (39%) than in the historical study (52%) [4].This suggests that the daily timing and the amount of reproductive isolation that it causes may vary geographically or across years in ECB.
Coupling of daily and seasonal mating time could occur due to selection against recombinant genotypes or pleiotropic alleles influencing both traits [3].We evaluated the hypothesis that genes within the circadian clock contribute to both seasonal and daily timing differences.Several circadian clock genes were significantly upregulated in one strain relative to the other during our nightly mating time course (Figures 4 and 5).The central regulator of the circadian clock is the CLK/CYC heterodimer which binds to E-box enhancers of various circadian genes, increasing their expression [75].At the end of photophase, we found levels of cyc, and to a lesser extent, clk, were higher in the E strain.Vri represses clk and its expression was lower in the E strain at the same time point [76].In our correlation analysis, we found that vri expression was significantly negatively correlated with clk, cyc, and per expression across multiple time points.The expression of slimb was positively correlated with cyc across time points and slimb was also more highly expressed in the E strain during photophase.SLIMB degrades PER, which should result in the decreased inhibition of CYC [77].This suggests that earlier mating in the E strain may be tied to differences in levels of cyc, clk, vri, and slimb at the end of photophase.Further, it suggests that the circadian clock of the E strain may have a shifted phase relative to that of the Z strain, with peaks in the expression of circadian genes occurring earlier in a given 24-h period, a hypothesis that could be tested by comparing patterns of activity and gene expression between strains over the entire 24-h cycle.Thus, changes in expression of these core circadian clock genes during late photophase may be triggering differences in the timing of pheromone biosynthesis and release, leading to differences in mating time.
Work in other species supports vri as an interesting candidate gene for shifting the daily timing of mating behavior.Selected early and late Drosophila eclosion chronotypes differ in the timing of the peak expression of vri by about 2.5 h during early scotophase [78].Both vri and slimb are known to be regulated by ecdysone, the insect steroid hormone that plays a key role in reproduction and gamete release [79][80][81].Finally, the pattern that we document here, of increased vri expression at the end of photophase being associated with later mating, was also found using RT-qPCR between corn and rice strains of the fall armyworm (Spodoptera frugiperda) that differ in mating time [35].Hanniger et al. (2017) found that vri was on the same chromosome as the QTL for mating time in armyworms, suggesting that it may be a causal factor for changes in mating time [35].Thus, vri and alterations in the expression of cyc and clk may be a general mechanism of regulating daily allochrony in moths, as corn borers and armyworms are distantly related moth species (Ostrinia are in the Pyraloidea superfamily; Spodoptera are in Noctuoidea).
The circadian clock pathway was hypothesized to be a shared genetic pathway between daily and seasonal allochrony.The circadian clock in the brain acts as the pacemaker that synchronizes the body's endogenous clock [30,82].We compared genes that were differentially expressed in adult female heads to those differentially expressed in larval heads during diapause break in ECB [54].We did not find any core circadian clock genes that differed in both time courses.Per is a candidate gene for diapause timing in ECB because it shows amino acids changes among E and Z strains and trends toward higher expression in the E strain during diapause break (days 1 and 7) [54].In addition, per alleles are linked to the QTL for diapause timing and oscillate in frequency across latitude with voltinism (generation number) [52,53].Pdp1 is also differentially expressed during diapause break [54] and has been associated with diapause timing in several other insects [83], but neither per nor pdp1 were differentially expressed across the first 4 h of the night.Further characterization of the expression of circadian clock genes is needed determine whether this result is robust across 24-h time periods and rule out that differential expression of per does not occur later in the night.Future work can also scan for naturally segregating genetic variation in regulatory or coding sequences between ECB strains in these important timing genes, as we would expect causal genes to contain mutations in cis-regulatory or protein-coding regions.Thus, while the circadian clock pathway might be involved in both seasonal and circadian timing, different genes within these pathways are associated with different types of timing.
Although associated circadian loci for daily and seasonal timing may be different, they may interact with common downstream physiological pathways that lead to phenotypic shifts in the biological rhythm.We did find evidence for shared differential expression of genes in the juvenile hormone pathway.Juvenile hormone (JH) and ecdysteroids play key roles in triggering developmental transitions and reproduction [84].In moths, JH stimulates female pheromone release [85] and is also involved in maintaining winter dormancy [86,87].Circadian clock proteins directly regulate a number of genes involved in hormone production, including the takeout family of JH binding proteins [88][89][90].Expression of takeout alters the timing of male courtship behavior in Drosophila [25,90].Six transcripts from the takeout family were significantly differentially expressed between strains at several night time points and four of these were also differentially expressed between strains during diapause break.For example, comp30593 is upregulated both in E strain when compared to Z strain at hour 1.3 (Figure S1) and in E strain as compared to Z on day 7 after diapause break [54].These results suggest that different changes in the circadian clock pathway activate the same downstream hormonal module to break diapause or induce mating [6].Such modularity is common among developmental genes [91] and may constrain the evolution of gene expression within modules [92].
Melatonin is also thought to be important in linking circadian rhythms to behavior, possibly by influencing the release of JH or ecdysteroids [84].aaNAT is a gene that is regulated by clk and cyc and stimulates melatonin production [93,94].In our data set, one aaNAT transcript was upregulated in the E strain (at photophase and hour 4).aaNAT is known to stimulate the release of prothoracicotropic hormone (PTTH).We observed that PTTH was also upregulated in the E strain during photophase.Expression differences of aaNAT and PTTH may represent another instance of daily and seasonal timing involving similar downstream pathways, as these genes are involved in seasonal diapause termination in the moth Antheraea pernyi [94,95].
We found little evidence for pleiotropy as a cause of coupling between daily and seasonal allochrony in ECB.Evidence for physical linkage among daily and seasonal differentially expressed candidate genes is also minimal.Per, cyc, and pdp1 are all located on the sex (Z) chromosome, while vri and slimb are autosomal [52,96,97].Per is over 40 centiMorgans (cM) away from pdp1 and cyc.This lack of pleiotropy or genetic linkage in ECB is consistent with the work in pitcher plant mosquitoes that found that seasonal and daily timing traits evolve independently under artificial selection in mosquitos [6,41].While such a genetic architecture suggests that future evolution of different types of allochrony could occur independently in ECB, the common activation of downstream hormone modules creates a susceptibility to pleiotropy and could produce constraints on the future evolution of these traits [8].
Given that coupling among different forms of allochronic isolation in ECB is not due to shared or physically linked genes, coupling is likely maintained by selection and occurs because recombinant individuals have low fitness [3].For example, reproductive output may be reduced in individuals who emerge from diapause in early spring, but reproduce later in the night (or individuals that emerge from diapause in the summer but reproduce early in the night) because they have difficulty finding mates.Future work using recombinant lines to decouple these phenotypes could test these hypothesized fitness effects as well as characterize their genetic architectures.Coupling of timing traits to other forms of reproductive isolation may further strengthen total reproductive isolation among sympatric ECB strains.Sexual isolation among E and Z strains occurs due to female pheromone production and male pheromone response, limiting ~79% of inter-strain mating [4].Timing loci are not genetically linked to the autosomal pheromone gene (pgFAR) and the major QTL for male response on the Z chromosome is far from per (11cM) and pdp1 (21 cM) [51,96,98].Thus, natural selection against hybrid phenotypes is likely to cause coupling of daily allochrony, seasonal allochrony, and sexual isolation in ECB, leading to almost complete reproductive isolation (95%) in geographic locations where populations differ in all three forms of isolation.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/9/4/180/s1, Figure S1: Heatmap of relative expression of downstream candidate transcripts in each contrast, Table S1: Summary statistics of transcriptome assembly in Trinity, Table S2: Overlapping differentially expressed genes between daily and seasonal time courses, Table S3: Summary of GO term enrichment among significantly upregulated genes within each strain at each time point.

Figure 1 .
Figure 1.Simplified pathway of the core Lepidopteran circadian clock.Arrows represent activation.Vertical lines represent inhibition.Gray lines indicate dimer formation.During the day, light activates the blue light photoreceptor CRYPTOCHROME1 (CRY1).CRY1 activation leads to SHAGGY (SGG) tagging TIMELESS (TIM) for degradation by JETLAG (JET).The PERIOD (PER) monomer is tagged for degradation by SLIMB when not bound to TIM.DOUBLETIME (DBT) regulates PER activity through phosphorylation.The PER/TIM heterodimer forms a complex with CRYPTOCHROME2 (CRY2) and enters the nucleus, where PER and CRY2 inhibit the activity of the CLOCK/CYCLE heterodimer.CLK/CYC binds to E-box enhancers to activate transcription.A second feedback loop is comprised of VRILLE (VRI) and PDP1, which bind to the V/P box of the clk promoter cyclically.VRI acts as repressor of clk transcription but PDP1 activates clk transcription.Modified from [29].

Figure 2 .
Figure 2. Schematic of specific contrasts tested in differential expression analysis.Within strain, comparisons were made between each pair of time points (photophase vs. 1.3 h of scotophase, photophase vs. 4 h of scotophase, and 1.3 vs. 4 h of scotophase).Between-strain, comparisons were made at each time point (BE vs. UZ).Numbers shown are total numbers of differentially expressed transcripts in each contrast.

Figure 3 .
Figure 3. Mating frequencies of the European corn borer E and Z pheromone strains during scotophase.Shown are the number of females mating for the first time during each hour of scotophase, for historical field populations with unknown seasonal timing [45] from Aurora, NY (E) and London, Ontario (Z), and contemporary colony populations: bivoltine E (Geneva, NY, first mating flight in May), univoltine Z (Bouckville, NY, first mating flight in July).Dashed vertical lines represent median mating times in contemporary populations (E = 1.3, Z = 4.0), and mean mating times for historical populations (E = 6.8,Z = 5.1).Contemporary populations show a significant difference in mean mating time between strains (p = 0.0038).When mating times are grouped into 1-h intervals, both of the strains showed significant differences between historical and contemporary mean mating times (E: historical = 7.33, contemporary = 2.81, p < 0.0001; Z: historical = 5.58, contemporary = 4.28, p = 0.014).

Figure 4 .
Figure 4. Heatmap of relative expression of genes in the circadian pathway and related photoreceptors in each contrast.Values are log FC.Photophase, 1.3 h, and 4 h contrasts are between strain with red indicating upregulation in BE (early daily and seasonal mating), and blue indicating upregulation in UZ (late daily and seasonal mating).For within strain/between time point contrasts the red color indicates upregulation at the first time point, and blue color indicates upregulation at the second time point.Log FC were capped to a minimum of −2 and a maximum of 2 for visualization purposes.Asterisks indicate genes that are significantly differentially expressed (q-value < 0.05) for a given contrast.

Figure 5 .
Figure 5. Trajectory of circadian transcript expression levels through time for transcripts that were significantly differentially expressed between strains.Each line represents a single component, with normalized counts at each time point averaged across libraries.Data collection time points were 1 h before scotophase, 1.3 h into scotophase, and 4 h into scotophase.Red indicates values for UZ, blue indicates BE.

Figure 6 .
Figure 6.Pairwise correlations among 22 transcripts for circadian genes.Only significant (p < 0.05) correlations shown.Blue indicates positive correlations, red indicates negative.Transcripts listed by gene symbol and comp number.