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.

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-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-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/d) 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 sperm-egg binding were isolated and partially characterized in Fucus [30-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].

Table 1

Strains used in the study

Lineage	Strain code or species	Short name	Strain no.		Origin
			Eca	CCAPb
4	Ec NZKU 1–3 m	Nze	-	CCAP 1310/56	Kaikoura, New Zealand
5a	Ec Cph 40–11 m	Cph	-	CCAP 1310/100	Copenhagen, Denmark
2d	Ec Tam 2b m	Tam	-	CCAP 1310/122	Tampa, Florida
5a	Ec Pen 2a m	Pen	-	CCAP 1310/111	Penikese Island, Massachusetts
2d	Ec PAr 18a m	Par	-	CCAP 1310/108	Port Aransas, Texas
5b	E. fasciculatus	Fas	Ec185	-	Perharidy, France
1a	E. siliculosus	Nap	Ec400	CCAP 1310/329	Naples, Italy
2c	E. crouaniorum	Cro	Ec477	-	Perharidy, France
1c	Ec 32 (genome strain)	Per	Ec32	CCAP1310/4	San Juan de Marcona, Peru
-	Scytosiphon “lomentaria”	-	ODC1349c		Brittany, France
-	Scytosiphon “lomentaria”	-	LT0114c		Australia

aMacroalgal Culture Collection in Roscoff, France; bCulture Collection of Algae and Protozoa, Oban, Scotland; cdeposited at Ghent University, Belgium

Table 2

Ectocarpus cross-strain fertility

		Male
		Nap1a	Per1c	Cro2c	Tam2d	PAr2d	NZ4	Cph5a	Pen5a	Fas5b
Female	Nap, 1a	F	post	post	z	z	post	z	pre	post
	Per, 1c	post	F
	Cro, 2c	post		F	z	z			pre	post
	Tam, 2d	z		z	F	z			prea
	PAr, 2d	pre		z	z	F			prea
	NZ, 4	post					F
	Cph, 5a	z
	Pen, 5a	pre		pre	prea	prea			F
	Fas, 5b	pre		pre						F

Data on cross-fertility summarized from [40, 41, 85–87]. Clade numbers corresponding to the phylogenetic position in are given with the strain names. F – full interfertility, z – zygote formation, no data on growth, pre – prezygotic barriers, no cell fusion, post – hybrids with reduced growth or non-functional reproductive structures. aNo 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

Fig. 1

Species tree and species sexual compatibility. Maximum likelihood tree based on the cox3 gene from nine Ectocarpus strains and two Scytosiphon lomentaria specimens. Strains used in this study are marked with blue bars. ML bootstrap (left) and Bayesian posterior probability (right) values are shown. Red lines define species with prezygotic barriers to fertilization; dotted lines describe species with uncertain gamete compatibility; grey shading indicates pairwise d /d  > 1 (PAML) (see Additional file 4)

Methods

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].

DNA extraction

DNA was extracted from 20–50 mg of culture material (Ectocarpus) or silica dried specimens (Scytosiphon). Samples were ground in liquid nitrogen and extracted using CTAB buffer (100 mM Tris–HCl, pH 8.0, 1.4 M NaCl, 20 mM EDTA, 1 % w/v CTAB, 1 % w/v PVP, 0,2 mg/ml proteinase K), followed by two extractions with chloroform:isoamylalcohol (24:1). DNA was precipitated with 80 % isopropanol and cleaned with Qiagen MagAttract Suspension G according to the manufacturer’s instructions.

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 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/d 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/d (ω) 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.

Results

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 cell-cell 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) and synonymous (d) 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/d) 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/d ratio >0.5, which could imply adaptive evolution (Table 3, d/d >0.5 marked in bold). It has been shown that the gene length may bias the estimations of d/d particularly for short genes with low divergence which are then preferably found to be under positive selection [71]. However, higher d/d 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.

Table 3

Parameter estimates under pairwise sequence analysis in CODEML, PAML 4 using Ectocarpus siliculosus lineage 1a vs Ectocarpus sp. lineage 1c

