Molecular Population Genetics of Accessory Gland Protein Genes and Testis-Expressed Genes in Drosophila mojavensis and D. arizonae

Abstract

Molecular population genetic investigation of Drosophila male reproductive genes has focused primarily on melanogaster subgroup accessory gland protein genes (Acp's). Consistent with observations from male reproductive genes of numerous taxa, Acp's evolve more rapidly than nonreproductive genes. However, within the Drosophila genus, large data sets from additional types of male reproductive genes and from different species groups are lacking. Here we report findings from a molecular population genetics analysis of male reproductive genes of the repleta group species, Drosophila arizonae and D. mojavensis. We find that Acp's have dramatically higher average pairwise K_a/K_s (0.93) than testis-enriched genes (0.19) and previously reported melanogaster subgroup Acp's (0.42). Overall, 10 of 19 Acp's have K_a/K_s > 1 either in nonpolarized analyses or in at least one lineage of polarized analyses. Of the nine Acp's for which outgroup data were available, average K_a/K_s was considerably higher in D. mojavensis (2.08) than in D. arizonae (0.87). Contrasts of polymorphism and divergence suggest that adaptive protein evolution at Acp's is more common in D. mojavensis than in D. arizonae.

MOLECULAR studies in a diverse array of animal taxa suggest that genes involved in reproduction evolve at an accelerated rate relative to other genes (reviewed in Swanson and Vacquier 2002). Positive selection has been inferred for some proteins (Swanson and Vacquier 1995; Metz and Palumbi 1996; Sutton and Wilkinson 1997; Wyckoff  et al. 2000; Torgerson  et al. 2002), although population genetic data are sufficiently sparse to leave unresolved the question of the relative importance of directional selection vs. genetic drift in reproduction-related proteins compared to other protein classes. In any case, rapid phenotypic/molecular evolution of reproductive characters/genes is consistent with the notion that male-male and male-female interactions may contribute to the rapid divergence between populations and the evolution of reproductive isolation (Eberhard 1996; Rice 1998).

Molecular evolutionary investigation of Drosophila reproduction has focused on male accessory gland protein genes (Acp's) of melanogaster subgroup species. The number of putative Acp's in these species is on the order of 83 (Swanson  et al. 2001), although <20 have extensive experimental support (Schäfer 1986; DiBenedetto  et al. 1987; Chen  et al. 1988; Monsma and Wolfner 1988; Wolfner  et al. 1997). Genetic analysis has shown that Acp's contribute to proper sperm storage (Neubaum and Wolfner 1999; Tram and Wolfner 1999; Chapman  et al. 2000), normal ovulation and oviposition (Herndon and Wolfner 1995; Heifetz  et al. 2000), increased egg-laying rates, and reduced female receptivity (Chen  et al. 1988; Aigaki  et al. 1991; Kalb  et al. 1993; Chapman  et al. 2003; Liu and Kubli 2003). Acp's show higher rates of protein divergence (Aguadé 1997, 1998, 1999; Tsaur and Wu 1997; Tsaur  et al. 1998; Begun  et al. 2000; Swanson  et al. 2001) and protein polymorphism (Coulthart and Singh 1988; Begun  et al. 2000) compared to “average” proteins in Drosophila melanogaster and D. simulans (e.g., Begun  et al. 2000). Less energy has been devoted to population genetic investigation of male reproductive genes primarily expressed in testes (but see Duvernell and Eanes 2000; Parsch  et al. 2001a). However, a few analyses suggest that Drosophila testis-expressed genes evolve quickly (Parsch  et al. 2001b; Meiklejohn  et al. 2004; Richards  et al. 2005) and may sometimes be associated with evolution of novel function (Long and Langley 1993; Nurminsky  et al. 1998; Betrán and Long 2003).

Because our current population genetic understanding of Drosophila is dominated by data from melanogaster subgroup species, we have no way of knowing whether the patterns of polymorphism and divergence or the functional biology of reproduction-related proteins will be similar in other Drosophila species (Wagstaff and Begun 2005). Given the hypothesis that the dynamics of certain male reproduction-related proteins may be driven by male-male and male-female postcopulatory interactions, one strategy for furthering our understanding of the evolution of these proteins is to investigate Drosophila species having different reproductive biology from D. melanogaster and D. simulans. D. arizonae and D. mojavensis are cactophilic fly species within the mulleri complex of the repleta group. As members of the subgenus Drosophila, these desert Drosophila are ∼40–60 million years diverged from D. melanogaster and other Sophophora subgenus flies (Powell and DeSalle 1995).

A major difference in the reproductive biology of desert Drosophila vs. D. melanogaster is that remating occurs more frequently in desert Drosophila. Within 24 hr of an initial mating, 95% of D. arizonae and D. mojavensis females tend to remate, while only 2% of D. melanogaster females remate in this same time period (Markow 1982, 1996). Frequent remating favors competition between male ejaculates, whereas infrequent remating would be more likely to favor genotypes successfully obtaining initial access to females (e.g., Markow 2002). Data from Drosophila species suggest that there is a positive correlation between high female remating rates and exaggerated ejaculates in the form of either sperm gigantism or excessive ejaculate donation to female tissues (Markow 2002). Although both desert Drosophila species discussed here contribute large ejaculate donations to ovaries, D. arizonae and D. mojavensis contribute small and large donations, respectively, to female somatic tissues (Pitnick  et al. 1997). Experiments in D. melanogaster revealed no detectable incorporation of ejaculate-derived material into female somatic or ovarian tissues (Pitnick  et al. 1997). While ejaculate donations are often perceived to be of nutritive value, a cost to remating has been observed in D. mojavensis females, suggesting the possibility of sexual conflict (Etges and Heed 1992). Another major difference in the reproductive biology of repleta group vs. melanogaster subgroup flies is that repleta group males require significantly more time to reach sexual maturity. For example, D. arizonae and D. mojavensis require 4–5 days posteclosion to reach maturity, compared to 1–2 days for D. melanogaster males (Pitnick  et al. 1995).

Data on natural variation in reproductive traits suggest a more dynamic postmating interaction between the sexes in desert Drosophila compared to melanogaster subgroup flies. Immediately after mating, a pronounced insemination reaction occurs in the female reproductive tract of D. arizonae and D. mojavensis (Patterson 1947; Patterson and Stone 1952) but is absent in D. melanogaster (Wheeler 1947; Markow and Ankney 1988). The reaction manifests itself as a large mass within the vaginal pouch and acts as a barrier that prevents remating for the several hours that it persists (Patterson 1947; Knowles and Markow 2001). Seminal fluid proteins may be the primary male contributor to this phenotype, as it is triggered in the absence of live spermatozoa (Patterson 1947). Comparisons between desert Drosophila species, as well as between different populations within species, show that postcopulatory male-female interactions change across short evolutionary time periods. For example, heterospecific matings between D. arizonae and D. mojavensis trigger an exaggerated insemination reaction that is both harder and longer lasting than that of the respective conspecific matings of either species (Patterson 1947). Moreover, both D. arizonae and D. mojavensis show larger and longer insemination reactions in interpopulation vs. intrapopulation crosses (Knowles and Markow 2001) within species. Further evidence of rapid evolution of reproductive traits comes from the observation that D. mojavensis shows significant among-population variation in the correlated traits of male sperm size and female sperm-storage organ length (Pitnick  et al. 2003).

These data support the idea that properties of ejaculates or ejaculate-female interactions evolve very quickly in desert Drosophila, possibly as a result of antagonistic coevolution between the sexes (Rice 1996, 1998) and/or cryptic female choice (Eberhard 1996). We should expect such elaboration of ejaculate characteristics to extend to the molecular level. The purpose of this study is to add a molecular framework to investigation of desert Drosophila reproduction. First, we report the composition of D. mojavensis male reproductive tract cDNA libraries relative to the gene annotations of D. melanogaster. Many of these data are presented as supplementary online material (http://www.genetics.org/supplemental). Second, we report results from molecular and evolutionary analyses of genes expressed in male reproductive tracts in D. mojavensis and D. arizonae and compare these results to those previously reported from D. melanogaster/D. simulans.

MATERIALS AND METHODS

D. mojavensis reproductive tract library:

Poly(A)-enriched mRNA was prepared with the MicroPoly(A)-Pure kit (Ambion, Austin, TX) from 50 whole reproductive tracts of adult male D. mojavensis flies. First-strand cDNA was reverse transcribed with the SMART PCR cDNA synthesis system reagents and protocol (CLONTECH, Palo Alto, CA). Second-strand product was produced with the Expand high-fidelity polymerase system (Roche Molecular Biochemicals, Indianapolis). Cycling parameters were programmed as instructed by the manufacturer, including a 4-min extension step for 10 total cycles. The second-strand product was cloned into the TOPO vector (Invitrogen, San Diego) and used for bacterial transformations according to the manufacturer's instructions. Colony PCR was carried out using cloning-vector-derived primers (M13 reverse and T7) on 480 colonies (i.e., five 96-well plates). The resulting PCR products were purified prior to sequencing with M13R and T7 primers on an Applied Biosystems (Foster City, CA) 377 automated sequencer (ABI, Columbia, MD). These sequences included 54 unique transcripts. Expressed sequence tags (ESTs) from this library can be found under accession nos. DR033184–DR033386 and DR033894–DR033895.

Preliminary expression analysis and D. mojavensis testis cDNA library production:

Dot blots prepared from PCR products of the 54 unique clones were hybridized separately to ³²P-labeled cDNAs derived from D. mojavensis accessory glands and testes. Hybridizations were carried out at 65° in a buffer consisting of 0.5 m NaPi (pH 7.2), 7% SDS, 1 mm EDTA. Filters were washed at 60° with buffer at 40 mm NaPi, 1% SDS, and 1 mm EDTA. Comparison of signal intensities from hybridization of labeled accessory gland vs. testis cDNA suggested that the majority of the clones represented accessory gland transcripts.

To increase the sample size of testis-enriched transcripts we made a testis cDNA library. This library was produced as described above for whole reproductive tracts, but with 50 D. mojavensis dissected testes as the source tissue. This library was sequenced to the point of producing 118 unique ESTs. ESTs from the testis library can be found under accession nos. DR033387–DR033542.

BLAST methodology and characterization of amino acid sequences:

All unique ESTs were compared to D. melanogaster through a pipeline of BLAST analyses to one or more FlyBase Release 3.1 databases (Altschul  et al. 1997). Default BLAST parameters were used except that the cutoff value for significance was set to E = 0.01. The pipeline started with BLASTp (protein to predicted D. melanogaster proteins) queries of all ESTs for which an open reading frame (ORF) was well established (as described below). ESTs that returned significant (E < 1e-8) D. melanogaster sequences were not queried further. The remaining ESTs were BLASTx (nucleotide to protein) queried to the same D. melanogaster database. Once again, ESTs that returned small E-values were not queried further. This pipeline continued through tBLASTx (nucleotide to nucleotide query, using all six possible protein translations of the sequences) and BLASTn (nucleotide to nucleotide) queries of predicted D. melanogaster genes and chromosome arms. For the ESTs that returned no D. melanogaster sequences at E < 0.0001, the NCBI whole-genome shotgun (wgs) database was tBLASTx queried with the same default parameters (Altschul  et al. 1997). The NCBI wgs database includes many complete genomes, including D. pseudoobscura and the mosquito, Anopheles gambiae. All D. mojavensis ESTs were also tBLASTx or BLASTn queried (BLASTn was used only if tBLASTx failed to return sequences of E < 0.0001) to the D. melanogaster dbEST database using default BLAST parameters and an E-score cutoff of 0.01. Finally, we queried the SignalP 3.0 (Nielsen and Krogh 1998; Bendtsen  et al. 2004) and NCBI CDD (Marchler-Bauer  et al. 2003) servers with amino acid sequences corresponding to ESTs with identifiable ORFs to identify the presence of signal peptides and conserved domains, respectively.

A subset of genes isolated from both libraries was scrutinized in greater detail to winnow candidates for population genetic analysis. Each clone sequence was subjected to an ORF analysis by the GeneJockey software program (Biosoft, Ferguson, MO). If a putative initiation codon followed by an ORF covering at least 70% of the EST could not be identified, we used RACE to gather additional cDNA sequence data.

Reproductively mature D. mojavensis adults of both sexes served as the tissue source for RACE-ready template. mRNA was isolated using the MicroPoly(A)-Pure kit (Ambion, Austin, TX). RACE-ready cDNA was prepared and target molecules were PCR amplified and isolated using the GeneRacer (Invitrogen) protocol, which preferentially selects full-length transcripts for first-strand cDNA synthesis. RACE products derived from such a library should provide high-quality information on the 5′ ends of transcripts. Several criteria were used to identify the set of ORFs ultimately used in molecular evolutionary analysis: (i) size and position of candidate ORFs within an EST, (ii) presence of a predicted signal peptide sequence for putative Acp's (Wolfner  et al. 1997; Swanson  et al. 2001), (iii) tBLASTx homology to genes in public databases (e.g., D. melanogaster genome release 3.1), and (iv) presence/absence of INDEL mutations and/or premature termination codons in polymorphism data from genomic DNA. Only strongly supported ORFs were used in evolutionary analysis.

Quantitative PCR evaluation of ESTs:

Genes targeted for population genetic analyses as accessory gland vs. testis-enriched in expression on the basis of dot blots were subjected to more rigorous quantification of transcript distribution and abundance by real-time quantitative PCR. For the subset of genes in which a related D. melanogaster gene was identified, quantitative PCR was also carried out in D. melanogaster to provide comparisons of expression between lineages. The purpose of this analysis was to assign genes to three expression classes: Acp, testis enriched, and other tissues. A total of 58 and 33 genes were investigated in D. mojavensis and D. melanogaster, respectively.

Tissue dissections consisted of 80 D. mojavensis and 40 D. melanogaster male flies. All flies were reproductively mature and were dissected in RNAlater (Ambion) into three tissue categories: accessory glands, testes, and carcasses without the reproductive tracts. Each collection of dissected tissues was divided equally into two replicate samples for RNA isolation. Likewise, whole, reproductively mature female flies from each species (n = 40) were evenly split into two replicate RNA preps. Total RNA was extracted using Trizol Reagent (Invitrogen), purified through RNeasy (QIAGEN) columns, and treated with DNase according to manufacturer's instructions (QIAGEN). RNAs were then reverse transcribed at a concentration of 20 ng/μl using TaqMan reverse transcription (RT) reagents (Applied Biosystems). These first-strand cDNAs served as the templates for quantitative PCR analysis.

Quantitative PCR was performed using an ABI Prism 7700 sequence detector and SYBR green PCR core reagents (Applied Biosystems). Amplification primers were designed with Primer Express (Applied Biosystems). For every 20-μl PCR reaction, 0.5 μl of first-strand cDNA was used. Quantitative PCR conditions were 94° for 10 min followed by 40 cycles of 94° for 20 sec, 59° for 30 sec, and 72° for 30 sec. A dissociation step was added to the end of the run to ensure that only a single amplicon was produced in each reaction. All primer pairs produced a single product. A total of 13 quantitative PCRs were processed for each gene. Three reactions were run for each of the four tissues: one for each of the two replicate RT reactions as well as a single minus-RT reaction derived by drawing equally from the minus-RT templates of paired replicates. The 13th reaction was a no-template control. We found no evidence of genomic contamination or primer-by-reagent interactions.

Quantitative PCR quantification:

Quantification followed the

methods of Livak and Schmittgen (2001). Quantitative PCR provides an estimate of C_T, the cycle at which the quantity of amplified product exceeds a predetermined threshold. Therefore, more abundant transcripts should yield lower C_T scores. To control for different first-strand cDNA concentrations across templates, as well as run and reagent effects, our ΔC_T scores were calculated by subtracting experimental gene C_T scores from housekeeping gene C_T scores derived from the same tissue and experimental microtiter plate. The housekeeping control for both species was the ribosomal protein gene CG7808, which was identified in the D. mojavensis reproductive tract cDNA library (moj12) and is highly conserved between D. mojavensis and D. melanogaster (96% protein similarity).

Our calculation of

reflects fold change in gene expression of the most abundant tissue template (lowest ΔC_T score) relative to the second most abundant tissue template for any given gene. There were two justifications for this approach. First, we observed several instances in which quantitative PCR product was detected in only two of the four templates. Second, compared to approaches estimating fold differences across all tissues, our approach minimizes fold difference values, thereby providing conservative lower-bound estimates for actual differences between tissue transcriptome profiles. The two replicate

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

scores for each gene were always independently calculated and then averaged for the reported values.

Quantitative PCR statistics:

Replicate

scores for every gene and for each of the four templates can be used to determine the amount of experimental error. A scatter plot of replicate ΔC_T scores for the most abundant tissue of each surveyed gene (n = 91, Figure 1) reveals a high degree of similarity between replicate pairs (R² = 0.979). The slope of this line (m = 0.985) is very close to m = 1, showing that the high repeatability of our assays holds across a large range of expression estimates.

We used our replicate

scores to determine threshold fold differences that are sufficiently disparate to represent significant differences. To approximate a gamma distribution, we calculated ratios of replicate pairs by dividing the higher

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

score by its counterpart and then subtracting one. A total of 91 replicate reaction pairs generated a distribution ranging from 0.0 to 18.24. We then used the x₀ value at which the area under the frequency distribution (0 ≤ x ≤ x₀) is equal to 0.95 to establish a critical threshold for significant differences between

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

scores. For the complete data set,

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

scores >7.84 represent significant differences between tissues (P < 0.05). This is a conservative critical threshold estimate because genes that are highly tissue specific (those with high

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

scores) are susceptible to larger error in terms of relative expression differences. This is a consequence of fold differences being derived by comparing the most abundant tissue (lowest ΔC_T) to the second most abundant tissue. Thus, fold difference for a gene that is highly tissue specific in expression is estimated by comparison to a tissue showing very low transcript abundance. In such cases, experimental error associated with the less abundant tissue expression will affect

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

scores of highly tissue-specific genes. Many of our genes have large

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

scores (see supplementary Table S2; http://www.genetics.org/supplemental), which indicate high tissue specificity. Restricting our statistical analysis to genes with

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

< 50 (n = 28), the critical threshold for significance is reduced to 3.25 (P < 0.05). Further narrowing the analysis to genes with

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

< 15 (n = 24) reduces the critical threshold to 2.10 (P < 0.05).

The different critical values for different subsets of the data support the idea that error variance of relative expression levels is greater for genes with the highest

scores. Therefore, we view the critical threshold of 2.10 as most informative because it is derived from the very data whose relative expression patterns are most in doubt. Even so, we choose a conservative critical threshold of

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}{=}5.0\)

(fivefold difference in relative expression for the most abundant vs. next most abundant tissue) for the purpose of categorizing genes as either Acp's or testis enriched. Though somewhat arbitrary, we note that categorization of genes would not be substantially altered by choosing a more conservative threshold. For example, a critical threshold of 18 would only recategorize three testis-enriched genes as genes showing no strong pattern of tissue enrichment.

D. mojavensis genomic library:

A genomic library was constructed to provide flanking data around gene sequences to help identify regions of homology between D. melanogaster and D. mojavensis (e.g., Wagstaff and Begun 2005; see supplementary material, http://www.genetics.org/supplemental). D. mojavensis genomic DNA was partially digested with Sau3A and size fractionated by electrophoresis through a 0.6% agarose gel. DNA fragments between 9 and 23 kb were selected via gel extraction (QIAGEN), ligated to λ-DASH II/BamHI vector (Stratagene, La Jolla, CA), and packaged using the Lambda DASH II/Gigapack II cloning kit (Stratagene). The resultant library consisted of ∼2.3 × 10⁶ plaque-forming units. Plaques were screened with ³²P-labeled D. mojavensis target DNA. Lambda DNA was purified from selected plaques and D. mojavensis genomic inserts were amplified using T3/T7 vector primers and LA-Taq long PCR polymerase (TaKaRa, Shiga, Japan). The resulting PCR products were sheared by sonication and the fragments were blunt ended using Klenow fragment of DNA polymerase and T4 DNA polymerase. Fragments of 1–2 kb were isolated from a low-melting agarose electrophoresis gel and cloned into the pUC18/SmaI/BAP vector with a Ready-to-Go kit (Amerisham Biosciences, Piscataway, NJ). Sequencing of the phage through ∼7× coverage was performed on an ABI Prism 3700 sequencer. Consensus sequences were assembled using the SeqMan program of the DNASTAR software package (Lasergene, Madison, WI).

Nomenclature:

Unique ESTs were assigned numbers (1–54 for reproductive tract library ESTs, 100–217 for testis library ESTs). Genes from the quantitative PCR analysis showing at least fivefold greater expression (⁠

⁠) in either accessory glands or testes were categorized as Acp's and testis enriched in expression (hereafter referred to as testis-enriched genes), respectively. Prefixes for numbered EST names were added according to these expression patterns, with Acp preceding accessory gland genes and Tes preceding testis-enriched genes. Those genes that did not exceed this threshold (moj9, moj29, moj30, moj32, moj137, and moj152) were given the moj prefix to avoid a connotation of tissue specificity. Four Acp's (Acp5, Acp16, Acp21, and Acp27) are members of recently duplicated gene families (B. J. Wagstaff, unpublished data) and are given an additional -a or -b suffix to differentiate between members. Five genes were named as Acp's (Acp4, Acp15, Acp17, Acp23, and Acp36) on the basis of very strong evidence from our dot blot data rather than from quantitative PCR experiments. The remaining ESTs were given the moj prefix, as no relative expression data were gathered for the associated genes.

Stocks and DNA sequencing:

A total of 15 fly stocks from the Drosophila Species Stock Center (Tucson, AZ) were used for collection of most population genetic data. D. arizonae (15081-1271.00, 15081-1271.04, 15081-1271.05, 15081-1271.08, 15081-1271.12, 15081-1271.13, and 15081-1271.14; various locations, mainland Mexico) and D. mojavensis were represented by seven lines each, while a single D. mulleri stock (15081-1371.00; Lake Travis, TX) provided outgroup data. Of the seven D. mojavensis stocks, four were D. mojavensis baja (15081-1351.03, 15081-1351.09, 15081-1351.12, and 15081-1351.14; various locations, Baja, Mexico) and three were D. mojavensis mojavensis (15081-1352.00, 15081-1352.01, and 15081-1352.02; various locations, southern California). Primers used for amplification of genomic DNA were designed from ESTs or from extended sequences identified by RACE analysis. Expand High-Fidelity polymerase (Roche Molecular Biochemicals) was used for PCR amplification. Single alleles for sequencing were isolated by cloning PCR products into the TOPO vector (Invitrogen) and selecting one bacterial colony for PCR amplification for each allele. Amplified colony-PCR products and their associated sequences were obtained using M13 reverse and T7 primers. A second collection of D. mojavensis mainland and Baja strains (kindly provided by W. Etges, University of Arkansas) was used for additional population sequencing of Acp7. PCR products from the Etges strains were directly sequenced. All sequencing was done on an Applied Biosystems 377 automated sequencer (ABI). Sequences were aligned and edited using the DNASTAR software package (Lasergene). Generally, the small, predicted size of most Acp's resulted in survey data for most codons. Compared to Acp's, testis-enriched genes, on average, provided lower coverage of codons on a per gene basis (see Table 1).

Statistical analysis of aligned sequences:

The DnaSP program (Rozas and Rozas 1999) was used for most of the population genetic analyses. Average levels of polymorphism or divergence for different groups of genes (e.g., Acp vs. testis enriched) refer to means weighted according to sequence length. For genes sampled for multiple alleles, replacement and synonymous divergence represent the average pairwise difference. Fixations for polarized McDonald-Kreitman tests were assigned using parsimony. Only codons with single mutations that could be clearly assigned to either the D. arizonae or D. mojavensis lineage were considered.

Lineage-specific synonymous and replacement divergences were estimated using the free-ratio maximum-likelihood model of the PAML computer program (Yang 1997). For most of these analyses we used one randomly selected allele from each of three species: D. arizonae, D. mojavensis, and D. mulleri. In some cases for which D. mulleri data were unavailable, we used a duplicated gene that predated the D. arizonae/D. mojavensis speciation event (B. J. Wagstaff, unpublished data). We used only duplicated genes showing synonymous divergence that was comparable to or less than the average D. mulleri synonymous divergence (see Table 3). Hypothesis testing was carried out using likelihood-ratio tests (Goldman and Yang 1994; Yang 1998). To determine whether or not K_a significantly exceeds K_s in a particular lineage, the likelihood value for the null hypothesis (K_a = K_s; i.e., the one-ratio model) was also calculated. Twice the log-likelihood difference between the two models is then compared to a χ²-distribution with one d.f. to determine the level of significance.

RESULTS

Content and characterization of D. mojavensis male reproductive tract libraries

The content and basic characteristics of the D. mojavensis male reproductive tract (ESTs 1–54) and testis (ESTs 100–217) libraries are listed in supplementary Table S1 (see http://www.genetics.org/supplemental). Genes with measured accessory gland or testis tissue enrichment are given the Acp and Tes prefixes, respectively. Six moj genes (moj9, moj29, moj30, moj32, moj137, and moj152; see supplementary Table S2, http://www.genetics.org/supplemental) are expressed in multiple tissues. No relative expression analyses were performed on genes corresponding to the remaining moj ESTs (see quantitative PCR section below).

Library content:

Minimal sequencing of the D. mojavensis male reproductive tract library revealed that most of the ESTs corresponded to just a few genes. Preliminary dot blot analysis of an initial set of clones revealed that most ESTs were accessory gland rather than testis derived. Of the first 139 sequenced clones, 35 corresponded to Acp1, 27 to Acp5, and 18 to Acp17. The 139 clones also contained 13 singletons and 10 transcripts represented by 2–9 clones each. The preponderance of Acp's in the reproductive tract library cannot be easily explained by size differences between accessory glands and testes, as D. mojavensis testes appear to be considerably larger than accessory glands (B. J. Wagstaff, personal observation). Thus, per unit of tissue, accessory glands likely produce much more mRNA than the testis. We conclude that the D. mojavensis accessory gland transcriptome has low complexity and high transcript abundance relative to that of the testis transcriptome. To increase the discovery rate of new transcripts, additional clones were screened by multiplexed PCR reactions that included primer pairs specific to Acp1, Acp5, and Acp17. Clones not corresponding to any of these three genes were then sequenced. This multiplex PCR strategy revealed 28 new ESTs from only 66 additional sequencing reactions. In total, 54 unique ESTs were revealed from the reproductive tract library. The average length of all 205 ESTs was 438 bp.

We constructed and screened a D. mojavensis testis cDNA library to increase our sample size of testis-expressed genes. The distribution of replicate ESTs differs dramatically from the original reproductive tract library (supplementary Table S1, http://www.genetics.org/supplemental). The testis library has a much higher complexity than the reproductive tract library, with 105 of 156 clones present as single-copy sequences. Similarly high complexity of a testis cDNA library was previously observed in D. melanogaster (Andrews  et al. 2000), suggesting that this might be a general property of the Drosophila testis transcriptome. In total, 156 sequencing reactions returned an average EST length of 451 bp and produced 118 unique ESTs.

The whole reproductive tract library contains a higher percentage of unique ESTs with potential signal peptide sequences, which is to be expected of a library derived primarily from accessory gland transcripts (Wolfner  et al. 1997; Swanson  et al. 2001). Of library sequences subjected to SignalP analysis, 64% (32/50) of whole reproductive tract-derived unique sequences and 10.3% (7/68) of testis-derived unique sequences contain putative signal sequences (those with hidden Markov model P > 0.75, supplementary Table S1, http://www.genetics.org/supplemental).

Library quality:

Completeness of 5′ cDNA ends was assessed by two methods on a total of 155 ESTs. First, for transcripts represented by more than one clone, we compared the similarity of 5′ ends among clones, with the assumption that the longest clone is likely to include the complete 5′ end of a gene. Second, several transcripts were subjected to 5′ RACE verification. RACE analysis showed that all 20 of the Tes100 ESTs were truncated products, each ∼113 bp shorter than the reference 5′ sequence. Thus, Tes100 clones appear to be outliers in terms of assessment of library quality. Using the multiple-clone method, we estimate 79.7% (63/79) of our ESTs to be complete at the 5′ end. For ESTs compared to a reference 5′ RACE sequence, ∼62.5% (60/96) contain the complete 5′ end. However, removing the Tes100 outliers increases the estimate to 78.9% (60/76), a ratio that is consistent with the multiple-clone estimate. Therefore, our estimates suggest that approximately four-fifths of cDNA clones were complete at the 5′ end.

BLAST analyses vs. D. melanogaster:

Results of D. mojavensis EST BLAST analyses to D. melanogaster databases, including the closest-matching genes and secondary E-scores, are listed in supplementary Table S1 (http://www.genetics.org/supplemental; E < 0.01 was the BLAST score threshold for inclusion). None of the ESTs that failed to match D. melanogaster sequences matched any other NCBI database sequences. Approximately 61% (33/54) and 58% (68/118) of whole reproductive tract and testis library unique ESTs, respectively, showed BLAST similarity to D. melanogaster sequences. However, there were major differences between accessory gland- vs. testis-derived sequences, with Acp's showing a much lower level of conservation between species than testis-enriched genes. Only 33% (8/24) of Acp's generated significant hits, compared to 82% (27/33) for testis-enriched genes. A 2 × 2 contingency table is significantly heterogeneous (P ≪ 0.01). Furthermore, the median E-value of the eight Acp's with E < 0.01 (E = 1e-3, a value too high to reliably indicate orthology) is much greater than the median for testis-enriched genes (E = 2e-21). The six genes that are more ubiquitously expressed on the basis of quantitative PCR data (moj9, moj29, moj30, moj32, moj137, and moj152) had highly significant BLAST matches to D. melanogaster sequences (median E = 5e-42). The remaining moj sequences are similar to the testis-enriched genes, with 55% (59/108) returning E < 0.01 vs. D. melanogaster (median E = 1e-27). This is not surprising, given that most moj sequences are from the testis cDNA library.

Twenty of the 27 D. mojavensis testis-enriched genes that appear to have D. melanogaster homologs have BLAST hits to the D. melanogaster testis EST collection (Andrews  et al. 2000), suggesting that testis expression patterns between species are generally conserved. Our quantitative PCR data from 6 of the remaining 7 genes isolated from the D. mojavensis testis library (see below) suggest that they too show testis-enriched expression in D. melanogaster in spite of their absence from the D. melanogaster testis EST collection, further supporting the notion for a generally conserved Drosophila testis transcriptome.

Certain biochemical functions, including proteases, protease inhibitors, and lipases, appear to be common in melanogaster subgroup Acp's, as inferred from sequence similarity to protein databases (Swanson  et al. 2001). This is in contrast to our results from the 54 unique D. mojavensis reproductive tract ESTs (most of which are likely Acp's), which revealed evidence for two protease inhibitors, Acp36 and Acp48, and a single lipase gene, moj37. None of the predicted 54 proteins contain putative protease domains. The proportion of known D. mojavensis Acp's (2/24) that contain any of these three domains is significantly different from the proportion (21/57) from the Swanson  et al. (2001) set of melanogaster subgroup Acp's (G-test, P = 0.026). This is suggestive of a fundamental, functional divergence in seminal fluid function in the two species, although more work, including direct biochemical assays, would be necessary to put this conclusion on firmer ground.

D. melanogaster-D. mojavensis orthology:

The existence of gene families and shared protein domains can yield small BLAST E-scores, yet obscure inferences regarding orthology between D. melanogaster and D. mojavensis. Alternatively, conserved intron-exon structure is expected for genes of shared ancestry (Meyer and Durbin 2004) but not for unrelated genes that share only a particular protein domain. For example, human-mouse orthologs have the same number of coding exons ∼86% of the time (Mouse  Genome  Sequencing  Consortium 2002). Thus, genes showing conserved intron-exon structure and large E-score differences (e.g., E > 1e-10) between primary and secondary BLAST hits are probably orthologs.

Comparison of genomic sequence from our population genetic data to our EST sequences allowed us to determine intron-exon structure for a subset of D. mojavensis genes (i.e., genes from Table 1). We used this information in concert with comparisons of primary vs. secondary BLAST E-values and protein size, to investigate putative D. melanogaster orthologs for many of our D. mojavensis genes (indicated by an asterisk, supplementary Table S1; http://www.genetics.org/supplemental). For the remaining ESTs we have data only on primary vs. secondary BLAST E-values, many of which are suggestive of orthology.

Acp's:

Of the eight Acp's that show BLAST similarity to D. melanogaster genes (E < 0.01), only Acp36 and CG16713 (supplementary Table S1; http://www.genetics.org/supplemental) are likely orthologs. Both consist of 82 residues and possess a Kunitz domain that covers 59 of those residues. The aligned predicted proteins are 57.3% identical (47/82) and require no gaps. Although Acp36 also returns a significant BLAST hit to another protein with a Kunitz domain (CG16712), its amino acid sequence is more similar to CG16713. Three additional Acp's (Acp1, Acp2, and Acp25) are part of a gene family and are clearly homologous to the Acp53 gene family (Holloway and Begun 2004) in D. melanogaster (supplementary Table S1; http://www.genetics.org/supplemental). However, a protein distance tree clusters the three D. mojavensis genes together, rather than generating the three interspecific pairs expected under one-to-one orthology and homogeneous rates of protein evolution. Thus, although the proteins appear homologous, orthology is uncertain. The remaining Acp's show no compelling evidence for orthology for several reasons, including poor BLAST scores, radically different protein lengths or intron-exon organization, or very different expression patterns between species (described below).

Testis-enriched and moj genes from the population genetics survey:

Most testis-enriched and all six moj genes from the population genetics survey have clear D. melanogaster orthologs on the basis of primary and secondary BLAST E-scores and gene organizations inferred from comparison of cDNA and genomic sequence (supplementary Table S1; http://www.genetics.org/supplemental). However, there are some exceptions. Tes33 and Tes104 are part of an SCP-related gene family and have no obvious orthologs among the many D. melanogaster SCP-related genes. Tes114, Tes120, and Tes123, are also part of gene families that obscure interspecific relationships. Finally, Tes101 and Tes109 are too dissimilar to their D. melanogaster primary BLAST hits (E = 6e-03 and E = 7e-04, respectively) to conclude that they represent orthologous pairs.

Of the 41 genes in our quantitative PCR analyses (see below) that return significant BLAST matches to D. melanogaster sequences, only Tes14 and Tes118 correspond to putative unannotated genes. This supports the observation that the D. melanogaster genome annotation is of high quality (Misra  et al. 2002; Drysdale 2003; Yandell  et al. 2005). Details regarding Tes14, Tes118, and other data bearing on D. mojavensis-D. melanogaster orthology are presented as supplementary material (http://www.genetics.org/supplemental).

Relative quantification of D. mojavensis gene expression:

Supplementary Table S2 (http://www.genetics.org/supplemental) summarizes the expression quantification results for all D. mojavensis genes surveyed, as well as several D. melanogaster genes that are discussed in the next section. Of the 58 total D. mojavensis genes selected for quantitative PCR, 19 are expressed primarily in the accessory glands, 33 are expressed primarily in the testis, and the remaining six (moj9, moj29, moj30, moj32, moj137, and moj152) are more evenly expressed, as indicated by

< 5. The vast majority of the 58 genes appear to be either tissue specific or highly tissue enriched in expression, with 46 of 58 genes being at least 50 times more abundant in one tissue than in any other. All 19 Acp's contain putative signal peptide sequences (supplementary Table S1; http://www.genetics.org/supplemental). Furthermore, the ΔC_T scores indicate that the six most abundantly expressed genes are Acp's. These data, as well as the preponderance of putative accessory gland transcripts in the D. mojavensis reproductive tract library, support the conclusion that Acp's are typically abundantly expressed, secreted peptides (Wolfner 1997).

Figure 2 depicts the relationship between ΔC_T and

scores. The highly significant negative correlation (R = −0.5, P = 0.0002) suggests that genes showing greater degrees of tissue specificity (high

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

scores) tend to have greater transcript abundance (lower ΔC_T). The

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

scores suggest that the 19 most tissue-specific genes are testis rather than accessory gland enriched. Although this could be genuine, we suspect that it is an artifact of trace accessory gland contamination of testis tissue dissections. The transparent and fragile nature of accessory gland tissue should lead to this type of contamination rather than the converse. However, low levels of this one-way contamination should not dramatically affect our conclusions. The fact that several putative Acp's clearly show very large fold differences suggests that this trace contamination is negligible. For example, Acp2 ranks as the most tissue-specific Acp, with transcript abundance in accessory glands estimated as 933 times greater than that in the testis (supplementary Table S2; http://www.genetics.org/supplemental). Conservatively assuming this gene is not transcribed in the testis, this would suggest that there are 933 parts accessory gland material in the accessory gland tissue preparation for every 1 part of contaminating accessory gland material in the testis tissue preparation. Thus, we would not conclude, for example, that Tes101 (⁠

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

= 36,656) is more tissue specific than Acp2 (⁠

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}{=}933\)

⁠). On the other hand, Acp2 is certainly more tissue specific than Acp25 (⁠

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}{=}51\)

⁠) since contamination would affect every Acp gene

\(2^{{-}\mathrm{{\Delta}{\Delta}}C_{\mathrm{T}}}\)

score in a similar manner.

Comparison of D. melanogaster and D. mojavensis expression patterns:

Our quantitative PCR data suggest that putative orthologs of D. mojavensis testis-enriched genes are also testis enriched in D. melanogaster. Nevertheless, the relative amount of testis specificity varies across genes. At the most extreme, D. melanogaster CG3708 is ∼164-fold more testis specific than Tes129. There are also large fold differences between Tes106/CG30334 (97-fold), Tes110/CG15219 (24-fold), and Tes127/CG10090 (53-fold). These comparisons reflect significant differences between D. mojavensis-D. melanogaster expression profiles at these genes (P < 0.05). Several additional testis-enriched genes are borderline significant with fold differences >5. These conclusions are all based on the premise that the housekeeping ribosomal protein gene used as an internal standard has not evolved substantial gene expression differences in D. melanogaster vs. D. mojavensis. Moreover, fold differences can be dramatically different between orthologous pairs solely because of regulatory changes in the secondary tissue and, as such, misrepresent actual differences between species in primary tissues. In this sense, ΔC_T scores are more revealing because they are correlated with absolute expression levels. Five moj and eight Tes orthologous gene pairs have sizable differences between ΔC_T scores (>4), with Tes124-CG14079 at the most extreme (11.35) (supplementary Table S2; http://www.genetics.org/supplemental). Because of the uncertainty associated with housekeeping gene regulation and primer efficiency (see materials  and  methods), strong, individual gene pairwise D. melanogaster vs. D. mojavensis conclusions are not warranted beyond the rank order of tissue enrichment. However, these data suggest that there have been gene regulation changes between lineages. Further discussion of expression differences between individual D. melanogaster-D. mojavensis pairs can be found in the supplementary material (http://www.genetics.org/supplemental).

Genome-wide assays of expression differences between melanogaster subgroup species suggest that male-biased genes show greater interspecific expression differences compared to other genes (Meiklejohn  et al. 2003; Ranz  et al. 2003; Rifkin  et al. 2003). Our data, although consistent with these reports, suggest that despite rapid evolution of male-biased expression, wholesale shifts in tissue specificity are uncommon. A potential caveat is that the apparent interspecific conservation of tissue specificity could be inflated by the fact that we focused on genes coding for more highly conserved proteins. If more highly conserved orthologs are less likely to change tissue specificity, then we have clearly underestimated the frequency of such changes. Comparative genomic analyses of more highly diverged orthologous genes will help address this question (Wagstaff and Begun 2005).

Molecular population genetics analysis

We surveyed a total of 56 genes (19 Acp's, 31 testis enriched, and 6 ubiquitously expressed) for our molecular population genetics analysis (see Table 1). Up to seven lines each of D. arizonae and D. mojavensis were analyzed for several genes. However, many genes are represented by only a single allele each from D. arizonae and D. mojavensis. A D. mulleri allele was sequenced whenever possible as an outgroup. An average of 9.3 alleles and 376 bp were sequenced for each gene surveyed.

Evidence of D. m. baja-D. m. mojavensis population substructure:

Our D. mojavensis data consist of up to four alleles of D. m. baja and three alleles of D. m. mojavensis from various locations of Baja, Mexico and southern California, respectively. Supplementary Table S3 (http://www.genetics.org/supplemental) shows our analysis of population substructure between D. m. baja and D. m. mojavensis. We use the fixation index, F_ST, to estimate genetic differentiation between subspecies. The small size of most surveyed regions and the small number of alleles make inferences from individual genes unreliable. A more accurate view of differentiation can be obtained by examining average F_ST-values, weighted according to sequence length. The average for all genes is 0.150, with the Acp subset of genes slightly higher at 0.168. These results are within the observed range for genetic differentiation between African and non-African D. melanogaster populations (Caracristi and Schlötterer 2003; Baudry  et al. 2004). Acp7 appeared to be something of an outlier with estimated F_ST of 0.864. Therefore, we included additional D. m. baja (n = 5) and D. m. mojavensis (n = 7) Acp7 alleles to the analysis. Our revised estimate showed that differentiation (F_ST = 0.429) at this locus, although at the high end compared to most loci, was not an obvious outlier.

We also investigate genetic differentiation by estimating divergence between subspecies (K_a and K_s) and comparing those values to nucleotide diversity (π) within subspecies (supplementary Table S3; http://www.genetics.org/supplemental). Since both measurements represent the probability that a particular nucleotide site drawn from two individuals is different, they can be directly compared. Again, our analysis shows some evidence of population substructure. Averaged across all genes, K_a (0.006) is higher than both replacement D. m. baja (0.005) and D. m. mojavensis (0.004) nucleotide diversities. However, there are no significant differences between sets of K_a and replacement π measurements (Mann-Whitney U-test, P = 0.77 and P = 0.41 for D. m. baja and D. m. mojavensis, respectively). There is also no evidence for differentiation at synonymous sites, with D. m. baja synonymous π at 0.016, K_s at 0.015, and D. m. mojavensis synonymous π at 0.013.

Given these results, we do not distinguish between D. m. baja and D. m. mojavensis alleles in our population genetics analyses. Although our estimates of polymorphism may be slightly inflated relative to those measured from a single population, our tests of adaptive evolution compare nucleotide substitution patterns at synonymous vs. replacement sites. Under neutrality, population substructure is expected to have little effect on rejecting the null in the direction of adaptive protein evolution.

Levels of synonymous and replacement polymorphism and divergence:

Summary statistics for heterozygosity and divergence for individual genes and for gene categories are presented in Tables 1–3. As suggested by previously published molecular population genetics data from these species (e.g., Begun 1997; Begun and Whitley 2002; Matzkin and Eanes 2003), they are highly variable (Table 1). Average synonymous heterozygosities for D. mojavensis and D. arizonae are 0.0181 and 0.0170, respectively (Table 2). Synonymous heterozygosity for D. mojavensis and D. arizonae is marginally lower for Acp's (0.0156 and 0.0135, respectively) compared to testis-enriched genes (0.0170 and 0.0175, respectively). Synonymous divergence between D. arizonae and D. mojavensis is similar across gene categories as well (Table 2, but see the polarized analysis below for between-species differences). Testis-enriched genes are the most divergent at 0.0682, followed by Acp's at 0.0643 and moj genes at 0.0518. None of the variation statistics are significantly different between gene classes or species by Mann-Whitney U-tests.

Patterns for replacement variation are quite different. First, mean replacement heterozygosity of Acp's in both species is greater than that of testis-enriched or moj genes (Table 2). This is especially striking for Acp vs. testis-enriched genes of D. mojavensis, with Acp's ∼3.7 times more variable than testis-enriched genes in D. mojavensis compared to 1.8 times more variable than testis-enriched genes in D. arizonae. D. mojavensis Acp's have the highest ratio of replacement to synonymous heterozygosity (0.5991), followed by D. arizonae Acp's at 0.4866 (Table 2). This observation is not attributable to population substructure in D. mojavensis, as the ratios of replacement to synonymous Acp nucleotide diversity (π) are 0.6429 and 0.6667 in D. m. baja and D. m. mojavensis, respectively (see supplementary Table S3; http://www.genetics.org/supplemental). The ratios of replacement to synonymous heterozygosity for testis-enriched genes (D. arizonae, 0.2095; D. mojavensis, 0.1476) and moj genes (D. arizonae, 0.1553; D. mojavensis, 0.1308) are much lower in both species. Average Acp replacement divergence between D. arizonae and D. mojavensis is also considerably higher (0.0595) than that observed at testis-enriched (0.0128) or moj genes (0.0060). The ratio of replacement to synonymous divergence for Acp's (0.9257) is 4.9 times greater than the corresponding Tes genes ratio (0.1873). Six genes, all Acp's, have K_a/K_s > 1 (Table 1). Several other pairwise Acp divergence estimates revealed unusually high K_a/K_s values (i.e., >0.5). In contrast, the highest K_a/K_s ratio among nonpolarized Tes and moj genes is 0.8992 for Tes109, with that of most genes being considerably lower (i.e., <0.5). A survey of Acp polymorphism and divergence in D. simulans and D. melanogaster also suggested that these genes evolve unusually quickly at replacement sites relative to other genes (Begun  et al. 2000). However, the relative amount of replacement to synonymous variation at Acp's in D. arizonae and D. mojavensis is much greater than that observed in D. simulans and D. melanogaster. For example, the ratio of replacement to synonymous polymorphism for desert Drosophila Acp's (0.5991 for D. mojavensis, 0.4866 for D. arizonae; Table 2) is about twofold greater than the corresponding ratio in D. simulans (0.2643). The same is true for replacement to synonymous divergence—the K_a/K_s ratio for desert Drosophila (0.9257) is more than twofold greater than the K_a/K_s ratio for D. melanogaster/D. simulans (0.4248). Thus, levels of both protein polymorphism and divergence are considerably greater at Acp's in D. arizonae/D. mojavensis than in D. melanogaster/D. simulans. Although ratios of replacement to silent Acp polymorphism appear to be heterogeneous across melanogaster subgroup species (Begun  et al. 2000; Kern  et al. 2004), we observed no such heterogeneity for D. mojavensis vs. D. arizonae Acp polymorphism (see Table 6; G-test, P = 0.574).

Polarized divergence:

The divergence estimates presented in Tables 1 and 2 result from pairwise comparisons and so provide no insight into evolution along the D. arizonae vs. D. mojavensis lineage. We investigated evolution along these two lineages using both parsimony and likelihood-based approaches. Table 3 shows the results for all genes for which an outgroup sequence was available. As one might expect from previous analyses, the rank order of K_a/K_s ratios is Acp > Tes > moj in each of the three lineages. Eight of nine Acp's have K_a/K_s > 1 in at least one of the three lineages in polarized analyses (Table 3). Tes genes contain just two examples of K_a/K_s > 1, Tes105 along the D. mojavensis lineage and Tes114 along the D. mulleri lineage. In each case, however, K_a/K_s > 1 is largely due to negligible K_s divergence (zero in both cases) rather than unusually rapid protein divergence. K_a/K_s ratios for polarized Acp's vs. Tes genes are highly significantly different (Mann-Whitney U-test, P ≪ 0.01).

The D. mojavensis lineage has a considerably greater average Acp K_a/K_s ratio than either the D. arizonae or D. mulleri lineage. Across all nine Acp's, the K_a/K_s ratio for D. mojavensis (2.0776) is 2.4 times greater than the ratio for D. arizonae (0.8715). Although Acp replacement divergence is higher in D. mojavensis (0.0273) than in D. arizonae (0.0220), the much lower K_s in D. mojavensis vs. D. arizonae Acp's makes a significant contribution to the higher D. mojavensis Acp K_a/K_s ratio. One possible reason for the low D. mojavensis K_s relative to the D. arizonae K_s could be different intensities of selection for codon bias between lineages. However, our estimates of effective number of codons (ENC) (Wright 1990) show no major differences between lineages. The average ENCs for D. mojavensis Acp's and testis-enriched genes, weighted according to size, are 51.8 and 50.8, respectively. The corresponding values for D. arizonae are 50.7 and 51.6, respectively. Thus, codon bias of D. mojavensis Acp's is actually slightly lower than that of D. arizonae Acp's, contrary to expectations if stronger selection at synonymous sites in D. mojavensis were contributing to the lower D. mojavensis K_s values.

Unfortunately, we have D. mulleri data from only five Acp's (Table 3). This limits our ability to directly compare Acp K_a/K_s across the three lineages in a comparable set of analyses. For these five genes the K_a/Ks average ratio is similar for D. arizonae and D. mulleri (0.8273 and 0.8484, respectively), while the D. mojavensis K_a/K_s ratio (1.7163) is roughly twofold greater. Note that the D. mulleri data are potentially biased because genes that are evolving more quickly would tend to be underrepresented as a result of PCR failure using primers designed from D. mojavensis sequence.

Two Acp's, Acp7 and Acp16a, have K_a/K_s significantly >1 in the D. mojavensis lineage, while neither gene is significant in the D. arizonae lineage. The significant K_a/K_s for D. mojavensis Acp7 reflects a contribution from low synonymous divergence (0.0000), as replacement divergence is similar in D. mojavensis (0.0275) to the Acp mean (0.0273) for the D. mojavensis lineage (Table 3). On the other hand, the high K_a/K_s ratio for D. mojavensis Acp16a is primarily attributable to the atypically high replacement divergence (0.1538) relative to the lineage mean (0.0273). D. mulleri provides a solitary example, Acp7, of K_a/K_s significantly >1 (P < 0.05).

Joint analysis of polymorphism and divergence:

The neutral theory of molecular evolution predicts that the ratio of replacement to synonymous substitutions should be similar to the ratio of replacement to synonymous polymorphisms (Kimura 1983). The McDonald-Kreitman (MK) test uses a 2 × 2 contingency table to detect differences in these ratios (McDonald and Kreitman 1991). Table 4 shows the polymorphism and fixation data for individual genes at synonymous and replacement sites. For cases in which an outgroup sequence was available (outgroups identical to those in Table 3), fixed differences between D. arizonae and D. mojavensis were polarized using parsimony. None of the 54 tests are significant after Bonferroni correction of critical values. The small sizes and large number of genes motivate analysis of pooled data (Table 5). The 2 × 2 table for Acp's is significantly heterogeneous in a direction consistent with adaptive protein evolution and remains marginally significant if Acp25 (the single Acp with P < 0.05) is removed from the analysis. Another individual gene that warrants mention is Acp48. With a total of 60 mutations to contribute to the 2 × 2 contingency table, one might speculate that it has a major effect on the overall conclusion. However, removing the Acp48 data increases the significance of the heterogeneity of the remaining Acp's. Overall, the analysis of pooled polymorphic and fixed mutations supports the notion that directional selection plays a role in accessory gland protein divergence. Data from testis-enriched and moj genes show no significant deviations from neutral expectation in 2 × 2 contingency tables.

Further evidence for different evolutionary processes among gene classes can be found in the ratios of replacement fixations to polymorphisms (Tables 4 and 5). While a total of seven Acp's have more replacement fixations than polymorphisms, no Tes or moj genes do, with the exception of Tes112, which has no replacement polymorphisms and just a single fixation. The ratio of fixed to polymorphic replacement mutations for Acp's (139:115) is highly significantly different from the ratio for testis-enriched genes (15:60; G-test, P ≪ 0.01), a result that cannot be explained by different neutral mutation rates for the two protein classes. The moj genes ratio (0:16) is more testis-like, although with so few data, strong conclusions are unwarranted.

Polarized McDonald-Kreitman analyses:

Investigation of polarized fixations provides more insight into evolutionary process in the D. arizonae vs. D. mojavensis lineages, although at a cost of reduced number of loci and substitutions included in the analysis (numbers of polarized vs. nonpolarized individual gene tests are 9:15, 13:17, and 2:2 for Acp's, Tes genes, and moj genes, respectively). The data for different gene classes, polarized using parsimony, are presented in Table 6. D. mojavensis Acp's show a highly significant (P = 0.004) deviation from neutral expectations. It is formally possible that the D. mojavensis data could be explained by too few silent fixations. However, such an explanation would require a change in silent neutral mutation at precisely the correct moment in time. Moreover, since the ratio of silent to replacement substitutions is similar in the two lineages in other gene categories, this explanation would require a bizarre perturbation of silent neutral mutation rate only in Acp's, which seems highly improbable. Thus, the D. mojavensis Acp data are more easily interpreted as a large excess of replacement fixations. Interestingly, however, the D. arizonae Acp data are not significantly heterogeneous (P = 0.181). The lineage differences in polarized MK tests, which are consistent with the greater K_a/K_s ratio in D. mojavensis Acp's noted earlier, support the idea that directional selection has greater effects on Acp divergence in D. mojavensis than in D. arizonae. Note that the number of fixed replacement vs. synonymous mutations (24:2) in D. mojavensis corresponds to a K_a/K_s ratio for fixed sites of ∼4 (assuming a ratio of replacement to silent sites of ∼3:1), providing additional support for the interpretation that the 2 × 2 table for D. mojavensis Acp's can plausibly be explained only by adaptive protein evolution. Polarized data from moj genes in both lineages and testis-enriched genes in D. mojavensis are not significantly heterogeneous. Data from D. arizonae testis-enriched genes are marginally significant (Fisher's exact test, P = 0.056; G-test, P = 0.026), but not in the direction of excess replacement fixations. Additional population genetic data will be required to investigate this pattern.

DISCUSSION

Population genetic investigation of accessory gland protein genes has previously focused on D. melanogaster and D. simulans (Aguadé 1997, 1998, 1999; Tsaur and Wu 1997; Tsaur  et al. 1998; Begun  et al. 2000; Swanson  et al. 2001; Kern  et al. 2004). Our study of Acp's and testis-enriched genes of desert Drosophila from the repleta group was motivated by our interest in understanding whether the highly diverged mating system of these flies (relative to D. melanogaster and D. simulans) is associated with different population genetic patterns and mechanisms for male reproduction-related genes. This question may be especially germane to the issue of Acp's (rather than testis-enriched genes).

Desert Drosophila from the repleta group remate much more frequently than do D. melanogaster or D. simulans, opening up the possibility for stronger or fundamentally different selection on male-male and male-female interactions in the repleta group. Previous results from within- and between-species matings of desert Drosophila (Patterson and Stone 1952; Knowles and Markow 2001; Pitnick  et al. 2003) support the idea of rapid evolution of ejaculate-female interactions. The fact that D. mojavensis males make detectable postmating donations to females whereas D. melanogaster and D. simulans do not (Markow and Ankney 1984; Pitnick  et al. 1997) is another interesting biological difference. If Acp's are major players in postcopulatory male-male and male-female interactions (Wolfner 1997, 2002; Chapman 2001), we might expect to observe different functions and patterns of evolution in desert Drosophila Acp's compared to melanogaster subgroup Acp's.

Our data do not directly address functional divergence of D. mojavensis/D. arizonae vs. D. melanogaster/D. simulans seminal fluid. However, our BLAST results to protein databases for D. mojavensis vs. D. melanogaster/D. simulans Acp's are suggestive of divergent functional biology (e.g., Wagstaff and Begun 2005), with D. mojavensis proteins enriched for unknown functions and depauperate of lipases, proteases, and protease inhibitors compared to those of D. melanogaster. Additional support for this inference and its possible connections to mating system variation await future investigation.

The population genetics of desert Drosophila Acp's showed some similarities and several important differences with respect to D. melanogaster/D. simulans. D. melanogaster and D. simulans Acp's are highly polymorphic and divergent at replacement sites compared to “typical” genes in these two species (Begun  et al. 2000; Swanson  et al. 2001). Acp's from D. arizonae and D. mojavensis showed a similar pattern in that they were much more polymorphic and divergent at replacement sites, at least compared to the non-Acp genes (mostly testis-enriched genes) surveyed here. However, D. arizonae/D. mojavensis Acp's are proportionally much more polymorphic and divergent in terms of protein variation compared to D. melanogaster/D. simulans Acp's (Table 2). One interpretation is that Acp's tend to be under less functional constraint in desert Drosophila compared to melanogaster subgroup flies. Alternatively, Acp's could be under stronger directional selection in desert Drosophila.

Two types of results support the idea that Acp's experience directional selection in desert Drosophila. First, the K_a/K_s ratio is significantly >1 for two of nine D. mojavensis Acp's. Given the small number of bases surveyed per gene and the fact that the K_a/K_s test is an extremely conservative test for directional selection, observing two of nine genes as individually significant is remarkable. The mean K_a/K_s for D. mojavensis Acp's is 2.078, an extremely high value for any class of genes. Second, the MK tests provide strong evidence for adaptive protein evolution in Acp's, but not in other genes.

Interestingly, Acp data strongly deviate from neutral expectations in D. mojavensis, but not in D. arizonae. Moreover, Table 4 suggests that the highly significant result from the pooled polymorphic and fixed mutations presented in Table 6 is attributable to a consistent excess of replacement fixations across most D. mojavensis Acp's rather than to unusual observations from one or two genes. In fact, almost all D. mojavensis Acp substitutions are amino acid changes. Note that polarized analyses of polymorphic and fixed, synonymous and replacement variation have not been carried out for the D. melanogaster/D. simulans comparison, as outgroup data are generally lacking. In this respect, the population genetic inferences for desert Drosophila are more incisive than those for D. melanogaster and D. simulans. These results support the notion that lineage differences in sexual selection may have detectable effects on patterns of protein evolution (Tsaur  et al. 2001; Dorus  et al. 2004).

Given their close evolutionary relationship and similar mating systems, the inference of directional selection on D. mojavensis Acp's and the lack of such an inference for D. arizonae are interesting. A notable distinction between mating systems is that the D. mojavensis ejaculate donation to female somatic tissues is three- to fourfold higher than that in D. arizonae, representing a far greater absolute difference than that observed for other sister species pairs from a large phylogenetic survey (Pitnick  et al. 1997). This suggests the possibility that this difference should be a focus of our attempts to understand effects of mating system variation on protein variation. Perhaps large somatic donations are correlated with more or stronger Acp-mediated postcopulatory male-female interactions. An intriguing possibility is that the somatic donation from the male is associated with mechanisms that provide males with direct access to the female soma, thereby allowing more direct manipulation of female physiology. In this sense, donation to the female soma could be a Trojan horse that exposes females to exploitation by males, thereby driving male-female conflict and associated Acp divergence. Data from other species pairs with differences in ejaculate donation will shed light on the role of variation in male somatic donations to females in Acp evolution.

An alternative explanation of the differences between D. arizonae and D. mojavensis Acp protein evolution is that our sampling of Acp loci has compromised our ability to make an unbiased comparison between lineages. Because our Acp's were isolated from a D. mojavensis accessory gland cDNA library, we are biased toward isolating genes that are more abundantly expressed in D. mojavensis than in D. arizonae. Therefore, a possible explanation for the differential importance of adaptive protein evolution in D. arizonae vs. D. mojavensis is that more abundantly expressed Acp's are under stronger directional selection. This possibility is easily addressed through additional quantitative analysis (for both expression and population genetics) of larger numbers of Acp's in both species.

There has been much speculation regarding the potential importance of adaptive protein evolution for male-reproduction-related genes. However, the data presented here are the first molecular population genetic analysis of a sample of Drosophila genes expressed primarily in testes. Our results show that in D. arizonae/D. mojavensis, testis-enriched genes evolve much more slowly than Acp's and show no evidence of adaptive protein divergence. Why might Acp's experience more directional selection than testis-enriched genes? Spermatogenesis requires several genes (Fuller 1993; Poccia 1994; Eddy 1998), many of which are unlikely to function directly in male-male and male-female postcopulatory interactions. This is in contrast to Acp's, which are more likely to regulate postcopulatory male-male and male-female interactions (Wolfner 1997, 2002; Chapman 2001; Birkhead and Pizzari 2002). Our contrasting population genetic data for Acp's vs. testis-enriched genes support the idea that proteins controlling postcopulatory, prefertilization phenotypes are more likely to be under directional selection compared to proteins controlling sperm phenotypes per se. However, we predict that proteins controlling sperm phenotypes directly involved in male-male or male-female interactions will show evolutionary patterns similar to those observed at Acp's. Our results suggest, not surprisingly, that the functional categorization of genes as male reproduction related or male biased (e.g., Zhang  et al. 2004) obscures a great deal of heterogeneity regarding mechanisms of evolution. More nuanced treatments of male reproduction-related genes with respect to expression and other aspects of biological annotation will likely add great additional insight into the factors explaining variance of protein evolution for such genes (e.g., Good and Nachman 2005).

Footnotes

Footnotes

Acknowledgement

We thank A. Long and two anonymous reviewers for useful comments. This work was supported by a National Institutes of Health grant (GM55298 to D.J.B.), a National Science Foundation grant (DEB-0327049 to D.J.B.), and a National Science Foundation doctoral dissertation improvement grant.

References