O R I G I NA L A RT I C L E doi:10.1111/evo.12352
INCIPIENT SPECIATION OF SEA STAR POPULATIONS BY ADAPTIVE GAMETE RECOGNITION COEVOLUTION Michael W. Hart,1,2 Jennifer M. Sunday,1 Iva Popovic,1 Kevin J. Learning,1 and Christine M. Konrad1 1
Department of Biological Sciences, Simon Fraser University, Burnaby, British Columbia, Canada 2
E-mail: [email protected]
Received September 3, 2013 Accepted December 16, 2013 Reproductive isolation—the key event in speciation—can evolve when sexual conflict causes selection favoring different combinations of male and female adaptations in different populations. Likely targets of such selection include genes that encode proteins on the surfaces of sperm and eggs, but no previous study has demonstrated intraspecific coevolution of interacting gamete recognition genes under selection. Here, we show that selection drives coevolution between an egg receptor for sperm (OBi1) and a sperm acrosomal protein (bindin) in diverging populations of a sea star (Patiria miniata). We found positive selection on OBi1 in an exon encoding part of its predicted substrate-binding protein domain, the ligand for which is found in bindin. Gene flow was zero for the parts of bindin and OBi1 in which selection for high rates of amino acid substitution was detected; higher gene flow for other parts of the genome indicated selection against immigrant alleles at bindin and OBi1. Populations differed in allele frequencies at two key positively selected sites (one in each gene), and differences at those sites predicted fertilization rate variation among male–female pairs. These patterns suggest adaptively evolving loci that influence reproductive isolation between populations. KEY WORDS:
bindin, coevolution, fertilization, sexual conflict, sexual selection.
Understanding the targets of selection for different locally adapted traits among populations of a single species, and the causes of reproductive isolation between them, is an important outstanding goal of evolutionary genetics (Coyne and Orr 2004; Schluter 2009; Ellegren et al. 2012; Feder et al. 2012; Nosil and Feder 2012). Sexual selection theory predicts that such targets of selection should include genes encoding molecules that mediate conflict between the sexes over interactions with different optimal values for males versus females, such as copulation rates between mates or encounter rates between gametes (Rice 1989; Gavrilets 2000; Levitan and Ferrell 2006; Palumbi 2009). Males often compete to mate with females and fertilize eggs, but entry of multiple sperm in a single egg (polyspermy) can disrupt mitosis and kill the zygote. Because sperm are abundant, small, and individually cheap compared to eggs, the costs of frequent mating between individuals or high fertilization rates between gametes
C
1
2014 The Author(s). Evolution
are typically lower for males and higher for females, and the consequences differ between the sexes: selection favors sperm that fertilize eggs first, but favors eggs that are fertilized once (Frank 2000; Rice 2000; Levitan and Ferrell 2006). This conflict between the reproductive interests of the two sexes can lead to an arms race between male adaptations (e.g., for maximizing sperm encounters with eggs) and female countermeasures (e.g., for repelling supernumerary sperm and avoiding polyspermy), and such a conflict may be resolved in different ways within different populations (Holland and Rice 1998; Frank 2000; Gavrilets 2000). This mechanism of population divergence in reproductive traits is of considerable interest because it can account for the rapid evolution of reproductive isolation among populations in similar habitats or environments (Mani and Clarke 1990; Schluter 2009; Nosil and Flaxman 2011), and because it can arise through conflictual interactions at the level of interacting
M . W. H A RT E T A L .
genes and molecules, such as sperm ligands and their egg coat receptors (Levitan and Ferrell 2006; Palumbi 2009; Vacquier and Swanson 2011; Vacquier 2012). Such a mechanism for speciation has strong theoretical (Gavrilets 2000; Gavrilets and Hayashi 2005) and experimental support (Rice 1989; Holland and Rice 1998; Civetta and Clark 2000; Stewart et al. 2005), but there is only limited evidence for this mechanism in the formation of incipient species in nature (Nosil and Flaxman 2011). Coyne and Orr (2004) suggested that additional evidence might come from analyzing both molecular divergence and sperm–egg compatibility among conspecific populations that have been geographically isolated for a long time, and have ample opportunity for coevolution of male and female traits under local selection. We followed this line of study in populations of a sea star (Patiria miniata) from northern and southern British Columbia that are separated by a relatively old (>200,000 years) phylogeographic break in the northeastern Pacific ocean (Keever et al. 2009; McGovern et al. 2010; Fig. 1). We recently reported significant population differences caused by population-specific responses to selection on the gene encoding the sperm acrosomal protein bindin (Sunday and Hart 2013). Here, we document molecular and functional coevolution between bindin and the egg surface receptor OBi1 (Hart and Foster 2013) in the same populations studied previously, including positive selection for amino acid differences between populations, low gene flow (and selection against immigrant alleles) for bindin and OBi1 compared to other loci, and correlations between bindin and OBi1 genotypes and fertilization rates. The congruence among these structural, functional, and genetic patterns strongly suggests that coevolution of this pair of male- and female-expressed gamete recognition genes driven by sexual conflict within sea star populations could serve as the basis for incipient species formation among diverging populations.
Materials and Methods BACKGROUND AND STUDY SYSTEM
Our analysis of population divergence in gamete recognition molecules and sperm–egg compatibility in P. miniata builds on previous studies of bindin evolution and fertilization ecology in sea urchins (Palumbi 2009; Lessios 2011), and biochemical or molecular genetic analyses of other gamete recognition proteins in sea urchins and sea stars (Hirohashi et al. 2008; Vacquier 2012). In sea urchins, bindin and its egg coat receptors mediate sperm binding to eggs, fertilization success within species, and reproductive isolation between species via molecular interactions between the bindin mass on the surface of the activated sperm head and glycosylated proteins in the egg extracellular matrix (Geyer
2
EVOLUTION 2014
and Palumbi 2003; Levitan and Ferrell 2006; Palumbi 2009; Vacquier 2012). In comparisons among some sea urchin populations and species, bindin variation shows high rates of nonsynonymous amino acid replacement substitutions (dN) relative to synonymous silent nucleotide substitutions (dS), with different codons and lineages under positive selection (ω = dN/dS > 1; e.g., Geyer and Palumbi 2003). These patterns of positive selection have been interpreted mainly in terms of the expected effects of sperm competition and sexual conflict, and variation in those effects among species with different mating system traits or among habitats with different ecological conditions for fertilization (Levitan 2008, 2012; Palumbi 2009; Lessios 2011; Vacquier and Swanson 2011). The key evidence (reviewed in the studies cited above) that implicates sexual conflict in bindin evolution and speciation in sea urchins (and analogous systems of fertilization genes in other organisms) comes from (1) experimental studies of fertilization rates under conditions of sperm competition and high polyspermy risk, and (2) codon model analyses of positive selection in the molecular coevolution of bindin and its egg coat receptor. Unfortunately, there are few case studies of coevolution between interacting male and female fertilization genes (Vacquier and Swanson 2011). Excellent examples come from studies of gastropod lysin and its egg receptor (Clark et al. 2009; Hellberg et al. 2012; Aagaard et al. 2013), but research in this area has emphasized comparisons between phenotypically divergent species at a relatively late stage in the speciation process, and does not directly address the early origins of population divergence leading to incipient species. In sea urchins, the interpretation of both experimental studies of fertilization rate and codon model analyses of selection has been complicated by obstacles to the direct study of variation and molecular evolution of the bindin receptor in the egg coat (Vacquier 2012). The sea urchin egg bindin receptor gene EBR1 (Kamei and Glabe 2003) has a large, complex, repetitive coding sequence structure (with many possible targets for selection distributed among >50 exons) that has been difficult to analyze using population genetic methods. Other genes encode proteins that interact with bindin, and might coevolve with bindin under selection, but have been more intensively studied by cell biologists than by evolutionary geneticists (Foltz et al. 1993; Foltz 1994; Hirohashi et al. 2008; Vacquier 2012). As a result (1) experimental studies of fertilization success among mates with different gamete recognition genotypes have so far used female bindin genotype as a proxy for the egg bindin receptor (Palumbi 1999; Levitan and Ferrell 2006); and (2) analyses of egg bindin receptor evolution have been limited to comparisons among full-length EBR1 coding sequences for a small number of species (Kamei and Glabe 2003) or codon model analyses of a small part of the EBR1 coding sequence for larger samples of species and populations (Pujolar and Pogson 2011).
G A M E T E R E C O G N I T I O N C O E VO L U T I O N I N S E A S TA R S
Figure 1.
Patiria miniata populations differ at positively selected codons in the OBi1 substrate-binding site for bindin. (A) Maximum
likelihood phylogeny of the single positively selected partition from an alignment of 12 OBi1 alleles obtained from transcriptomes (RNASeq data) of northern (cyan) and southern (green) females. The gene tree is drawn using midpoint rooting and HKY genetic distances (scale bar shows 0.008 substitutions per nucleotide site). Labels indicate individual (1–3), population (N, S), and allele (a, b). Positive selection (ω >> 1) was inferred for two codons along the foreground branch shown as the red line. Populations were fixed for Ser-Thr or Ile-Ser differences at codons 658 and 659. (B) Locations of study sites (asterisks) in diverging northern (Haida Gwaii, cyan) and southern (Vancouver Island, green) populations in the northeastern Pacific Ocean. (C) Structural model prediction for OBi1 inferred using SWISSMODEL and rendered in Chimera. A southern allele is shown (northern and southern alleles had nearly identical structures). Backbone structure is rendered as ribbons. Two binding sites identified by SMART are shown in blue: the substrate-binding domain (ribbons), and the ATP-binding domain (spheres). Positively selected codons 658 and 659 are shown in red, between two of the alpha helices predicted to form a lid over the substrate-binding groove.
To extend this line of study to genes encoding interacting protein products from diverging conspecific populations, we first characterized population differences in sperm–egg compatibility. We designed our fertilization experiment specifically to test for local adaptation in fertilization success within a southern and a northern P. miniata population. Because we did not know the genotypes of individual sea stars before running the fertilization experiment, the males and females in this experiment included an unbiased sample of bindin and OBi1 (and other) variation. We then characterized population divergence in bindin (Pati˜no et al. 2009; Sunday and Hart 2013) using the same samples of males and females from the fertilization experiments. After discovering strong population divergence under positive selection at specific bindin sites, we then used RNA-Seq methods (Hart 2013; Hart and Foster 2013) to characterize coding sequence variation in P. miniata homologs of egg coat genes that encode proteins known or likely to interact with sperm (and specifically with bindin) on
the basis of previous studies in sea urchins. We then reanalyzed the fertilization data to test for an interactive effect of variation at positively selected sites in two of these genes (bindin and OBi1) on fertilization success. In the sections below, we describe the population genetic methods and results first, followed by the population differences in fertilization rates and the effects of individual male and female genotypes on sperm–egg compatibility, although this was not the same order in which we analyzed those different types of data. CODON MODELS OF POSITIVE SELECTION
We analyzed coding sequences of OBi1 and five other genes that encode proteins expressed in the egg extracellular matrix, including EBR1, three members of the ARIS gene family that induce the acrosome reaction in sperm, and Rendezvin. We compared those patterns to five other genes that are expressed in eggs or ovaries but are not implicated in gamete recognition, including Dual
EVOLUTION 2014
3
M . W. H A RT E T A L .
oxidase 1, Histone deacetylase 1, Mothers against decapentaplegic 1, Vasa, and Vitellogenin (Appendix S1). The data come from six ovary transcriptomes (three from northern females, three from southern females). All methods for generating, assembling, and analyzing those RNA-Seq data were described previously and are given in full detail elsewhere (Hart 2013; Hart and Foster 2013). Illumina RNA-Seq reads are available from the NCBI Sequence Read Archive (BioProject PRJNA175319); contig assemblies of the transcriptome data can be downloaded from the Dryad digital repository. Alignments (Dataset S1) for each egg-expressed gene were analyzed for evidence of recombination (Kosakovsky Pond et al. 2006; Delport et al. 2010) and divided into two or more partitions corresponding to nonrecombining blocks of codons (see Hart 2013). Input data for each subsequent analysis was an alignment of 12 haplotypes (two gene copies per individual, three individuals per population) for one data partition. Individual partitions for each gene were analyzed using codon models: we used the MEME method (Murell et al. 2012) to identify codons under positive selection (on some unspecified branches in the gene tree), and the branch-site random effects likelihood model (BSR; Pond et al. 2011) to identify lineages of alleles under positive selection (at some unspecified codons in the sequence alignment). Both methods use genetic algorithms that test hypotheses about among-site and among-lineage variation in ω, while avoiding false positives associated with the older branch sites model of positive selection (Yang 2007). For one data partition in which we found congruent evidence of positive selection in both the MEME and BSR models, we used the branch sites model to test the specific hypothesis of divergence between northern and southern populations. We compared the results from these analyses of several egg-expressed genes from ovary transcriptomes (Appendix S1, Table S1) to our recently published analyses of bindin (Sunday and Hart 2013), which were obtained from a larger sample of individuals from the same populations by polymerase chain reaction (PCR) and Sanger sequencing of cloned amplicons. ISOLATION-WITH-MIGRATION MODELS OF GENE FLOW
For genes that were inferred to be under positive selection (in the codon model analyses described above), we estimated rates of gene flow between the northern and southern P. miniata populations, and compared those results to gene flow estimated for loci that were found not to be under positive selection. We used these analyses to test the prediction that positively selected loci would have lower rates of gene flow between populations that reflected selection against immigrant alleles, compared to loci without evidence of population differences caused by selection. To obtain the necessary power for this analysis, we grouped loci into four multilocus datasets based on our codon model results (Table S1):
4
EVOLUTION 2014
one partition of bindin and one partition of OBi1 in which positive selection for population differences had been detected (group A), all other partitions from these same genes for which we found no evidence of positive selection (group B), all partitions of egg coat genes with no evidence of positive selection (group C), and all partitions of five other egg-expressed genes (group D). We used isolation-with-migration models of population demographic history to characterize patterns of gene flow between northern and southern populations using each of those four datasets A through D. We used data from our analyses of eggexpressed genes and our previously published bindin data. We treated partitions as loci and carried out multilocus analyses in IMa (Hey and Nielsen 2007; e.g., see Puritz et al. 2012; Keever et al. 2013). We used all partitions of all genes except Rendezvin (Table S1): for this gene, we were unable to infer all recombination sites simultaneously in one analysis; instead, we divided the Rendezvin alignment into three subalignments (with similar numbers of repetitive protein-coding domains), analyzed each of those subalignments for recombination, and included the longest partition from each as three loci in the IMa dataset. We assumed a generation time (average age of breeding adults) of five years for large, long-lived sea stars. This approach is similar to the IMa-based method of Sousa et al. (2013) that assigns loci to groups with or without evidence of selection against gene flow based on amonglocus differences in isolation-with-migration model parameters, but has the advantage of prior knowledge of among-locus variation in selection (based on the results of codon model analyses). Input files (including sequence alignments) and IMa output are available as Datasets S1 and S2. To standardize the gene flow estimates from different parts of the genome for comparison, we included an additional locus with a known mutation rate in each of the four IMa datasets. We used population samples of mitochondrial DNA haplotypes from the same Sandspit and Bamfield populations analyzed in our previous phylogeographic study of P. miniata populations (Keever et al. 2009). We used the same mtDNA mutation rate calibration as in our previous studies: the lowest observed divergence in mtDNA between geminate pairs of sea urchin species (4.57% for Diadema species) assuming a three million year divergence time (Lessios et al. 2001; Hickerson et al. 2006; McGovern et al. 2010; Puritz et al. 2012). Including a locus with an empirical mutation rate calibration such as this in the IMa analysis allowed us to convert the IMa model parameters estimated from each dataset (all in units scaled by the mutation rate specific to each set of data) to population parameters in the same demographic units (individuals, years) shared among all four datasets. A corollary of this approach is that, because the same single mtDNA coalescent pattern contributed to each set of multilocus population parameter estimates for the four groups of loci (Table S1; Dataset S1), significant differences among the four results are probably
G A M E T E R E C O G N I T I O N C O E VO L U T I O N I N S E A S TA R S
conservative estimates of real differences in demographic history caused by differences among those four parts of the P. minata genome. ADDITIONAL MOLECULAR GENETIC METHODS
Methods for characterizing predicted protein domains of coding sequences from transcriptomes, including the substrate- and ATP-binding domain of the hsp70-like bindin receptor, are given elsewhere (Hart 2013; Hart and Foster 2013). We used the Limbo method (Van Durme et al. 2009) to find domains in a full-length P. miniata bindin coding sequence (GenBank FJ439659.1) that are predicted to be substrates for the OBi1 substrate-binding site. The OBi1 structural model prediction based on the SWISS-MODEL method is given in Dataset S3. To obtain larger sample sizes of allele frequencies for one part of the OBi1 coding sequence that was under positive selection in our analyses of small samples of haplotypes from RNA-Seq data (and implicated in prezygotic reproductive isolation), we amplified and sequenced short (100–200 bp) portions of the OBi1 exon containing the positively selected codon 659. That PCR product was amplified from genomic DNA of all females and some males used in our fertilization experiments (see below). Custom primers (5 -tgtaaaacgacggccagtGTTGAACAAACACAGGGTCAG-3 , including the lower-case sequencing primer at the 5 tail; and 5 -GTTTTCYGTCTTGTTCAGGAG-3 ) were designed from RNA-Seq coding sequences aligned to an OBi1 genomic sequence from the draft P. miniata genome assembly (GenBank AKZP01045554.1). Amplicons varied in length due to an indel in the adjacent intron. We purified and direct-sequenced each amplicon using the M13 universal primer, and scored the three-nucleotide genotype of each individual at codon 659 directly from the chromatograms. FERTILIZATION EXPERIMENTS
We conducted full-factorial fertilization crosses in sperm-limiting noncompetitive conditions. Adult sea stars were collected at low tide from a northern population at Sandspit, Haida Gwaii, and from a southern population at Bamfield, Vancouver Island (see Fig. 1) and held at Simon Fraser University in approximately 1 L plastic containers in large recirculating UV-sterilized seawater tanks under constant low light conditions at 11°C to 12°C and fed regularly with mussels. We acclimated all adults to these laboratory conditions for at least two weeks before using them in fertilization experiments. We conducted 10 independent fertilization crosses on separate days. Each cross included a male and a female individual from each of the two populations, crossed reciprocally. Each individual was removed from seawater, brought to room temperature, and given an intracoelomic injection of 1 to 3 mL 100 mM 1-methyl-adenine, which induces maturation of oocytes (in females) and spawning (in both sexes). Eggs
were collected and washed in filtered (0.45 μm) seawater, then allowed to settle in a glass beaker at approximately 12°C. Concentrated sperm were collected directly from the gonopore and held on ice. Concentrated sperm were diluted 100-fold in filtered seawater. We fixed a subsample of this stock in 5% formalin in seawater, and measured absorbance using a Nanodrop 2000c spectrophotometer (Thermoscientific) to estimate cell concentration (absorbance at 600 vs. 280 nm as recommended by the manufacturer) relative to a predetermined standard curve based on hemocytometer counts of fixed sperm in a dilution series. The experimental sperm stock was then independently diluted to 104 , 105 , and 106 mL−1 in 20 mL of filtered seawater. Aliquots of sperm dilutions used in experimental treatments were also fixed in 5% formalin and counted on a hemocytometer after each experiment to confirm sperm concentrations. Four milliliters of each sperm treatment was immediately dispensed into two 5 mm sterile plastic Petri dishes. One drop of settled concentrated eggs from each female (approximately 0.1 mL) was added to each dish containing sperm (and to a filtered seawater control), and the inseminated eggs were cultured at approximately 12°C. For each dish, we scored fertilization success for the first 100 eggs observed under a dissecting microscope by (1) presence or absence of fertilization envelopes at 30 min postinsemination, (2) presence or absence of cleavage after 2 h, and (3) development to the blastula stage after 24 h. The three measures of fertilization success were highly correlated, and we report proportion of blastula embryos after 24 h because this measure integrates the effects of both successful fertilization encounters (between an egg and single sperm) and unsuccessful encounters (between an egg and multiple sperm leading to polyspermy and embryo death before blastula formation). We tested for localized coevolution of sperm–egg compatibility using a generalized mixed effects linear model of fertilization rate with a binomial error distribution for proportion data (Crawley 2007); sperm concentration, population source of males, population source of females, and the interaction of male and female population source were fixed effects. With this approach, a positive coefficient for the interaction term would indicate a disproportionately greater rate of fertilization success between males and females from the same population. To account for the effect of crossing each male and female to two different mates (one from the same population source, and one from the other population), we included both male and female identity as crossed random effects. We tested that model assumptions were met by plotting model residuals against each variable and fitted values (Zuur et al. 2009). Because fertilization success was close to saturation (100%) in most crosses at the high sperm concentration, inclusion of that treatment introduced a nonlinearity that led to
EVOLUTION 2014
5
M . W. H A RT E T A L .
violation of the model assumption of equal variances, and we therefore excluded the high treatment from these models. For the two positively selected codons that showed strong population differences in allele frequencies (one in bindin, one in OBi1), we fitted a similar linear model to test effects of female OBi1 genotype and male bindin genotype (and their interaction) on fertilization success. Individual genotypes were categorized by the presence or absence of the less common allele at the positively selected male- or female-expressed gene. The linear model included population and sperm concentration as fixed effects. Using this approach, a significant interaction between male and female genotype would indicate that particular pairings of genotypes have greater fertilization success, consistent with a hypothesis of selection on both genes leading to their coevolution. To avoid having to account for repeated sampling of individuals, we used only within-population mating crosses and thus treated the 10 crosses within each population as independent replicates.
Results CODON MODELS OF POSITIVE SELECTION
We found little evidence of positive selection acting on EBR1 coding sequences in P. miniata and no evidence of population divergence under selection (Hart 2013). However, we found strong evidence for population divergence among coding sequences for the other bindin receptor OBi1, the sea star homolog (Hart and Foster 2013) of the gene that encodes the hsp70-like sperm receptor in sea urchin eggs (Foltz et al. 1993; Vacquier 2012). One partition (codons 613–700; Dataset S1) of OBi1 consisted of monophyletic clades of northern and southern alleles (Fig. 1A). We found one site (a serine-isoleucine polymorphism at codon 659, P = 0.0076; Table S1) under positive selection in a MEME analysis. Evidence for positive selection at that codon is especially strong because the serine-isoleucine polymorphism includes nonsynonymous substitutions at both the first and second codon positions, and implies a two-step change in amino acid sequence via a phenylalanine or a threonine intermediate that was not observed in any of our population samples. Analyses focused on variation among lineages revealed complementary evidence of selection that was consistent with this MEME result. We found one lineage (red in Fig. 1, P = 0.01) under positive selection for high ω in a BSR analysis. A branch sites model focused on that single basal lineage estimated a high posterior probability (P = 0.99) of strong positive selection on codon 659 (ω = 291; i.e., a high but imprecisely estimated value that was much greater than 1), and weaker evidence of positive selection at the adjacent codon 658 (P = 0.90). Among all of our codon model analyses of egg-expressed genes, this was the only partition in any gene that produced such a genealogy of recip-
6
EVOLUTION 2014
rocally monophyletic population samples (Table S1; Appendix S1), and the only locus showing positive selection for population differences. A similar pattern of positive selection for population differences was previously found among the same population samples in bindin (Sunday and Hart 2013). In that previous analysis, several different codons were found to be under positive selection within the southern population and a different set of codons under positive selection in the northern population, with just a single codon (site 842 in the coding sequence alignment) under positive selection within both the northern and the southern population samples. Most remarkably, the region of the OBi1 gene under positive selection for population differences is also the part of the modeled OBi1 protein structure predicted to bind the sperm ligand bindin (Dataset S3). OBi1 codon 659 encodes part of a small loop between alpha helices in the predicted substrate-binding domain of this hsp70-like bindin receptor (Fig. 1C). In structural and functional studies of bacterial and animal hsp70 proteins, these alpha helices together form a molecular lid that interacts with the substrate peptide within the binding site groove (Mayer et al. 2000), and the lid includes several amino acid sites in the region corresponding to P. miniata OBi1 codon 659 that mediate the strength of interaction between the substrate peptide and the substrate-binding site (e.g., Aponte et al. 2010). Using a predictive model (Van Durme et al. 2009) for identifying hsp70 binding site substrates, we found that P. miniata bindin contains a single substrate for the predicted OBi1 binding site: the LGLLLRHLRHH amino acid motif in the bindin core domain (also called the B18 domain of sea urchin bindin; Ulrich et al. 1998; Vacquier 2012). This domain has not previously been identified as an OBi1 binding site substrate, but it is strongly conserved among bindin coding sequences of sea urchins and sea stars, and mediates sperm–egg fusion at fertilization in sea urchins (Foltz 1994). ISOLATION-WITH-MIGRATION MODELS OF GENE FLOW
Our analysis of gene flow in partitions of egg receptor genes and bindin under positive selection showed a clear signal of population-specific selection against immigrant alleles. Group A partitions of bindin and OBi1 under positive selection had remarkably low gene flow (m, the proportion of gene copies that are new immigrants into a population) relative to other parts of the genome: m < 1 × 10−7 per generation from the northern to the southern population and m = 4.8 × 10−6 from the southern to the northern population, with posterior probability distributions that were strongly skewed toward m = 0 (Fig. 2A). Likelihood ratio testing of nested population models with gene flow fixed at m = 0 in either direction gave nonsignificant changes in model likelihood scores, consistent with zero gene flow (into the south P = 0.94; into the north P = 0.10; Dataset S2). However, a model
G A M E T E R E C O G N I T I O N C O E VO L U T I O N I N S E A S TA R S
Figure 2.
Gene flow is zero for positively selected parts of bindin and OBi1. Posterior probability distributions for gene flow (m) into the
northern (cyan) or into the southern population (green) from isolation-with-migration models fitted to data from (A) the two positively selected partitions of bindin plus OBi1; (B) three other partitions of bindin plus OBi1; (C) five other genes encoding egg coat proteins; or (D) five other highly expressed genes in egg transcriptomes (see Table S1 and Appendix S1 for codon model results; see Datasets S1 and S2 for IMa input and output files). Comparisons of the shapes of the posterior distributions to a null hypothesis with m = 0 (using likelihood ratio tests) are used to distinguish zero from nonzero gene flow. Note the wider prior distribution (m = 0–0.002 generation−1 ) needed to characterize gene flow posterior distributions in (C) and (D). Inferences about the magnitude of gene flow based on likelihood ratio tests of nested models are indicated as m = 0 or m > 0 in each panel.
with both gene flow parameters fixed to zero was a significantly poorer fit to these data (P < 0.001), indicating that gene flow in one direction (probably into the northern population) may be low but not zero. We obtained similar insights into gene flow patterns (but with less well-resolved posterior probability distributions for the parameter estimates) when we dropped one of the positively selected partitions from the analysis (results not shown). This consistency across genes suggests that both bindin and OBi1 experience population-specific selection against immigrant alleles. Analyses of the partitions of OBi1 plus bindin not subject to positive selection (group B), five other gamete recognition genes (group C), and a comparison set of five other egg-expressed genes (group D) all resulted in gene flow estimates that were one to two orders of magnitude higher (m10−5 to 10−4 ), and in all cases the likelihood ratio test indicated that immigration in both directions was significantly greater than zero (P
INCIPIENT SPECIATION OF SEA STAR POPULATIONS BY ADAPTIVE GAMETE RECOGNITION COEVOLUTION Michael W. Hart,1,2 Jennifer M. Sunday,1 Iva Popovic,1 Kevin J. Learning,1 and Christine M. Konrad1 1
Department of Biological Sciences, Simon Fraser University, Burnaby, British Columbia, Canada 2
E-mail: [email protected]
Received September 3, 2013 Accepted December 16, 2013 Reproductive isolation—the key event in speciation—can evolve when sexual conflict causes selection favoring different combinations of male and female adaptations in different populations. Likely targets of such selection include genes that encode proteins on the surfaces of sperm and eggs, but no previous study has demonstrated intraspecific coevolution of interacting gamete recognition genes under selection. Here, we show that selection drives coevolution between an egg receptor for sperm (OBi1) and a sperm acrosomal protein (bindin) in diverging populations of a sea star (Patiria miniata). We found positive selection on OBi1 in an exon encoding part of its predicted substrate-binding protein domain, the ligand for which is found in bindin. Gene flow was zero for the parts of bindin and OBi1 in which selection for high rates of amino acid substitution was detected; higher gene flow for other parts of the genome indicated selection against immigrant alleles at bindin and OBi1. Populations differed in allele frequencies at two key positively selected sites (one in each gene), and differences at those sites predicted fertilization rate variation among male–female pairs. These patterns suggest adaptively evolving loci that influence reproductive isolation between populations. KEY WORDS:
bindin, coevolution, fertilization, sexual conflict, sexual selection.
Understanding the targets of selection for different locally adapted traits among populations of a single species, and the causes of reproductive isolation between them, is an important outstanding goal of evolutionary genetics (Coyne and Orr 2004; Schluter 2009; Ellegren et al. 2012; Feder et al. 2012; Nosil and Feder 2012). Sexual selection theory predicts that such targets of selection should include genes encoding molecules that mediate conflict between the sexes over interactions with different optimal values for males versus females, such as copulation rates between mates or encounter rates between gametes (Rice 1989; Gavrilets 2000; Levitan and Ferrell 2006; Palumbi 2009). Males often compete to mate with females and fertilize eggs, but entry of multiple sperm in a single egg (polyspermy) can disrupt mitosis and kill the zygote. Because sperm are abundant, small, and individually cheap compared to eggs, the costs of frequent mating between individuals or high fertilization rates between gametes
C
1
2014 The Author(s). Evolution
are typically lower for males and higher for females, and the consequences differ between the sexes: selection favors sperm that fertilize eggs first, but favors eggs that are fertilized once (Frank 2000; Rice 2000; Levitan and Ferrell 2006). This conflict between the reproductive interests of the two sexes can lead to an arms race between male adaptations (e.g., for maximizing sperm encounters with eggs) and female countermeasures (e.g., for repelling supernumerary sperm and avoiding polyspermy), and such a conflict may be resolved in different ways within different populations (Holland and Rice 1998; Frank 2000; Gavrilets 2000). This mechanism of population divergence in reproductive traits is of considerable interest because it can account for the rapid evolution of reproductive isolation among populations in similar habitats or environments (Mani and Clarke 1990; Schluter 2009; Nosil and Flaxman 2011), and because it can arise through conflictual interactions at the level of interacting
M . W. H A RT E T A L .
genes and molecules, such as sperm ligands and their egg coat receptors (Levitan and Ferrell 2006; Palumbi 2009; Vacquier and Swanson 2011; Vacquier 2012). Such a mechanism for speciation has strong theoretical (Gavrilets 2000; Gavrilets and Hayashi 2005) and experimental support (Rice 1989; Holland and Rice 1998; Civetta and Clark 2000; Stewart et al. 2005), but there is only limited evidence for this mechanism in the formation of incipient species in nature (Nosil and Flaxman 2011). Coyne and Orr (2004) suggested that additional evidence might come from analyzing both molecular divergence and sperm–egg compatibility among conspecific populations that have been geographically isolated for a long time, and have ample opportunity for coevolution of male and female traits under local selection. We followed this line of study in populations of a sea star (Patiria miniata) from northern and southern British Columbia that are separated by a relatively old (>200,000 years) phylogeographic break in the northeastern Pacific ocean (Keever et al. 2009; McGovern et al. 2010; Fig. 1). We recently reported significant population differences caused by population-specific responses to selection on the gene encoding the sperm acrosomal protein bindin (Sunday and Hart 2013). Here, we document molecular and functional coevolution between bindin and the egg surface receptor OBi1 (Hart and Foster 2013) in the same populations studied previously, including positive selection for amino acid differences between populations, low gene flow (and selection against immigrant alleles) for bindin and OBi1 compared to other loci, and correlations between bindin and OBi1 genotypes and fertilization rates. The congruence among these structural, functional, and genetic patterns strongly suggests that coevolution of this pair of male- and female-expressed gamete recognition genes driven by sexual conflict within sea star populations could serve as the basis for incipient species formation among diverging populations.
Materials and Methods BACKGROUND AND STUDY SYSTEM
Our analysis of population divergence in gamete recognition molecules and sperm–egg compatibility in P. miniata builds on previous studies of bindin evolution and fertilization ecology in sea urchins (Palumbi 2009; Lessios 2011), and biochemical or molecular genetic analyses of other gamete recognition proteins in sea urchins and sea stars (Hirohashi et al. 2008; Vacquier 2012). In sea urchins, bindin and its egg coat receptors mediate sperm binding to eggs, fertilization success within species, and reproductive isolation between species via molecular interactions between the bindin mass on the surface of the activated sperm head and glycosylated proteins in the egg extracellular matrix (Geyer
2
EVOLUTION 2014
and Palumbi 2003; Levitan and Ferrell 2006; Palumbi 2009; Vacquier 2012). In comparisons among some sea urchin populations and species, bindin variation shows high rates of nonsynonymous amino acid replacement substitutions (dN) relative to synonymous silent nucleotide substitutions (dS), with different codons and lineages under positive selection (ω = dN/dS > 1; e.g., Geyer and Palumbi 2003). These patterns of positive selection have been interpreted mainly in terms of the expected effects of sperm competition and sexual conflict, and variation in those effects among species with different mating system traits or among habitats with different ecological conditions for fertilization (Levitan 2008, 2012; Palumbi 2009; Lessios 2011; Vacquier and Swanson 2011). The key evidence (reviewed in the studies cited above) that implicates sexual conflict in bindin evolution and speciation in sea urchins (and analogous systems of fertilization genes in other organisms) comes from (1) experimental studies of fertilization rates under conditions of sperm competition and high polyspermy risk, and (2) codon model analyses of positive selection in the molecular coevolution of bindin and its egg coat receptor. Unfortunately, there are few case studies of coevolution between interacting male and female fertilization genes (Vacquier and Swanson 2011). Excellent examples come from studies of gastropod lysin and its egg receptor (Clark et al. 2009; Hellberg et al. 2012; Aagaard et al. 2013), but research in this area has emphasized comparisons between phenotypically divergent species at a relatively late stage in the speciation process, and does not directly address the early origins of population divergence leading to incipient species. In sea urchins, the interpretation of both experimental studies of fertilization rate and codon model analyses of selection has been complicated by obstacles to the direct study of variation and molecular evolution of the bindin receptor in the egg coat (Vacquier 2012). The sea urchin egg bindin receptor gene EBR1 (Kamei and Glabe 2003) has a large, complex, repetitive coding sequence structure (with many possible targets for selection distributed among >50 exons) that has been difficult to analyze using population genetic methods. Other genes encode proteins that interact with bindin, and might coevolve with bindin under selection, but have been more intensively studied by cell biologists than by evolutionary geneticists (Foltz et al. 1993; Foltz 1994; Hirohashi et al. 2008; Vacquier 2012). As a result (1) experimental studies of fertilization success among mates with different gamete recognition genotypes have so far used female bindin genotype as a proxy for the egg bindin receptor (Palumbi 1999; Levitan and Ferrell 2006); and (2) analyses of egg bindin receptor evolution have been limited to comparisons among full-length EBR1 coding sequences for a small number of species (Kamei and Glabe 2003) or codon model analyses of a small part of the EBR1 coding sequence for larger samples of species and populations (Pujolar and Pogson 2011).
G A M E T E R E C O G N I T I O N C O E VO L U T I O N I N S E A S TA R S
Figure 1.
Patiria miniata populations differ at positively selected codons in the OBi1 substrate-binding site for bindin. (A) Maximum
likelihood phylogeny of the single positively selected partition from an alignment of 12 OBi1 alleles obtained from transcriptomes (RNASeq data) of northern (cyan) and southern (green) females. The gene tree is drawn using midpoint rooting and HKY genetic distances (scale bar shows 0.008 substitutions per nucleotide site). Labels indicate individual (1–3), population (N, S), and allele (a, b). Positive selection (ω >> 1) was inferred for two codons along the foreground branch shown as the red line. Populations were fixed for Ser-Thr or Ile-Ser differences at codons 658 and 659. (B) Locations of study sites (asterisks) in diverging northern (Haida Gwaii, cyan) and southern (Vancouver Island, green) populations in the northeastern Pacific Ocean. (C) Structural model prediction for OBi1 inferred using SWISSMODEL and rendered in Chimera. A southern allele is shown (northern and southern alleles had nearly identical structures). Backbone structure is rendered as ribbons. Two binding sites identified by SMART are shown in blue: the substrate-binding domain (ribbons), and the ATP-binding domain (spheres). Positively selected codons 658 and 659 are shown in red, between two of the alpha helices predicted to form a lid over the substrate-binding groove.
To extend this line of study to genes encoding interacting protein products from diverging conspecific populations, we first characterized population differences in sperm–egg compatibility. We designed our fertilization experiment specifically to test for local adaptation in fertilization success within a southern and a northern P. miniata population. Because we did not know the genotypes of individual sea stars before running the fertilization experiment, the males and females in this experiment included an unbiased sample of bindin and OBi1 (and other) variation. We then characterized population divergence in bindin (Pati˜no et al. 2009; Sunday and Hart 2013) using the same samples of males and females from the fertilization experiments. After discovering strong population divergence under positive selection at specific bindin sites, we then used RNA-Seq methods (Hart 2013; Hart and Foster 2013) to characterize coding sequence variation in P. miniata homologs of egg coat genes that encode proteins known or likely to interact with sperm (and specifically with bindin) on
the basis of previous studies in sea urchins. We then reanalyzed the fertilization data to test for an interactive effect of variation at positively selected sites in two of these genes (bindin and OBi1) on fertilization success. In the sections below, we describe the population genetic methods and results first, followed by the population differences in fertilization rates and the effects of individual male and female genotypes on sperm–egg compatibility, although this was not the same order in which we analyzed those different types of data. CODON MODELS OF POSITIVE SELECTION
We analyzed coding sequences of OBi1 and five other genes that encode proteins expressed in the egg extracellular matrix, including EBR1, three members of the ARIS gene family that induce the acrosome reaction in sperm, and Rendezvin. We compared those patterns to five other genes that are expressed in eggs or ovaries but are not implicated in gamete recognition, including Dual
EVOLUTION 2014
3
M . W. H A RT E T A L .
oxidase 1, Histone deacetylase 1, Mothers against decapentaplegic 1, Vasa, and Vitellogenin (Appendix S1). The data come from six ovary transcriptomes (three from northern females, three from southern females). All methods for generating, assembling, and analyzing those RNA-Seq data were described previously and are given in full detail elsewhere (Hart 2013; Hart and Foster 2013). Illumina RNA-Seq reads are available from the NCBI Sequence Read Archive (BioProject PRJNA175319); contig assemblies of the transcriptome data can be downloaded from the Dryad digital repository. Alignments (Dataset S1) for each egg-expressed gene were analyzed for evidence of recombination (Kosakovsky Pond et al. 2006; Delport et al. 2010) and divided into two or more partitions corresponding to nonrecombining blocks of codons (see Hart 2013). Input data for each subsequent analysis was an alignment of 12 haplotypes (two gene copies per individual, three individuals per population) for one data partition. Individual partitions for each gene were analyzed using codon models: we used the MEME method (Murell et al. 2012) to identify codons under positive selection (on some unspecified branches in the gene tree), and the branch-site random effects likelihood model (BSR; Pond et al. 2011) to identify lineages of alleles under positive selection (at some unspecified codons in the sequence alignment). Both methods use genetic algorithms that test hypotheses about among-site and among-lineage variation in ω, while avoiding false positives associated with the older branch sites model of positive selection (Yang 2007). For one data partition in which we found congruent evidence of positive selection in both the MEME and BSR models, we used the branch sites model to test the specific hypothesis of divergence between northern and southern populations. We compared the results from these analyses of several egg-expressed genes from ovary transcriptomes (Appendix S1, Table S1) to our recently published analyses of bindin (Sunday and Hart 2013), which were obtained from a larger sample of individuals from the same populations by polymerase chain reaction (PCR) and Sanger sequencing of cloned amplicons. ISOLATION-WITH-MIGRATION MODELS OF GENE FLOW
For genes that were inferred to be under positive selection (in the codon model analyses described above), we estimated rates of gene flow between the northern and southern P. miniata populations, and compared those results to gene flow estimated for loci that were found not to be under positive selection. We used these analyses to test the prediction that positively selected loci would have lower rates of gene flow between populations that reflected selection against immigrant alleles, compared to loci without evidence of population differences caused by selection. To obtain the necessary power for this analysis, we grouped loci into four multilocus datasets based on our codon model results (Table S1):
4
EVOLUTION 2014
one partition of bindin and one partition of OBi1 in which positive selection for population differences had been detected (group A), all other partitions from these same genes for which we found no evidence of positive selection (group B), all partitions of egg coat genes with no evidence of positive selection (group C), and all partitions of five other egg-expressed genes (group D). We used isolation-with-migration models of population demographic history to characterize patterns of gene flow between northern and southern populations using each of those four datasets A through D. We used data from our analyses of eggexpressed genes and our previously published bindin data. We treated partitions as loci and carried out multilocus analyses in IMa (Hey and Nielsen 2007; e.g., see Puritz et al. 2012; Keever et al. 2013). We used all partitions of all genes except Rendezvin (Table S1): for this gene, we were unable to infer all recombination sites simultaneously in one analysis; instead, we divided the Rendezvin alignment into three subalignments (with similar numbers of repetitive protein-coding domains), analyzed each of those subalignments for recombination, and included the longest partition from each as three loci in the IMa dataset. We assumed a generation time (average age of breeding adults) of five years for large, long-lived sea stars. This approach is similar to the IMa-based method of Sousa et al. (2013) that assigns loci to groups with or without evidence of selection against gene flow based on amonglocus differences in isolation-with-migration model parameters, but has the advantage of prior knowledge of among-locus variation in selection (based on the results of codon model analyses). Input files (including sequence alignments) and IMa output are available as Datasets S1 and S2. To standardize the gene flow estimates from different parts of the genome for comparison, we included an additional locus with a known mutation rate in each of the four IMa datasets. We used population samples of mitochondrial DNA haplotypes from the same Sandspit and Bamfield populations analyzed in our previous phylogeographic study of P. miniata populations (Keever et al. 2009). We used the same mtDNA mutation rate calibration as in our previous studies: the lowest observed divergence in mtDNA between geminate pairs of sea urchin species (4.57% for Diadema species) assuming a three million year divergence time (Lessios et al. 2001; Hickerson et al. 2006; McGovern et al. 2010; Puritz et al. 2012). Including a locus with an empirical mutation rate calibration such as this in the IMa analysis allowed us to convert the IMa model parameters estimated from each dataset (all in units scaled by the mutation rate specific to each set of data) to population parameters in the same demographic units (individuals, years) shared among all four datasets. A corollary of this approach is that, because the same single mtDNA coalescent pattern contributed to each set of multilocus population parameter estimates for the four groups of loci (Table S1; Dataset S1), significant differences among the four results are probably
G A M E T E R E C O G N I T I O N C O E VO L U T I O N I N S E A S TA R S
conservative estimates of real differences in demographic history caused by differences among those four parts of the P. minata genome. ADDITIONAL MOLECULAR GENETIC METHODS
Methods for characterizing predicted protein domains of coding sequences from transcriptomes, including the substrate- and ATP-binding domain of the hsp70-like bindin receptor, are given elsewhere (Hart 2013; Hart and Foster 2013). We used the Limbo method (Van Durme et al. 2009) to find domains in a full-length P. miniata bindin coding sequence (GenBank FJ439659.1) that are predicted to be substrates for the OBi1 substrate-binding site. The OBi1 structural model prediction based on the SWISS-MODEL method is given in Dataset S3. To obtain larger sample sizes of allele frequencies for one part of the OBi1 coding sequence that was under positive selection in our analyses of small samples of haplotypes from RNA-Seq data (and implicated in prezygotic reproductive isolation), we amplified and sequenced short (100–200 bp) portions of the OBi1 exon containing the positively selected codon 659. That PCR product was amplified from genomic DNA of all females and some males used in our fertilization experiments (see below). Custom primers (5 -tgtaaaacgacggccagtGTTGAACAAACACAGGGTCAG-3 , including the lower-case sequencing primer at the 5 tail; and 5 -GTTTTCYGTCTTGTTCAGGAG-3 ) were designed from RNA-Seq coding sequences aligned to an OBi1 genomic sequence from the draft P. miniata genome assembly (GenBank AKZP01045554.1). Amplicons varied in length due to an indel in the adjacent intron. We purified and direct-sequenced each amplicon using the M13 universal primer, and scored the three-nucleotide genotype of each individual at codon 659 directly from the chromatograms. FERTILIZATION EXPERIMENTS
We conducted full-factorial fertilization crosses in sperm-limiting noncompetitive conditions. Adult sea stars were collected at low tide from a northern population at Sandspit, Haida Gwaii, and from a southern population at Bamfield, Vancouver Island (see Fig. 1) and held at Simon Fraser University in approximately 1 L plastic containers in large recirculating UV-sterilized seawater tanks under constant low light conditions at 11°C to 12°C and fed regularly with mussels. We acclimated all adults to these laboratory conditions for at least two weeks before using them in fertilization experiments. We conducted 10 independent fertilization crosses on separate days. Each cross included a male and a female individual from each of the two populations, crossed reciprocally. Each individual was removed from seawater, brought to room temperature, and given an intracoelomic injection of 1 to 3 mL 100 mM 1-methyl-adenine, which induces maturation of oocytes (in females) and spawning (in both sexes). Eggs
were collected and washed in filtered (0.45 μm) seawater, then allowed to settle in a glass beaker at approximately 12°C. Concentrated sperm were collected directly from the gonopore and held on ice. Concentrated sperm were diluted 100-fold in filtered seawater. We fixed a subsample of this stock in 5% formalin in seawater, and measured absorbance using a Nanodrop 2000c spectrophotometer (Thermoscientific) to estimate cell concentration (absorbance at 600 vs. 280 nm as recommended by the manufacturer) relative to a predetermined standard curve based on hemocytometer counts of fixed sperm in a dilution series. The experimental sperm stock was then independently diluted to 104 , 105 , and 106 mL−1 in 20 mL of filtered seawater. Aliquots of sperm dilutions used in experimental treatments were also fixed in 5% formalin and counted on a hemocytometer after each experiment to confirm sperm concentrations. Four milliliters of each sperm treatment was immediately dispensed into two 5 mm sterile plastic Petri dishes. One drop of settled concentrated eggs from each female (approximately 0.1 mL) was added to each dish containing sperm (and to a filtered seawater control), and the inseminated eggs were cultured at approximately 12°C. For each dish, we scored fertilization success for the first 100 eggs observed under a dissecting microscope by (1) presence or absence of fertilization envelopes at 30 min postinsemination, (2) presence or absence of cleavage after 2 h, and (3) development to the blastula stage after 24 h. The three measures of fertilization success were highly correlated, and we report proportion of blastula embryos after 24 h because this measure integrates the effects of both successful fertilization encounters (between an egg and single sperm) and unsuccessful encounters (between an egg and multiple sperm leading to polyspermy and embryo death before blastula formation). We tested for localized coevolution of sperm–egg compatibility using a generalized mixed effects linear model of fertilization rate with a binomial error distribution for proportion data (Crawley 2007); sperm concentration, population source of males, population source of females, and the interaction of male and female population source were fixed effects. With this approach, a positive coefficient for the interaction term would indicate a disproportionately greater rate of fertilization success between males and females from the same population. To account for the effect of crossing each male and female to two different mates (one from the same population source, and one from the other population), we included both male and female identity as crossed random effects. We tested that model assumptions were met by plotting model residuals against each variable and fitted values (Zuur et al. 2009). Because fertilization success was close to saturation (100%) in most crosses at the high sperm concentration, inclusion of that treatment introduced a nonlinearity that led to
EVOLUTION 2014
5
M . W. H A RT E T A L .
violation of the model assumption of equal variances, and we therefore excluded the high treatment from these models. For the two positively selected codons that showed strong population differences in allele frequencies (one in bindin, one in OBi1), we fitted a similar linear model to test effects of female OBi1 genotype and male bindin genotype (and their interaction) on fertilization success. Individual genotypes were categorized by the presence or absence of the less common allele at the positively selected male- or female-expressed gene. The linear model included population and sperm concentration as fixed effects. Using this approach, a significant interaction between male and female genotype would indicate that particular pairings of genotypes have greater fertilization success, consistent with a hypothesis of selection on both genes leading to their coevolution. To avoid having to account for repeated sampling of individuals, we used only within-population mating crosses and thus treated the 10 crosses within each population as independent replicates.
Results CODON MODELS OF POSITIVE SELECTION
We found little evidence of positive selection acting on EBR1 coding sequences in P. miniata and no evidence of population divergence under selection (Hart 2013). However, we found strong evidence for population divergence among coding sequences for the other bindin receptor OBi1, the sea star homolog (Hart and Foster 2013) of the gene that encodes the hsp70-like sperm receptor in sea urchin eggs (Foltz et al. 1993; Vacquier 2012). One partition (codons 613–700; Dataset S1) of OBi1 consisted of monophyletic clades of northern and southern alleles (Fig. 1A). We found one site (a serine-isoleucine polymorphism at codon 659, P = 0.0076; Table S1) under positive selection in a MEME analysis. Evidence for positive selection at that codon is especially strong because the serine-isoleucine polymorphism includes nonsynonymous substitutions at both the first and second codon positions, and implies a two-step change in amino acid sequence via a phenylalanine or a threonine intermediate that was not observed in any of our population samples. Analyses focused on variation among lineages revealed complementary evidence of selection that was consistent with this MEME result. We found one lineage (red in Fig. 1, P = 0.01) under positive selection for high ω in a BSR analysis. A branch sites model focused on that single basal lineage estimated a high posterior probability (P = 0.99) of strong positive selection on codon 659 (ω = 291; i.e., a high but imprecisely estimated value that was much greater than 1), and weaker evidence of positive selection at the adjacent codon 658 (P = 0.90). Among all of our codon model analyses of egg-expressed genes, this was the only partition in any gene that produced such a genealogy of recip-
6
EVOLUTION 2014
rocally monophyletic population samples (Table S1; Appendix S1), and the only locus showing positive selection for population differences. A similar pattern of positive selection for population differences was previously found among the same population samples in bindin (Sunday and Hart 2013). In that previous analysis, several different codons were found to be under positive selection within the southern population and a different set of codons under positive selection in the northern population, with just a single codon (site 842 in the coding sequence alignment) under positive selection within both the northern and the southern population samples. Most remarkably, the region of the OBi1 gene under positive selection for population differences is also the part of the modeled OBi1 protein structure predicted to bind the sperm ligand bindin (Dataset S3). OBi1 codon 659 encodes part of a small loop between alpha helices in the predicted substrate-binding domain of this hsp70-like bindin receptor (Fig. 1C). In structural and functional studies of bacterial and animal hsp70 proteins, these alpha helices together form a molecular lid that interacts with the substrate peptide within the binding site groove (Mayer et al. 2000), and the lid includes several amino acid sites in the region corresponding to P. miniata OBi1 codon 659 that mediate the strength of interaction between the substrate peptide and the substrate-binding site (e.g., Aponte et al. 2010). Using a predictive model (Van Durme et al. 2009) for identifying hsp70 binding site substrates, we found that P. miniata bindin contains a single substrate for the predicted OBi1 binding site: the LGLLLRHLRHH amino acid motif in the bindin core domain (also called the B18 domain of sea urchin bindin; Ulrich et al. 1998; Vacquier 2012). This domain has not previously been identified as an OBi1 binding site substrate, but it is strongly conserved among bindin coding sequences of sea urchins and sea stars, and mediates sperm–egg fusion at fertilization in sea urchins (Foltz 1994). ISOLATION-WITH-MIGRATION MODELS OF GENE FLOW
Our analysis of gene flow in partitions of egg receptor genes and bindin under positive selection showed a clear signal of population-specific selection against immigrant alleles. Group A partitions of bindin and OBi1 under positive selection had remarkably low gene flow (m, the proportion of gene copies that are new immigrants into a population) relative to other parts of the genome: m < 1 × 10−7 per generation from the northern to the southern population and m = 4.8 × 10−6 from the southern to the northern population, with posterior probability distributions that were strongly skewed toward m = 0 (Fig. 2A). Likelihood ratio testing of nested population models with gene flow fixed at m = 0 in either direction gave nonsignificant changes in model likelihood scores, consistent with zero gene flow (into the south P = 0.94; into the north P = 0.10; Dataset S2). However, a model
G A M E T E R E C O G N I T I O N C O E VO L U T I O N I N S E A S TA R S
Figure 2.
Gene flow is zero for positively selected parts of bindin and OBi1. Posterior probability distributions for gene flow (m) into the
northern (cyan) or into the southern population (green) from isolation-with-migration models fitted to data from (A) the two positively selected partitions of bindin plus OBi1; (B) three other partitions of bindin plus OBi1; (C) five other genes encoding egg coat proteins; or (D) five other highly expressed genes in egg transcriptomes (see Table S1 and Appendix S1 for codon model results; see Datasets S1 and S2 for IMa input and output files). Comparisons of the shapes of the posterior distributions to a null hypothesis with m = 0 (using likelihood ratio tests) are used to distinguish zero from nonzero gene flow. Note the wider prior distribution (m = 0–0.002 generation−1 ) needed to characterize gene flow posterior distributions in (C) and (D). Inferences about the magnitude of gene flow based on likelihood ratio tests of nested models are indicated as m = 0 or m > 0 in each panel.
with both gene flow parameters fixed to zero was a significantly poorer fit to these data (P < 0.001), indicating that gene flow in one direction (probably into the northern population) may be low but not zero. We obtained similar insights into gene flow patterns (but with less well-resolved posterior probability distributions for the parameter estimates) when we dropped one of the positively selected partitions from the analysis (results not shown). This consistency across genes suggests that both bindin and OBi1 experience population-specific selection against immigrant alleles. Analyses of the partitions of OBi1 plus bindin not subject to positive selection (group B), five other gamete recognition genes (group C), and a comparison set of five other egg-expressed genes (group D) all resulted in gene flow estimates that were one to two orders of magnitude higher (m10−5 to 10−4 ), and in all cases the likelihood ratio test indicated that immigration in both directions was significantly greater than zero (P
