Evolution of multiple phosphodiesterase isoforms in stickleback involved in cAMP signal transduction pathway

Phylogeny and synteny among stickleback PDE1C genes

The maximum likelihood (ML)/Bayesian molecular phylogeny of chordate PDE1A and PDE1C (Figure 1) showed that seven of the duplicated PDE1C genes in stickleback (PDE1Cb1-b7) arose from repeated duplications that occurred after the divergence between stickleback and pufferfish (see Figure 1B). In addition, the teleost (zebrafish, medaka, stickleback, and pufferfish) PDE1C genes were grouped within two major clades, PDE1Ca and PDE1Cb, with relatively high support values (see black circle in Figure 1A; LR-ELW [the expected-likelihood weights applied to local rearrangements of tree topology] [24] edge support = 90%; Bayesian posterior probability = 95%). This branching pattern suggests that the PDE1C gene was duplicated in a teleost ancestor after its split from the tetrapods (see double circle in Figure 1A), possibly through the 3R-WGD [18, 19]. This possibility was examined by investigating the synteny of the surrounding genes of the PDE1C loci in teleosts and tetrapods as described below.

To clarify the genomic events that generated the two paralogous genes in teleosts (PDE1Ca and PDE1Cb) and multiple PDE1Cb genes in stickleback, we investigated the genomic regions around the PDE1C loci. We found conserved synteny between the PDE1C locus in tetrapods (human, frog, and chicken) and the PDE1Ca locus in medaka, stickleback, and pufferfish (described as "conserved synteny [CS]-1," Figure 2). One more homologous region corresponding to the CS-1 was found in medaka, pufferfish, and stickleback (described as "CS-2," Figure 2). CS-1 and CS-2 were considered doubly conserved syntenies derived from the 3R-WGD, which corresponds to a single tetrapod chromosome segment. However, CS-2 appeared to have lost the putative region that contains PDE1Cb. Zebrafish PDE1Ca, pufferfish PDE1Cb, and stickleback PDE1Cb1-b7 were found within other chromosome segments, implying secondary translocation(s) of these PDE1C genes after the 3R-WGD. The pufferfish PDE1Cb region partially corresponded to human chromosomes 1 (94.6–99.6 Mb) and 10 (16.6–17.3 Mb) [see Additional file 1: Table S1]. However, we were unable to find human chromosome regions corresponding to the zebrafish PDE1Ca and stickleback PDE1Cb1-b7 regions. Among the stickleback PDE1Cb1-b7 regions, the surrounding genes were very similar (see "Teleost PDE1Cb" in Figure 2), implying that repeated segmental duplications generated the multiple PDE1Cb regions in stickleback. We found one more ninth PDE1C locus in stickleback (provisionally named PDE1Cbx) [see Additional file 1: Table S1], whose sequence was 5'-and 3'-truncated and partial (246 residues in total whereas the others have >395 residues) probably due to incompleteness of the sequence contig where the locus was located (scaffold_809: only 11,237 base pairs [bp] in total). We could not conclude whether this is pseudogene or not, and thus we did not include this locus in our analysis.

Our overall results from the phylogenetic and synteny analyses clearly revealed the evolutionary relationships among PDE1C genes in the bony vertebrates and the evolutionary origin of multiple PDE1C genes in stickleback. This provided the basis for our subsequent analysis on the molecular evolution of the multiple stickleback PDE1C genes.

Molecular evolution of multiple PDE1C genes

The multiple PDE1C genes of stickleback were analyzed by ML-estimation of the nonsynonymous to synonymous substitution rates (d_N/d_S = ω) during evolution, which is a possible indicator of adaptation at the protein sequence level [25–27]. In this analysis, we tested whether some portion of the PDE1Cb sequences shows ω >1, which is the signature of adaptive amino acid changes, by comparing maximum-likelihood values of simple evolution model having fewer ω parameter(s) (M0, M1, and M7 in Table 1) with those of more complex model having more ω parameters, some of which were allowed to be >1 (M2, M3, and M8 in Table 1; for details, see Methods). The likelihood ratio test (LRT)-1, the comparison between M0 with M3, implied that the stickleback PDE1Cb genes (PDE1Cb1-b7) were under positive diversifying selection with regard to their protein sequences (ω₂ = 17.83). However, the LRT-2, LRT-3, and "branch-site test", the comparison between M1 and M2, M7 and M8, and M2 and M3, respectively, did not support this implication (ω₂ = 0.15, ω₂ = 0.45, and ω₂ = 0.08, respectively; see Table 1). The LRT-1 (M3 model) detected five codon sites under positive selection, including sites 76–78, 82, and 83 (indicated by stars in Figure 3). On the other hand, the LRT-2, LRT-3, and branch-site test (ω >1) did not identify these individual sites as being under positive selection.

