Evolutionary Analysis of the Highly Conserved Insect Odorant Coreceptor (Orco) Revealed a Positive Selection Mode, Implying Functional Flexibility.

Odorant coreceptor (Orco) represents one of the essential genes in the insect olfactory system, which facilitates signal transduction and heterodimerization with different odorant receptors (Ors) in the insect antennal dendritic membrane. Evolutionary analysis by detecting positive selection is important to examine the functional flexibility of Orco that potentially supports insect survival. The maximum likelihood codon substitution model was applied using CODEML program as implemented in PAML ver 4.9e package across 59 Orco codon sequences available from GenBank. These sequences represented five major insect orders and two reproductive systems (holometabola and nonholometabola). In the site model that identified common ω values for Orco, it was clearly shown that Orco was under strong purifying selection, indicated by the ω value that was far from 1 (ω: 0.03). However, in to the branch model, positive selection was detected to be acting on Dipteran lineages, whereas in the branch-site model, several sites were under significant positive selection occurring in the following four clades: Coleoptera, Diptera, Lepidoptera, and Psocodea. The typical evolutionary mode acting on Orco was consistent with the entropy value [H(x)], confirming that 48.9% of the Orco site was under conservation (H(x) < 0.5), whereas 26.9% of the Orco sites was under high variation (H(x) > 1). These findings confirmed that Orco genes are generally highly conserved and can possibly be used for the manipulation of insect pest control programs. However, positive selection that acts on certain lineages suggested future adaptive evolutionary ability of Orco to anticipate flexible functions for successful olfactory processes.

One of the important physiological elements for any living organism is the olfactory system. Insects, which represent the most abundant group on earth, extremely depend on the olfactory system, through which they can detect and discriminate a diverse array of chemicals, which result in their behavior of finding food, mates, oviposition sites, and prey (Gross 2006, Hallem and Carlson 2006, Hallem et al. 2006, Ha and Smith 2008). The insect olfactory system primarily deals with signal transduction, which involves the transformation of chemical molecules into an electrical signal, wherein an extracellular signaling molecule activates a cell surface receptor, which in turn alters the intracellular molecules, thus generating a specific response (Ruebenbauer 2006). A critical step in the insect olfactory system occurs in the interaction of odorant molecules with odorant receptors (Ors) that heterodimerizes with an odorant coreceptor (Orco). Orco is hypothetically considered as conserved and plays an essential role during signal transduction (Vosshall 2003, Ha and Smith 2008, Sato et al. 2008, Wicher et al. 2008, Pellegrino and Nakagawa 2009, Leal 2013, Stengl and Funk 2013).

Orco in most insect species has a unique atypical membrane topology; the N-terminal end is located in the cytoplasm, whereas the C-terminal end faces outside of the cell. These experimentally confirmed features imply that Orco proteins differ from the members of the mammalian seven-transmembrane protein family or G protein-coupled receptors (Benton et al. 2006, Tsitoura et al. 2010). Another unique characteristic feature of Orco is the C-terminal conserved region (Pitts et al. 2004, Qi et al. 2010, Yang et al. 2012, Dong et al. 2013). It has been suggested that the conserved amino acid (AA) residues at the C-termini attempted to maintain the primary function of Orco during the 250 million years of evolution since the divergence of Lepidopteran and Dipteran lineages, which mediated the functional interactions of Orco and Ors (Pitts et al. 2004, Jones et al. 2005, Yang et al. 2012, Butterwick et al. 2018). Current real structure Orco analysis by Cryo-electron microscopy technique had deepened the understanding of Orco architecture which is a tetrameric channel, and had confirmed previously reported Orco structure with more detail possible structure–function relation (Butterwick et al. 2018).

