Expression and evolutionary divergence of the non-conventional olfactory receptor in four species of fig wasp associated with one species of fig

Background The interactions of fig wasps and their host figs provide a model for investigating co-evolution. Fig wasps have specialized morphological characters and lifestyles thought to be adaptations to living in the fig's syconium. Although these aspects of natural history are well documented, the genetic mechanism(s) underlying these changes remain(s) unknown. Fig wasp olfaction is the key to host-specificity. The Or83b gene class, an unusual member of olfactory receptor family, plays a critical role in enabling the function of conventional olfactory receptors. Four Or83b orthologous genes from one pollinator (PFW) (Ceratosolen solmsi) and three non-pollinator fig wasps (NPFWs) (Apocrypta bakeri, Philotrypesis pilosa and Philotrypesis sp.) associated with one species of fig (Ficus hispida) can be used to better understand the molecular mechanism underlying the fig wasp's adaptation to its host. We made a comparison of spatial tissue-specific expression patterns and substitution rates of one orthologous gene in these fig wasps and sought evidence for selection pressures. Results A newly identified Or83b orthologous gene was named Or2. Expressions of Or2 were restricted to the heads of all wingless male fig wasps, which usually live in the dark cavity of a fig throughout their life cycle. However, expressions were widely detected in the antennae, legs and abdomens of all female fig wasps that fly from one fig to another for oviposition, and secondarily pollination. Weak expression was also observed in the thorax of PFWs. Compared with NPFWs, the Or2 gene in C. solmsi had an elevated rate of substitutions and lower codon usage. Analyses using Tajima's D, Fu and Li's D* and F* tests indicated a non-neutral pattern of nucleotide variation in all fig wasps. Unlike in NPFWs, this non-neutral pattern was also observed for synonymous sites of Or2 within PFWs. Conclusion The sex- and species-specific expression patterns of Or2 genes detected beyond the known primary olfactory tissues indicates the location of cryptic olfactory inputs. The specialized ecological niche of these wasps explains the unique habits and adaptive evolution of Or2 genes. The Or2 gene in C. solmsi is evolving very rapidly. Negative deviation from the neutral model of evolution reflects possible selection pressures acting on Or2 sequences of fig wasp, particularly on PFWs who are more host-specific to figs.