The positively selected sites inferred by the LRT-1 described above were not located on the known active sites of the enzyme or on specific domains or motifs of the PDE proteins (Figure 3). The active enzyme sites that are generally conserved in PDE proteins (indicated by black shading in Figure 3) [28] were also conserved among the stickleback PDE1Cb genes, with exceptions at sites 238 in PDE1Cb4, 417 in PDE1Cb4, and 443 in PDE1Cb5. Moreover, the known PDE-specific protein domains were successfully detected in all stickleback PDE1Cb genes (indicated as solid boxes in Figure 3) via queries to protein domain databases Pfam [29] and PROSITE [30], suggesting that they conserve PDE function. On the other hand, the sequence region containing the inferred positively-selected sites (sites 66–88; indicated as dashed boxes in Figure 3) yielded no hits in these databases. Those positively-selected sites were distributed on the opposite side of a substrate-binding pocket of this enzyme [see Additional file 1: Figure S2]. In addition, this region did not correspond to the important variant region of the three functional splicing variants reported in mouse PDE1C [21].

Spatial expression of the multiple PDE1C genes

Spatial expression patterns across tissues were investigated and compared among the multiple PDE1C and G(olf) genes in adult stickleback (Figure 4). Reverse transcription polymerase chain reaction (RT-PCR) analysis showed that seven of the eight PDE1C genes were co-expressed with G(olf) in the nose of stickleback. Among these seven genes, five (PDE1Ca, PDE1Cb3, b4, b5, and b7) were strongly expressed, and two (PDE1Cb2 and b6) were weakly expressed in the nose. In addition, those stickleback PDE1Cs were roughly divided into two groups in terms of their overall expression patterns. One group included the PDE1Cs that were expressed strongly in a few tissues ("group I"; 0–3 tissues; PDE1Ca, PDE1Cb1, b2, b6, and b7) and the other group included the PDE1Cs that were strongly expressed in several tissues ("group II"; 4–9 tissues; PDE1Cb3, b4, and b5). Among the "group I" genes, the expression of PDE1Ca was clearly detected in the nose, brain, and intestine and was similar to that of G(olf), which was strongly expressed in the nose, gill raker, and intestine. On the other hand, PDE1Cb1, b2, b6, and b7 seemed to have lost the strong co-expression with G(olf) in the nose (PDE1Cb6), intestine (PDE1Cb7), or both tissues (PDE1Cb1 and b2). Among the "group II" genes, PDE1Cb4 was strongly expressed in the nose, brain, skin, and skeletal muscle, and PDE1Cb3 and b5 were strongly expressed in almost all tissues.

In silico simulations for the multiple PDE1C genes

We examined the possible effects of increased gene/product dosage from the multiple PDE1C genes in stickleback on the output (depolarization) of the OT pathway using in silico simulations. Figure 5 depicts a schematic diagram of the OT simulation model constructed according to the KEGG pathway database (panel A) [31], and a representative result of the simulation (panels B and C; the case of "single-PDE1C [threshold = 1]" model). Figure 5B shows that the odorant and olfactory receptor (OR) produced OR-odorant complex (violet line) that stimulated the G(olf) (blue line), led to depolarization occurring on the simulation time scale of 0–10 (orange line). The depolarization is blocked by activation of the PDE1C. Figure 5C shows that the activated G(olf) increased the concentration of Ca²⁺ (green line), which stimulates PDE1C (red line) via the mediation of calmodulin (CaM; see Figure 5A). When the activity level of the PDE1C has gone over its firing threshold (set to 1 here; see Figure 5C), the activated PDE1C decomposed cAMP (see Figure 5A), and finally terminated the depolarization (see Figure 5C). Although this preliminary model approach may be simple, we thought that this model would provide a basic framework to analyze the effect of an increased dosage of a specific element; that is, the PDE1C gene. We compared the depolarization signals of the "single-PDE1C" model and "multiple-PDE1Cs" model. Since five PDE1Cs were strongly expressed, and two were weakly expressed in the nose (see Figure 4), we set the number of PDE1C circuits equal to six in the "multiple-PDE1Cs" model as an approximate representation of the estimated gene dosage in the OT.

