Molecular evolution of candidate male reproductive genes in the brown algal model Ectocarpus

Background Evolutionary studies of genes that mediate recognition between sperm and egg contribute to our understanding of reproductive isolation and speciation. Surface receptors involved in fertilization are targets of sexual selection, reinforcement, and other evolutionary forces including positive selection. This observation was made across different lineages of the eukaryotic tree from land plants to mammals, and is particularly evident in free-spawning animals. Here we use the brown algal model species Ectocarpus (Phaeophyceae) to investigate the evolution of candidate gamete recognition proteins in a distant major phylogenetic group of eukaryotes. Results Male gamete specific genes were identified by comparing transcriptome data covering different stages of the Ectocarpus life cycle and screened for characteristics expected from gamete recognition receptors. Selected genes were sequenced in a representative number of strains from distant geographical locations and varying stages of reproductive isolation, to search for signatures of adaptive evolution. One of the genes (Esi0130_0068) showed evidence of selective pressure. Interestingly, that gene displayed domain similarities to the receptor for egg jelly (REJ) protein involved in sperm-egg recognition in sea urchins. Conclusions We have identified a male gamete specific gene with similarity to known gamete recognition receptors and signatures of adaptation. Altogether, this gene could contribute to gamete interaction during reproduction as well as reproductive isolation in Ectocarpus and is therefore a good candidate for further functional evaluation. Electronic supplementary material The online version of this article (doi:10.1186/s12862-015-0577-9) contains supplementary material, which is available to authorized users.