Although the characterization of Orco has received increasing attention, there is still lack of a detailed study on the molecular evolutionary process of Orco. To measure the rate of evolution of Orco, a substitution model at the codon level was implemented by estimating the ratio of nonsynonymous-to-synonymous substitutions as ω = dN/dS, with ω >1 indicating positive (adaptive or diversifying) selection, ω = 1 indicating neutral evolution, and ω <1 indicating negative (purifying) selection. The strength and mode of natural selection acting on Orco were measured using this ratio, which also specifically identified the potential site (codon) or branch that may have undergone adaptive evolution (Jeffares et al. 2015). Up to date there were 176 Orcos identified from different insect species (Butterwick et al. 2018). Codon-based substitution model, which includes the site, branch, and branch-site models, was used to indicate the possibilities of positive selection acting on Orco under different sites and lineages. This study is aimed at revealing the Orco evolution pattern at codon level, which would probably provide an insight into their adaptive olfactory capabilities in the environment.

Material and Methods

Sequence Data Set

Sequence data set of Orco was collected manually from the literatures, and from GenBank using keyword searches in the Geneious ver. 11.03 software. However, only 59 sequences identified can be used for the codon-based substitution analysis using CODEML program in PAML package (listed in Table 1), because the program requires both nucleotide (Nt) and amino acid (AA) sequences for each sample. Trimming of nucleotide sequences was performed to retain only the ORF region, which corresponded to their AA residues. Three files were generated to run the CODEML program, including the control file (ctl file) for each model, paml file containing the codon alignment, and tree file, where a preselected branch was coded for the branch and branch-site models. Codon alignment was generated by aligning the selected 59 Orcos AA sequences using Clustal Omega as implemented in Geneious software. Web-based PAL2NAL service was used to convert the aligned AA sequences to their corresponding codon nucleotides. All the gaps and poor alignment regions of the codon were truncated and adjusted by the eye to ensure minimum codon deletion and maintenance of structural motifs, and the final codon alignment was saved as.paml extension. Tree files were constructed using MEGA 6 software by Clustal Omega AA alignment (as implemented in Geneious software) and run by the maximum likelihood method with the best fit model of LG+G according to the lowest BIC resulting from 48 possible models run in MEGA 6 data best fit analysis. The final phylogenetic tree was used for all the codon-based substitution models.

Table 1.

Summary of topological features and the conserved Cam-binding site across insect orders (represented as percentage of total species number for each order level)

Order	N	Topologya		Cam-bindingb
		7TM	in	VRSAIKYWV	ARSAIKYWV
Coleoptera	10	100	70	100	0
Diptera	15	80	80	100	0
Hemiptera	9	77	100	100	0
Hymenoptera	2	100	100	100	0
Lepidoptera	22	100	100	4	96
Psocodea	1	100	100	100	0

For topology, the data are presented as the percentage of total number (N) of sequences; ‘7TM’ indicates the seven-transmembrane protein, while ‘in’ indicates the direction of the N-terminus toward the cytoplasm.

Cam-binding site as predicted by http://calcium.uhnres.utoronto.ca/ctdb/ctdb/sequence.html, resulting in two groups of AA sequences: VRSAIKYWV and ARSAIKYWV.

Topology and Entropy Value Analysis

AA variation in all the Orco genes across different insects was measured by calculating the entropy value [H(x)] at all the positions using BioEdit (ver. 7.1.9) software with the aligned AA data set. The entropy value starts from 0, which indicates no AA variation at a given position; an increase in the entropy value indicates increased variation. The resulting entropy values were further divided into four (arbitrary) categories according to increasing entropy as follows: H(x): 0, 0.1–0.5, 0.5–1, and >1. Conserved regions were assumed based on the entropy values of 0 and <0.5.