Our OT model simulations implied that the increase in the number of PDE1C circuits affects the longevity of the depolarization signal (Figure 6). In the case of the single- and multiple-PDE1C models with a PDE1C threshold = 0, the resultant depolarization signals were weak (intensity = 0.05–0.2) and short-lived (longevity = 5) when the number of PDE1C circuits was one or six (Figure 6A and 6B). In these simulation models, PDE1C seems to work rapidly after odorant stimulation, leading to an instant block of depolarization by the PDE1C circuit regardless of whether multiple PDE1Cs were involved. When the PDE1C threshold = 1, on the other hand, the resultant depolarization signals became more intense (intensity = 0.3–0.5; Figure 6C and 6D), and the "multiple-PDE1C [threshold = 1]" model yielded an elongated signal of depolarization (longevity = 17–25; Figure 6D) compared to the "single-PDE1C [threshold = 1]" model (longevity = 8–14; Figure 6C).

We also found that the depolarization longevity of the "multiple-PDE1C [threshold = 1]" model was further extended by limiting the availability of Ca²⁺, which is an upstream regulator of the PDE1C circuit. Ca²⁺ dosage limitation affected the single- and multiple-PDE1C [threshold = 1] models similarly in terms of their depolarization intensities, which were reduced to 0.2–0.4, compared to the model with unlimited Ca²⁺ (Figure 6E and 6F). In the "multiple-PDE1C [threshold = 1] + Ca²⁺-limited" model, the longevity of the depolarization signal was greatly elongated (longevity = 25–30; Figure 6F), compared to the "single-PDE1C [threshold = 1] + Ca²⁺-limited" model (longevity = 10–15; Figure 6E).

Since the PDE1C is activated by Ca²⁺-activated CaM, it is possible that the finite Ca²⁺ dosage invoked competition among increased PDE1Cs, resulting in a delay in blocking of the depolarization. This situation likely led to a positive-feedback circuit (the processes of 21–25 shown in Figure 5A), also activated by upstream Ca²⁺, to be more prevalent in the OT system. Therefore, the depolarization longevity was extended in the "multiple-PDE1C" models (see Figure 6D and 6F). This possibility was supported by an in silico mutation analysis (Figure 6G and 6H), in which we knocked out the positive feedback circuit (processes 21–25 and the related elements) in the simulation model. In this "multiple-PDE1C [threshold = 1] + positive-feedback loop defeated" model, the resultant depolarization signals were not greatly elongated (longevity = 13–17; Figure 6H), and the longevity was similar to that of the "single-PDE1C [threshold = 1]" model (Figure 6C).

Evolutionary origins of the multiple PDE1C genes in stickleback