The host-specificity behaviour of insects relies heavily on olfaction [24], which has evolved to a level of extreme sensitivity and specificity [25]. This highly discriminative sensory mode is facilitated by an odour-activated seventransmembrane-domain G protein-coupled receptor (GPCR) that signals cascades [26].  [30][31][32]. Different species of wasps appear to have unique volatile profiles which could account for host-specificity [28]. Such chemical constraints likely reduce host switching. Currently, investigations on adaptive changes in fig wasps are limited to morphological characters associated with an ecological niche [14,[33][34][35][36]. For example, a pollinator must enter the syconium to deposit eggs. Thus, their antennae are easily broken to ease crawling through the ostiole. Unlike PFWs, most NPFWs cannot enter syconia. They usually possess a long ovipositor that is inserted through the syconial wall for egg laying [37][38][39]. Furthermore, both eyes and wings in most adult male PFWs exhibit highly vestigial traits, a correlate of living in the dark fig cavity all their life [40,41]. Unlike these attributes, very little is known about genetic changes in fig wasp's olfactory system.
As members of the GPCR superfamily, olfactory receptors (ORs) for odorous compounds play the critical role in the olfactory process [42]. The process consists of several linked systems ranging from stereo chemical recognition to the generation of an odour code in the brain. A characteristic trait of conventional odour ligand-binding OR types is the tremendous diversity of their sequences, often exhibiting only ~20% identities to each other [43][44][45]. A highly conserved, non-conventional member in the insect OR family is known as Or83b [46]. Orthologs have been identified from Drosophila melanogaster (DOr83b) [47][48][49], Anopheles gambiae (AgOr7) [50], Heliothis virescens (HvirR2) [51], Apis mellifera (AmelR2) [52] and others [50,53,54]. Apparently, Or83b does not directly respond to odorants but rather acts as a chaperone receptor to form heterodimers with other odorant and pheromone receptors, thus ensuring dendritic localization [55][56][57].
While spatial expression patterns of conventional ORs are restricted in small subpopulations of olfactory sensory neurons (OSNs), Or83b is co-expressed with conventional ORs in most, if not all, neurons [52,55]. Tissue-specific expression patterns of putative ORs have been observed in ecologically distinct species. For example, in the hematophagous mosquitoes An. gambiae and Aedes aegypti, Or83b orthologs are expressed in the antennae, legs and proboscis (i.e., general gustatory organs). As such, Or83b orthologs may be involved in locating human hosts for blood feeding [50]. Sex-specific expression patterns occur in both mosquitoes and moths, for blood ingestion and mate searching, respectively [58][59][60]. Insects employ ORs to recognize and discriminate various quantitative or special odours in their ecological niche. Therefore, host specialization could reflect selection acting on ORs [61,62]. For example, bee-specific rapid expansion of the OR family presumably underlies their remarkable olfactory abilities, including perception of several pheromone blends, kin recognition signals, and diverse floral odours [45]. In Drosophila sechellia, rapid evolution and lack-of-function mutations in olfactory and gustatory receptor genes following a host shift reflect positive selection and/or relaxed constraints associated with an altered ecological niche [62].
In view of prior studies in mosquitoes and moths, different tissue-specific expression patterns of Or83b orthologs are expected to occur between resident male and hostsearching female fig wasps, or between PFWs that deposit eggs within the figs, and NPFWs that do not. If the tissuespecific expression patterns occur together with the crucial role of Or83b orthologs in locating a host, we can infer selection pressure on Or83b orthologous genes. The different degree of host specificity between pollinator and nonpollinators [13,[18][19][20][21][22][23] likely makes Or83b orthologous genes subject to different magnitudes of selection. Our experimental group consists of one species of PFW (Ceratosolen solmsi) and three species of NPFWs (Apocrypta bakeri, Philotrypesis pilosa and Philotrypesis sp.), all associated with Ficus hispida. Herein, we identify Or83b orthologous genes in these four species of fig wasps, compare spatial expression patterns and substitution rate of these genes in fig wasps, examine possible evolutionary forces, and explore molecular mechanisms involved in co-evolution.

cDNA cloning
The newly identified genes of C.solmsi, A.bakeri, P. pilosa and P. sp. were named as CsmOr2, AbOr2, PpOr2 and PsOr2, respectively, and the sequences were deposited in GenBank [Accession numbers: EU281848, EU281849, EU281850, EU281851]. The coding regions of these data had the same length in all fig wasp species (1422 bp when the stop codon was not considered). The length was same as in Nasonia vitripennis, but differed from Microplitis mediator by a one-codon indel (3 bp) and from Ap. mellifera by two indels. The amino acid (aa) sequences were extraordinarily conservative among our fig wasps, with > 90.3% identity (the percentage of identical matches between the two sequences over the reported aligned region) and > 96.4% similarity (the percentage of matches between the two sequences over the reported aligned region where the scoring matrix value is greater or equal to 0.0). This suggested that the four species of fig wasp are closely related relative to other species. Sequence identities ranged from 75.3% to 91.4% when compared with other hymenopteran orthologs including N. vitripennis, M. mediator and Ap. mellifera. Alignment of this non-conventional receptor and the other insects showed that the 19 protein sequences shared greater than 60% identity and 72% similarity. Remarkably, extreme conservation was discovered in the final 164 aa of the C-terminal, where the 19 protein sequences shared nearly 90% identity.
We discovered the typical membrane topological structure of insect ORs in the four fig wasp Or2 protein sequences for seven putative transmembrane (TM) domains ( Figure  1). The N-terminal of these receptors, located intracellularly [56,63], occurred in all protein sequences. The putative membrane spanning domains were inferred to occur at similar relative positions in CsmOr2, AbOr2, PpOr2 and PsOr2. Remarkably, the N-terminal between TM 1 and TM 2, and from TM 3 towards TM 4, contained the most variable regions. The variable second intracellular loop (IC2) connecting TM 4 and TM 5 was much longer than the other five loops.

Phylogeny of the Or2 genes
The phylogenetic analyses using 1st and 2nd codon positions showed that four fig wasp plus N. vitripennis Or2 genes formed a well-supported group (Figure 2; 100% bootstrap support). A sister group relationship of Or2 genes was detected for P. pilosa and P. sp. (99% bootstrap support). The sequence of N. vitripennis was more closely related to Or2 genes of P. pilosa and P. sp. than to A. bakeri, although this association was weakly supported (60% bootstrap). The Or2 gene of the pollinator was dissimilar to those of the three NPFWs. The hymenopterans Ap. mellifera and M. mediator clustered with the fig wasps plus N. vitripennis (100% bootstrap support). All dipteran and lepidopteron sequences formed strongly supported groups (97% and 100%, respectively). Removal of the sequence of Ceratitis capitata from the analysis did not affect tree structure. Thus, the sequence was retained in subsequent analyses. Use of aa sequences instead of 1st and 2nd codon positions yielded a nearly identical tree topology except that N. vitripennis first clustered with A. bakeri (47% bootstrap support) and then formed a group with P. pilosa and P. sp. (98% bootstrap support) (tree not show). When we used a Sankoff (step matrix) in PAUP* [64] to force the minimum number of mutations (steps) required to transform from one aa to another, we obtained the same trees found without using the matrix (tree not show). The difference in tree topologies likely resulted from a difference in potentially parsimony informative characters between the aa data (210 sites) and 1st and 2nd nucleotide codon positions (355 sites).

Tissue specificity of expression
Following PCR for tissue specific-expression (Figure 3), all bands were the size expected from the primer design (330 bp). Bands of similar intensities obtained with primers specific to the actin control indicated the integrity of the cDNA preparations. The intensity of the products (relative to an internal control) indicated variable levels of expression in different tissues. Strong expressions of this ORtype occurred in male heads (including antennae and maxillary palps) and in female antennae, in PFW and all types of NPFW. Lower levels of expression were discovered in the abdomens and legs of females in all species of fig wasp. In contrast to female NPFWs, a band of lower intensity was also obtained with cDNA from the thorax of female C. solmsi. No transcripts were detected within nonolfactory tissues (e.g., thorax, abdomen and legs) of all males. In all cases, actin amplifications were more robust for tissue templates, reflecting higher template quantities in the parallel control reactions. This further demonstrated the absence of detectable expression of ORs in non-olfactory tissues. Finally, to confirm this pattern, an additional 10 cycles of PCR were added to the protocols. Even under these extremely sensitive conditions, RT-PCR products were not detected in all male non-olfactory tissues and in female NPFW thoraxes. Genomic contamination of cDNA templates was clearly distinguishable from cDNA products by primers that spanned predicted introns (data not shown). To verify their specificity, the RT-PCR products from each tissue cDNA were sequenced, revealing that a specific product had indeed been obtained in each instance.

Selection pressure on individual amino acid sites
Because adaptive evolution usually occurs on a small region of a gene's sequence and can even be restricted to a single aa site [65][66][67], we sought aa sites that were subjected to selection pressure. Selection pressure can be inferred from the ratio of nonsynonymous to synonymous changes, Ka/Ks (also known as d N /d S or ω). When Ka/Ks = 1, neutral selection is indicated and when Ka/Ks > 1, positive selection is implied. Alternatively, when Ka/Ks < 1, purifying selection is suggested [68]. In the present study, two thresholds for significance (0.1 and 0.2) were taken into account in order to identify residues potentially involved in ligand-binding activities. In our analysis, Ka/ Ks = 0.0844161; no positively selected nucleotide sites were detected at the P < 0.1 in the SLAC analysis. However, 380 of 474 aa sites were negatively selected. Or2-type receptors appeared to be primarily negatively selected. When significance level was adjusted to P < 0.2, one positively selected site was discovered at the 21st aa (P = 0.157371). Serine (S) occurred at this position in the PFWs, but all NPFWs had Glycine (G). The other species of insect displayed various amino acids at this site, including Isoleucine (I), Alanine (A), Tyrosine (Y), Leucine (L), Methionine (M) and G, but not S.

Analyses of substitution rate
Substitution rates were estimated for synonymous and nonsynonymous sites within the 1365 bp Or2 sequences. The same data [Accession numbers: FJ606763] were collected from a closely related species of C. solmsi, Ceratosolen cornutus, which is the pollinator of Ficus auriculata. This Alignment of Or83b orthologous amino acids in insects Figure 1 Alignment of Or83b orthologous amino acids in insects. Conservation of amino acid sequence is displayed as a sequence logo. The relative frequency with which an amino acid appears at a given position is reflected by the height of its oneletter amino acid code, with the total height at a given position proportional to the level of sequence conservation. Transmembrane domains (TM 1-7) and extracellular (EC) and intracellular domains (IC) are numbered and indicated. Threonine (T) and tyrosine (Y) residues that are sites of potential phosphorylation are enclosed in heavy yellow boxes.
addition increased the statistical power of the comparison of substitution rates [69]. Estimations of synonymous substitutions per synonymous site (Ks) between the dipteran species and C. solmsi ranged from 2.1785 to 3.1533 (Table 1). However, the estimates between the dipterans and C. cornutus were much lower, ranging from 1.3811 to 1.5638 only. In contrast, the rather similar estimates of substitutions occurred both between the dipterans and P. pilosa (1.1044-1.4665) and between the dipterans and P. sp. (1.0173-1.4974). Estimations of the number of nonsynonymous substitutions per nonsynonymous site (Ka) were slightly higher between the dipterans and C. solmsi (0.3064-0.3511) than between the dipterans and C. cornutus (0.2995-0.3373). The elevated substitution rate of the Or2 gene in C. solmsi also occurred when the analyses involved other hymenopterans (N. vitripennis and Ap. mellifera) instead of dipterans. Between-species comparisons in pollinators (Ceratosolen) yielded higher values of both Ks and Ka than in non-pollinators (Philotrypesis) ( Table 1).
Tajima's relative rate test [70] detected a significant difference in the rate of synonymous substitutions between C.  solmsi and C. cornutus (Table 2). In contrast, these rates of substitutions were very similar between P. pilosa and P. sp. Significantly elevated substitution rates were not found when C. cornutus was compared with most of the NPFWs (except P. sp.). However, significantly elevated substitution rates were obtained when C. solmsi was compared with all NPFWs. No significant differences in the rate of nonsynonymous substitutions were detected between any lineages, although the number of substitutions at the 1st and 2nd codon positions was higher in C. solmsi than in C. cornutus and NPFWs.
The extent of the non-random substitutions at synonymous codons in different species was measured by three algorithms (Table 3). Congruent with the rate of synonymous substitutions (Ks) being inversely related to codon usage bias [71,72], C. solmsi showed the lowest codon bias. The NPFWs exhibited a much higher codon bias. C. solmsi had the lowest G+C content in their codons. Because the G+C content at second codon positions was similar in all species, the lower G+C content at (synonymous) third coding positions may have been the main cause of the lower codon bias of C. solmsi.

Nucleotide diversity and neutrality tests for populations of fig wasps
Sequences from the N-terminal of the Or2 gene were collected from 20 individuals of each species in order to test the null hypotheses of strictly neutral evolution [73,74]. . Assuming a randomly mating population at equilibrium, the average number of pairwise nucleotide differences between sequences (π) is expected to equal to average number of nucleotides segregating per site (θ). However, π was always smaller than θ (Table 4)  For all sites of N-terminal, the majority of tests were significantly less than zero, suggesting a departure from neutrality [74]. There were two notable exceptions to the departure from neutrality. D* tests in P. pilosa and P. sp. did not indicate significant departure from neutrality, although both the values were close to significance. However, both D T and F* tests for the two species of Philotrypesis were significantly less than zero. Similar to tests of neutrality for all sites, most of tests for nonsynonymous sites significantly rejected the model of strict neutrality, suggesting selective forces. Unlike neutrality test for all sites and nonsynonymous sites, none of the tests of synonymous sites in the NPFWs Or2 data was significant; these data were consistent with a neutral model of evolution. However, significantly negative D T , D* and F* tests were observed for synonymous sites of Or2 data in C. solmsi, indicating too many rare nucleotide polymorphisms with respect to predictions of the neutral theory [76].

Identification of Or2 genes
We identified Or83b orthologs for the first time and from four species of fig wasps associated with F. hispida. As expected, their primary aa sequences with the Or83b subfamily were highly conserved relative to other insects [43,[47][48][49][50][51][52][53][54]. They had > 60% identity and > 72% similar- Ma and Mb are the numbers of substitutions at the 3 rd codon positions in the lineages leading to species A and B, respectively. Na and Nb are the number of substitutions at the 1st and 2nd codon positions in the lineages leading to species A and B, respectively. Species names are as in Table 1. P < 0.05 was used to reject the null hypothesis of equal rates between lineages. * 0.05 > P > 0.01; ** 0.01 > P > 0.001; *** P < 0.001.
ity, respectively, suggesting that the genes were orthologs in the different species. The extremely conservative C-terminal, especially at the last 164 aa, are involved in G-protein binding/activation required for downstream signal transduction [77]. The high identity and similarity of the aa sequences suggested that the four species of fig wasps are closely related species. This high degree of conservatism has not been observed in conventional ORs, even when compared to all 170 candidate ORs of Ap. mellifera [45], to 62 candidate ORs of D. melanogaster [44] and to 79 candidate ORs of An. gambiae [43]. This level of con-servatism reflected a strong selective pressure on the aa sequence and was consistent with the important role of non-conventional receptors in the olfactory process [42]. The aa sequences of the four fig wasps included all the seven transmembrane domains of the G-protein coupled receptor. The group of threonines and tyrosines were extremely conserved (aa 465 and 466) in the region of TM 7 (Figure 1). They constituted candidate phosphorylation sites that may be important for regulating protein function [78]. An unusually long second intracellular loop was detected between TM 4 and TM 5. The function of this structural specialization remained unknown, even though the second extracellular loop of certain types of mammalian G protein-coupled receptors may be critical for ligand binding and affinity [79,80]. Recent bioinformatics and experimental investigations revealed that the membrane topology of ORs in Drosophila was, in fact, the inverse of mammalian GPCRs, with the N-terminal of these receptors located intracellularly [56,63]. Thus, we speculate that the long second intracellular loop of the non-conventional receptor in insects likely plays an important role in binding conventional ligand-binding ORs.

Adaptation of spatial expression patterns of Or2 genes
We characterized the spatial expression patterns of Or2 genes in four species of fig wasp. RT-PCR experiments demonstrated that receptors were expressed only in adult male heads (including antennae and maxillary palps), but more widely in various tissues from adult females, including the antennae, legs and abdomens. The current view is that Or83b gene does not directly bind odorant ligands but rather it acts to form heterodimers with conventional ORs; this ensures appropriate dendritic localization and function [55][56][57]. The broad spatial expression patterns in OSNs further support the essential role of Or83b orthologs for a sense of smell [52,55]. Given the different habits of female and male fig wasps, the sex-tissue-specific expression of Or2-type receptors may indicate a role in host searching and oviposition in adult females.
Host plants are distributed patchily. Therefore, adult females must be able to distinguish the particular odours emitted by their host plant from the myriad of other volatile compounds. Such remarkable sensitivity and specificity is likely achieved by ORs expressed in the legs of adult female fig wasps. Likewise, Or83b orthologs are also expressed in legs of An. gambiae and Ae. Aegypti. Mosquito legs are known only as having a gustatory function [50,53]. Gustatory receptors (GRs) are the only chemosensory receptors whose expression has been detected in the legs of D. melanogaster [81,82], the location of gustatory sensilla. It seems likely that in addition to their olfactory function, Or83b orthologs in legs of fig wasps might function in a contact chemosensory pathway. Because they function in heterodimerization, perhaps Or83b   pulling females out of their galls, and, most importantly, for chewing an exit tunnel for newly transformed female wasps to leave the syconium [85][86][87]. In the male's enclosed environment, the broad expression of ORs may not be required. Thus, it is not surprising that female wasps express ORs in a greater number of tissues than do male wasps. In support of this hypothesis, few sensory hairs occur on a male's body.
We can not exclude the possibility that several conventional ORs independently express in some OSNs of other tissues, beside the head of males. Or83b orthologs do not always co-express with conventional ORs in all OSNs [48][49][50]. Indeed, the lack of AgOR7 (mosquito orthologs) expression in grooved peg sensillum of An. gambiae suggests the presence of an alternative pathway for olfactory signal transduction that is independent of OR7 function [50]. Unexpectedly, a weak band was obtained with cDNA prepared from the thorax of C. solmsi. Most NPFWs oviposit into ovaries of female flowers by inserting their ovipositor through the syconium wall (while not entering the syconium) [37][38][39] [88]. The primary function of legs and the abdomen are movement and reproduction, respectively. The role of olfaction in these organs is a secondary function. The diversity of expression patterns may be important for determining species-specific olfactory profiles, such as in the detection of fruit odours by fruit flies [62] and human host odours by mosquitoes [50,53]. The patterns in fig wasps may reflect species-specific adaptations to ecology, habitat and physiology. Further study into the function and characteristics of Or2 and other conventional ORs will facilitate our understanding of co-evolution in this model system.

Strong purifying selection for most orthologous amino acids
The evaluation of selection pressure used the Ka/Ks ratio, based on 19 orthologous Or2 genes. A very low average Ka/Ks ratio indicated that Or2 genes were mainly subjected to purifying pressure. This finding differed drastically from the properties of the ORs family, most often characterized by rapid evolution and highly species-specific gene repertoires [44,45,62]. A higher level of functional constraint on protein-coding exon sequences should lead to lower level of nonsynonymous variation which usually generates low average nonsynonymous:synonymous substitution ratios. Thus, the low average Ka/Ks ratio suggests that the Or2-type receptor plays an important role in the olfactory process of insects and is subject to a higher level of functional constraint. The 21st aa of Or2 appeared to be subjected to positive selection, although the P-value was not very significant. The aa in PFWs (S) differed from that in all NPFWs (G). Suzuki and Gojobori (1999) demonstrated that two amino acid sites of the human leukocyte antigen (HLA) gene undergoing positive selection might be involved in antigen recognition [89]. Other studies also showed that positive selection is focused mainly on the binding site and the distinct DNA-binding properties are determined by one or a few critical amino acids [90][91][92][93][94]. If the 21st aa performed a binding function, such divergence could have reflected ligand binding specificity between PFWs and NPFWs. Species-specific adaptive divergence might have been driven by the environment of the host [62]. Ecological niches differ between PFWs and NPFWs [37][38][39]. The composition and concentration of bouquets could vary between the inside and outside the syconium, although no data document this to occur. PFWs that enter the syconium are likely to smell extra syconian odours. Perhaps the change of aa at site 21 helped the Or2 of pollinators bind conventional ORs that respond to odours in the syconium.

Elevated rate of substitution in pollinators
All PFWs and most NPFWs are currently classified as being part of the same chalcid family, Agaonidae [95][96][97]. However, recent molecular studies suggest that this family is paraphyletic and all PFWs form a monophyletic group [17,98]. An evaluation of closely related species should help to control for possible differences in mutation rate [69,99]. Thus, sequence data for Or2 genes were also collected from C. cornutus, a pollinator of F. auriculata and closely related to C. solmsi, to increase the power of our analysis of substitution rates.
Between-species comparisons of substitution rates involved two species of PFWs and three NPFWs. An elevated substitution rate was not detected in Philotrypesis, but it was in C. solmsi. The rapid evolution of Or2 genes in C. solmsi became more apparent when the closely related species in the genus Ceratosolen were compared. Low effective population sizes may be accompanied by relatively weak selection compared to drift yet usually drive an increase in the frequency of slightly deleterious substitutions [62,[100][101][102]. However, no evidence suggests that PFWs have smaller population sizes than NPFWs. Typically, pollinators are almost always the dominant wasp (Herre EA & Machado CA, personal communication), though NPFWs might affect pollinator numbers by successfully invading the fig-pollinator mutualism system [103]. Thus, the hypothesis of a low effective population size can not be employed to explain these observations. Machado [104] found that species of Ceratosolen exhibited accelerated mutation in their mitochondrial genome relative to species in other genera. It is possible that this increased rate has affected the nuclear genome as well through cytonuclear associations [105][106][107][108]. For example, Rand et al. (2004) demonstrated that the association between nuclear and mitochondrial substitutions drives the evolutionary divergence [105]. Mishamar et al. (2006) showed that mitochondrial DNA and nuclear DNA complex I genes may have co-evolved [109].
Our results are consistent with the proposal that the rate of synonymous substitutions is inversely related to codon usage bias [71,72]. Codon usage bias was much lower in C. solmsi than in the NPFWs. The changed codon bias indicated that mutations were not strictly neutral [101,110]. Weak selection on synonymous mutations caused a codon usage bias in bacteria, yeast and flies [111][112][113]. In Drosophila, preferred codons also correspond to the more abundant tRNAs [114,115]. Because the most abundant tRNA translates their corresponding codon more rapidly, the preferred codon speeds up translation [116,117]. However, there are many counterexamples of highly expressed genes with little or no codon bias [118,119]. Therefore, we cannot unequivocally state that Or2 genes in PFWs that have lower codon usage bias are expressed at a lower level than in NPFWs. Conversely, the broader tissue expression of Or2 in PFWs than in NPFWs, together with the higher degree of host specificity of PFWs than NPFWS [13,[18][19][20][21][22][23], implies that Or2 genes in PFWs are expected to be expressed at a higher level than in NPFWs.

Non-neutral patterns of nucleotide variation
Olfaction genes of insects are widely assumed to experience selection [43,45,62,[120][121][122]. Our investigation of Or2 polymorphism in four species of fig wasp sought evidence for selection, and the signature of selection can be detected by various tests [74,75,123]. These tests also tend to implicate the corresponding evolutionary forces. The significantly negative D T and F* values for Or2 data from populations of C. solmsi, A. bakeri, P. pilosa and P. sp. indicate a higher-than-expected number of low-frequency variations. This is consistent with purifying selection or directional selection reducing deleterious mutations in Or2 genes [74]. Once deleterious mutations appear in population, they will be maintained at a relatively low frequency due to selection pressures. And low-frequency variation will be increased more than in neutral conditions. For the DNA data, θ will be increased, and the test value will be significantly negative. With D*, Or2 data from C. solmsi and A. bakeri reject the neutrality hypothesis, while those from P. pilosa and P. sp. do not, although the values are close to significance (0.10 > P > 0.05) [75]. This result likely reflects either short length of sequence or a small sample size [124]. Increasing the length of the sequence and/or sample size could yield a more accurate result. Because both the D T and F* tests support a non-neutral model of evolution, most likely the Or2 data in P. pilosa and P. sp. also reflect selection.
Most tests on nonsynonymous sites of Or2 also show a significantly negative departure from neutrality, again suggesting possible purifying selection. Certainly, Or2 plays an important role in fig wasps and variations on nonsynonymous sites are maintained at a relatively low proportion. An alternative explanation for negative test values is the occurrence of directional selection. It eliminates deleterious mutations from the population and pro-motes the fixation of advantageous mutations that optimize function and lead to adaptation to an ecological niche. The D T test, together with D* and F* tests, clearly rejects the null hypothesis that the synonymoussites of Or2 gene in C. solmsi are evolving under a strictly neutral model of molecular evolution [74]. Selection pressures seem to also act against slightly deleterious mutations at synonymous sites within C. solmsi. Regarding synonymous sites in NPFWs, none of tests significantly rejected the null hypothesis. Thus, these data are consistent with a neutral model of evolution [74,75].
The use of several tests appears to be a more powerful means of inferring patterns of selection affecting nucleotide variance. Combined with test results on nonsynonymous sites, Or2 in C. solmsi appears to experience more effective selection relative to that of NPFWs. This might be related to the greater host specificity of PFWs relative to NPFWs, and/or it may also be related to the larger effective population size of PFWs. It remains to be determined whether this trend occurs just for Or2 or throughout the entire genome of C. solmsi. Data from additional loci are required to evaluate this possibility.

Conclusion
Or83b orthologous genes obtained from one pollinator and three non-pollinator species of fig wasps associated with F. hispida were evaluated. We examined spatial expression patterns, evolutionary rate and selective forces acting on these genes in fig wasps.
Or2 genes were expressed beyond the known primary olfactory tissues and potentially this has functional implications [52]. Or2 plays an essential role in the localization and function of co-expressed OR proteins [55,57]. Thus, the expressions of Or2 genes in non-olfactory tissues of female fig wasps strongly indicated the presence of cryptic olfactory inputs in these tissues. Olfactory responses obtained from non-olfactory tissues may have indicated a long-term adaptation to figs. Future investigations can determine whether these genes have a chemosensory function expressed in non-olfactory tissues, or not.
Based on the initial results, four additional degenerate primers were designed and used to obtain the remaining N-terminal coding region sequences: FW5'3: 5'-ATGAT-GAARWYNAAGCAWCARGG-3'; FW3'3: 5'-TTGCTRTA-DATNCCWCGNASRTC-3'; FW5'4: 5'ATGAARWYNAAGCAWCARGGNYTRRTNGCSGA-3'; FW3'4: 5'TGRTCNGCRCTKCCRGCCTTGAA-3'. A reaction volume of 50 μl contained 2 mM MgCl 2 , 0.2 mM of each dNTP, 0.2 μM of each primer, and 2.5 U of EasyTaq DNA polymerase (TransGen Biotech). PCR program involved 5 min at 94°C, then 35 cycles with 94°C for 30 s, 52°C for 40 s and 72°C for 1 min, followed by incubation for 10 min at 72°C. The amplified DNA products were purified and automated DNA sequencing was performed on an ABI3730 with an ABI PRISM BigDye terminator cycle sequencing ready reaction kit (Perkin-Elmer Biosystems). The same N-terminal data were obtained from 20 individuals per species of PFWs and NPFWs associated with F. hispida. These taxa were randomly collected from different trees and different crops. Thus, we were able to assess levels of intraspecific variation and test the neutrality model. Sequences from N-terminal and C-terminal were assembled as contiguous fragments using ContigExpress of Vector NTI Advance 9 (Invitrogen). The four full-length sequences of Or83b orthologous genes were obtained from four species of wasps associated with F. hispida. We also obtained a full-length sequence from C. cornutus. Thus, we were able to compare substitution rates in closely related species. The TMHMM v2. . A significantly heterogeneous pattern was detected at 3rd codon positions in most comparisons (P < 0.05). In contrast to 3rd codon position, the 1st + 2nd nucleotide positions and the amino acid sequences were found to be more homogeneous, except for 1st + 2nd codon position of C. capitata (data not show). Therefore, we used the 1st + 2nd codon positions and the amino acid sequences in the phylogenetic analysis. All positions containing gaps and missing data were eliminated from the dataset (complete deletion option). The final dataset here contained 864 aligned nucleotide positions and 430 aligned amino acid positions. The best tree was selected using the maximum parsimony (MP) criterion. Nodal stability was assessed using bootstrap proportions (1000 replicates) and the values are shown next to the branches. The MP trees were obtained using the Close-Neighbour-Interchange algorithm with search level 3, in which the initial trees were obtained with the random addition of sequences (10 replicates). Phylogenetic analyses were conducted in MEGA 4.0 [132] and PAUP* 4.0b10 [64].

RNA expression
Antennae, thoraxes, abdomens, legs and heads (male) of four species of fig wasp associated with F. hispida were dissected in Sample Protector (TaKaRa), and total RNA of each tissue was isolated by using TRIzol procedures according to the manufacturer's instructions (Transgen Biotech). A series of non-quantitative RT-PCR experiments were performed by using cDNA preparations from various tissues of all male or female fig wasps and the same primer pairs were used for each type of tissue. In order to control for genomic DNA contamination, primers spanning predicted introns were designed for subsequence NEST-PCR reactions: FW5'5: 5'-AGTGCBATCAARTAYTGGGTNGA-3'; FW5'6: 5'-CTNGC-NTACCARGCNAC NAA-3'; FW3'5: 5'-TTACTTCAGCT-GCACCARN ACCATRAA-3'. All RT-PCR reactions were replicated at least three times. We also amplified the actin gene [133] from each tissue as a control for cDNA integrity by use of the following primers: β-actin F: 5'-ATGT-GCAAGGCHGGHTTCGC-3'; β-actin R: 5'-CRTGGATRCCGCA VGAYTCC-3'. PCR products were purified and directly sequenced as described above.

Selective pressure on individual codons
Selective pressure (positive selection and purifying selection) on individual codons (sites) within the coding region of the 19 amino acid sequences of Or83b orthologs were inferred using the Single Likelihood Ancestor Counting (SLAC) package http://www.datamonkey.org [89]. This codon-based maximum likelihood method did not assume equal synonymous substitution rates throughout the sequence and it chose the most appropriate model for nucleotide substitution.

Substitution rate
Evolutionary divergence at the gene-coding region was estimated by the number of synonymous substitutions per synonymous site (Ks) and the number of nonsynonymous substitutions per nonsynonymous site (Ka) [68]. Analyses were conducted using the modified Nei-Gojobori (Jukes-Cantor) method (assumed transition/transversion bias = 1.021) in MEGA4 [132,134,135]. The results were based on the sequences alignment of seven hymenopteran species (C. solmsi, C. cornutus, A. bakeri, P. pilosa, P. sp., N. vitripennisl and Ap. mellifera) and four dipteran species (Ae. aegypti, C. quinquefasciatus, An. gambiae and D. melanogaster). All positions containing alignment gaps and missing data were eliminated only in pairwise sequence comparisons (pairwise deletion option).
Equality of evolutionary rate between two lineages of fig wasps was tested using an outgroup in Tajima's relative rate test in MEGA4 [70,132].P < 0.05 was used to reject the null hypothesis of equal rates between lineages.  [113] methods were implemented in DnaSP program version 4.50.3 [138]. These methods estimated the codon bias present at the gene-coding region. The value of ENC ranges from 20 (only one codon is used for each amino acid; i.e., the codon bias is maximum) to 61 (all synonymous codons for each amino acid are equally used; i.e., no codon bias). CBI values range from 0 (uniform use of synonymous codons) to 1 (maximum codon bias). SChi2 measured the difference between the observed number of codons and those expected from equal usage of codons. A higher SChi2 value indicated a stronger deviation from the random use of synonymous codons.

Neutrality test
The signature of selection can be detected by various tests [74,75,123], but it is not clear which is most powerful. Thus, we exploited Tajima's D statistic [74], and Fu and Li's D* and F* tests [75] to estimate deviations from neutral expectations. These tests were implemented in DnaSP version 4.50.3 [138]. A negative Tajima's D signified an excess of low frequency polymorphisms, indicating population size expansion, purifying selection, recent directional selection and/or background selection of deleterious mutation. A positive Tajima's D signified low levels of both low and high frequency polymorphisms, indicating a decrease in population size and/or balancing selection [74]. The same thing held for Fu and Li's tests.

Authors' contributions
BL carried out the majority of the molecular work, designed and conceived the study and drafted the manuscript. NW designed some primers and carried out a portion of the molecular work. JX participated in the design of the study. YX coordinated the study. RWM contributed to data analysis, interpretations of phylogenetic inference and participated in the writing. DH participated in the design of the study, directed the research and helped to draft the manuscript. All authors read and approved the final manuscript.