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

Abstract

Background: Duplicate genes are considered to have evolved through the partitioning of ancestral functions among duplicates (subfunctionalization) and/or the acquisition of novel functions from a beneficial mutation (neofunctionalization). Additionally, an increase in gene dosage resulting from duplication may also confer an advantageous effect, as has been suggested for histone, tRNA, and rRNA genes. Currently, there is little understanding of the effect of increased gene dosage on subcellular networks like signal transduction pathways. Addressing this issue may provide further insights into the evolution by gene duplication.

Results: We analyzed the evolution of multiple stickleback phosphodiesterase (PDE, EC: 3.1.4.17) 1C genes involved in the cyclic nucleotide signaling pathway. Stickleback has 8-9 copies of this gene, whereas only one or two loci exist in other model vertebrates. Our phylogenetic and synteny analyses suggested that the multiple PDE1C genes in stickleback were generated by repeated duplications of >100-kbp chromosome segments. Sequence evolution analysis did not provide strong evidence for neofunctionalization in the coding sequences of stickleback PDE1C isoforms. On the other hand, gene expression analysis suggested that the derived isoforms acquired expression in new organs, implying their neofunctionalization in terms of expression patterns. In addition, at least seven isoforms of the stickleback PDE1C were co-expressed with olfactory-type G-proteins in the nose, suggesting that PDE1C dosage is increased in the stickleback olfactory transduction (OT) pathway. In silico simulations of OT implied that the increased PDE1C dosage extends the longevity of the depolarization signals of the olfactory receptor neuron.

Conclusion: The predicted effect of the increase in PDE1C products on the OT pathway may play an important role in stickleback behavior and ecology. However, this possibility should be empirically examined. Our analyses imply that an increase in gene product sometimes has a significant, yet unexpected, effect on the functions of subcellular networks.

Figure 1

Molecular phylogeny of vertebrate PDE1C. (A) Maximum-likelihood (ML) tree of the PDE1C and PDE1A genes in four teleosts, three tetrapods, and an ascidian, constructed under the GTR + I + Γ model with 930 base pairs (bp) of the coding region. Numbers indicate support values (percentages) from 1,000 LR-ELW edge support tests (left) and percent posterior probabilities from the Bayesian method (right). Single numbers indicate the LR-ELW edge support for the nodes, for which the Bayesian tree inference resulted in a different branching pattern. (B) ML tree of the teleost PDE1C genes constructed under the TrN + Γ model with 1248 bp of the coding region. Numbers indicate the LR-ELW edge support values (1,000 replications).

Figure 2

Conserved synteny around the PDE1C locus (loci) in tetrapods and teleost fishes. Triangles indicate gene loci and their direction of transcription. Doubly conserved synteny, which was derived from the teleost-specific genome duplication [18,19], is indicated by yellow shading and labeled "CS (conserved synteny)-1" and "CS-2." Orthologous/paralogous relationships among PDE1C genes are shown by solid magenta lines. The dashed lines indicate putative orthologous relationships predicted in the Ensembl genome database [40]. The solid yellow and green lines show phylogenetic relationships of neighboring "unknown" genes around stickleback PDE1Cb loci, which are estimated in the present study [see Additional file 1, Figure S1]. The PDE1Cx (Ensembl ID: ENSGACP00000001336) [see Additional file 1, Table S1] of the stickleback is located alone in a small contig (Scaffold 809; 12 kbp), and therefore has no synteny information.

Figure 3

Protein sequence alignment of multiple stickleback PDE1Cb and human PDE4B. The PDE tertiary structure and active enzyme sites (black shading) are reported for the human PDE4B (site numbers of active sites are according to ref. [28]). The signature domains of PDE protein are designated by solid boxes, and the sequence region that contains the ML-inferred positively-selected sites is designated by a dashed box. The stars above and gray shading indicate the inferred positively selected sites.

Figure 4

Spatial expression patterns of PDE1C and G(olf) in the stickleback. Semi-quantitative RT-PCR analysis was performed to assess expression levels and patterns of G-protein subunit alpha olfactory type (G [olf]) and multiple PDE1C genes in stickleback. Plus (+) and minus (-) signs indicate PCR assays using reverse-transcribed cDNA for each tissue type and assays using total RNA without reverse transcription (negative controls), respectively. The glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was amplified as a positive control. The overall expression patterns of the genes were essentially similar between the two stickleback individuals investigated.

Figure 5

Schematic view of a simulation model of the olfactory transduction (OT) pathway and its simulation performance. (A) OT simulation model constructed under the KEGG pathway database [31]. (B) Observed resultant oscillations of the odorant, OR, OR-odorant complex, G(olf), and depolarization under the "single-PDE1C [threshold = 1]" model. The x-axis indicates the simulation time scale, and the y-axis indicates the concentration of the odorant and/or activity intensities of involved proteins. (C) Observed oscillations of the G(olf), depolarization, Ca²⁺, and PDE1C, the latter two are the key molecules of the negative feedback circuit of the OT, which finally blocks the depolarization.

Figure 6

Comparison of the simulated depolarization signals between single- and multiple-PDE1C models. Simulation results of single- and multiple-PDE1C models are shown in the left and right panels, respectively. The x-axis indicates the simulation time scale, and the y-axis indicates the intensity of depolarization of the olfactory receptor neutron. The values of "threshold" indicate the firing threshold of PDE1C in terms of their activity levels (see Figure 5C) set in the respective simulations. Depolarization signals were obtained using 50 replications of the respective simulation. (A) Results under the single-PDE1C model in which the firing threshold of PDE1C was set to 0 (threshold = 0). (B) Results under the multiple-PDE1C model (threshold = 0). (C) Results under the single-PDE1C model (threshold = 1). (D) Results under the multiple-PDE1C model (threshold = 1). (E) Results under the model of single-PDE1C (threshold = 1) and limited Ca²⁺ availability. (F) Results under the model of multiple-PDE1C (threshold = 1) and limited Ca²⁺ availability. (G) Results under the model of single-PDE1C (threshold = 1) and positive feedback circuit knocked out. (H) Results under the model of multiple-PDE1C (threshold = 1) and positive feedback circuit knocked out.