Molecular phylogenies (Figure 1) and conserved synteny (Figure 2) analyses showed that the teleost PDE1C was duplicated through the 3R-WGD, generating two isoforms, PDE1Ca and PDE1Cb, and that PDE1Cb was repeatedly duplicated in the stickleback lineage 6–7 times through chromosome segment duplications. Such multiple gene copies could be generated by fewer than six duplications especially in tandemly located genes; however, this seems not to be the case of the stickleback PDE1Cb genes, which are located on different scaffolds respectively. The stickleback PDE1Cb loci and the orthologous pufferfish PDE1Cb locus appears to have translocated from CS-2, which is one of the pair of conserved chromosome regions derived from 3R-WGD. The corresponding PDE1Cb locus in medaka and zebrafish has been lost (Figure 2). Since no synteny was found between stickleback PDE1Cb loci and the pufferfish PDE1Cb locus, the translocation of the PDE1Cb gene may have occurred independently in the two lineages. Alternatively, the PDE1Cb translocation may have occurred once in an earlier ancestor, and either the stickleback or pufferfish experienced an additional PDE1Cb translocation. After the translocation(s), the stickleback PDE1Cb gene along with the neighboring tandemly repeated "unknown" genes (denoted by Ensembl transcript IDs in Sca. 37, 94, 134, 154, 188, 215, and 223 in Figure 2) appeared to have undergone lineage-specific expansion (LSE) [32] by chromosome segment duplications of >100 kbp. Such plasticity in the evolution of the 3R-WGD-derived duplicates, such as an independent loss or persistence among teleost lineages, is also shown in a recent study on androgen receptors [33]. The persistence of these receptors in various Percomorph fishes is suggested to be responsible for their neofunctionalization in ligand-binding activities that may be associated with plasticity of sex determination in those teleosts. It is also interesting to explore how the multiple stickleback PDE1Cb loci presented here have been retained.

The various stickleback chromosome segments that contain PDE1Cb1-b7 (see Figure 2) may have been retained because of advantages related to the LSE of PDE1Cb and/or the "unknown" genes mentioned above. LSEs are thought to have played an important role in the proteome evolution of multicellular eukaryotes [32], particularly in the evolution of proteins involved in ligand recognition, pathogen resistance, etc [34–36]. The function of the "unknown" genes is not annotated in the Ensembl stickleback genome. However, at least a part of their corresponding human genes is suggested to function as an angiotensin II/arginine vasopressin (AII/AVP) receptor-like protein responsible for an autoimmune disease (cold autoinflammatory syndrome) (found by a BLASTP search; data not shown). The "unknown" stickleback genes may also have immune-related functions. The repeated segmental duplication of these "unknown" genes (see Figure 2) might be advantageous in stickleback responses to diverse pathogens or other ligands. This possibility should be examined further. On the other hand, the stickleback PDE1Cb genes are intact and do not have stop codons or frameshifts in their coding sequences. Thus, they should also have been allowed to be highly duplicated. Accordingly, the PDE1Cb genes in stickleback may also have important phenotypic consequences; for this reason, we analyzed their molecular evolution, gene expression, and possible gene dosage effects, as discussed below.

Sequence and expression evolution of the multiple stickleback PDE1C genes

The sequence evolution analysis of the multiple PDE1Cb genes in stickleback did not provide strong evidence for neofunctionalization [1, 2] in their coding proteins. The known active sites and specific domains/motifs of PDE enzymes were highly conserved among the multiple PDE1Cb genes (Figure 3). In addition, although the ML-estimation of the d_N/d_S ratio implied adaptive sequence evolution in the case of LRT-1 (Table 1), the estimated positively-selected amino acid sites are located on the opposite side of a substrate-binding pocket of the PDE1 [see Additional file 1: Figure S2] and these positively-selected sites are included in the sequence region for which the possible function is not reported (see Results). From these observations, we can hardly conclude that the multiple PDE1Cb genes were maintained through the acquisition of new enzyme functions. This implies that evolutionary mechanisms other than neofunctionalization in coding sequence may have retained the multiple stickleback PDE1C genes.

Gene expression analysis (Figure 4) showed that the PDE1Ca and relatively "basal" PDE1Cb genes (see Figure 1B; designated as "group I" in Figure 4) were strongly expressed in a few (0–3) tissues, and their expression patterns were similar to G(olf). Since PDE1C and G(olf) are involved in the OT [20, 21, 37], expression patterns of these basal ("group I") PDE1Cb genes seem to represent ancestral states of the teleost PDE1Cb. This might be verified by further expression analysis on the PDE1Cb from many other teleosts in the future. On the other hand, the "derived" PDE1Cb genes (designated as "group II" in Figure 4) were strongly expressed in many more tissues, implying that they acquired expression in new organs. These PDE1Cb genes, which seem to be neofunctionalized in terms of expression patterns, may be under functional adaptation to their new subcellular environments. Furthermore, they might experience neofunctionalization through the acquisition of new mutations in the future.