Gene	Typea	Alignment length (AA)	ω = d N /d S	d N	d S	κ
Esi0018_0186	GOI	513	0.154	0.016	0.107	5.66
Esi0026_0077	GOI	310	0.596	0.034	0.057	2.45
Esi0030_0082	GOI	975	0.435	0.041	0.094	4.33
Esi0033_0091	GOI	184	0.106	0.004	0.040	9.69
Esi0043_0064	GOI	679	0.387	0.029	0.074	2.89
Esi0101_0018	GOI	727	0.025	0.005	0.190	2.82
Esi0116_0079	GOI	375	0.049	0.006	0.123	7.02
Esi0130_0068	GOI	1349	0.529	0.022	0.042	2.09
Esi0132_0081	GOI	1464	0.577	0.048	0.083	2.64
Esi0146_0068	GOI	1537	0.705	0.041	0.058	2.86
Esi0180_0035	GOI	410	0.499	0.026	0.052	2.79
Esi0183_0054	GOI	428	0.338	0.024	0.071	2.76
Esi0186_0062	GOI	687	0.255	0.025	0.098	4.69
Esi0188_0041	GOI	429	0.137	0.011	0.081	2.34
Esi0660_0004	GOI	745	0.533	0.038	0.071	2.05
Esi0008_0135	HKG	533	0.214	0.021	0.098	2.565
Esi0010_0097	HKG	316	0.136	0.0099	0.0733	4.154
Esi0010_0133	HKG	363	0.157	0.0122	0.0777	3.494
Esi0021_0112	HKG	228	0.192	0.0125	0.065	1.986
Esi0044_0085	HKG	319	0.035	0.0037	0.1052	5.886
Esi0069_0059	HKG	609	0.108	0.0106	0.0981	3.173
Esi0116_0065	HKG	92	0	0	0.2612	2.829
Esi0138_0009	HKG	642	0.002	0.0006	0.3315	1.414
Esi0159_0021	HKG	375	0.079	0.0065	0.0814	5.385
Esi0197_0055	HKG	89	0.479	0.0231	0.0483	0.645
Esi0289_0026	HKG	372	0	0	0.0777	2.165
Esi0387_0021	HKG	435	0.016	0.0017	0.1097	1.215

aGOI - gene of interest; HKG - housekeeping gene

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).

Table 4

Positively selected sites identified by the site-prediction methods in PAML 4 and HyPhy (DATAMONKEY)

		Pr > 90 %		p value <0,1
Gene	ω = d N /d S a	M1a-M2a	M7-M8	MEME
Esi0130_0068	0.7	155	155	224, 230, 303, 626, 819, 820, 823
Esi0101_0018	0.1	none	none	none
Esi0289_0026	0	none	none	none
Esi0069_0059	0	none	none	none
Esi0138_0009	0	none	none	none

aEstimate of d /d assuming no rate heterogeneity (model M0, CODEML, PAML4)

Table 5

Likelihood ratio statistics (2delta L) for Esi0130_0068

Comparison	2delL	df	Chi-squared 5 %	Chi-squared 1 %
M0 (one ratio) vs. M3 (discrete)	15.229618	4	3753	13,28
M1a (nearly neutral) vs. M2a (positive selection)	6.603632	2	5,99	9,21
M7 (beta) vs. M8 (beta&w)	6.723490	2	5,99	9,21

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+, Ca2+ channel activated by Ca2+) 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/d >1 (Fig. 2). This result is consistent with an evolutionary history characterized by frequent episodes of positive selection. All three models (M3, M2a, M8) suggest ~7 % of sites under positive selection with ω2 = 4.07 (Table 6). The codons inferred to be under positive selection by PAML 4 with posterior probability >90 % lie within the REJ domain in the extracellular region of the protein. Sites identified by HyPhy are adjacent to the N-terminal site of the REJ domain (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].

Fig. 2

Posterior means of ω and probabilities for site classes in Esi0130_0068. a Posterior probabilities for site classes in Esi0130_0068 as calculated by M8 (beta&ω) model in PAML. The ω ratios are 0.00000, 0.00005, 0.00618, 0.12717, 0.64633, 0.96054, 0.99825, 0.99997, 1.00000, 1.00000 and 3.55267. Each of the first 10 categories has a proportion of 0.08717 and the last category has a proportion of 0.12826. Categories are grouped by ω < 0.65, 1 ≥ ω > 0.96, ω = 3.55. b Posterior means of ω, calculated as the average over 11 site classes. The amino acid sequence represents the genome strain (Ec32), the site numbers correspond to the position in the coding sequence

Table 6

Parameter estimates and log-likelihood values under models of variable w ratios among sites for Esi0130_0068

Model	p	Parameters	l	d N /d S	Positively selected sites BEB (Pr > 90 %)
M0: one ratio	1	ω = 0.65860	–1812.54	=ω	None
M1a: nearly neutral	2	p0 = 0.45184, ω0 = 0.08242, p1 = 0.54816, ω1 = 1.00000	–1808.23	0.5854	Not allowed
M2a: positive selection	4	p0 = 0.93134, ω0 = 0.45793, p1 = 0.00000, ω1 = 1.00000, p2 = 0.06866, ω2 = 4.06616	–1804.93	0.7057	155
M3: discrete	5	p0 = 0.88395, ω0 = 0.45793, p1 = 0.04738, ω1 = 0.45796, p2 = 0.06866, ω2 = 4.06615	–1804.93	0.7057	155, 626, 821
M7: beta	2	p = 0.02684 q = 0.01590	–1808.29	0.6092	Not allowed
M8: beta& ω	4	p0 = 0.93179 p = 83.98499q = 99.00000, (p1 = 0.06821)ω = 4.07807	–1804.93	0.7058	155

Fig. 3

Topology of Esi0130_0068 as predicted by Phyre2 and Blast2GO. The extracellular and cytoplasmic sides of the membrane are labeled and the beginning and end of each transmembrane helix are indicated by the residue index. Positively selected sites identified by PAML 4 are marked with orange arrows, sites identified by HyPhy are marked in blue

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–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–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.

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/d 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/d 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 pre-zygotic 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 sex-related genes in Ectocarpus.