Introduction {#s1}
============

USP/RXR, a group II nuclear receptor and transcription factor which belongs to the steroid receptor superfamily, is found across a diversity of Metazoan species ranging from sponges to mammals [@pone.0023416-Wiens1]. First characterized in *Drosophila melanogaster*, USP is the insect homolog of vertebrate RXR [@pone.0023416-Oro1], [@pone.0023416-Mangelsdorf1]. In insects USP/RXR is involved in an array of functions including metamorphosis, reproduction, growth and development (see [@pone.0023416-Hall1], [@pone.0023416-Barchuk1], [@pone.0023416-Martn1], [@pone.0023416-Tan1], [@pone.0023416-Xu1] for examples). USP/RXR is known primarily as the partner of the ecdysone receptor (EcR). In response to the binding of 20-hydroxyecdysone (20E) to EcR, heterodimerization between USP/RXR and EcR occurs and the complex then binds DNA to transactivate genes [@pone.0023416-Thomas1], [@pone.0023416-Yao1], [@pone.0023416-Yao2]. However, it is the controversy over whether or not USP/RXR binds juvenile hormone (JH) in insects, and is activated independent of EcR, which has led to the most debate regarding the function and evolution this receptor [@pone.0023416-Jones1], [@pone.0023416-Xu2], [@pone.0023416-Hayward1].

In vertebrates, 9-*cis* retinoic acid (9cRA) was identified as the high affinity ligand of the RXR receptor [@pone.0023416-Levin1], [@pone.0023416-Heyman1], [@pone.0023416-Mangelsdorf2]. However, in insects no natural ligand has been conclusively identified, and USP/RXR remains an orphan receptor. JH has been suggested as a candidate ligand [@pone.0023416-Jones1], but only experimental evidence from the dipteran *D. melanogaster* supports this hypothesis. In cell lines expressing *Drosophila* USP, the application of JH III induces the transcription of a transfected promoter, suggesting that JH binds USP resulting in a functional outcome [@pone.0023416-Xu2], [@pone.0023416-Wozniak1], [@pone.0023416-Fang1]. Fluorescence-binding assays have shown that *Drosophila* USP binds not only JH III but also the JH precursors farnesol, farnesoic acid, and methyl farnesoate [@pone.0023416-Jones2]. However, JH does not appear to directly bind with USP/RXR in less derived insects such as the holometabolous *Tribolium castaneum* or the hemimetabolous *Locusta migratoria* [@pone.0023416-Hayward1], [@pone.0023416-Iwema1]. There is also little evidence to suggest a high affinity for retinoids in insects [@pone.0023416-Oro1], [@pone.0023416-Iwema1]. However, recent displacement binding experiments have shown that USP/RXR may bind 9cRA and all-*trans* RA at the high nanomolar range in *L. migratoria* [@pone.0023416-Nowickyj1]. Given the small quantity of RA found in *Locusta*, the physiological significance of this *in vitro* finding is unclear.

Lineage specific variation is also evident in the structure of USP/RXR. A comparative analysis of structural data demonstrates key differences in the ligand-binding pocket (LBP) of USP/RXR in insects. Crystallography data have shown that Dipteran and Lepidopteran USP share a conserved ligand-binding domain (LBD) with a large hydrophobic cavity capable of accepting a natural ligand [@pone.0023416-Clayton1], [@pone.0023416-Billas1], [@pone.0023416-Billas2]. In both the Diptera and Lepidoptera, USP copurified from the bacterial expression system with a phospholipid occupying the LBP. However, the identity of the endogenous ligand is unknown. In contrast, crystal structures of USP/RXR from two less derived insects, the red flour beetle *T. castaneum* and the sweet potato whitefly *Bemisia tabaci*, reveal a collapsed LBP [@pone.0023416-Iwema1], [@pone.0023416-Carmichael1]. Given the taxonomic relationships of insects from which USP/RXR has been crystallized, this data implies a derived open pocket in higher insects. Despite the open LBP, Dipteran and Lepidopteran USP may not function independently as ligand-dependent transcription factors. In these species the ligand-dependent activation function (AF-2) domain in α-helix H12, is locked in an antagonistic conformation [@pone.0023416-Clayton1], [@pone.0023416-Billas1], [@pone.0023416-Billas2]. As a consequence of interactions between H3, H11--H12 and a highly conserved insertion in the loop connecting helices H1 and H3 (L1--3), residues in H12 occupy the coactivator groove (H3, loop L3--4 and H4), precluding the agonist conformation and preventing the binding of transcriptional coactivators [@pone.0023416-Clayton1], [@pone.0023416-Billas1], [@pone.0023416-Bonneton1].

These striking lineage-specific differences in both structure and function have sparked an interest in exploring the molecular evolution of insect USP/RXR, in part to understand the evolution of endocrine regulation, but also to aid in the design of selective pesticides which target the USP/RXR-EcR complex. The use of relative-rate tests revealed high divergence rates in the Lepidoptera and Diptera for both USP/RXR and EcR, particularly in the LBD of USP/RXR [@pone.0023416-Iwema1], [@pone.0023416-Bonneton1]. These results suggested that the heterodimer may have coevolved to become functionally divergent in the Holometabola. It was later shown that this event was unique to the Mecopterida, a suborder of insects which includes the Diptera, Trichoptera, Mecoptera, Siphonaptera and Lepidoptera [@pone.0023416-Bonneton2], [@pone.0023416-Kristensen1]. More recently, this work has been extended to the heterodimer interface between USP/RXR and EcR. By reconstructing and modeling the ancestral Mecopterida heterodimer, Iwema *et al.* [@pone.0023416-Iwema2] demonstrated that an expanded dimerization surface was common to this branch. Furthermore, the data suggested that this enlarged surface was the result of torsion in the structure of USP/RXR caused by the position of loop L1--3 in these insects.

Based on experimental data from *Drosophila*, crystallography work from *Heliothis virescens*, and molecular evolutionary studies, a model for the functional evolution of USP/RXR has been proposed. Iwema *et al.* [@pone.0023416-Iwema1], [@pone.0023416-Iwema2] and Tocchini-Valentini *et al.* [@pone.0023416-TocchiniValentini1] suggest the following three groups of USP/RXR: 1. The RXR type which evolved retinoid-binding early, a function retained in the Cnidaria, Mollusca, and Chordata; 2. The Non-Mecopterida type which lost ligand-binding function along an ancient arthropod lineage, a state found in insect orders such as the Hymenoptera, Coleoptera, Orthoptera and Dictyoptera; 3. Finally, the Mecopterida USP type which may have gained a novel ligand-binding function and an enlarged dimerization interface with EcR during the evolution of some higher insects.

To investigate these hypotheses of USP/RXR evolution, we implemented codon-based maximum likelihood methods to estimate the ratio of non-synonymous (*d* ~N~) to synonymous (*d* ~S~) substitutions across an insect phylogeny and independently along the Mecopterida lineage. This ratio serves as an indicator of selective constraint and such methods have been successfully used to detect positive selection, the signature of functional gain, in a variety of insect gene families and vertebrate nuclear receptors [@pone.0023416-Kelleher1], [@pone.0023416-Wu1], [@pone.0023416-Smadja1], [@pone.0023416-Bulmer1], [@pone.0023416-Chen1]. We also examined among-site variation in *d* ~N~/*d* ~S~ across a dataset of Mecopterida and Non-Mecopterida insects using random-sites models to investigate changes in site-specific evolutionary rates across key structural and functional motifs. Additionally, several of these analyses were repeated with a dataset of EcR sequences to examine evolutionary constraint in the heterodimer. Finally, we employed newly developed methods which allow the independent estimation of *d* ~N~ and *d* ~S~ in order to determine the effect among-site variation in *d* ~S~ on the patterns of site-specific evolutionary rates observed under random-sites analyses.

Methods {#s2}
=======

Sequence data collection and dataset assembly {#s2a}
---------------------------------------------

USP/RXR and EcR sequences were collected from literature and GenBank using a combination of BLAST and keyword searches ([Table S1](#pone.0023416.s005){ref-type="supplementary-material"}). Where possible, EcR-A isoform data were used in order to incorporate a larger portion of the protein in subsequent analyses. Both USP/RXR and EcR sequences were not available for all species, and in such cases sequences from the most closely related insects were used instead. Amino acid multiple sequence alignments for USP/RXR and EcR were constructed using ClustalW [@pone.0023416-Thompson1] as implemented in MEGA 4 [@pone.0023416-Tamura1] and adjusted by eye to ensure structural motifs were maintained. Poorly aligned regions and major gaps were deleted (see supporting [Figure S1](#pone.0023416.s001){ref-type="supplementary-material"}, [S2](#pone.0023416.s002){ref-type="supplementary-material"}).

To explore selective constraint across the entire phylogeny, both USP/RXR and EcR alignments were truncated such that only the well conserved LBD was included for use in branch and branch-sites analyses (termed USP/RXR LBD and EcR LBD datasets). To examine the variation in evolutionary rates across sites and across functional domains within each group, the USP/RXR and EcR alignments were each split into Mecopterida and Non-Mecopterida datasets (termed USP/RXR A/B-LBD and EcR A/B-LBD datasets for Mecopterida and Non-Mecopterida respectively).

Estimation of evolutionary rates {#s2b}
--------------------------------

In order to test for positive selection and examine changes in site-specific evolutionary rates codon-based likelihood methods were used to estimate *d* ~N~/*d* ~S~ (ω) ratios across the USP/RXR gene. Under no selective pressure, sequences evolve neutrally, and this is indicated by ω = 1, whereas ω\<1 indicates purifying selection, and ω\>1 positive selection [@pone.0023416-Kimura1], [@pone.0023416-Yang1]. To estimate ω several branch-specific [@pone.0023416-Yang2], [@pone.0023416-Yang3], branch-site [@pone.0023416-Yang4], [@pone.0023416-Zhang1], and random-sites [@pone.0023416-Nielsen1], [@pone.0023416-Yang5] models were implemented using the codeml program of the PAML software package (version 4.2b; [@pone.0023416-Yang6]). Random-sites models which allow for variation in *d* ~S~ [@pone.0023416-Pond1] were also implemented using the HyPhy software package (version 1.0; [@pone.0023416-Pond2]). In addition, many of these analyses were also carried out for EcR to examine evolutionary constraint in the USP/RXR-EcR heterodimer.

The amino acid alignments described above for each of the six datasets were converted into nucleotide data. A tree reflecting current understanding among major insect lineages was used for both USP/RXR and EcR ([Figure 1](#pone-0023416-g001){ref-type="fig"}; [@pone.0023416-Whiting1], [@pone.0023416-Wheeler1] Insecta; [@pone.0023416-Hunt1] Coleoptera; [@pone.0023416-Weller1] Lepidoptera; [@pone.0023416-Yeates1] Diptera; [@pone.0023416-Brothers1], [@pone.0023416-Castro1] Hymenoptera). Gene trees generated with our USP/RXR and EcR alignments, using both maximum-likelihood [@pone.0023416-Guindon1] and neighbor-joining [@pone.0023416-Tamura1], [@pone.0023416-Saitou1] methods, were generally consistent with these known inter-species relationships (see supporting [Figure S3](#pone.0023416.s003){ref-type="supplementary-material"}). However, low branch support for many of the ordinal relationships resulted in trees with poor resolution of the Non-Mecopterida, therefore all analyses were only performed using the species trees. For analyses carried out on the USP/RXR and EcR A/B-LBD datasets, the trees were modified to include only Mecopterida or Non-Mecopterida taxa. However, this resulted in a tree bifurcated at the root with two clades when only the Mecopterida were considered, so *T. castaneum* was added as an outgroup to form a tripartition tree. Correspondingly, the *T. castaneum* sequence was added to both USP/RXR and EcR A/B-LBD Mecopterida datasets.

![Phylogeny of insect species used in our analysis.\
The topology of the trees used for the USP/RXR (left) and EcR (right) datasets are based on known taxonomic relationships. Species contained in the Diptera (Dip.) are shaded green, Lepidoptera (Lep.) blue, Trichoptera (Tri.) yellow, Hymenoptera (Hym.) orange, Coleoptera (Col.) red, Hemiptera (Hem.) purple, Dictyoptera (Dic.) and Orthoptera (Ort.) grey. The origin of holometabolous insects is indicated by an encircled bold H. The highlighted branch is the foreground lineage for detecting positive selection, where a separate *d* ~N~/*d* ~S~ ratio was estimated for the Mecopterida (ω~m~).](pone.0023416.g001){#pone-0023416-g001}

It is important to note that our dataset was composed of insect taxa which diverged over 140 million years ago [@pone.0023416-Gao1], and that the saturation of synonymous substitution rates over such large distances may result in the estimation of *d* ~N~/*d* ~S~ ratios higher than the actual value [@pone.0023416-Kimura2]. However, this would affect all estimations across all datasets, and our primary objective was to examine relative differences and patterns in substitution rates.

Branch and branch-site models {#s2c}
-----------------------------

To investigate the gain in function hypothesized to have occurred in Mecopterida USP/RXR, branch models were implemented which allowed for an additional ω parameter along the lineage leading to the Mecopterida (Two-ratios), or alternatively along the lineage leading to the ancestor of Mecopterida and Hymenoptera. Codon frequencies were estimated using the F3x4 matrix and models were run from varying starting ω values ranging above and below 1, or κ (relative transition to transversion rate ratio) values ranging from 0 to 5 in order to test for convergence. Likelihood ratio tests (LRTs; [@pone.0023416-Felsenstein1]) were then conducted to determine statistically significant differences between nested models. Models with the additional ω parameter were compared against M0 where only one ω value was estimated across the phylogeny. Since coevolution of USP/RXR and EcR has been suggested, the above branch models were also applied to the EcR LBD dataset.

Given the elevation in ω identified by branch models, variation in selective pressure among codon sites in Mecopterida USP/RXR was examined using branch-site models applied to the USP/RXR LBD dataset. Branch-site models allow the use of a Bayes Empirical Bayes (BEB) analysis [@pone.0023416-Yang7] to identify specific positively selected sites within the gene along a given branch. A model with Mecopterida designated as foreground, Model A alt, was compared against a stringent model for positive selection, Model A null, and a less stringent model, M1a. Positive selection is detected if 'Model A alt' is a better fit than 'Model A null' where classes of positively selected sites have ω set to 1, whereas relaxed purifying selection (possibly suggestive of positive selection) is indicated if 'Model A alt' is a better fit than M1a which does not allow a class of sites with ω\>1 (referred to as test 2 and test 1 of positive selection by [@pone.0023416-Zhang1]). Sites identified as under positive selection by BEB analysis with a posterior probability (P)\>0.95 were subsequently mapped onto the 1.65 Å crystal structure of *H. virescens* USP (PBD 1G2N) or the 2.90 Å crystal structure of the *H. virescens* USP-EcR heterodimer (PBD 1R1K) using the VMD software package [@pone.0023416-Humphrey1].

Random-sites models {#s2d}
-------------------

In order to compare among-site variation in *d* ~N~/*d* ~S~ across functional domains of USP/RXR and EcR between Mecopterida and Non-Mecopterida insects, random-sites models were implemented using PAML [@pone.0023416-Nielsen1], [@pone.0023416-Yang5]. Codon frequencies and initial starting values were the same as those used in the above branch and branch-site models. Comparison between M3, with three discrete classes of ω, and null M0 tests for variation in selection pressure among sites, but does not explicitly incorporate positive selection. However, comparisons between both M2a and null M1a (discrete distribution), and M8 and null M7 (where ω is drawn from a beta distribution), test for positive selection by allowing for an additional class of sites with ω\>1. Posterior probabilities for site classes under sites models can be calculated either by the Naïve Empirical Bayes (NEB) approach [@pone.0023416-Nielsen1], [@pone.0023416-Yang5] or by the BEB approach [@pone.0023416-Yang7] which considers sampling errors in the estimation of ω. For models such as M8 where BEB data was available, the posterior mean ω at each site was plotted using site class assignments as calculated by BEB, not NEB.

Random-sites models that incorporate variation in *d* ~S~ {#s2e}
---------------------------------------------------------

To account for the effect of variation of synonymous substitution rates across sites, we implemented several codon selection analyses using the dNdSRateAnalysis program of the HyPhy software package [@pone.0023416-Pond2]. Similar to M3 in PAML the variable nonsynonymous rates model, or 'Nonsynonymous', implemented in HyPhy assumes a constant *d* ~S~ (α~S~ = 1), but samples *d* ~N~ (β~S~) values from a given rate distribution. However, the dual variable rates model, or 'Dual', estimates *d* ~S~ (α~S~) and *d* ~N~ (β~S~) independently, sampling both from a given rate distribution [@pone.0023416-Pond1]. In this study, both models were run using the MG94×REV core rate matrix with the GDD (general discrete distribution) rate distribution method using the 'Independent Discrete' setting. For each model three synonymous and nonsynonymous rate classes (3×3) were specified. All models were run several times using the randomized initial value option to find the global optimum. To test for site-to-site variation in *d* ~S~ across both USP/RXR and EcR A/B-LBD datasets for the Mecopterida and Non-Mecopterida, nested Dual and Nonsynonymous models were compared using LRTs.

Results {#s3}
=======

Branch models {#s3a}
-------------

The current theory of USP/RXR functional evolution proposes that the ligand-binding function was lost in an ancient arthropod lineage followed by a subsequent gain in function along the Mecopterida lineage. Additionally, it has been suggested that Mecopterida USP/RXR acquired an expanded dimerization interface with EcR. We tested these hypotheses of functional gain by estimating evolutionary rates across a dataset composed of insect USP/RXR ligand-binding domain sequences (USP/RXR LBD dataset), using codon-based models of substitution ([Figure 1](#pone-0023416-g001){ref-type="fig"}). Likelihood scores and ω (*d* ~N~/*d* ~S~) values, as calculated by PAML, are shown in [Table 1](#pone-0023416-t001){ref-type="table"}. Branch models implemented in PAML allow ω to be freely estimated along specified foreground branches while all other background branches are constrained to the same ω across the phylogeny (Two-ratios model). A branch model analysis, for which the Mecopterida lineage was set as foreground, demonstrated an elevated value (ω~m~ = 0.166) compared to the background (ω = 0.040), but cannot distinguish between relaxed constraint or positive selection at a subset of sites ([Table 1](#pone-0023416-t001){ref-type="table"}). However, the likelihood ratio test (LRT) indicated that the additional parameter did not yield a significantly better fit for the data than the null model M0 where one ω is estimated across the entire phylogeny (p = 0.347).

10.1371/journal.pone.0023416.t001

###### Parameter estimates for branch and branch-site models (LBD).

![](pone.0023416.t001){#pone-0023416-t001-1}

                                                                                                                                                               LRT                          
  ---------------------------------------------- ---- ----------- ------- ------------------------------------ -------------------------- --------------------------------------------- --- -------------------------------------------------
  **USP/RXR**                                                                                                                                                                               
  M0: one-ratio                                   53   −11110.02   1.537               ω = 0.041                           \-                                                               
  Branch-specific models:                                                                                                                                                                   
  Two-ratios                                      54   −11109.58   1.533       ω~0~ = 0.040, ω~1~ = 0.166                  \-                                  M0                        1                        0.347
  Branch-site models:                                                                                                                                                                       
  M1a[a](#nt102){ref-type="table-fn"}             54   −11097.30   1.544       ω~0~ = 0.040, ω~1~ = 1.000                  \-                                                               
                                                                            *p* ~0~ = 0.990, *p* ~1~ = 0.010                                                                                
  Model A null[b](#nt103){ref-type="table-fn"}    55   −11082.55   1.535       ω~0~ = 0.039, ω~1~ = 1.000                  \-                                                               
                                                                              ω~2a~ = 1.000, ω~2b~ = 1.000                                                                                  
                                                                            *p* ~0~ = 0.814, *p* ~1~ = 0.009                                                                                
                                                                           *p* ~2a~ = 0.175, *p* ~2b~ = 0.002                                                                               
  Model A alt                                     56   −11079.65   1.536       ω~0~ = 0.040, ω~1~ = 1.000       30 sites: 222, 230, 231,                       M1a                       2   2.174×10^−08^ [\*](#nt104){ref-type="table-fn"}
                                                                              ω~2a~ = 8.122, ω~2b~ = 8.122      252, 253,272, 302, 403,    Br-s A null[b](#nt103){ref-type="table-fn"}   1       0.016[\*](#nt104){ref-type="table-fn"}
                                                                            *p* ~0~ = 0.796, *p* ~1~ = 0.009        411 (at P\>0.95)                                                        
                                                                           *p* ~2a~ = 0.194, *p* ~2b~ = 0.002      296, 301, 324, 398                                                       
                                                                                                                      (at P\>0.99)                                                          
  **EcR**                                                                                                                                                                                   
  M0: one-ratio                                   53   −9701.18    1.661               ω = 0.029                                                                                            
  Branch-specific models:                                                                                                                                                                   
  Two-ratios                                      54   −9700.11    1.664       ω~0~ = 0.029, ω~1~ = 0.008                                                      M0                        1                        0.143

NOTE -- np is the number of parameters for each model, df is the degrees of freedom, while *p* is the proportion of sites in a given site class. Sites under positive selection are listed according to the site numbering of the *H. virescens* reference sequence, accession number AX383958.

Null model for test 1 of positive selection from Zhang *et al.* [@pone.0023416-Zhang1].

Null model for test 2 from Zhang *et al.* [@pone.0023416-Zhang1], a more stringent test of positive selection.

\*Significant (p-value\<5%).

Since coevolution of USP/RXR and EcR has been suggested, a branch model for which the Mecopterida lineage was set as foreground was also fit to a dataset composed of EcR ligand-binding domain sequences (EcR LBD dataset). Unlike USP/RXR, the foreground lineage demonstrated a decreased evolutionary rate (ω~m~ = 0.008) suggestive of stronger purifying selection ([Table 1](#pone-0023416-t001){ref-type="table"}). However, the LRT again indicated that the additional parameter did not yield a significantly better fit for the data than M0 (p = 0.143).

In a separate analysis, the lineage leading to the ancestor of Hymenoptera and Mecopterida was also freely estimated across both USP/RXR and EcR LBD sequences to determine if ω was elevated in the branch preceding the Mecopterida. However, in both datasets, ω was found to be below background (ω = 0.002 and ω = 0.0001, respectively; data not shown). The added parameter along that branch also did not yield a statistically better fit than M0 in either dataset (p = 0.072 and p = 0.084, respectively; data not shown). Overall, the branch model results, which only allow for a constant ω parameter across sites, suggest a weak trend where evolutionary rates might be elevated in USP/RXR along the Mecopterida lineage.

Branch-sites models {#s3b}
-------------------

The insignificance of the branch models tests may have been a consequence of a lack of positive selection, or positive selection acting only on a few sites within the gene along a given branch. Therefore branch-site models were implemented to estimate ω for each site in the USP/RXR LBD dataset along the specified foreground lineage. Branch-site Model A alt demonstrated that a proportion of sites in the USP/RXR gene have ω\>1 in the Mecopterida lineage ([Table 1](#pone-0023416-t001){ref-type="table"}). When compared to both M1a and branch-site Model A null, a more stringent test for positive selection, the LRT showed that the result was statistically significant (p = 2.174×10^−8^ and p = 0.016, respectively; [Table 1](#pone-0023416-t001){ref-type="table"}).

Using a Bayes Empirical Bayes (BEB) analysis in Model A alt, a class of sites with ω\>1 was identified ([Table S2](#pone.0023416.s006){ref-type="supplementary-material"}). Nine sites showed ω\>1 at a posterior probability (P) 0.99\>P\>0.95, and four at P\>0.99 ([Table 1](#pone-0023416-t001){ref-type="table"}). These sites mapped onto important structural and functional regions of the *H. virescens* USP [@pone.0023416-Billas1], [@pone.0023416-Billas2] crystal structure ([Figure 2A](#pone-0023416-g002){ref-type="fig"}). Of the 13 sites, one directly interacted with the ligand, and three others lay near ligand-binding sites ([Figure 2B](#pone-0023416-g002){ref-type="fig"}). Several sites lay within loop L1--3, and were involved in a hydrogen bond network with loop L11--12 and H3 which stabilizes the position of structural elements in Mecopterida USP/RXR ([Figure 2C](#pone-0023416-g002){ref-type="fig"}). Although only one site is known to form direct contact with EcR in the crystal structure of the heterodimer, several others fell within H9, a component of the dimerization core ([Figure 2D](#pone-0023416-g002){ref-type="fig"}). One site was located in the coactivator groove in loop L3--4, and two others lay in the region immediately adjacent to loop L5-S1. In addition, several of the sites with P\<0.95 also participated in interactions with EcR or were located beside sites which do so (supporting [Table S2](#pone.0023416.s006){ref-type="supplementary-material"}). Overall, these results indicate significantly elevated ω values along the Mecopterida lineage, indicating positive selection in regions of the protein important for ligand-binding, structural stability and dimerization.

![Location of putative positively selected sites in Mecopterida USP/RXR.\
(A) The distribution of sites inferred by PAML to be under positive selection (P\>0.95) along the Mecopterida branch (red spheres) across the crystal structure of *H. virescens* USP (green) with the ligand shown in grey. (B) Positively selected sites (red spheres) located near ligand binding sites (white sidechains). Only L230 directly interacts with the phospholipid ligand (grey) via van der Waals (red dashes). (C) The involvement of sites L203 and R231 in a hydrogen bond network between H3, loop L1--3 and H12. Hydrogen bonds are indicated by dotted blue lines and waters by small red dots (adapted from [@pone.0023416-Billas1]). The backbone and sidechains of residues not under positive selection are shown in white, with blue and red indicating nitrogen and oxygen respectively. (D) Polar interactions (dashed black lines) between the positively selected site E411 in USP loop L9--10 (green) and two arginine residues in H9 of the EcR (cyan). The images were created using PDB 1G2N and 1R1K, site numbering according to *H. virescens* USP (AX383958) and EcR (Y09009).](pone.0023416.g002){#pone-0023416-g002}

Random-sites models {#s3c}
-------------------

In order to compare among-site variation in selection pressure across USP/RXR functional domains between Mecopterida and Non-Mecopterida insects, random-sites models were conducted using PAML with datasets for each group which spanned all domains (USP/RXR A/B-LBD datasets) ([Table 2](#pone-0023416-t002){ref-type="table"}). M3 was the best fit for the Mecopterida and Non-Mecopterida datasets (p = 2.231×10^−127^ and p = 5.063×10^−47^, respectively), suggestive of significant among-site variation in ω. However, M3 is generally a better fit than the null M0 for most proteins as a single ω value across all sites is unlikely. For both the Mecopterida and Non-Mecopterida USP/RXR A/B-LBD datasets, tests of positive selection were not significant. Neither M2a nor M8 were a better fit than the respective null models M1a and M7 (p = 1.0 in both cases). Although no positively selected sites were identified in either Mecopterida or Non-Mecopterida USP/RXR, ω was in fact, elevated at particular regions of the gene ([Figure 3A, B](#pone-0023416-g003){ref-type="fig"}). Site class assignments were generally consistent across models, with peaks in ω occurring at the same codon sites, so only the results for M8 are presented (for M3 results see supporting [Figure S4](#pone.0023416.s004){ref-type="supplementary-material"}). Overall, more sites with elevated ω values were observed in the Mecopterida compared to the Non-Mecopterida.

![Posterior mean ω at each amino acid site across the USP/RXR gene.\
The values of ω as estimated by M8 in PAML (black line) and the Dual model of HyPhy (shaded blue) are shown for each codon site across the gene for the Mecopterida (A) and Non-Mecopterida (B). All sites inferred to be under positive selection by branch-site analysis the along the Mecopterida branch are shown as red boxes on the x-axis in of the Mecopterida plot. A schematic of USP/RXR secondary structure is shown above both plots to illustrate the position of each functional domain (A/B, C, and D) as well as the helices (H1--H12) and β sheets (S1--S2) of the ligand-binding domain. The schematic for the Mecopterida and Non-Mecopterida genes are based of the crystal structures of *H. virescens* and *B. tabaci* USP/RXR, respectively. Site numbering is based on the USP/RXR alignment given in supporting [figure S1](#pone.0023416.s001){ref-type="supplementary-material"}.](pone.0023416.g003){#pone-0023416-g003}

10.1371/journal.pone.0023416.t002

###### Parameter estimates for USP/RXR gene (A/B-LBD).

![](pone.0023416.t002){#pone-0023416-t002-2}

                                                                                                          LRT      
  ------------------------ ---- ----------- ------- ---------------------------------------------------- ----- --- --------------------------------------------------
  **Mecopterida:**                                                                                                 
  M0: one-ratio             35   −10594.82   1.532                       ω = 0.035                                 
  M1a: neutral              36   −10531.81   1.622               ω~0~ = 0.029, ω~1~ = 1.000                        
                                                              *p* ~0~ = 0.950, *p* ~1~ = 0.050                     
  M2a: selection            38   −10531.81   1.622        ω~0~ = 0.029, ω~1~ = 1.000, ω~2~ = 1.000        M1a   2                        1.000
                                                     *p* ~0~ = 0.950, *p* ~1~ = 0.002, *p* ~2~ = 0.048             
  M3: discrete (*K* = 3)    39   −10297.50   1.564        ω~0~ = 0.002, ω~1~ = 0.046, ω~2~ = 0.194        M0    4   2.231×10^−127^ [\*](#nt106){ref-type="table-fn"}
                                                     *p* ~0~ = 0.457, *p* ~1~ = 0. 450, *p* ~2~ = 0.092            
  M7: beta                  36   −10299.37   1.557                *p* = 0.417, *q* = 9.184                         
  M8: beta& ω               38   −10299.37   1.557                *p* = 0.417, *q* = 9.185                         
                                                           *p* ~0~ = 0.99999, (*p* ~1~ = 0.00001)         M7    2                        1.000
                                                                         ω = 1.000                                 
  **Non-Mecopterida:**                                                                                             
  M0: one-ratio             19   −5537.53    1.625                       ω = 0.016                                 
  M1a: neutral              20   −5496.97    1.789               ω~0~ = 0.013, ω~1~ = 1.000                        
                                                              *p* ~0~ = 0.967, *p* ~1~ = 0.033                     
  M2a: selection            22   −5496.97    1.789        ω~0~ = 0.013, ω~1~ = 1.000, ω~2~ = 1.000        M1a   2                        1.000
                                                     *p* ~0~ = 0.967, *p* ~1~ = 0.019, *p* ~2~ = 0.014             
  M3: discrete (*K* = 3)    23   −5426.21    1.735        ω~0~ = 0.001, ω~1~ = 0.020, ω~2~ = 0.119        M0    4   5.063×10^−47^ [\*](#nt106){ref-type="table-fn"}
                                                     *p* ~0~ = 0.521, *p* ~1~ = 0.394, *p* ~2~ = 0.085             
  M7: beta                  20   −5429.44    1.728               *p* = 0.300, *q* = 13.893                         
  M8: beta& ω               22   −5429.44    1.728               *p* = 0.300, *q* = 13.895                         
                                                          *p* ~0~ = 0. 99999, (*p* ~1~ = 0. 00001)        M7    2                        1.000
                                                                         ω = 1.000                                 

NOTE -- np is the number of parameters for each model, df is degrees of freedom, while *p* ~0~--*p* ~2~ is the proportion of sites in a given site class, and *p*, *q* describe the beta distribution.

\*Significant (p-value\<5%).

In the Mecopterida, there were a dramatic number of sites with increased ω in the LBD and D domain compared to the A/B and C domains ([Figure 3A](#pone-0023416-g003){ref-type="fig"}). These sites clustered within H1, the amino-terminal of loop L1--3, loop L3--4, H4, H7, H9 and H12 of the LBD (see supporting [Table S3](#pone.0023416.s007){ref-type="supplementary-material"} for details). Two sites were located beside ligand-binding sites and six other sites lay near residues which make dimerization contacts with EcR. Three sites in H12 lay near residues that interact with loop L1--3, H3 and the coactivator binding groove. Four sites lay in the coactivator binding groove, a region blocked in the Mecopterida as a consequence of contacts with H12 and loop L1--3 [@pone.0023416-Clayton1], [@pone.0023416-Billas1]. Several of the above sites also overlapped with those identified to be under positive selection along the Mecopterida lineage by branch-site analysis ([Figure 3A](#pone-0023416-g003){ref-type="fig"}, red squares). However, many sites differed between these types of analyses as branch-sites models compare the specified foreground branch to the background lineages, whereas random-sites models detect among-site rate variation within a clade. Unlike the Mecopterida, only a handful of sites with an elevated ω were observed in the Non-Mecopterida USP/RXR dataset ([Figure 3B](#pone-0023416-g003){ref-type="fig"}; supporting [Table S3](#pone.0023416.s007){ref-type="supplementary-material"}). These were primarily located in loop L1--3, a region with extensive sequence variation among Non-Mecopterida species [@pone.0023416-Hayward1], [@pone.0023416-Hayward2], [@pone.0023416-Maestro1]. There were also a few sites with elevated ω in H4, H9, H12 and the carboxy-terminal DNA-binding domain (C domain). Overall, among-site variation in the USP/RXR suggests that regions involved in dimerization and the blockage of the coactivator binding groove may be under relaxed constraint within the Mecopterida clade.

We also sought to examine among-site variation across EcR functional domains (EcR A/B-LBD datasets), as shifts in USP/RXR structure and function may also have affected the molecular evolution of EcR ([Table 3](#pone-0023416-t003){ref-type="table"}). As with USP/RXR, M3 also yielded the best fit for both Mecopterida and Non-Mecopterida datasets (p = 3.320×10^−97^ and 9.273×10^−32^, respectively). No positively selected sites were identified and neither M2a nor M8 were a better fit than the null models M1a and M7 for either dataset (p = 1.0 in both cases). As in the USP/RXR analyses, *d* ~N~/*d* ~S~ site-profiles were plotted for EcR based on the M8 result ([Figure 4A, B](#pone-0023416-g004){ref-type="fig"}; see supporting [Figure S4](#pone.0023416.s004){ref-type="supplementary-material"} for M3 results). Again, more sites with elevated ω values were observed in Mecopterida compared to the Non-Mecopterida.

![Posterior mean ω at each amino acid site across the EcR gene.\
The values of ω as estimated by M8 in PAML (black line) and the Dual model of HyPhy (shaded blue) are shown for each codon site across the gene for the Mecopterida (A) and Non-Mecopterida (B). A schematic of EcR secondary structure is shown above both plots to illustrate the position of each functional domain (A/B, C, and D) as well as the helices (H1--H12) and β sheets (S1--S3) of the ligand-binding domain. The schematic for the Mecopterida and Non-Mecopterida genes are based of the crystal structures of *H. virescens* and *B. tabaci* EcR, respectively. Site numbering is based on the EcR alignment given in supporting [figure S2](#pone.0023416.s002){ref-type="supplementary-material"}.](pone.0023416.g004){#pone-0023416-g004}

10.1371/journal.pone.0023416.t003

###### Parameter estimates for EcR gene (A/B-LBD).

![](pone.0023416.t003){#pone-0023416-t003-3}

                                                                                                         LRT      
  ------------------------ ---- ----------- ------- --------------------------------------------------- ----- --- -------------------------------------------------
  **Mecopterida:**                                                                                                
  M0: one-ratio             35   −10215.15   1.690                       ω = 0.029                                
  M1a: neutral              36   −10198.25   1.737              ω~0~ = 0.028, ω~1~ = 1.000                        
                                                             *p* ~0~ = 0.988, *p* ~1~ = 0.012                     
  M2a: selection            38   −10198.25   1.737       ω~0~ = 0.028, ω~1~ = 1.000, ω~2~ = 1.000        M1a   2                        1.000
                                                     *p* ~0~ = 0.988, *p* ~1~ = 0.008, *p* ~2~ = 0.003            
  M3: discrete (*K* = 3)    39   −9987.57    1.721       ω~0~ = 0.000, ω~1~ = 0.023, ω~2~ = 0.107        M0    4   3.320×10^−97^ [\*](#nt108){ref-type="table-fn"}
                                                     *p* ~0~ = 0.404, *p* ~1~ = 0.393, *p* ~2~ = 0.203            
  M7: beta                  36   −9988.17    1.725               *p* = 0.357, *q* = 10.575                        
  M8: beta& ω               38   −9988.17    1.719               *p* = 0.357, *q* = 10.575                        
                                                         *p* ~0~ = 0. 99999, (*p* ~1~ = 0. 00001)        M7    2                        1.000
                                                                         ω = 1.000                                
  **Non-Mecopterida:**                                                                                            
  M0: one-ratio             19   −6357.75    2.022                       ω = 0.016                                
  M1a: neutral              20   −6357.74    2.029              ω~0~ = 0.016, ω~1~ = 1.000                        
                                                             *p* ~0~ = 0.999, *p* ~1~ = 0.001                     
  M2a: selection            22   −6357.74    2.029       ω~0~ = 0.016, ω~1~ = 1.000, ω~2~ = 1.000        M1a   2                        1.000
                                                     *p* ~0~ = 0.999, *p* ~1~ = 0.001, *p* ~2~ = 0.000            
  M3: discrete (*K* = 3)    23   −6281.96    2.085       ω~0~ = 0.000, ω~1~ = 0.013, ω~2~ = 0.0737       M0    4   9.273×10^−32^ [\*](#nt108){ref-type="table-fn"}
                                                     *p* ~0~ = 0.412, *p* ~1~ = 0.435, *p* ~2~ = 0.153            
  M7: beta                  20   −6282.43    2.085               *p* = 0.326, *q* = 17.365                        
  M8: beta& ω               22   −6282.43    1.906               *p* = 0.326, *q* = 17.365                        
                                                          *p* ~0~ = 0.99999, (*p* ~1~ = 0.00001)         M7    2                        1.000
                                                                         ω = 1.000                                

NOTE -- np is the number of parameters for each model, df is the degrees of freedom, while *p* ~0~--*p* ~2~ is the proportion of sites in a given site class, and *p*, *q* describe the beta distribution.

\*Significant (p-value\<5%).

In the Mecopterida there was a concentration of sites with elevated ω in carboxy-terminal D domain and the dimerization core ([Figure 4A](#pone-0023416-g004){ref-type="fig"}; supporting [Table S4](#pone.0023416.s008){ref-type="supplementary-material"}). Several sites were located in loop L2--3, a region which varies in sequence and is not well modeled in insect crystal structures [@pone.0023416-Iwema1], [@pone.0023416-Billas2], [@pone.0023416-Carmichael1]. One site in the β sheet region was located between several ligand-binding sites. Six sites in H7, H9 and H10 were located near residues which make contact with USP/RXR. Three sites at the amino-terminus of H10 fell into a region associated with dimerization and ligand-binding specificity. Substitutions in this region confer *Aedes aegypti*-like sensitivity to ecdysone to *D. melanogaster* EcR, which is normally only responsive to 20E [@pone.0023416-Wang1]. As in USP/RXR, there was far less among-site variation in ω across Non-Mecopterida EcR ([Figure 4B](#pone-0023416-g004){ref-type="fig"}; supporting [Table S4](#pone.0023416.s008){ref-type="supplementary-material"}). Three sites in loop L9--10 and H10 align with the aforementioned region in *A. aegypti* which affects sensitivity to ecdysone. In addition, three other sites were located near ligand-binding sites. Overall, these results indicate that the dimerization core of EcR, like that of USP/RXR, may also be under relaxed constraint in the Mecopterida clade.

Random-sites models that incorporate variation in *d* ~S~ {#s3d}
---------------------------------------------------------

Traditional *d* ~N~/*d* ~S~ methods estimate *d* ~N~ and *d* ~S~ as a ratio, and assume a constant *d* ~S~ [@pone.0023416-Nielsen1], [@pone.0023416-Yang5]. However, this may be an unrealistic assumption for some datasets, and synonymous substitutions are thought to be under selection in flies, nematodes and yeast [@pone.0023416-Sharp1]. Recently, models have been developed to allow for the independent estimation of *d* ~N~ and *d* ~S~. To examine the effect of among-site variation in *d* ~S~ on the estimation of ω in our datasets, we implemented several such models in the HyPhy software package (see supporting [Tables S5](#pone.0023416.s009){ref-type="supplementary-material"}, [S6](#pone.0023416.s010){ref-type="supplementary-material"}, [S7](#pone.0023416.s011){ref-type="supplementary-material"} for parameter estimates). Overall, the results of HyPhy are similar to those of PAML ([Figure 3](#pone-0023416-g003){ref-type="fig"} and [4](#pone-0023416-g004){ref-type="fig"}, shaded blue). For all datasets the regions with elevated ω were consistent between both methods, with some variation in peak height. The baseline of the *d* ~N~/*d* ~S~ site-profile plots appears elevated for the PAML M8 results compared to the HyPhy results because the former employs BEB methods, which account for sampling error in the maximum likelihood *d* ~N~/*d* ~S~ (and proportion) estimates, whereas the latter does not.

For both USP/RXR and EcR A/B-LBD datasets in the Mecopterida site-to-site variation in *d* ~S~ was significant when the Dual model compared to the Nonsynonymous null model where *d* ~S~ is fixed at all sites (p = 0.037 and 0.025 respectively; [Table 4](#pone-0023416-t004){ref-type="table"}). However, this was not the case with either dataset for the Non-Mecopterida in which significant variation in *d* ~S~ among sites was not detected. Given the small size of the Non-Mecopterida datasets and the limited taxonomic sampling, we may lack the statistical power to estimate *d* ~N~ and *d* ~S~ accurately in this group. Overall the HyPhy results indicate that despite significant among-site variation in *d* ~S~ within some datasets, the patterns in substitution rates across both the USP/RXR and EcR A/B-LBD datasets are the same when *d* ~N~ and *d* ~S~ are estimated independently.

10.1371/journal.pone.0023416.t004

###### Synonymous rate variation among sites (HyPhy).

![](pone.0023416.t004){#pone-0023416-t004-4}

                                       MG94×REV Nonsynonymous GDD 3   MG94×REV Dual GDD 3×3   LRT                  
  ----------------------------------- ------------------------------ ----------------------- ----- ----------- --- ----------------------------------------
  Mecopterida USP/RXR (A/B-LBD)                     43                      −10237.10         47    −10232.01   4   0.037[\*](#nt110){ref-type="table-fn"}
  Non-Mecopterida USP/RXR (A/B-LBD)                 27                      −5415.51          31    −5410.81    4                   0.052
  Mecopterida EcR (A/B-LBD)                         43                      −9946.36          47    −9940.77    4   0.025[\*](#nt110){ref-type="table-fn"}
  Non-Mecopterida EcR (A/B-LBD)                     27                      −6269.43          31    −6269.40    4                   0.999

NOTE -- np is the number of parameters for each model, Log *L* is the logarithm of the maximum likelihood value, and df is the degrees of freedom used in the LRT calculation when comparing Dual versus Nonsynonymous models.

\*Significant (p-value\<5%).

Discussion {#s4}
==========

Our results indicate that USP/RXR is under positive selection along the branch leading to the Mecopterida, with positively selected sites tending to be located in regions involved in ligand-binding, interactions with loop L1--3 and dimerization ([Figure 2](#pone-0023416-g002){ref-type="fig"}). Furthermore, random-sites analyses showed that ω was elevated across the Mecopterida clade compared to the Non-Mecopterida for both USP/RXR and EcR ([Figure 3](#pone-0023416-g003){ref-type="fig"}, [4](#pone-0023416-g004){ref-type="fig"}). Sites with elevated ω tended to be concentrated in components of the dimerization interface, suggesting relaxed constraint in this region among the Mecopterida ([Figure 3A](#pone-0023416-g003){ref-type="fig"}, [4A](#pone-0023416-g004){ref-type="fig"}). There were also several sites with elevated ω in both proteins located near sites important for ligand-binding. Overall, these results are consistent with the acquisition of functional gains with respect to ligand-binding and dimerization in USP/RXR in the Mecopterida.

Although higher relative substitution rates have been shown in the Mecopterida, our study represents the first instance of codon-based likelihood phylogenetic methods being used to estimate *d* ~N~/*d* ~S~ ratios across the USP/RXR gene [@pone.0023416-Iwema1], [@pone.0023416-Bonneton1], [@pone.0023416-Bonneton2]. The elevated *d* ~N~/*d* ~S~ values we observed along the Mecopterida lineage are consistent the with high amino acid substitution rates reported previously [@pone.0023416-Bonneton2]. Strikingly, the regions identified as under either positive selection or relaxed constraint play a major role in defining the structure of USP/RXR in the Mecopterida. In order to function as ligand-dependent transcription factors, nuclear receptors must undergo a shift in the position of H12 upon ligand binding to generate an interface for coactivators [@pone.0023416-Wurtz1]. Residues in the AF-2 domain and the coactivator binding groove become accessible, allowing the binding of transcriptional coactivators [@pone.0023416-Wurtz1], [@pone.0023416-Bourguet1], [@pone.0023416-Bourguet2]. However, the position of loop L1--3 in Mecopterida USP/RXR prevents shifts in the H12 position, regardless of ligand binding, thereby locking the receptor in an antagonist conformation [@pone.0023416-Clayton1], [@pone.0023416-Billas1].

These conformational changes are evident in the evolutionary history of USP/RXR. Studies using reconstruction and homology modeling of the ancestral Mecopterida USP/RXR have investigated the origin of an expanded interface with EcR; these novel contacts were a result of intramolecular epistasis whereby the position of loop L1--3 created torsion in the protein, shifting the position of loop L8--9 closer to EcR [@pone.0023416-Iwema2]. We identified sites under positive selection along the Mecopterida lineage in the loop L1--3 region which interact with the ligand and contribute to a hydrogen bond network with H3 and loop L11--12 ([Figure 2](#pone-0023416-g002){ref-type="fig"}). Site R385 (according to *H. virescens*) in H9, a Mecopterida-specific contact with EcR, was also shown to be under positive selection along the Mecopterida lineage by branch-site analysis. Furthermore, sites surrounding R385 appear to be under relaxed constraint within the Mecopterida in our random-sites analysis ([Figure 3A](#pone-0023416-g003){ref-type="fig"}). This was also the case at sites near S447 (according to *H. virescens*) in EcR H7, which establishes contact with USP/RXR R385 ([Figure 4A](#pone-0023416-g004){ref-type="fig"}). This is consistent with the acquisition of an expanded dimer interface between USP/RXR and EcR in the Mecopterida.

It has been postulated that the expansion of this interface may reflect a strengthened inter-dependency between USP/RXR and EcR for activation [@pone.0023416-Iwema2]. Based on structural studies, Diptera and Lepidoptera USP/RXR cannot bind the canonical NR-box LXXLL motif of coactivators because H12 sits within the binding site [@pone.0023416-Clayton1], [@pone.0023416-Billas1], [@pone.0023416-Darimont1], [@pone.0023416-Nolte1]. Our random-sites analyses revealed sites with elevated ω within both H12 and the coactivator binding groove (loop L3--4, H4) of Mecopterida USP/RXR ([Figure 3A](#pone-0023416-g003){ref-type="fig"}). These results are consistent with a change in the function of these regions, leading to a relaxation of evolutionary constraint at these sites. Without a novel interaction interface, such shifts may render USP/RXR more dependent on partner proteins for coactivator recruitment [@pone.0023416-Clayton1], [@pone.0023416-Moras1]. Interestingly, similar shifts in cofactor interfaces may also have occurred in EcR. An elevation in ω was observed in the carboxy-terminal D domain of Mecopterida EcR, a region implicated in corepressor binding in the mammalian thyroid hormone receptor [@pone.0023416-Hrlein1]. The biological significance of an increased dependence on partner proteins remains unclear, but there are some possibilities.

One mechanism which may have led to increased evolutionary rates in the Mecopterida dimer interface is the increased recruitment of USP/RXR as a protein hub. In *D. melanogaster* and *A. aegypti*, USP/RXR is known to interact with many other proteins such as Seven-up (Svp), hormone receptors 38 and 78 (HR38, 78), as well as the *Methoprene-tolerant* gene product (MET) [@pone.0023416-Zhu1], [@pone.0023416-Hirai1], [@pone.0023416-Baker1], [@pone.0023416-Li1]. Hubs at major branch points in protein-protein interaction networks tend to display evidence of elevated ω rates and positive selection [@pone.0023416-AlvarezPonce1]. The strong positive selection we observed could be congruent with such a role for USP/RXR. Indeed network level analyses suggest that several proteins in the ecdysone cascade, including USP/RXR and EcR, experienced an increased evolutionary rate at the base of the Mecopterida [@pone.0023416-Bonneton3]. However, it is unclear if this is due to functional divergence or simply parallel accelerations across the network in order to maintain interactions. Similarly, the elevated evolutionary rates we detected within the dimerization interface of the Mecopterida clade may be the result of a combination of strong positive selection and purifying selection, not simply relaxed constraint. Substitutions in this region may affect both the affinity and specificity of protein-protein interaction, but might not necessarily improve both of these attributes. In fact, there may be some trade-offs, whereby some substitutions confer a tighter affinity for one protein, but reduce the specificity of interactions with other partner proteins. However, at this point the residues involved in the interface between USP/RXR and alternate partner proteins remain unknown. Furthermore, the network of USP/RXR protein interactions is unknown in lower insects, and a lack of sequence data from outside of the Holometabola prevented the analysis of alternate USP/RXR partners in our study. Therefore, further work is needed to clarify the protein-protein interactions of USP/RXR among insect lineages.

The link between shifts in USP/RXR molecular evolution and shifts in physiology at the organismal level has been difficult to establish. Studies have shown that the role of USP/RXR in molting and metamorphosis seems conserved across insect lineages [@pone.0023416-Martn1], [@pone.0023416-Tan1]. Differences in the reproductive biology of insect lineages have been largely unexplored, but may hold more promise [@pone.0023416-Bonneton2]. Although there are some exceptions, JH is the major hormone regulating reproductive events such as vitellogenesis in most insect groups [@pone.0023416-Raikhel1]. In the Lepidoptera, there is a decreased dependence on JH for vitellogenesis, whereas in the Diptera, vitellogenesis is largely regulated by ecdysteroids. In many Diptera, JH is instead required during previtellogenesis in order for tissues to acquire competence to 20E [@pone.0023416-Zhu2]. The site of vitellogenin production also varies between insect lineages. In many basal lineages, the fat body appears to be the sole site of synthesis, in others, both the ovary and fat body contribute, whereas in some Diptera, the ovary may be the major source [@pone.0023416-Valle1], [@pone.0023416-Huebner1]. Thus, a stronger relationship between USP/RXR and EcR for heterodimerization is also consistent with the increased dependence on ecdysteroids for the regulation of reproductive events in some higher insects.

USP/RXR is involved in many reproductive events in higher insects such as egg chamber formation and chorion gene expression in *Drosophila* and regulation of the cyclicity of vitellogenesis in *A. aegypti* [@pone.0023416-Zhu1], [@pone.0023416-Carney1], [@pone.0023416-Bernardi1]. Ovarian morphology and structure also varies among insects groups, particularly in lineages between the basal Non-Mecopterida and Mecopterida, such as the Paraneoptera, Coleoptera and Neuropterida [@pone.0023416-Bning1], [@pone.0023416-Trauner1]. The genes which control reproductive events differ among ovary types as well. For example, the set of genes responsible for chorion formation in meroistic ovaries differs from those in the panoistic type ovary, found in the basal Non-Mecopterida [@pone.0023416-Bning1], [@pone.0023416-Irles1]. In addition, USP/RXR and EcR are highly expressed in nurse cells, a cell type unique to the meroistic ovary [@pone.0023416-Carney1], [@pone.0023416-Oro2]. Thus, shifts in ovarian morphology, development and the endocrine control of reproduction may have also required plasticity in nuclear receptors. However, the role of USP/RXR in the reproductive events of basal insects is not known, as is its role in highly derived insects contained within the Mecopterida.

Finally, it is important to consider caveats of *d* ~N~/*d* ~S~ based methods. In recent years a debate has emerged regarding the use of these methods to reliably detect positive selection because of potential false positives identified in some genes [@pone.0023416-Nozawa1], [@pone.0023416-Yang8], [@pone.0023416-Yang9]. However, in our dataset the positively selected sites identified along the Mecopterida lineage by branch-site analysis were consistent with the structural and experimental work of other researchers. Furthermore, we used patterns of *d* ~N~/*d* ~S~ variation across sites, not merely positive selection, to draw conclusion regarding USP/RXR evolution. Finally, we utilised new models which allow for among-site variation in *d* ~S~ and observed strikingly similar patterns in site-specific evolutionary rates even across these different methods of estimating *d* ~N~/*d* ~S~.

While computational methods, such as those used here, are useful in identifying patterns of selective constraint acting on USP/RXR across insect lineages, the functional changes generated by the substitutions identified here remain speculative. For example, we cannot clarify whether or not sites under positive selection near the ligand-binding pocket are related to the potential acquisition of JH binding. Further experimental evidence is required to understand the physiological implications of these evolutionary shifts. In particular, examining the role of USP/RXR in the reproductive system of more basal insects, for example, with the use of techniques such as RNAi, may shed light on these issues. Additionally, the apparent increase in the rigidity of Mecopterida USP/RXR makes this an excellent system for the use molecular dynamics studies to investigate subtle shifts in structure. Future work in this area may yield more insights into the role of the positively selected sites identified in this study which impact the conformation of the receptor complex. Overall, our results highlight the need for more comparative physiological work in the basal Non-Mecopterida and provide critical site-specific information for the design of future site directed mutagenesis studies on USP/RXR.

Supporting Information {#s5}
======================

###### 

**Alignment of insect USP/RXR sequences.** The green and orange bars indicate the Mecopterida and Non-Mecopterida taxa, respectively, and the schematic below the alignment shows the USP/RXR domain structure. For ease of viewing one alignment is shown. However, the complete alignment was not used for analysis as some regions do not align (e.g. D domain). Only the carboxy-terminal E/F (\*), or ligand-binding domain, was used for the branch and branch-sites analyses reported in [table 1](#pone-0023416-t001){ref-type="table"}. Arrows indicate where poorly aligned regions and major gaps were deleted. Full length sequences were used for the random-sites and HyPhy analyses where the larger dataset was split into Mecopterida and Non-Mecopterida only datasets in order to compare evolutionary rates between the two groups. For clarity, site numbering for the full length Mecopterida (green) and Non-Mecopterida (orange) datasets is shown above and below the alignment, respectively. Note that species names have been abbreviated to six characters, complete names can be found in supporting [table S1](#pone.0023416.s005){ref-type="supplementary-material"}.

(PDF)

###### 

Click here for additional data file.

###### 

**Alignment of insect EcR sequences.** The green and orange bars indicate the Mecopterida and Non-Mecopterida taxa, respectively, and the schematic below the alignment shows the EcR domain structure. For ease of viewing one alignment is shown. However, the complete alignment was not used for analysis as some regions do not align (e.g. D domain). Only the carboxy-terminal E (\*), or ligand-binding domain, was used for the branch analyses reported in [table 1](#pone-0023416-t001){ref-type="table"}. Arrows indicate where poorly aligned regions and major gaps were deleted. Full length sequences were used for the random-sites and HyPhy analyses where the larger dataset was split into Mecopterida and Non-Mecopterida only datasets in order to compare evolutionary rates between the two groups. For clarity, site numbering for the full length Mecopterida (green) and Non-Mecopterida (orange) datasets is shown above and below the alignment, respectively. Note that species names have been abbreviated to six characters, complete names can be found in supporting [table S1](#pone.0023416.s005){ref-type="supplementary-material"}.

(PDF)

###### 

Click here for additional data file.

###### 

**USP/RXR and EcR gene trees.** Gene trees for USP/RXR and EcR were generated using the alignment of LBD sequences given in supporting [figures S1](#pone.0023416.s001){ref-type="supplementary-material"} and [S2](#pone.0023416.s002){ref-type="supplementary-material"}. Maximum-likelihood trees for USP/RXR (A) and EcR (C) were constructed in PhyML [@pone.0023416-Guindon1] using the WAG substitution model, with four rate categories to estimate the gamma parameter shape. Neighbor-joining [@pone.0023416-Saitou1] trees for USP/RXR (B) and EcR (D) were constructed in MEGA 4 [@pone.0023416-Tamura1] using the Poisson correction model, with the pair-wise deletion of gaps. For all analyses 100 bootstrap replicates were performed, and nodes with values less than 60 were later collapsed. Each tree was then rooted along the branch leading to *B. germanica* and *L. migratoria*. Note that species names have been abbreviated, see supporting [table S1](#pone.0023416.s005){ref-type="supplementary-material"}.

(PDF)

###### 

Click here for additional data file.

###### 

***d*** **~N~/** ***d*** **~S~ site-profile plots for PAML model M3.** The values of ω as estimated by M3 in PAML using the NEB method are shown for each codon site across Mecopterida USP/RXR (A), Non-Mecopterida USP/RXR (B), Mecopterida EcR (C) and Non-Mecopterida EcR (D). A schematic of USP/RXR and EcR secondary structure is shown above each plot to illustrate the position of each functional domain (A/B, C, and D) as well as the helices (H1--H12) and β sheets of the ligand-binding domain. These schematics are based on the crystal structure of each gene in *H. virescens* and *B. tabaci*. Site numbering is the same as [figures 3](#pone-0023416-g003){ref-type="fig"} and [4](#pone-0023416-g004){ref-type="fig"}.

(PDF)

###### 

Click here for additional data file.

###### 

**Sequence data information.**

(DOC)

###### 

Click here for additional data file.

###### 

**Putative positively selected sites.**

(DOC)

###### 

Click here for additional data file.

###### 

**Sites with elevated ω in the USP/RXR gene.**

(DOC)

###### 

Click here for additional data file.

###### 

**Sites with elevated ω in the EcR gene.**

(DOC)

###### 

Click here for additional data file.

###### 

**Statistics for rate distributions inferred by HyPhy.**

(DOC)

###### 

Click here for additional data file.

###### 

**Rate class assignment for MG94×REV Nonsynonymous GDD 3.**

(DOC)

###### 

Click here for additional data file.

###### 

**Rate class assignment for MG94×REV Dual GDD 3×3.**

(DOC)

###### 

Click here for additional data file.

We thank Jingjing Du and Cameron J. Weadick of the department of Ecology and Evolutionary Biology at the University of Toronto for their insightful comments, technical expertise and critiques of earlier drafts.

**Competing Interests:**The authors have declared that no competing interests exist.

**Funding:**This study was supported by Natural Sciences and Engineering Council of Canada (NSERC) Discovery grants (SST and BSWC) and an Early Researcher Award to BSWC. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

[^1]: Conceived and designed the experiments: EFH SST BSWC. Performed the experiments: EFH BSWC. Analyzed the data: EFH SST BSWC. Contributed reagents/materials/analysis tools: SST BSWC. Wrote the paper: EFH SST BSWC.