To examine the principal characteristic of Orco orthologs across different insect Orco, a web-based transmembrane prediction was performed using MEMSAT2 (http://www.sacs.ucsf.edu/cgi-bin/memsat.py). Particular attention was given to any AA site that had a high entropy value (indicating high variation) and calmodulin (Cam)-binding site (http://calcium.uhnres.utoronto.ca/ctdb/ctdb/sequence.html), which is thought to be the most conserved region across Orco AA sites. TOPO2 software (Web-based service; http://www.sacs.ucsf.edu/cgi-bin/open-topo2.py) was used to visualize the seven-transmembrane structure with Dmel (Drosophila melanogaster (Meigen) (Diptera:Drosophilidae)) Orco used as a model. A representative homology 3D modeling of Dmel Orco protein (AAT71306) was generated with known Orco protein structure from Apocrypta bakeri (Joseph) (Hymenoptera: Pteromalidae) (PDB ID: 6C70) using MODELLER software (Webb and Sali 2014) (Fig. 2B).

Fig. 2.

(A) Seven-transmembrane (7TM) topology of representative Orco (Dmel Orco), as predicted using MEMSAT2. Visualization was performed using the TOPO software. Entropy value of the aligned Orco was mapped with different colors, including white, green, orange, and blue, according to increased AA variations. Calmodulin-binding site is indicated in red, which represents the most conserved region across Orco from different insects. (B) 3D protein structure of representative Dmel Orco, as generated using MODELLER software with known Orco structure from Apocrypta bakeri (PDB ID: 6C70) and visualized using pyMOL. The calmodulin-binding site, N-terminus, and C-terminus are indicated in arrows.

Evolutionary Analysis Using Codon-Based Substitution Model

To examine positive selection acting on Orco, a statistical test for evaluating adaptive evolution among the codons was conducted using CODEML program as implemented in PAML software package (version 4.2b) (Jeffares et al. 2015). In codon analysis, it was assumed that there are mutations that can be classified as synonymous (S-silent) and nonsynonymous (N-AA replacing). By estimating the dN/dS ratio (ω), the strength and mode of natural selection acting on protein-coding genes can be measured (Jeffares et al. 2015). The ratio of dN/dS (ω) can be <1, =1, and >1, which indicate purifying (negative) selection, neutral evolution, and positive selection (Darwinian selection), respectively. The CODEML program has three primary models, including the site, branch, and branch-site models. All these models were executed to detect the possibility of positive selection acting at a particular site, specific lineage, and on few sites along particular lineages (known as foreground branches) (Yang and Nielsen 2002).

Site Model

Under the site model, different models were implemented to obtain the best fit model for detecting positive selection (Yang et al. 2000). The first model was the null model one-ratio (M0 in CODEML), which calculates the average ω value over all the branches and sites of Orco in the phylogeny. Model M0 did not allow ώ variation among the site and branches; therefore, only single ώ value was generated. Although the M0 model is useful to obtain a general view of the selection mode occurring in the entire data set, this model has been considered as unrealistic to detect adaptive evolution. Therefore, two alternative models were generated; one was the nearly neutral model (M1a) that classified Orco codon sites into the following two groups: one codon group that has codons subjected to purifying selection (ω < 1), and other that has codons under neutral evolution (ω = 1). The second alternative model was the positive selection model (M2a), which was similar to the M1a model, with the exception of addition of a third codon classification that considered positive selection (ω > 1) (Jeffares et al. 2015). The hypothesis in the site model was that either classification of the codons into two groups (M1a) improves the fitness data compared with that by the M0 model. Followed by testing of addition of a third class of sites (M2a) with ω >1 (adaptive evolution), the data were fit better than those by the model with two classes [ω < 1 (M1a)]. Model 2a also calculated the Bayes empirical Bayes (BEB) data and posterior mean ω at each site, which was found to be significantly different when higher than either 95% or 99% (Jeffares et al. 2015).

Branch and Branch-Site Models

To investigate the positive selection of Orco occurring in different lineages leading to different insect order phylogeny, including Coleoptera, Diptera, Hemiptera, Hymenoptera, Lepidoptera, and Psocodea, branch models were implemented, which allowed for an additional ω parameter along the specific lineage, leading to different insect order lineages (two ratios). The hypothesis was that the two-ratio model was significantly different compared with the one-ratio model (null, M0 one-ratio model), where only one ω value was estimated across the phylogeny.

The branch-site model was used for detecting positive selection occurring in individual codons along specific lineages. Particularly, the branch-site model examined the likelihood ratio test (LRT), which compares the tested model (allows positive selection on each of the foreground lineages) with the null model that does not allow such positive selection (Zhang et al. 2005). The branch-site model provided BEB analysis to identify specific positively selected sites in a particular lineage (foreground), with a posterior probability (P) of either 0.95 or 0.99.

Both branch and branch-site model analyses were conducted to examine the six primary lineages (foreground) indicated on the optimized maximum likelihood phylogenetic tree (Fig. 1) by incorporating the notation of #1 adjacent to the internal node of each primary lineage in the newick format tree file assisted by TreeView software (Jeffares et al. 2015).

Fig. 1.

Average AA frequencies for all Orco sequences, representing five major insect orders: Coleoptera, Lepidoptera, Hemiptera, Hymenoptera, and Diptera. The AAs are grouped according to their biochemical properties, with the following percentages for each: 42, 27, 13, 11, and 7.4% for nonpolar, polar, aromatic, positively charged, and negatively charged AAs, respectively.

Likelihood ratio test

The CODEML program resulted in parameters that were best fit to the data and a likelihood value for each model. The natural log of likelihood value (lnL) was then used to determine whether positive selection had occurred in a particular site or lineage by determining whether the tested model (where ω >1 at some sites or lineages) fits the data better than that which does not include positive selection (null model; no sites with ω >1). This was achieved by calculating LRTs (2Δl). The null model can be rejected if LRT (2Δl) is greater than the chi-square critical value according to the k degrees of freedom, with a significance level of 0.01 or 0.001 (Yoder et al. 2014, Jeffares et al. 2015).

Results

The typical structure of Orco was examined for the seven-transmembrane structure prior to evolutionary analysis. The seven-transmembrane structure is a typical structure of any odorant receptor, but the orientation of the N-terminus toward cytoplasm is specific to insects. Across the 59 Orco sequences corresponding to 59 insect species that were grouped into Coleoptera, Diptera, Hemiptera, Hymenoptera, Lepidoptera, and Psocodea, only Diptera and Hemiptera were found to have variations in the seven-transmembrane structure. However, the N-terminus orientation also varied particularly in Coleoptera and Diptera. Interestingly, the Cam-binding site was detected in all the species, except those of the Lepidoptera order, in which most species had one mutation from V to A (Table 1).

Regarding the polarity of AAs, nonpolar AAs were found to dominate the structure of Orco protein, followed by polar uncharged AAs. The AAs possessing aromatic, negative, and positive charges were found to exhibit a similar distribution across Orco (Fig. 1).

The entropy value variation was visualized by the seven-transmembrane topology of Orco using Dmel Orco as a representative gene (Fig. 2A). 3D protein homology modeling (DmelOrco) was generated as in Fig. 2B. As shown by the visual topology, most conserved region existed in the region toward the C-terminus (white or green color); this is particularly interesting in terms of molecular behavior when the protein is exposed to the extracellular region. The total 501 aligned AA sequences exhibited variation in entropy values, including 23.95, 24.95, 24.15, and 26.95 for H(x) = 0, 0.1–0.5, 0.5–1, and >1, respectively.

Phylogenetic tree of Orco was constructed using the maximum likelihood statistics with the LG+G best fit model (Fig. 3). Of the 59 Orco sequences, the tree diverged into five primary clades in parallel with the common insect classification (morphological based) at the order level, including Coleoptera, Diptera, Hemiptera, Hymenoptera, Lepidoptera, and Psocodea; the tree diverged similarly when all the available Orco sequences were included (146 sequences).

Fig. 3.

Phylogenetic tree of 59 Orco genes constructed using the maximum likelihood method with the TG+G best fit model. The tree was rooted to Hemipteran lineages. Each primary lineage is indicated in different color, including the internal nodes (in circles), where the specific lineages were analyzed using the branch and branch-site models.

Site Model

Positive selection was first examined using the site model (Yang et al. 2000). Three models were implemented, including M1a (nearly neutral), M2a (positive selection), and null (M0;one-ratio). LRT analysis confirmed that the M1a model was a significantly better fit than the M0 model, whereas the M2a model was not significantly different compared with the M1a model (Table 2).

Table 2.

Positive selection analysis using the site model on 59 insect Orco sequences

Model	np	InL	K	Parameter estimates	2Δl
Site model
M0: One-ratio (null)	118	−32518.44	1.73	ω: 0.03315
M1a: Nearly neutral	119	−32448.059	1.79	p0 = 0.98; p1 = 0.02	140.76a
				ω0 = 0.03; ω1 = 1
M2a: Positive selection	121	−32448.059	1.79	p0 = 0.98; p1 = 0; p2 = 0.02;	0
				ω0 = 0.03; ω1 = 1; ω2 = 1;
M7	119	−31635.722	1.69	p = 0.46019; q = 8.38195
M8	121	−31635.726	1.69	p0 = 0.99; p = 0.46; q = 8.38	0.0083
				(p1 = 0.00001); ω = 1.05

Significantly different according to the χ2 value with specific k degrees of freedom at a significance level of 0.001.

Branch and Branch-Site Models

Both the branch and branch-site models used each primary insect order lineage as the foreground. The foreground ω varied depending on the selected lineage; Coleoptera and Hemiptera showed an increasing ω value compared with their background value, including the slight increases for Psocodea and Lepidoptera, but none of them were significantly different in LRT analysis. Hymenoptera and Diptera showed decreasing ω values, and significant LRT results were observed for the Dipteran lineage. In this branch model analysis, no positive selection was detected (Table 3).

Table 3.

Positive selection analysis using branch and branch-site models on 59 insect Orco sequences

Model	np	InL	K	Parameter estimates	2Δla	Positively selected sites (BEB)b
Branch model
M0: One-ratio (null)	118	−32518.4	1.728	ω: 0.033		Not allowed
Coleoptera	119	−32517.9	1.723	ω0: 0.033, ω1:10.503	1.16	Not allowed
Diptera	119	−32512.6	1.735	ω0: 0.033, ω1:0.0008	11.6**	Not allowed
Hemiptera	119	−32518.4	1.728	ω0: 0.033, ω1: 15.292	0.01	Not allowed
Hymenoptera	119	−32518.3	1.731	ω0: 0.033, ω1:0.016	0.34	Not allowed
Lepidoptera	119	−32518.4	1.728	ω0: 0.033, ω1: 0.043	0.02	Not allowed
Psocodea	119	−32518.4	1.728	ω0: 0.0332, ω1:0.037	0.01	Not allowed
Branch-site model
Coleoptera	121	−32427.4	1.794	p0 = 0.91; p1 = 0.02; p2a = 0.07; p2b = 0	9.59*	16V, 20K, 99S, 322V
				ω0 = 0.03; ω1 = 1; ω2a = 32.04; ω2b = 32.04
Diptera	121	−32429.4	1.796	p0 = 0.92; p1 = 0.02; p2a = 0.06; p2b = 0	14.5**	252I*, 269S
				ω0 = 0.03; ω1 = 1; ω2a = 999; ω2b = 999
Hemiptera	121	−32423.9	1.791	p0 = 0.89; p1 = 0.02; p2a = 0.09; p2b = 0	17.73**
				ω0 = 0.03; ω1 = 1; ω2a = 999; ω2b = 999
Hymenoptera	121	−32432.1	1.796	p0 = 0.91; p1 = 0.02; p2a = 0.06; p2b = 0	4.61
				ω0 = 0.03; ω1 = 1; ω2a = 55.14; ω2b = 55.14
Lepidoptera	121	−32427.6	1.785	p0 = 0.89; p1 = 0.02; p2a = 0.09; p2b = 0	2.62	117S
				ω0 = 0.03; ω1 = 1; ω2a = 6.71; ω2b = 6.71
Psocodea	121	−32457.6	1.779	p0 = 0; p1 = 0; p2a = 0.98; p2b = 0.02; ω0 = 0.03; ω1 = 1; ω2a = 1; ω2b = 1	42.04**	18F*,33Q*,36F*,50Q,77S,81T*,87N,97S,109M*,1 30I*,137V*, 67G*, 170L*,171G*,209Q, 229V*, 238A*, 249V*, 251E*,325F, 328I

Significantly different values according to the χ2 value with specific k degrees of freedom at a significance level of 0.001 (**) or 0.01 (*).

Significantly different values with P > 95% (*) and P > 99% (**).

To examine which site of a particular lineage may be involved in positive selection, the branch-site model was implemented for each lineage. Comparing the null model and tested model for the same tree file clearly demonstrated that most tested models (through LRT analysis) showed better fit than that of the null models, except the Hymenopteran and Lepidopteran lineages. Positive selection was detected for tested models that showed significant LRT results, particularly for Coleoptera, Diptera, and Psocodea.

Discussion

In molecular evolution, under no selective pressure, sequences most often evolve neutrally, and most genes are under strong purifying selection, due to which the occurrence of positive selection is rare. Orco has been considered as one of the most conserved genes in insects, whose primary function is to receive chemical molecules from the environment and translate them into electrical signals, known as signal transduction, for response processing by the nervous system (Vosshall 2003, Ha and Smith 2008, Sato et al. 2008, Jeffares et al. 2015). The Orco conservation had been reported previously, and in parallel with our finding through entropy value [H(x)] method, where C terminus had more conserved proportion. Recent report on Orco architecture and testing with VUAA1 agonist had confirmed that conservation is related to the possible adaptation to diverse chemical stimulus (Butterwick et al. 2018).

Although Orco has been hypothesized to have undergone a 250-million-year conservation (Jones et al. 2005), there is limited detailed information to support this hypothesis. A total of 176 Orco genes have been characterized recently across different insect orders (Butterwick et al. 2018); however, most of them have been intended for potential utilization in the field of pest management, as indicated by recent reports (Soffan et al. 2016, Zhang et al. 2016, Liu et al. 2017). Little attention has been given to determine the pattern of evolutionary of Orco. Structural topology analysis of insect Orco proteins revealed that in most cases, the N-termini of Orco proteins were oriented toward the cytoplasmic region; this feature is particularly specific for insect Orco proteins (Table 1; Fig. 2A). The detail of Orco architecture even revealed that indeed Orco had a tetrameric structure with specific pore structure, potentially involved in the binding with Ors (Butterwick et al. 2018). In term of polarity, Orco had a common feature of any transmembrane protein, where majority of Orco portion are nonpolar (Fig. 1) (Engelman et al. 1986). The presence of multiple charged amino acid polarity across the structure has been suggested to stabilize the transmembrane domain in each sub unit (Butterwick et al. 2018). The entropy value analysis clearly showed that most sites were generally conserved, reaching 48.9% of the sequence (H(x) < 0.5), which was bias to C terminus (Fig. 2A). However, minor but high variations indeed occurred as shown in the visual topology of Orco. These findings are in accordance with the general knowledge regarding Orco conservation in insects (Jones et al. 2005). This typical conservation pattern had been hypothesized related to the Orco flexibility to bind with diverse array of Ors, particularly when it was related to less conserved region in N terminus (Butterwick et al. 2018).

The phylogenetic tree of Orco constructed using the maximum likelihood method revealed the consistent organization of insects based on the classical insect tree according to order level as well as reproductive system. All holometabolous insects (ants, bees, beetles, moths, and flies) were located in the same clade, separated from hemimetabolous insects (aphids, bugs, lice, and locusts). Based on the similarity between the Orco gene tree and morphological taxonomy, it has been speculated that the olfactory system is indeed a product of evolution, and that Orco may have occurred even before Ors (Missbach et al. 2014).

Orco has been suggested as a conserved gene, which was clearly verified using the site model that resulted in an ω-value of 0.03. However, some positive selection was detected through the branch-site model in parallel with the significant better fit LRT model compared with the null model. The lineages (insect order) which undergo Orco positive selection include Coleoptera, Diptera, Lepidoptera, and Psocodea. Similar study was conducted to examine the positive selection for each lineage in Vomeronasal receptor genes (Yoder et al. 2014).

Positive selection reflected by the situation of allele, which indicated the potential to promote AA substitutions, caused parallel AA changes, increased the rate of insertion or deletion substitutions, accelerated gene loss, and enhanced gene expression noise (Zhang 2010). Although the functional effect of this positive selection needs to be further studied, it indicated that a gene in a particular lineage has the potential to promote the beneficial allele (Zhang 2008). Since current understanding had confirmed that possible mutation in Orco indeed could be rescued by Ors through heterodimer system (Butterwick et al. 2018). Therefore, the positive selection presence possibly had relation with potential flexibility mechanism of Orco to complement and compensate for Ors diversity, which in fact might enable Orco–Ors heterodimer to adapt with diverse chemical stimulus in environment.

Click here for additional data file.