In addition, the number of PDE1C loci involved in the stickleback OT seems to be increased, as suggested by the fact that at least seven PDE1C isoforms were co-expressed with G(olf) in the stickleback nose (Figure 4). Such an increase in PDE1C gene products may have some functional significance in the performance of the OT. However, directly addressing this question is difficult within the scope of this study, because it would require gene knockout, neurophysiological, and behavioral analyses. Furthermore, gene expression levels of the PDE1C would be assessed more quantitatively by a real-time PCR approach. Such an approach will also useful to explore the existence of gene number variation of PDE1C within and between populations of stickleback, which includes a variety of ecomorphs (e.g., marine, anadromous, and freshwater populations) [22, 23] in the future. To obtain and describe the primary predictions of the phenotypic effects from increased PDE1C dosage, we used in silico pathway simulations, as discussed below.

Possible effect of increased PDE1C products on olfactory transduction

The effect of increased PDE1C dosage on stickleback OT was surveyed using in silico pathway simulation. The simulation was based on limited information and knowledge of OT, and the predictions resulting from the simulation should be empirically evaluated. Regardless, such approaches may provide insight into the evolutionary significance of the multiple PDE1C genes in the stickleback.

According to the results of the OT simulation, the increased PDE1C dosage extends the longevity of the depolarization signal of the olfactory receptor neuron (Figure 6). PDE1C is involved in the negative feedback circuit of the OT and decomposes cyclic adenylic acid (cAMP) and eventually terminates the depolarization depending on upstream Ca²⁺/CaM [20, 21, 38] (see Figure 5). Accordingly, the increased PDE1C products may compete with each other for binding to these upstream Ca²⁺/CaM molecules, and consequently delay the termination of depolarization (Figure 6C and 6D). This proposed delay mechanism seems to be supported from our additional simulation results that Ca²⁺/CaM limitation further extends the longevity of the depolarization signal (Figure 6E and 6F). Competition among the increased PDE1Cs involved in the negative feedback circuit may result in positive feedback in the OT system to be dominant. This situation would also elongate the depolarization signal. Our in silico mutation-simulation analysis supported this hypothesis (Figure 6G and 6H). That is, a defect in the positive feedback circuit antagonized the effect of the increased PDE1C dosage, suppressing the depolarization elongation. To summarize, our simulation analyses suggested that increased PDE1C dosage induces competition among PDE1Cs and causes the positive feedback loop to be dominant in the OT system, resulting in a delay in the termination of the output depolarization of the olfactory receptor neuron.

It is proposed that an extension in the duration of olfactory signals is associated with the territorial ecology of the house mouse Mus domesticus [39]. In M. domesticus, the male scent mark contains lipocalin proteins called major urinary proteins (MUPs). The MUPs bind with the semiochemical molecules and release them gradually, which eventually extends the longevity of the odor signal. This makes it difficult for other male mice to tell whether the odor signals come from scent marks or a territorial male, and the other males are hesitant to approach the territorial zone because the scent-mark odor has an aggressive message. This is thought to be evolutionarily advantageous for both the territorial and the other individual because potential male invaders reduce their risk of damage or death due to conflict [39]. Although this phenomenon and the underlying mechanisms in house mice are different from those proposed for the stickleback, the multiple PDE1C genes may play an important role in stickleback ecology and behavior. This speculation may be appealing when considering that sticklebacks also hold territories where they build nests and reproduce.

Of course, these hypotheses should be empirically examined. Additionally, the application of the vertebrate OT system described in the KEGG [31] to this stickleback study should be verified in future research. However, we propose that the evolutionary significance of multiple gene duplicates may be evaluated more comprehensively using available biological information and analytical tools such as whole genome sequences, pathway data, and in silico simulation software, as was attempted here. Such a comprehensive approach would be particularly favorable for questions that are not entirely addressed using the molecular evolutionary analysis of a particular gene/protein. For example, this approach could be used in testing the effects of increased gene dosage in signal transduction pathways. With improvements in the pathway models and their parameters, the in silico pathway simulation, which can perform a synthetic analysis of molecular dynamics for multiple gene products and other biomolecules, will become one of the most powerful approaches in understanding complex macro-phenotypic evolution.