Background
Sexual selection and the evolution of mating specificity are a central focus to evolutionary biology [1]. Recent advancement in molecular tools allowed for identification of specific genes in several organisms controlling particular stages of male-female interactions at fertilization (e.g., [2,3]) and with the aid of next-generation sequencing methods we are now closer to understand the evolution of mating specificity [4]. In particular, evolutionary studies of genes involved in recognition between male and female gametes (GRPs) provide important insights into the evolution of reproductive isolation and speciation. It has been shown, that fertilization genes generally evolve faster and are targets of sexual selection, reinforcement, and other evolutionary forces (reviewed in [5][6][7][8]. Strength of selection on reproductive genes increases with the potential risk of cross-species hybridization or polyspermy within species, and may be correlated with the rate of adaptive evolution. If progressing genetic isolation would be manifested by reduced fitness or sterility of hybrids, selection acting on the genes involved in fertilization may result in the formation of prezygotic reproductive barriers and consequently enhance probabilities of speciation [9][10][11][12]. Signatures of rapid evolution of reproductive genes have been detected across different lineages of the eukaryotic tree, including land-plants [13], insects [14], invertebrates [15,16] and mammals [17]. For example, genes underlaying fertilization or mating traits in Drosophila species showed a lack of evolutionary constraints, suggesting that they were most likely shaped during the early stages of speciation by directional sexual selection [18]. However, complicated mating systems as found in copulating animals or animal-pollinated plants involve synergy of many genetic traits, and make it difficult to disentangle the genetic components under sexual selection. In contrast, marine free-spawning species offer a relatively simple model in which male-female interactions are determined by spatial and temporal factors such as the time and place of gamete release; and sexual selection is narrowed down to the gamete interactions [3]. The latter involve specific recognition processes mediated by peptides or glycoproteins expressed on the surface of sperm or egg (reviewed in [19]) which have been shown to be under strong directional selection (reviewed in [3]). Research conducted on sea urchin bindin, an acrosomal protein binding in a species-specific manner to a receptor on the egg plasma membrane prior to the fusion, is exemplar in this respect (e.g., [16,20,21]). Similarly, studies of the sperm receptor lysin and its egg ligand VERL in Haliotis revealed a strong excess of nonsynonymous to synonymous divergence (d N /d S ) between species driven by positive selection and concerted evolution [22,23], making it one of the fastest evolving metazoan proteins known. In addition, variable selective pressure exerted on the lysin gene in closely related (ω >1) and more diverged species (ω <1) of abalone suggested, that diversifying selection acts on closely related sympatric species, whereas distantly related species are already relieved from it [15,24]. This could be indicative of reinforcement [25] where sex involved genes would be under selective pressure to establish barriers to reproduction in reunited populations.
Contrary to the rich literature on GRPs in marine invertebrates, studies identifying sperm-egg recognition proteins and their putative adaptive evolution are relatively scarce in other branches of the eukaryotic tree. Here, we focus on the evolution of genes hypothesized to be involved in sexual reproduction and gamete recognition in the brown algal model Ectocarpus sp. (Phaeophyceae, Stramenopiles) [26,27]. Brown algae present an opportunity to investigate gamete receptors in a lineage, that has been evolving independently from land plants and animals for over a billion years [28,29]. Similarly to sea urchin and abalone, brown algae represent free spawning species with gamete interaction limited to pheromone signalling and surface recognition. Although proteins involved in fertilization in brown seaweeds have not been described so far, a number of proteins involved in spermegg binding were isolated and partially characterized in Fucus [30][31][32][33][34], but the underlying genes were never identified. In addition to the studies focusing on Fucus, a Sexually Induced Gene 1 (Sig1) in a diatom Thalassiosira spp. was shown to have high divergence, both within and between species [35]. Transcription of the gene is upregulated during mating [36], however its exact function is not known. Nevertheless, its extreme divergence indicates a possible function as a barrier to hybridization between geographically distant strains [37,38]. Homologs of the family of Sig genes have also been found in other Stramenopiles, including Ectocarpus [26,39].
Ectocarpus is a cosmopolitan genus composed of three recognized species: E. siliculosus, E. fasciculatus and E. crouaniorum. Nonetheless, cross-fertilization experiments imply that more species exist in reality [40,41]. Previous experiments suggested that gamete recognition is mediated by N-acetylglucosamine (GlcNAc) residues exposed on the plasma membrane of the female gametes and a lectin-like receptor at the tip of male anterior flagella [42,43]. However, detailed structural information on the proteins involved and gene sequences for the corresponding receptors are lacking.
Here, we used large scale genomic data covering different life stages of Ectocarpus [44,45] to identify genes with expression limited to the male gametic stage, expecting surface receptors mediating fertilization to be expressed uniquely in gametes. These genes (Additional file 1) provided a subset to test for divergence and positive selection at the amino acid level using Ectocarpus species of known sexual compatibility (Tables 1 and 2, Fig. 1). We found divergence-based evidence of selective pressure acting on at least one of the investigated genes. Interestingly, that gene displayed domain similarities to the receptor for egg jelly (REJ) protein involved in sperm-egg recognition in sea urchins [46].

Selection of genes of interest
Transcriptome data covering male and female gametes [45] as well as gametophyte and sporophyte life stages [44] were compared to obtain a list of male gamete specific genes. An expression filter of RPKM > 1 was applied to remove background noise. Selected genes were further screened for the presence of transmembrane domains using TMHMM [47,48] and signal peptides targeting the outer membrane using HECTAR [49].

Sequence data collection
To perform pairwise divergence analysis, we used genomic data available for the genome strain Ectocarpus sp. linage 1c [26,50] and Ectocarpus siliculosus lineage 1a (S. Coelho, unpublished). The two lineages are considered separate 'species' with post-zygotic barriers to reproduction (A. Peters, personal communication). Additionally, cultures of nine Ectocarpus strains representing clades 1a, 1c, 2c, 2d, 4, 5a and 5b as described by Stache Crain et al. [41] were retrieved from the Culture Collection of Algae and Protozoa (Oban, Scotland) or macroalgal culture collection (Roscoff, France) for DNA extraction. Scytosiphon "lomentaria" (specimens ODC1349 and LT0114 deposited at Ghent University, Belgium) was used as an outgroup (Fig. 1, Table 1). Primers were designed based on the Ectocarpus genome sequence [26,50] using Primer3 [51,52] (Additional file 2). Gene fragments were amplified from genomic DNA using a touch-down PCR procedure with initial denaturation for 3 min at 95°C followed by 10 cycles of denaturation at 95°C for 30s, 30s annealing at 65°C decreasing 1°C per cycle and elongation at 72°C for 1 min and then 25 cycles of denaturation at 95°C for 30s, 30s annealing at 55°C and elongation at 72°C for 1 min, with a final elongation step of 10 min at 72°C. Amplicons were sequenced using Sanger sequencing. DNA chromatograms were edited and checked using BioNumerics (Applied Maths). Sequences can be retrieved under accession numbers LN901218 -LN901253. We refer to the Additional file 3 for full details.

Phylogenetic analyses
Sequences of the mitochondrial cytochrome oxidase subunit 3 (cox3) were generated for phylogenetic analysis. Obtained sequences were compared with Ectocarpus sequences in GenBank using blastn [53]. Since Bayesian inference and maximum likelihood methods may produce conflicting results, gene sequences were analyzed using both methodologies. Maximum likelihood analyses were carried out with RAxML version 7.7.1 [54] on the RAxML   [40,41,[85][86][87]. Clade numbers corresponding to the phylogenetic position in are given with the strain names. Ffull interfertility, zzygote formation, no data on growth, preprezygotic barriers, no cell fusion, posthybrids with reduced growth or non-functional reproductive structures. a No data on the actual Penikese strain (Pen) are available; fertilization data were inferred from Woods Hole, Massachusetts strain of similar restriction in the mating pattern with other strains, but completely interfertile with the Penikese strain blackbox server (http://phylobench.vital-it.ch/raxml-bb/) using the CAT model [55]. Searches were started from 200 distinct randomized maximum parsimony starting trees and branch support was assessed with the classic bootstrapping algorithm (1000 replicates). Bayesian phylogenetic inference was carried out with MrBayes version 3.2 [56]. Two independent runs, each consisting of four incrementally heated chains, were run for 3 million generations using default priors and other settings. Trees were sampled every thousand generations. Convergence of likelihood and parameter values were assessed with Tracer version 1.5 [57] and a suitable burn-in value was chosen (burnin = 500). Bayesian posterior probabilities for clades were computed from the post burn-in sample of trees and indicated on the Maximum Likelihood (ML) tree (Fig. 1).

Divergence and positive selection analyses
To estimate the rates of evolution of putative male reproductive genes (receptors), we performed a pairwise d N /d S analysis using coding sequences of Ectocarpus sp. genome strain lineage 1c and Ectocarpus siliculosus lineage 1a [41]. Orthologous sequences were aligned using ClustalW implemented in MEGA6 [58,59] and manually curated. A pairwise d N /d S (ω) analysis was performed using Phylogenetic Analysis by Maximum Likelihood (PAML, CODEML, F3x4 model, runmode = −2) [60]. Positive selection was estimated by the maximum likelihood method available in CODEML, PAML 4 using the F3X4 model of codon frequencies and paired nested site models of sequence evolution (M0, M3; M1a, M2a; M7, M8) [24,61,62]. Nested models were compared using the likelihood ratio test (LRT) with either 2 (M0, M3) or 4 degrees of freedom (M1a,M2a; M7, M8). Empirical Bayes methods allowed for identification of positively selected sites a posteriori [60,62]. Codon-based nucleotide alignments were used in conjunction with the species tree (see above). Individual exon sequences representing the same gene for a given strain were concatenated using FaBox [63].
Since signatures of positive selection might be overlooked due to its transient or episodic nature, we performed, in parallel, an analysis of occasional selection at individual sites using MEME (Mixed Effects Model of Evolution), implemented in HyPhy [64] available at Datamonkey.org server [65,66]. The F81 codon substitution model was used in HyPhy [67].

Protein functional characteristics prediction
Protein structure prediction was done using Phyre2 (Protein Homology/Analogy Recognition Engine) [68]. Transmembrane helices and their topology were inferred from memsat-svm implemented in Phyre2 or from the TMHMM server v. 2.0 [47,48]. To determine functional domains we performed a Gene Ontology (GO) and protein domain search (InterPro database) using Blast2GO v. 2.6.6 [69] with an E-value Hit Filter set to 1.0e-6.

Phylogeny
Maximum likelihood and Bayesian analyses of the mitochondrial cox3 gene dataset yielded essentially the same phylogenetic trees (Fig. 1), which divided Ectocarpus into 4 well-supported clades. The relationships among these clades, however, remained unresolved. Each clade is subdivided into subclades, again with near full support. While some of these subclades bear formal taxonomic names (e.g., E. crouaniorum, E. fasciculatus and E. siliculosus), others are known by informal identifiers only. Major clades are reproductively isolated by prezygotic barriers. For subclades such data are either not available or no prezygotic barriers were identified by means of no-choice experiments, however, post-zygotic barriers may exist (see Table 2).

Screening for candidate reproductive genes
Comparative analysis of transcriptomic data of male and female gametes [45], male and female immature gametophytes and sporophyte of Ectocarpus [44] identified 109 male gamete specific genes. These were investigated further for the presence of signal peptides, transmembrane helixes and functional domains potentially involved in cellcell recognition. Fertilization experiments in Ectocarpus have shown the presence of lectin-like receptors in male gametes localized on the anterior flagellum [42]. Therefore, we restricted the search to male expressed genes coding for extracellular or cell surface proteins with potential receptor activity. With these criteria, twelve out of 109 genes were identified as putative gamete recognition genes and were selected for evolutionary analyses together with twelve house-keeping genes as controls (Additional file 4). In addition, we included three genes belonging to the Sig family (Sexually induced genes), previously described as potential recognition genes in diatoms [35]. The Ectocarpus genome codes for three Sig-like proteins, two of which showed no expression in gametes (Esi0033_0091 and Esi0116_0079) and one homolog of Sig1 (Esi0101_0018) which was found abundant in both male and female gametes. It has to be noted, that during the initial phase of gamete release both male and female gametes are flagellated in Ectocarpus. Therefore, the presence of Sig1 in both sexes may indicate a structural role rather than the involvement in gamete recognition as hypothesized by Honda et al. [39] and Blackman et al. [70].

Evolution of putative male receptors
To test for evolutionary divergence of candidate male receptor genes, we calculated levels of nonsynonymous (d N ) and synonymous (d S ) substitution using pairwise comparisons between Ectocarpus sp. (lineage 1c) and its sister species Ectocarpus siliculosus (lineage 1a). Overall putative receptor genes showed significantly faster evolutionary rates (d N /d S ) compared to housekeeping genes (U test, p = 0.006), with Sig genes showing the highest sequence conservation among putative receptors. Six out of twelve selected male specific genes showed d N /d S ratio >0.5, which could imply adaptive evolution (Table 3, d N /d S >0.5 marked in bold). It has been shown that the gene length may bias the estimations of d N /d S particularly for short genes with low divergence which are then preferably found to be under positive selection [71]. However, higher d N /d S in this study were associated with GOIs, which were on average substantially longer than housekeeping genes (720 bp vs 358 bp, respectively). We therefore believe that the gene length does not have a major influence on the evolutionary rates of the genes analyzed here.
Esi0130_0068 was particularly interesting due to the presence of a REJ-like domain (IPR002859). REJ-like domains are found in the sperm proteins of sea urchins, where they mediate egg-sperm binding [21,46]. The remaining genes had either no functional domains (Esi0180_0035), could be involved in substrate transport (Esi0026_0077), lipid metabolism (Esi0132_0081, Esi0146_0068) or nucleic acid hydrolysis (Esi0660_0004) (Additional file 4). We therefore selected 130_0068, sig1 and a couple of house-keeping genes for sequencing in the 9 Ectocarpus strains to search for signatures of adaptive evolution using maximum likelihood method implemented in CODEML, PAML 4 [60] and HyPhy [64]. With this approach we obtained evidence that Esi0130_0068 is evolving under positive selection (Tables 4 and 5).

Interspecies polymorphism and evidence for positive selection in Esil0130_0068
The Esi0130_0068 protein consists of 1427 amino acids and presumably contains 7 (TMHMM algorithm) or 6 (memsat-svm algorithm) transmembrane domains (Fig. 2). The N-terminal region, composed of 976 amino acids, is predicted to be located extracellularly. A functional domain scan for this part of the protein identified a GPS domain (E-values = 3.78e-05), a REJ domain (E-value = 3.98e-30) and a polycystin cation channel (Na + , K + , Ca 2+ channel activated by Ca 2+ ) within the REJ domain (E-value = 3.01e-05). The N-terminal fragment was targeted for resequencing in representative Ectocarpus strains. Esi0130_ 0068 showed statistical evidence for adaptive evolution in the PAML 4 and MEME (HyPhy) analysis. However, the latter analysis presented different sites under selection depending on the model used (Table 4). All models allowing individual sites to evolve under positive selection (M3, M2a, M8) gave a significantly better fit to the Esi0130_0068 data (Table 5) and identified a substantial proportion of sites with d N /d S >1 (Fig. 2). This result is consistent with an evolutionary history characterized by frequent episodes of positive selection. All three models   (Fig. 3) and are also indicated by the Empirical Bayes analysis in PAML 4 with lower probability. Noteworthy, sequence analysis of sperm PKDREJ in primates also revealed several positively selected sites in the REJ domain and its flanking extracellular regions [72].

Discussion
Adaptive evolution has commonly been observed in proteins responsible for egg-sperm interactions (for a review see [5]) with a particularly large proportion of positively selected sites found in the sperm-egg binding moieties [6,[73][74][75]. This indicates that gamete recognition might be subject to selective pressure and sexual selection may operate at the gamete level [5]. Although proteins responsible for egg-sperm interactions in brown algae have not been identified so far, Sig1 was originally hypothesized to play a role in gamete adhesion in the diatom Thalassiosira [35,36]. Sig1 shows high nonsynonymous sequence divergence between closely related species of Thalassiosira [35]. Although Ectocarpus possesses one homolog of Sig1 with similar length and domain architecture, no evidence was obtained for positive selection in the C-terminal region (amino acids 514-721, exons 5-6). Taken into account that Sig1 is expressed in both types of gametes, these findings support the structural role as mastigoneme proteins rather than their involvement in gamete recognition. It should be noted that the remaining exons (1)(2)(3)(4), that might bear the positively evolving sites, could not be amplified by PCR from the genomic DNA and thus are missing in this analysis.  98  101  104  107  110  113  116  119  122  125  128  131  134  137  140  143  146  149  152  155  158  213  216  219  222  225  228  231  234  237  240  285  288  291  294  297  300  303  306  309  312  315  619  622  625  628  631  634  637  640  643  646  649  652  655  658  661  664  667  670  673  676  679  682  685  688  691  696  699  702  705  708  821  824  827  830  833  836  839  842  845  One of the genes expressed specifically in male gametes (Esi0130_0068) revealed significant variation in selective pressure acting on different amino acids. The changes in the evolutionary rates could be a result of weaker selective constrains acting on Esi0130_0068 and/or positive selection. The former could be caused by gene expression bias, if deleterious alleles are less effectively removed when expressed only by one sex in the population [76]. Therefore, gene expression bias could result in relaxation of purifying selection on protein-coding genes. This phenomenon is common in genes with sex-biased expression ( [77,78]; reviewed in [79]) and has been observed also in Ectocarpus [44]. The relaxed selective constrain may also facilitate specialized adaptation [76,80]. To test for signatures of adaptive evolution in the genes of interest, we performed a maximum likelihood test using PAML [60] and HyPhy [66]. Both methods found that divergence in Esi0130_0068 sequence is promoted by adaptive evolution. Analysis of d N /d S ratio among individual sites identified particular amino acids with good statistical support of positive selection (Table 6), indicative of their putative importance for the function of Esi0130_0068.
Pairwise d N /d S analysis of Esi0130_0068 showed highest values (ω > 1) in comparisons between closely related strains, which concerns clades 5a, 5b and 2d as well as clades 1a and 1c (Additional file 5, Fig. 1). While prezygotic barriers to fertilization already exist between clades 5a and 2d, clades 1a and 1b are able to form zygotes which are later arrested during development ( Table 2). In this case, excess of nonsynonymous to synonymous substitutions could be a sign of diversifying selection acting to reinforce reproductive barriers between strains 1a and 1c. Similar scenario could account for the high diversity between clades 5a and 5b, however compatibility studies are lacking in this group.
Interestingly, Esi0130_0068 resembles the topology of the egg recognition protein in sea urchin sperm [46]. Although the exact function of the REJ domain containing proteins is unknown in brown algae the unique localization and expression pattern of REJ proteins in sea urchins and humans (hPKDREJ) suggest a central role in fertilization [46,81]. Sea urchin sperm contains several members of the REJ protein family (named SpREJ1-10), out of which SpREJ1 binds a fucose sulfate polymer of the egg jelly, triggering the acrosome reaction and transforming the sperm into a fusogenic cell [21,46,82]. Ectocarpus Esi0130_0068 protein does not contain the lectin domain upstream of REJ, which probably interacts with the female  gamete surface glycoprotein in sea urchin, however this domain is also not found in the human hPKDREJ. In sea urchins, the carbohydrate recognition domains that are subjected to positive selection [75] in contrast to the sperm PKDREJ in primates were positively selected sites are found within and around the REJ domain [72]. It is unclear at present whether the REJ domain itself can take part in the interaction with the egg glycoprotein coat The REJ domain has the potential role in regulating the ion channel which triggers the acrosome reaction [81,83]. Therefore, the REJ domain may be important for triggering the species-specific recognition cascade through control of the ion flux. Additionally, almost all of the members of the SpREJ and human PKD families possess a G-protein coupled receptor cleavage site (GPS) upstream of the first transmembrane helix [82]. REJ domains are predicted to affect the cleavage at the GPS site [84], which may be another way of influencing the fertilization process. The presence of a REJ domain in combination with a GPS motif and possible cation channel function make Esi0130_0068 an appealing candidate for in situ evaluation.

Conclusions
This study focused on genes for which there is evidence that expression is limited only to male gametes of Ectocarpus and possibly subjected to adaptive evolution. By extrapolation, the observed positive selection may pinpoint the genes that are directly involved in male reproduction, which would be an important step towards understanding the molecular basis of gamete interaction during reproduction in Ectocarpus. In particular, one male gamete specific gene (Esi0130_0068) appears to be a good candidate due to signatures of positive selection and its similarity to the sperm-egg recognition protein in sea urchin. However, nucleotide sequences used in this study represent only a partial coding sequence of selected genes in a limited number of strains. Future work would require a larger sample size and complete gene sequences for a better estimation of evolution over time and forces shaping the divergence in the sexrelated genes in Ectocarpus.