Phylogenetic analysis

To analyze the evolutionary origins and relationships of the multiple PDE1C genes in the stickleback, we performed a phylogenetic analysis of PDE1C and its closely related PDE1A genes from eight chordate species (human, chicken, frog, pufferfish, medaka, stickleback, zebrafish, and ascidian) using available whole-genome sequence data. The primary sequences of the PDE1 genes were gathered via queries to the Ensembl genome database [40] and its Ortholog Predictions section. We confirmed that no additional PDE1A and PDE1C genes existed in these genomic sequences using BLAST searches (E-value cut-off of < 10^-3). When a partial sequence was detected in the Ensembl database, we predicted the full-length coding sequence from the genomic sequence using WISE2 [41]. The corresponding PDE1 of the sea lamprey Petromyzon marinus was searched using BLAST against the UCSC Genome Browser Database [42]. However, we were unable to find the full length sequence, which can be used as an outgroup for the phylogenetic analysis. The species names and Ensembl IDs of the analyzed PDE1A and PDE1C genes are provided in supplementary table S2 [see Additional file 1].

The nucleotide sequences of the PDE1A and PDE1C genes from the seven vertebrates and ascidian (outgroup) were aligned using ClustalW [43]. The alignment was manually adjusted according to the amino acid sequences using MacClade ver. 4.06 [44]. After removing the gaps, 930 bp of the PDE1 coding region were phylogenetically analyzed using maximum likelihood (ML) and Bayesian methods in TREEFINDER (version June, 2007) [45, 46] and MrBayes (version 3.0b4) [47] under the GTR + I + Γ model [48], which was selected as the best-fitting model for nucleotide substitution by hierarchical LRT (hLRT) [49]. The ML analysis was assessed using 1000 replications of the LR-ELW edge support tests [24]. Bayesian posterior probabilities of the phylogeny and its branches were determined from 9901 trees. The re-aligned teleost PDE1C genes (1248 bp), excluding the partial zebrafish PDE1Ca (1185 bp), were analyzed using the ML method under the TrN + Γ model of nucleotide substitution [50], which was chosen by the hLRT. Bayesian method was not applied to this analysis, because MrBayes does not allow to use the TrN model.

Synteny analysis

To investigate the chromosomal/genomic events that generated multiple PDE1C genes in the stickleback, genomic regions around the stickleback PDE1C loci were investigated and compared to those of human, chicken, frog, pufferfish, medaka, and zebrafish. Physical mapping data nearby each PDE1C locus were obtained from the Ensembl database [40]. An orthology of the neighboring genes [see Additional file 1: Table S1] within each species was examined according to descriptions in the Orthologous Prediction section of Ensembl database. Phylogenetic relationships of a part of neighboring genes of the stickleback PDE1Cb loci were analyzed by ML method according to the procedure described above [see Additional file 1: Figure S1]. The genomic location data of the genes near the PDE1C genes were used to rebuild the synteny maps.

Molecular evolutionary analysis

To examine whether the multiple PDE1C genes in the stickleback were subjected to diversifying selection in terms of their amino acid sequences, we analyzed the d_N/d_S ratio (ω) using ML inferences of the ω values in codeml [25]. The re-aligned teleost PDE1C genes (1248 bp; 416 codons) were analyzed, excluding the partial zebrafish PDE1Ca (1185 bp) gene. An ML tree of these teleost PDE1C genes (shown in Figure 1B) was used as a reference tree. In this analysis, a simple ML model was compared to a more complex model having more parameters to obtain an adequate ω estimate. Statistical significance of the comparisons was assessed by LRT with degrees of freedom equal to the differences in the number of free parameters between the two models. First, the one-ratio model (M0), which assumes a one ω (ω₀:0 < ω₀<1) for all codon sites, was compared to the discrete model (M3), which assumes three ω (ω₀, ω₁, and ω₂) for three site classes with proportions p₀, p₁, and p₂. Second, the neutral model (M1), which assumes conserved sites (0<ω₀<1) and neutral sites (ω₁ = 1) with proportions p₀ and p₁ (p₁ = 1-p₀) was compared to the selection model (M2), which assumes the conserved, neutral, and selection sites with proportions p₀, p₁, and p₂. The ω value of the selection sites (ω₂) was allowed to be greater than 1. Third, the beta model (M7), which assumes that the ω varies according to the beta distribution β (p, q), was compared to the beta and ω model (M8), which assumes an additional site class with a ω value >1.

In addition, we performed a "branch-site test" [26] to detect positive selection at individual codon sites, if it exists, along respective branches leading to the multiple PDE1C genes of the stickleback. For this purpose, we set branches connecting the multiple PDE1C genes of the stickleback as "foreground" branches. The other branches leading to the pufferfish and medaka PDE1C genes were considered "background" branches. To obtain the adequate proportion estimates of site classes and their ω values, the selection model (M2) and discrete model (M3) were compared on the basis of their log-likelihood scores (l) estimated using codeml. Individual codon sites were assessed in terms of their posterior probability to belong to the site class for which the ω value was allowed to be >1.

RT-PCR based co-expression analysis

For the RT-PCR experiment, we used two adults of the anadromous form of the three-spine stickleback. Live specimens were collected at Akkeshi Lake, Hokkaido, Japan, in May 2008, and were treated according to the ethical recommendations of the Ichthyological Society of Japan and the University of Tokyo. Total RNA was extracted from the lip, nose, gill raker, brain, heart, liver, intestine, skin, and skeletal muscle of fresh stickleback samples, using 1 ml TRIZOL reagent (Invitrogen). Residual genomic DNA was removed using DNase I (Takara), and 168 ng of the repurified total RNA from each tissue were reverse-transcribed into first-strand cDNA with oligo-dT adaptor primer using TaKaRa RNA PCR kit ver. 3.0 (Takara). Genomic DNA was also extracted from a piece approximately 5 mg in size of the caudal fins using the AquaPure DNA extraction kit (BioRad).

To assess the expression patterns of the PDE1C and G(olf) genes across tissues, reverse-transcribed cDNA from each tissue was subjected to PCR reactions with the GSPs (Table 2). The thermal-cycle profile was as follows: 1 cycle at 94°C for 2 min; 35 cycles at 94°C for 30 sec, 55°C for 30 sec, and 72°C for 30 sec; followed by 1 cycle at 72°C for 3 min. As a positive control for gene expression, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was amplified using a primer pair designed by Aoki et al. [51]. As a negative control, PCR amplification was also conducted for each RNA sample without a reverse-transcribed reaction. The amplified DNA fragments were separated by electrophoresis on a 2% LO3 agarose gel (TaKaRa; 35 min at 50 V; constant voltage setting), stained with ethidium bromide, and visualized under UV light. GeneRuler 100-bp DNA Ladder Plus (MBI Fermentas) was used as a size marker.

Pathway simulation

Since the exact parameter settings in modeling biological pathways are generally difficult because of the limited amount of available experimental data [54], we took a simple approach in constructing the OT model. In this approach, almost all parameters in the model were set to the default values (threshold = 0 and no priority) apart from some exceptions explained in the caption of Table 3. After confirming that the odorants successfully elicit depolarization of the olfactory receptor neuron modeled in the OT simulation, we assessed the intensity and longevity of the depolarization signals in the following situations: (i) PDE1C circuit is single (single-PDE1C model); (ii) the number of PDE1C circuits was increased according to the number of OT-involving PDE1C genes in the stickleback estimated by the RT-PCR-based analysis of co-expression with G(olf) (multiple-PDE1C model); and (iii) Ca²⁺ ion is limited by adding the degradation process of Ca²⁺ (Ca²⁺-limited model). To further examine the predicted effect of the multiple PDE1C genes on the negative feedback loops and depolarization signals, we performed an in silico mutation analysis, where the elements and processes involved in the positive feedback loop were knocked out. The cell system markup language (CSML) files of the OT models constructed and used in the analyses are available online [see Additional file 2].