Additive transcriptomic variation associated with
reproductive traits suggest local adaptation in a recently
settled population of the Pacific oyster, Crassostrea gigas

Rossana Sussarellu, Arnaud Huvet, Sylvie Lapègue, Virgile Quillen, Christophe Lelong, Florence Cornette, Lasse Fast Jensen, Nicolas Bierne, Pierre Boudry

► To cite this version:

Rossana Sussarellu, Arnaud Huvet, Sylvie Lapègue, Virgile Quillen, Christophe Lelong, et al. Additive transcriptomic variation associated with reproductive traits suggest local adaptation in a recently settled population of the Pacific oyster, Crassostrea gigas. BMC Genomics, 2015, 16, pp.808. 10.1186/s12864-015-1972-8 . hal-01258474

HAL Id: hal-01258474

https://hal.sorbonne-universite.fr/hal-01258474

Submitted on 19 Jan 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Additive transcriptomic variation associated with reproductive traits suggest local adaptation in a recently settled population of the Pacific oyster, *Crassostrea gigas*

Rossana Sussarellu<sup>1,7*</sup>, Arnaud Huvet<sup>1</sup>, Sylvie Lapègue<sup>2</sup>, Virgile Quillen<sup>1</sup>, Christophe Lelong<sup>3</sup>, Florence Cornette<sup>2</sup>, Lasse Fast Jensen<sup>4</sup>, Nicolas Bierne<sup>5,6</sup> and Pierre Boudry<sup>1</sup>

**Abstract**

**Background:** Originating from Northeast Asia, the Pacific oyster *Crassostrea gigas* has been introduced into a large number of countries for aquaculture purpose. Following introduction, the Pacific oyster has turned into an invasive species in an increasing number of coastal areas, notably recently in Northern Europe.

**Methods:** To explore potential adaptation of reproductive traits in populations with different histories, we set up a common garden experiment based on the comparison of progenies from two populations of Pacific oyster sampled in France and Denmark and their hybrids. Sex ratio, condition index and microarray gene expression in gonads, were analyzed in each progeny (n = 60).

**Results:** A female-biased sex-ratio and a higher condition index were observed in the Danish progeny, possibly reflecting an evolutionary reproductive strategy to increase the potential success of natural recruitment in recently settled population. Using multifarious statistical approaches and accounting for sex differences we identified several transcripts differentially expressed between the Danish and French progenies, for which additive genetic basis is suspected (showing intermediate expression levels in hybrids, and therefore additivity). Candidate transcripts included mRNA coding for sperm quality and insulin metabolism, known to be implicated in coordinated control and success of reproduction.

**Conclusions:** Observed differences suggest that adaptation of invasive populations might have occurred during expansion acting on reproductive traits, and in particular on a female-biased sex-ratio, gamete quality and fertility.

**Keywords:** *Crassostrea gigas*, Adaptation, Invasiveness, Fertility, Sex-ratio, Transcriptome, Microarray

**Background**

Invasive species are privileged models to analyse the evolution and adaptation of life history traits in new environments [1, 2] While an increasing number of studies have documented adaptation during invasion in terrestrial species [3, 4], few studies have been conducted in marine species to date [5–7]. Many marine species have a benthic sessile adult phase and disperse *via* a planktonic

larval stage. In these organisms, high gene flow between populations can be ensured by high larval dispersal, high fecundity and huge population sizes [8, 9]. Nevertheless, paradoxical levels of genetic differentiation has been observed in marine populations [10, 11] originating from either local adaptation through natural selection or purely random processes such as genetic drift. In this context, marine invasive species offer the opportunity to study evolution and adaptation of life history traits upon introduction into new habitats.

Originating from Northeast Asia, the Pacific oyster *Crassostrea gigas* has been introduced and translocated worldwide, mainly for aquaculture purposes. Self-

<footnote>* Correspondence: Rossana.sussarellu@ifremer.fr</footnote>

<footnote n="1">¹Ifremer, Laboratoire des Sciences de l'Environnement Marin UMR 6539 (UBO/CNRS/IRD/Ifremer), Plouzané, France</footnote>

<footnote n="7">⁷Present address: Ifremer, Laboratoire d'Ecotoxicologie, Nantes, France
Full list of author information is available at the end of the article</footnote>
sustaining populations have today been recorded in at least 17 of them [12]. Although highly variable, the invasiveness of *C. gigas* has been demonstrated in several countries and this species is considered as pest or noxious in an increasing number of coastal areas [13]. In European waters, *C. gigas* is cultured from Norway to Portugal as well as in the Mediterranean Sea. This species rapidly settled along the Atlantic coasts of France following its massive introduction at the end of the 1960s [14]. More recently, feral populations of *C. gigas* have been reported in northern Europe [15, 16] as far north as Sweden, where dense populations of settled oysters can now be observed in several shallow water sites [17]. This expansion of Pacific oyster in the North Sea occurred much later after their first introduction than along the French Atlantic coast, suggesting that increase in population size may have been retarded by false or irregular recruitment depending on water temperature [15–18]. Irregular recruitment has been reported to coincide with above-average summer temperatures in late summer [19]. Indeed, *C. gigas* now reproduce and settle in Scandinavian waters as far as 60° N, while in the beginning of the 70’s, attempts of culture of this species failed at higher latitudes [18]. In this context, the success of the species and especially its northward expansion might be explained by climate change [20, 21] but also local adaptation as well as phenotypic plasticity, or both [22].

History of initial introductions and later transfers, natural connectivity resulting from natural dispersal and their resulting genetic structure are needed to assess the potential local adaptation of newly settled populations (for review see [6]). Moehler et al. [23] proposed that the genetic population structure of *C. gigas* in the Wadden sea could have been shaped by aquaculture practices. The presence of two separate genetic groups, one in the southern and one in the northern Wadden sea was suggested to be the result of two independent invasions. Rohfritsch et al. [24] investigated the possible effect of adaptation with a genome-scan approach. As previously reported [25, 26], no significant genetic structure was noted when comparing the population from Japan to European populations sampled along the French coasts and between populations of southern Europe. However, a significant genetic structure was observed in Europe among northern populations (located in Germany, Denmark and Sweden) and between northern and southern populations, together with lower genetic diversity in the north. Recently, Lallias et al. [27] highlighted how the number of oyster introduction events, aquaculture practices, genetic bottlenecks followed by genetic drift and natural dispersal could shape the genetic diversity and structure of introduced populations. Furthermore the analysis of $F_{ST}$ outliers revealed 6 candidate loci for adaptation which could either reflect (i) parallel adaptation to similar environmental pressures (fjord-like

environment) within each of the two groups or (ii) a footprint of a secondary introduction of an alternative genomic background, maintained by multifarious isolation factors [24]. This suggests that adaptation or reshuffling of pre-existing genetic backgrounds could have occurred during the invasion. In general, however, genome scans with an insufficient marker density have proved unsatisfactory to identify adaptation during marine invasion [24, 28], and the study of phenotypic differentiation has been proposed as a more promising approach [29, 30].

Colonization of new environments may promote rapid population divergence as a by-product of local adaptation to differential selective pressures [31]. Local adaptation in populations can be examined by analysing phenotypic traits that are likely to be differentially selected in the wild, but genetic bases of trait variation are necessary to allow response to natural selection. In most cases, the heritability estimates of these traits are low and their study requires the analysis under common conditions of progenies of wild genitors (rather than their direct study in the field), allowing disentangling environmental and genetic effects. In this context, gene expression studies can reveal adaptive mechanism, and the genetic basis of traits affecting fitness [32]. Furthermore, comparing the extent of quantitative genetic differences of phenotypic traits among populations could be assessed by the degree of differentiation in quantitative traits $Q_{ST}$ [33] that has been widely used to assess the relative contributions of selection to phenotypic traits divergence [34, 35] and gene transcription profiles [36]. Transcriptomic scan ($eQ_{ST}$) can indeed reveal adaptation, as outlier genes, showing the highest levels of differentiation between populations, will represent those that are most likely to evolve under directional selection [37]. Conversely, $P_{ST}$ (differentiation phenotypic traits) and $eP_{ST}$ scan can reveal selective selection [38]. In a context of climate change facilitating the establishment of invasive species, transcriptomic studies are promising in understanding basis of variation in phenotypic plasticity [39].

In this study we set out to investigate the evolutionary processes associated with the establishment and invasion into new environments of Pacific oyster population. Specifically, we raised progenies originating from Northern and Southern Europe along with their hybrids under common environmental conditions. A number of phenotypic traits, including gene expression patterns were analysed to search for evidence of differentiation in reproductive traits that might result from local adaptation.

## Methods

### Biological material

Wild oysters were collected in October 2009 at the Ile de Ré (France) and Limfjord (Denmark), transported to the Ifremer's facilities in La Tremblade (Charente-Maritime, France) and kept in common controlled conditions. The
sampling site in France is located 30 km from the location "MAR" that was studied by Rohfritsch et al. [24], while the sampling site in Limfjorden is identical to the location "LIM" [24]. Maturation was carried in shared conditions suitable for germ cell maturation according to Fabioux et al. [40] in Ifremer's facilities in La Tremblade. In May 2010, crosses within and between populations were produced, using 40 parental individuals (20 females and 20 males) from each of the two populations. For the hybrid progeny, the same individuals were employed as parental oysters as for the within site ones, crossing French males with Danish females. Gametes were collected as described in Huvet et al. [41]. Briefly, sperm and oocytes were collected in seawater by stripping the gonads. For the two groups of females, oocytes were counted and equally pooled. Oocytes were then distributed in fertilization beakers and then fertilized separately by each male at a ratio of 100 spermatozoa/oocyte.

Progenies were reared under standard and common hatchery and nursery conditions. In October 2010, when about 10 mm large, juveniles were transferred to the field in the Marennes-Oléron basin until sampling in June 2011. For each progeny ("French", "Danish" and "hybrid"), 100 oysters were randomly collected (in 4 replicated baskets) and weighed (total and wet flesh weights). The condition index, as specified by the French norm AFNOR was calculated using the following equation: (wet flesh weight/total weight) × 100 [42]. Gonads were immediately dissected from each oyster, a transversal section of the gonadic area was made for histological examination and the rest of the gonad was frozen in liquid nitrogen. Gonad tissues were crushed to a fine powder at -196 °C with an oscillating mill mixer and stored in liquid nitrogen until RNA extraction and biochemical analyses.

## Histological analysis

For each sample, a 3 mm cross-section of the visceral mass was excised in front of the pericardic region and immediately fixed in modified Davidson's solution [43] at 4 °C for 48 h. Sex-ratio (i.e. the ratio males/males + females) and gonad developmental stage were determined by histological methods (see details in Fabioux et al. [40]) according to the reproductive scale of Steele and Mulcahy [44]. Sixty individuals for each progeny were chosen (for a total of 180) on the basis of their maturation stage (stage 3 of gametogenesis, ripeness) and sex for transcriptomic analysis.

## RNA extraction, amplification, labeling and microarray hybridization

Total RNA was isolated using 1.5 mL of Extract-all Reagent® (Eurobio AbCys, Courtaboeuf, France) per 50 mg of gonad powder, according to the manufacturer's

instructions. RNA concentrations were determined using a ND-1000 spectrophotometer (Thermo Scientific, Waltham MA, USA) at 260 nm, using the conversion factor 1 OD = 40 µg/mL RNA. RNA integrity was assessed on an Agilent bioanalyzer using RNA 6000 Nano kits® (Agilent Technologies, Santa Clara, CA, USA), according to manufacturer's instructions. Over the 180 samples extracted, the RNA Integrity Number (RIN), obtained by setting the threshold "Unexpected Ribosomal Ratio" to 2 in the Agilent 2100 Bioanalyzer software considering the co-migration of 28S and 18S rRNA fragments in bivalves [45], varied from 7.9 to 9.6 (mean = 9.1 ± 0.3). For microarray hybridizations, 200 ng of total RNA were indirectly labelled with Cy3 using the Low Input Quick Amp Labeling kit One-Color® (Agilent Technologies), according to the manufacturer's instructions. Qiagen RNeasy® (Qiagen, Venlo, Netherlands), mini spin columns were used for purifying amplified RNA samples. After purification, RNA amplification and dye incorporation rates were verified using a ND-1000 spectrophotometer (Thermo Scientific) and shown to lie between 100 and 200 ng/ µL (RNA concentration) and between 1 and 5 pmol/µL RNA (dye incorporation). Hybridization was performed using the Agilent Gene expression hybridization kit® (5188-5242), as described by the manufacturer, with 1.65 µg of labelled RNA at 65 °C for 16 h. The employed arrays were Agilent 60-mer 4x44 K custom microarrays, containing 31,918 *C. gigas* transcripts, designed by Dheilly et al. [45]. Samples were randomly hybridized onto 48 different slides, which were subsequently treated with Gene expression wash buffer solution® (5188-5327; Agilent Technologies), Stabilization and Drying solutions® (5185-5979; Agilent Technologies). Slides were scanned on an Agilent Technologies G2565AA Microarray Scanner system® at 5 µm resolution, using default parameters. Features were extracted using the Agilent Feature Extraction software 6.1 (Agilent Technologies).

## Pre-processing and microarray data analysis

Microarray data were processed and analysed using the language R/BioConductor (R Development Core Team 2008 [46]). Quantile normalization was performed on background-corrected features with the limma package [47]. Arrays having more than 3 % of not uniform features were eliminated for subsequent analysis. Filtering step was performed according to Agilent Feature Extraction software results on spot quality and spot intensity reliability. Negative filtered features were excluded from subsequent analyses. Missing values were imputed with impute package [48]. Raw and normalized hybridization values are deposited in the gene expression omnibus (GEO) repository with the accession number GSE66103.
We first applied a redundancy analysis (RDA), a con-
strained ordination method implemented in the Vegan
package [49], to obtain a global view of the extent to
which the explanatory variables “sex” (male, female) and
“progeny” (Denmark, France, hybrid) influenced the ex-
pression levels. Three RDA were performed. The first
one was a full model with sex and progeny as explana-
tory variables (Sex + Progeny). Then, two partial models
were applied in order to partitioning explicable variance
on sex and progeny (respectively: Sex + Condition(Pro-
geny) and Progeny + Condition(Sex)). Contributions of
genes to RDA axis 2 were retrieved in order to identify
individual transcripts that contribute the most. Secondly,
differentially expressed transcripts between French and
Danish progenies were also identified with analysis of
variance. Fixed factors for the two-way ANOVA were
sex and progeny. Multiple testing *p*-values were adjusted
using Benjamini-Hochberg (BH) correction. Finally, a
scan of phenotypic differentiation of expression levels
was performed (eP<sub>ST</sub> scan). Additivity on all mRNA ex-
pression levels was specifically evaluated by assessing if
hybrid expression was not different from the theoretical
mid-parent level (tested by t-test at *p*-value < 0.01) as
proposed by Hedgecock et al. [50]. Then eP<sub>ST</sub> estimates
were calculated using the equation: $\sigma_{GB}^2/(\sigma_{GB}^2 + 2\sigma_{GW}^2)$,
where $\sigma_{GB}^2$ and $2\sigma_{GW}^2$ represent among- and within-
population components of the genetic variance for
quantitative traits respectively [33], and tested by per-
mutations (N = 5000). Transcripts for which the eP<sub>ST</sub>
values exceeded the 0.999 quantile of the permuta-
tion distribution were retained as outliers. The script
used for RDA, ANOVA and eP<sub>ST</sub> calculations is in
Additional file 1.

Transcripts putative annotations were identified using ngKlast blast (Korilog) against a protein data base (E-value 1.0 × 10⁻⁵) obtained from the *C. gigas* sequenced genome and transcriptome deposited on Genbank [51]. GO terms were obtained using ngKlast against the Swissprot database (E-value 1.0 × 10⁻⁵). GO terms enrichment analysis were performed using the Fisher's Exact test on Blast2Go [52]. Hierarchical clustering was performed using Ward method and 1-correlation as dissimilarity matrix.

Statistical analyses of sex-ratio and condition index

Data for sex-ratio, condition index and biochemical ana-
lyses were processed and analysed using the language R/
BioConductor [46]. Comparisons of sex distributions be-
tween progenies were made using Chi-square tests (Χ²).
The condition index, was tested using a two-way
ANOVA with sex and progeny as fixed factors, post-hoc
test (LSD test) was used to determine which groups
were significantly different. Normality was checked using

Shapiro-Wilk test and homogeneity of variances matri-
ces with Bartlett test.

Results

Sex-ratio and condition index

Sex-ratio differed significantly (p <0.01) between French and Danish progenies with a *female*-biased sex-ratio observed in the Danish progeny (40 %) compared to the French progeny (73 %). Hybrids presented an intermediate but not significantly different sex-ratio (58 %) from French and Danish progenies. The two-way ANOVA for the condition index (CI) showed a progeny effect (p <0.003). The condition index of the 3 progenies clustered into two significantly different groups, hybrid and Danish progenies showing significant higher values (17.5 %) than French one (16.4 %, Fig. 1). Histological examination showed that all individuals of the three progenies were in stage 3, corresponding to full ripeness.

Gene expression analysis

RDA analysis on all mRNA and individuals clearly dis-
criminated males from females on its first component
RDA1 (Fig. 2). From the two following partial RDA, sex
accounted for 95.5 % and progeny for 4.3 % of the ex-
plained variance, respectively. The two explanatory vari-
able vectors were orthogonal on the RDA plan, and
hybrids were intermediate between French and Danish
progenies on second component RDA2. Scores of con-
tribution to second component RDA2 allowed us to ob-
tain a first list of 64 genes whose contribution accounted
in the 0.1 and 99.9 percentile of the distribution. Gen-
bank access, putative annotations and E values of these

Condition index

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Boxplot of the condition index. Condition index by progeny and sex, small letters indicate progeny effects found by the LSM post-hoc test. DAN = Danish, FRA = French, HYB = Hybrids, F = Females, M = Males (N = 27-31)</figcaption></figure>
<figure><img src="image_4.png" /></figure>

transcripts are presented in Table 1 and their mean expression *per* progeny in the heatmap in Fig. 3. The enriched GO biological processes attributed to those genes are presented in Additional file 2.

The two-way ANOVA performed on French and Danish progenies confirmed the strong sex effect on gene expression already identified with the RDA, with 11,472 differentially expressed transcripts between sexes, and a progeny effect with 80 transcripts differentially expressed (BH correction, cut-off *p*-value of $10^{-4}$ presented in Additional file 3). Twenty of these 80 genes were in common with the candidates identified with the contribution to RDA2.

Due to the strong differentiation in gene expression between males and females, the $eP_{ST}$ were estimated separately by sexes on all transcripts. Only 953 and 980 transcripts respectively (corresponding to 6 % of the total 31,918 *C. gigas* transcripts present on the microarray), showed an mRNA level intermediate in hybrids between the French and Danish progenies. Average $eP_{ST}$ values were 0.04 for both males and females. Transcripts considered as outliers had $eP_{ST}$ values higher than 0.15 for females and 0.17 for males. The number of outlier transcripts was 53 on males and 69 on females, and they are presented in Additional file 4. Ten transcripts were in common between all analyses, two of these having

known putative annotations (A-kinase anchor protein 7 isoform gamma and Galectin-4).

## Discussion

So far, colonization of *C. gigas* on European coasts has been investigated by ecological or marker-based population genetic approaches, but the variability of phenotypic traits potentially associated with this range expansion has only been recently investigated [7]. Despite strong difference in expression levels between sexes, as previously reported [53], we nonetheless detected a clear difference in expression levels between progenies of oysters originated from France and Denmark, and their hybrids. From the RDA approach on microarrays we found that 4.3 % of the explained variance is ascribed by the progeny in a mainly orthogonal way to sex differences, which means that despite the huge transcriptomic differences between sexes, the more modest differences between progenies proved to be mainly independent of the sex. Only characters being heritable could allow adaptation because significant genetic bases of trait variation (rather than only phenotypic plasticity) are necessary to allow response to eventual natural selection. In our study, the use of hybrid progeny allowed to identify transcripts presumed to show additive genetic variation. In
<table><caption>**Table 1** RDA sex significant transcripts 64 transcripts whose contribution accounted for the 0.1 and 99.9 percentile of the RDA2 axis distribution scores; Genbank accession number, putative annotations, E values, cluster and presence on these transcripts in the other analyses (a = ANOVA, m = eP<sub>ST</sub> males, f = eP<sub>ST</sub> females)</caption><thead><tr><th>Genbank</th><th>Hit definition</th><th>E value</th><th>Cluster</th></tr></thead><tbody><tr><td>CU998588</td><td></td><td></td><td>1</td></tr><tr><td>AM861480</td><td>hypothetical protein CGI_10017906</td><td>1.00E-17</td><td>1</td></tr><tr><td>FP007310</td><td></td><td></td><td>1</td></tr><tr><td>FP002060</td><td></td><td></td><td>1</td></tr><tr><td>FP007067</td><td></td><td></td><td>1</td></tr><tr><td>CU998537</td><td></td><td></td><td>1</td></tr><tr><td>FP001074</td><td></td><td></td><td>1</td></tr><tr><td>AM855692</td><td>hypothetical protein CGI_10011429</td><td>6.00E-53</td><td>1</td></tr><tr><td>AM866859</td><td>A-kinase anchor protein 7 isoform gamma</td><td>2.00E-70</td><td>1</td></tr><tr><td>CU995043</td><td>Cadherin-23</td><td>2.00E-25</td><td>1</td></tr><tr><td>AM863477</td><td>Barrier-to-autointegration factor</td><td>2.00E-08</td><td>1</td></tr><tr><td>FP005183</td><td></td><td></td><td>1</td></tr><tr><td>AM854599</td><td></td><td></td><td>1</td></tr><tr><td>CU999500</td><td></td><td></td><td>1</td></tr><tr><td>FP004770</td><td>Zinc finger protein 26</td><td>3.00E-12</td><td>1</td></tr><tr><td>FP000261</td><td>CD63 antigen</td><td>3.00E-106</td><td>1</td></tr><tr><td>FP011598</td><td></td><td></td><td>1</td></tr><tr><td>FP004638</td><td></td><td></td><td>1</td></tr><tr><td>AM859278</td><td></td><td></td><td>1</td></tr><tr><td>AM858298</td><td>RING-box protein 1</td><td>6.00E-42</td><td>1</td></tr><tr><td>CX069237</td><td>Actin</td><td>8.00E-38</td><td>1</td></tr><tr><td>AM856393</td><td>hypothetical protein CGI_10011429</td><td>1.00E-31</td><td>1</td></tr><tr><td>CU986854</td><td></td><td></td><td>1</td></tr><tr><td>AM856074</td><td>hypothetical protein CGI_10023286</td><td>4.00E-113</td><td>1</td></tr><tr><td>ES789933</td><td>Coatomer subunit gamma-2</td><td>3.00E-10</td><td>1</td></tr><tr><td>AM865168</td><td></td><td></td><td>1</td></tr><tr><td>AM858762</td><td>hypothetical protein CGI_10024669</td><td>6.00E-35</td><td>1</td></tr><tr><td>AM863236</td><td>Transaldolase</td><td>3.00E-86</td><td>1</td></tr><tr><td>FP003575</td><td></td><td></td><td>1</td></tr><tr><td>FP003052</td><td></td><td></td><td>1</td></tr><tr><td>AJ565528</td><td></td><td></td><td>1</td></tr><tr><td>CU999266</td><td>hypothetical protein CGI_10005871</td><td>1.00E-42</td><td>1</td></tr><tr><td>FP006374</td><td>Interferon-induced, double-stranded RNA-activated protein kinase</td><td>5.00E-10</td><td>2</td></tr><tr><td>CU682388</td><td>Ephrin type-A receptor 2</td><td>3.00E-25</td><td>2</td></tr><tr><td>FP000015</td><td>Glutathione reductase, mitochondrial</td><td>2.00E-12</td><td>2</td></tr><tr><td>AM854329</td><td></td><td></td><td>2</td></tr><tr><td>FP005519</td><td></td><td></td><td>2</td></tr><tr><td>CU685625</td><td>Galectin-4</td><td>2.00E-43</td><td>2</td></tr><tr><td>AM866437</td><td></td><td></td><td>2</td></tr><tr><td>BQ427216</td><td>tRNA 2'-phosphotransferase 1</td><td>1.00E-06</td><td>2</td></tr><tr><td>AM867714</td><td></td><td></td><td>2</td></tr><tr><td>AM867735</td><td>Ribosomal RNA-processing protein 8</td><td>1.00E-39</td><td>2</td></tr></tbody></table>
<table><caption>**Table 1** RDA sex significant transcripts 64 transcripts whose contribution accounted for the 0.1 and 99.9 percentile of the RDA2 axis distribution scores; Genbank accession number, putative annotations, E values, cluster and presence on these transcripts in the other analyses (a = ANOVA, m = eP<sub>ST</sub> males, f = eP<sub>ST</sub> females) (Continued)</caption><tr><td>AM858003</td><td>Mortality factor 4-like protein 1</td><td>6.00E-101</td><td>2</td></tr><tr><td>EE677761</td><td>L-rhamnose-binding lectin CSL3, partial</td><td>6.00E-46</td><td>2</td></tr><tr><td>CU681535</td><td></td><td></td><td>2</td></tr><tr><td>AM860599</td><td>hypothetical protein CGI_10018932</td><td>3.00E-40</td><td>2</td></tr><tr><td>CU984193</td><td>hypothetical protein CGI_10008458</td><td>7.00E-15</td><td>2</td></tr><tr><td>AM858265</td><td></td><td></td><td>2</td></tr><tr><td>AM860282</td><td></td><td></td><td>2</td></tr><tr><td>AM859564</td><td></td><td></td><td>2</td></tr><tr><td>DW713995</td><td>Dedicator of cytokinesis protein 6</td><td>2.00E-32</td><td>2</td></tr><tr><td>CB617483</td><td>hypothetical protein CGI_10016110</td><td>8.00E-11</td><td>2</td></tr><tr><td>CU988889</td><td>DNA polymerase epsilon subunit 4</td><td>1.00E-58</td><td>2</td></tr><tr><td>CX068783</td><td>hypothetical protein CGI_10013277</td><td>4.00E-26</td><td>2</td></tr><tr><td>FP001190</td><td></td><td></td><td>2</td></tr><tr><td>AM868733</td><td>Anaphase-promoting complex subunit 13</td><td>4.00E-39</td><td>2</td></tr><tr><td>ES789380</td><td>hypothetical protein CGI_10020002</td><td>8.00E-72</td><td>2</td></tr><tr><td>FP003774</td><td>Multimerin-1</td><td>3.00E-10</td><td>2</td></tr><tr><td>AM866600</td><td>DNA topoisomerase 1</td><td>3.00E-54</td><td>2</td></tr><tr><td>FP007621</td><td></td><td></td><td>2</td></tr><tr><td>FP002393</td><td></td><td></td><td>2</td></tr><tr><td>AM867802</td><td></td><td></td><td>2</td></tr><tr><td>AM860621</td><td></td><td></td><td>2</td></tr><tr><td>AM861439</td><td></td><td></td><td>2</td></tr></table>

the RDA the hybrid progeny appeared intermediate between French and Danish, and 50 % of the transcripts contributing the most to RDA2 axis proved to behave mainly additively. Our analyses allowed us to obtain a list of candidate genes, the expression of which might have reflected adaptation during invasion or were already genetically differentiated in the founding populations. By the way overall microarray features, only a very limited proportion of the total genes showed intermediate values in hybrids. The extensive non additivity of the transcriptome has been observed in Drosophila (Gibson et al. <2 %, [54]) and maize (Auger et al., ~30 % [55]) in contrast with the classical assumption in quantitative genetics of predominately additive genetics effects. Furthermore, a low proportion (2 %) of additive patterns of gene expression was previously observed by transcriptionic analysis in larvae of partially inbred Pacific oyster populations [50]. Finally, another study on oyster showed the non-additive nature of genetic variance for fitness-related traits and that the non-additive genetic component of yield is often the largest [56]. However, by focusing on a population likely to have adapted to a new environment during invasion (Denmark), we have access to the adaptation filter on phenotypic evolution and uncover additively

behaving traits (intermediate expression in hybrids) in the subset of transcripts that are the most differentially expressed. Interestingly, Wendling and Wegner [7] recently investigated the adaptive potential of North Sea *C. gigas* populations to local Vibrio spp., proposing that dominantly inherited resistance could facilitate fast adaptation.

In our study, ANOVA and eP<sub>ST</sub> analyses gave complementary results on phenotypic traits differentiation between the Danish and the French progenies. The majority of the mRNA differentially expressed in the ANOVA analysis was not in common with eP<sub>ST</sub> analysis. This could be explained by the fact that eP<sub>ST</sub> estimates were only performed on transcripts showing additive patterns and that most of the eP<sub>ST</sub> outliers did not have a high p-value in the ANOVA analysis. Roberge et al. [37], who find similar discordant results between ANOVA and Q<sub>ST</sub> estimates on salmon, suggested evolution as consequence of selection but without an additive genetic basis for these traits. eP<sub>ST</sub> estimates performed separately on males and females, demonstrates that considering sexes separately, when sex effect is so strong on gene expression patterns like in in gonads, could help in highlighting new candidates. Furthermore, this suggests that sex-dependent adaptation might be involved in the
<figure><img src="image_5.png" /></figure>

observed genetically-based phenotypic and transcriptional variations of reproductive traits, supported by an increasing evidence of a role for sex differentiated effects in the architecture of complex traits [57]. Sex interaction effects are common in model organisms for a wide range of traits, and can often explain a substantial part of the genetic basis of phenotypic variation [58].

In males, functional annotation of outlier eP<sub>ST</sub> estimates pinpointed mRNA encoding proteins involved in

sperm motility. The A-kinase anchor protein 7 isoform gamma [Genbank:AM866859], identified as outlier eP<sub>ST</sub> both in males and females and more expressed in Danish progeny, is a A-kinase anchoring protein (AKAPs) promoting the selective sequestration of intracellular cAMP-dependent protein kinase (PKA) involved in epithelial sodium channel regulation [59–61]. This protein has been suggested to have additional and perhaps unique function in spermatozoa as a scaffolding protein for the Rho-
GTPase pathway that regulates sperm motility. Furthermore, Kington et al. [62] found high expression levels of proteins having an important role in sperm motility of *C. gigas* involved in Rho signalling pathways. Sperm motility is commonly used as criteria of male gamete quality linked with the fertilization success and therefore potentially favouring colonization of new habitats. Furthermore, Filamin-A [Genbank:CU686267], more expressed in Danish progeny, is an Actin-binding protein that participates in the anchoring of membrane proteins for the actin cytoskeleton and is involved in sperm morphogenesis in mammals [63]. Finally, as cells protection against oxidative injury in marine invertebrates [64], Alternative oxidase (AOX, [Genbank:BQ426710]), more expressed in Danish progeny, may have a role in oxidative protection in gonads, and potentially in gamete quality. In females, the Molluscan insulin-related peptide 5 [Genbank:CU987248] was more expressed in Danish progeny. In mammals, the role of insulin pathway in fertility is well known [65] and in fish the insulin pathway has been positively associated to gamete quality [66]. In *C. gigas*, an Insulin-related peptide receptor has previously been identified in oysters by Gricourt et al. [67] as well as several factors of the insulin signalling pathway. Jouaux et al. [68] found that insulin pathway elements can modulate germinal cells proliferation during food deprivation in the first stages of gametogenesis with expected consequences on fertility.

A comparison of divergence in neutral markers, F<sub>ST</sub>, to divergence in phenotypic traits, Q<sub>ST</sub>, has widely been used as a method to assess the relative strength of genetic drift and selection [34]. The Q<sub>ST</sub> levels typically exceed that observed in F<sub>ST</sub> suggesting an important role of natural selection on quantitative traits [34, 69]. In our study, we did not have access to Q<sub>ST</sub> values but we restricted P<sub>ST</sub> values to additively behaving transcript with intermediate expression levels in hybrid. We found that mean eP<sub>ST</sub> (0.04) was similar to the estimated F<sub>ST</sub> values [24] on the same populations, while the outliers eP<sub>ST</sub> had values greater than 0.15. Our results therefore suggest that diversifying selection has most probably acted on outlier gene expressions. Finally, gene expression can also be controlled by epigenetic mechanisms, potentially in a transgenerational manner, meaning in a way, by genetic mechanisms non-DNA dependents. There are yet only a few papers studying epigenetics in oysters [70, 71] and this is clearly an emerging topic. Moreover, few studies until now focused on epigenetic-mediated adaptation in invasive species [72].

## Sex-ratio and population expansion

In our study, we observed a female biased sex-ratio in the Danish progeny. Sex-ratio biased to female in

invasive species has also been observed in an estuarine shrimp, *Palaemon macrodactylus* [73] and was proposed as a good descriptor to detect invading populations in signal crayfish, *Pacifastacus leniusculus* [74]. Although sex determinism in *C. gigas*, an alternative and irregular protandrous hermaphrodite, is complex, it seems to be determined both by environmental and genetic factors [75]. In our study, environmental effects were minimized by rearing progenies in common garden conditions, indicating that the observed differences in sex-ratio are genetically based. However interactions with environmental conditions have been suggested to result from regulatory pathways involved in sex determination [76]. Two genetic models have been proposed for sex determinism in *C. gigas*. The first suggests the presence of 2-genotypes, a dominant male M allele and a protandric recessive F allele [77]. The second propose a 3-genotypes model, FF for true female oysters, MM for true male oysters and FM for individuals that may mature as females or males [75]. Furthermore, an energy-mediated sex determinism was proposed [78], in which sex-ratio could be a possible way to select faster-growing populations when more females are produced in the first year of the life. The greater condition index observed in the Danish progeny could translate a greater reproductive effort in this new-expanded Danish population. This index is strictly correlated to ripeness and gonadal occupation production during gametogenesis in oysters [79]. In species like *C. gigas*, having an “r” demographical strategy, characterized by high fecundity, reproductive success is greatly dependent on the quantity of gametes produced, especially oocytes, as well as their quality. Furthermore, Cardoso et al. [80] found that in Northern European locations, oysters produce smaller eggs in larger quantities, suggesting an increasing reproductive output. The authors proposed that, since smaller oocytes are thought to have a longer development time, the environmental conditions along the Northern European coasts may result in increased larval dispersal and possibly in further population expansion. In this context, a greater reproductive effort together with a female-biased sex-ratio in *C. gigas* could favour a rapid colonization of new habitats.

## Conclusions

To conclude, genetic differentiation previously reported between *C. gigas* Southern Europe populations and population north to the Wadden Sea is corroborated by phenotypic differentiation based on transcriptomic data, biased sex-ratio and condition index. Overall, our results suggest a population expansion strategy of the studied Pacific oyster Danish population, potentially relying on a females-biased sex-ratio, a greater reproduction effort and gamete quality, noticeable on molecular signatures.
## Data availability

Microarray data are deposited in the gene expression omnibus (GEO) repository with the accession number GSE66103 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66103).

## Ethics statement

All the experiments were conducted according to the regulations of local and central government, and the study protocol was conducted in accordance with institutional guidelines

## Additional files

**Additional file 1: Script used for RDA, ANOVA and ePST calculations.**
(R 10 kb)

**Additional file 2: GO terms of RDA sex significant transcripts.** GO associated terms for transcripts whose contribution accounted for the 0.1 and 99.9 percentile of the RDA2 axis distribution scores. (XLS 179 kb)

**Additional file 3: ANOVA sex significant transcripts.** 80 mRNA differentially expressed in the two-way ANOVA between French and Danish progeny, Genbank accession number, putative annotations E values and GO associated terms. (XLS 132 kb)

**Additional file 4: ePST outliers for male and females oysters.** 53 and 69 outliers ePST mRNAs for male and female oysters respectively, Genbank accession number, putative annotations, E-values, fold changes and GO associated terms. (XLS 160 KB)

## Competing interests

The authors declare that they have no competing interests.

## Authors' contributions

RS VQ and CL carried out the molecular lab work; RS AH VQ and NB participated in data analysis; RS carried out the statistical analyses; AH SL VQ and FC participate in animal sampling; RS AH NB SL LFJ and PB participated in the design of the study; RS AH NB SL and PB conceived of the study, coordinated the study and helped draft the manuscript. All authors gave final approval for publication.

## Acknowledgements

The authors are grateful to P. Favrel for his support during the course of this work and as the supervisor of the "GametoGenes" project. We thank all the staff of the LGPMM Ifremer La Tremblade, particularly L. Degremont, conditioning oysters and producing progenies, and J-B. Lamy for revising R code. This study was conducted as part of the ANR "Hi-Flo" (contract n° ANR-08-BLAN-0334-01) and ANR "GametoGenes" (contract n° ANR-08-GENM-041). LFJ was supported by a grant (grant number 2009-7.40.01/56103-0003) from the Heritage Agency of Denmark.

## Author details

¹Ifremer, Laboratoire des Sciences de l'Environnement Marin UMR 6539 (UBO/CNRS/IRD/Ifremer), Plouzané, France. ²Ifremer, SG2M-LGPMM, Laboratoire de Génétique et Pathologie des Mollusques Marins, 17390 La Tremblade, France. ³UNICAEN, UMR BOREA MNHN, UPMC, UNICAEN, CNRS-7208, IRD207, F-14032 Caen, France. ⁴Fisheries and Maritime Museum, DK-6710 Esbjerg V, Denmark. ⁵Université Montpellier 2, Montpellier, France. ⁶CNRS - Institut des Sciences de l'Evolution, UMR5554, Station Méditerranéenne de l'Environnement Littoral, Sète, France. ⁷Present address: Ifremer, Laboratoire d'Ecotoxicologie, Nantes, France.

Received: 2 April 2015 Accepted: 3 October 2015

Published online: 19 October 2015

## References

1. Façon B, Genton BJ, Shykoff J, Jarne P, Estoup A, David P. A General Eco-Evolutionary Framework for Understanding Bioinvasions. 2006;21:130-135. doi:10.1016/j.tree.2005.10.012.

2. Sax DF, Stachowicz JJ, Brown JH, Bruno JF, Dawson MN, Gaines SD, et al. Ecological and Evolutionary Insights from Species Invasions. 2007;22:465-471. doi:10.1016/j.tree.2007.06.009

3. Prentis PJ, Wilson JRU, Dormontt EE, Richardson DM, Lowe AJ. Adaptive Evolution in Invasive Species. 2008;13:288-94. doi:10.1016/j.tplants.2008.03.004.

4. Bock DG, Caseys C, Cousens RD, Hahn MA, Heredia SM, Hübner S, et al. What We Still Don't Know about Invasion Genetics. 2014. doi:10.1111/mec.13032.

5. Lockwood BL, Somero GN. Transcriptomic Responses to Salinity Stress in Invasive and Native Blue Mussels (genus Mytilus). 2011;20:517-529. doi:10.1111/j.1365-294X.2010.04973.x.

6. Rius M, Turon X, Bernardi G, Volckaert FAM, Viard F. Marine Invasion Genetics: From Spatio-Temporal Patterns to Evolutionary Outcomes. 2014. doi:10.1007/s10530-014-0792-0.

7. Wendling CC, Wegner KM, Wendling CC. Adaptation to Enemy Shifts: Rapid Resistance Evolution to Local Vibrio Spp. in Invasive Pacific Oysters. 2015;282:2014.2244. doi:10.1098/rspb.2014.2244

8. Strathmann RR. Feeding and Nonfeeding Larval Development and Life-History Evolution in Marine Invertebrates. 1985;16:339-361. doi:10.2307/2097052.

9. Palumbi SR. Marine Speciation on a Small Planet. 1992;7:114-118. doi:10.1016/0169-5347(92)90144-Z.

10. Ward RD, Skibinski DOF, Woodwark M. Protein heterozygosity, protein structure, and taxonomic differentiation. In Hecht MK, Wallace B, Macintyre RJ (editors) Evolutionary Biology. Springer US; 1992:73-159. (Evolutionary Biology)

11. Small KS, Brudno M, Hill MM, Sidow A. Extreme Genomic Variation in a Natural Population. 2007;104:5698-5703. doi:10.1073/pnas.07008909104.

12. Ruesink JL, Lenihan HS, Trimble AC, Heiman KW, Micheli F, Byers JE, Kay MC. Introduction of Non-Native Oysters: Ecosystem Effects and Restoration Implications. 2005;36:643-689. doi:10.1146/annurev.ecolsys.36.102003.152638.

13. Orensanz JM, Schwindt E, Pastorino G, Bortolus A, Casas G, Darrigran G, et al. No Longer the Pristine Confines of the World Ocean: A Survey of Exotic Marine Species in the Southwestern Atlantic. 2002;4:115-143. doi:10.1023/A:1020596916153

14. Grizel H, Héral M: Introduction into France of the Japanese Oyster (*Crassostrea Gigas*). 1991;47:399-403. doi:10.1093/icesjms/47.3.399

15. Reise K. Pacific Oysters Invade Mussel Beds in the European Wadden Sea. 1998;28:167-175. doi:10.1007/BF03043147.

16. Nehls G, Diederich S, Thieltgès DW, Strasser M: Wadden Sea Mussel Beds Invaded by Oysters and Slipper Limpets: Competition or Climate Control? 2006;60:135-143. doi:http://dx.doi.org/10.1007/s10152-006-0032-9

17. Strand A, Blanda E, Bodvin T, Davids J, Jensen L, Hom-Hansen T, et al. Impact of an Icy Winter on the Pacific Oyster (*Crassostrea Gigas*) Thunberg, 1793 Populations in Scandinavia. 2012;7:433-440. doi:10.3391/ai.2012.7.3.014.

18. Wrange A-L, Valero J, Harkestad LS, Strand Ø, Lindegarth S, Christensen HT, et al. Massive Settlements of the Pacific Oyster, *Crassostrea Gigas*, in Scandinavia. 2010;12:1145-1152. doi:10.1007/s10530-009-9535-z.

19. Diederich S, Nehls G, Beusekom JEE van, Reise K: Introduced Pacific Oysters (*Crassostrea Gigas*) in the Northern Wadden Sea: Invasion Accelerated by Warm Summers? 2005;59:97-106. doi:10.1007/s10152-004-0195-1

20. Troost K. Causes and Effects of a Highly Successful Marine Invasion: Case-Study of the Introduced Pacific Oyster Crassostrea Gigas in Continental NW European Estuaries. 2010;64:145-165. doi:10.1016/j.seares.2010.02.004.

21. Jones MC, Dye SR, Pinnegar Jk, Warren R, Cheung WWL: Applying Distribution Model Projections for an Uncertain Future: The Case of the Pacific Oyster in UK Waters. 2013;23:710-722. doi:10.1002/aqc.2364.

22. Lande R. Evolution of Phenotypic Plasticity in Colonizing Species. 2015. doi:10.1111/mec.13037.

23. Moehler J, Wegner KM, Reise K, Jacobsen S. Invasion Genetics of Pacific Oyster *Crassostrea Gigas* Shaped by Aquaculture Stocking Practices. 2011;66:256-262. doi:10.1016/j.seares.2011.08.004

24. Rohfritsch A, Bierre N, Boudry P, Heurtebise S, Cornette F, Lapègue S: Population Genomics Shed Light on the Demographic and Adaptive Histories of European Invasion in the Pacific Oyster, *Crassostrea Gigas*. 2013;6:1064-1078. doi:10.1111/eva.12086.

25. Huet A, Boudry P, Ohresser M, Delsert C, Bonhomme F. Variable Microsatellites in the Pacific Oyster *Crassostrea Gigas* and Other Cupped Oyster Species. 2000;31:71-72. doi:10.1111/j.1365-2052.2000.579-5x.
26. Meistertzheim A-L, Arnaud-Haond S, Boudry P, Thébaut M-T. Genetic Structure of Wild European Populations of the Invasive Pacific Oyster *Crassostrea Gigas* due to Aquaculture Practices. 2013;160:453-463. doi:10.1007/s00227-012-2102-7

27. Lallias D, Boudry P, Batista FM, Beaumont A, King JW, Turner JR, Lapègue S: Invasion Genetics of the Pacific Oyster *Crassostrea Gigas* in the British Isles Inferred from Microsatellite and Mitochondrial Markers. 2015. doi:10.1007/s10530-015-0896-1

28. Riquet F, Daguin-Thiébaut C, Ballenghien M, Bierne N, Viard F. Contrasting Patterns of Genome-Wide Polymorphism in the Native and Invasive Range of the Marine Mollusc *Crepidula Fomicata*. 2013;22:1003–18. doi:10.1111/mec.12161.

29. Mayrose M, Kane NC, Mayrose I, Dlugosch KM, Riesenberg LH. Increased Growth in Sunflower Correlates with Reduced Defences and Altered Gene Expression in Response to Biotic and Abiotic Stress. 2011;20:4683–94. doi:10.1111/j.1365-294X.2011.05301.x

30. Hodgins KA, Lai Z, Nurkowski K, Huang J, Riesenberg LH. The Molecular Basis of Invasiveness: Differences in Gene Expression of Native and Introduced Common Ragweed (*Ambrosia Artemisiifolia*) in Stressful and Benign Environments. 2013;22:2496–510. doi:10.1111/mec.12179.

31. Schlüter D. Ecological Character Displacement in Adaptive Radiation. 2000;156:54–S16. doi:10.1086/303412.

32. Wheat CW, Fescemyer HW, Kvist J, TAS EVA, Vera JC, Frilander MJ, Hanski I, Marden JH: Functional Genomics of Life History Variation in a Butterfly Metapopulation. 2011;20:1813–1828. doi:10.1111/j.1365-294X.2011.05062.x.

33. Spitze K: Population Structure in *Daphnia Obtusa*: Quantitative Genetic and Allozymic Variation. 1993;135:367–374.

34. Merilä J, Crnokrak P. Comparison of Genetic Differentiation at Marker Loci and Quantitative Traits. 2001;14:892–903. doi:10.1046/j.1420-9011.2001.00348.x.

35. Koskinen MT, Haugen TO, Primmer CR. Contemporary Fisherian Life-History Evolution in Small Salmonid Populations. 2002;419:826–830. doi:10.1038/nature01029.

36. Gibson G, Weir B. The Quantitative Genetics of Transcription. 2005;21:616–623. doi:10.1016/j.tig.2005.08.010.

37. Roberge C, Guderley H, Bernatchez L. Genomewide Identification of Genes Under Directional Selection: Gene Transcription QST Scan in Diverging Atlantic Salmon Subpopulations. 2007;177:1011–1022. doi:10.1534/genetics.107.073759.

38. Raeymaekers JAM, Van Houdt JKJ, Larmuseau MHD, Geldof S, Volckaert FAM. Divergent Selection as Revealed by P(ST) and QTL-Based F(ST) in Three-Spined Stickleback (*Gasterosteus Aculeatus*) Populations along a Coastal-Inland Gradient. 2007;16:891–905. doi:10.1111/j.1365-294X.2006.03190.x

39. Chown SL, Hodgins KA, Griffin PC, Oakeshott JG, Byrne M, Hoffmann AA: Biological Invasions, Climate Change and Genomics. 2014;8:23–46. doi:10.1111/eva.12234.

40. Fabioux C, Huvet A, Le Souchu P, Le Pennec M, Pouvreau S. Temperature and Photoperiod Drive Crassostrea Gigas Reproductive Internal Clock. 2005;250:458–470. doi:10.1016/j.aquaculture.2005.02.038.

41. Huvet A, Gérard A, Ledu C, Phélipot P, Heurtebise S, Boudry P. Is Fertility of Hybrids Enough to Conclude That the Two Oysters Crassostrea Gigas and Crassostrea Angulata Are the Same Species? 2002;15:45–52. doi:10.1016/S0990-7440(01)01148-2

42. Bodoy A, Prou J, Berthome J-P. Etude comparative de différents indices de condition chez l'huitre creuse (*Crassostrea gigas*). 1986;15:173-182.

43. Latendresse JR, Warbritton AR, Jonassen H, Creasy DM. Fixation of Testes and Eyes Using a Modified Davidson's Fluid: Comparison with Bouin's Fluid and Conventional Davidson's Fluid. 2002;30:524-533. doi:10.1080/01926230290105721.

44. Steele S, Mulcahy M F. Gametogenesis of the Oyster Crassostrea Gigas in Southern Ireland. 1999;79:673-686.

45. Dheilly NM, Lelong C, Huvet A, Favrel P. Development of a Pacific Oyster (*Crassostrea Gigas*) 31,918-Feature Microarray: Identification of Reference Genes and Tissue-Enriched Expression Patterns. BioMed Central Ltd; 2011;12:468. doi:10.1186/1471-2164-12-468

46. Gentry J, Hornik K, Hothorn T, Huber W, Iacus S, Irizarry R, et al. Bioconductor: Open Software Development for Computational Biology and Bioinformatics. 2004;5:R80. doi:10.1186/gb-2004-5-10-r80

47. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. 2015;gkv007-. doi:10.1093/nar/gkv007

48. impute: impute: Imputation for microarray data. R package version 1.42.0.

49. vegan: Community Ecology Package. R package version 2.2-1 [http://cran.r-project.org/package=vegan]

50. Hedgecock D, Lin J-ZZ, DeCola S, Haudenschild CD, Meyer E, Manahan DT, Bowen B: Transcriptomic Analysis of Growth Heterosis in Larval Pacific Oysters (*Crassostrea Gigas*). 2007;104:2313–2318. doi:10.1073/pnas.0610880104

51. Zhang G, Fang X, Guo X, Li L, Luo R, Xu F, et al.: The Oyster Genome Reveals Stress Adaptation and Complexity of Shell Formation. 2012;490:49–54. doi:10.1038/nature11413.

52. Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: A Universal Tool for Annotation, Visualization and Analysis in Functional Genomics Research. 2005;21:3674–3676. doi:10.1093 bioinformatics/bti610.

53. Dheilly NM, Lelong C, Huvet A, Kellner K, Dubos M-P, Riviere G, et al. Gametogenesis in the Pacific Oyster *Crassostrea Gigas*: A Microarrays-Based Analysis Identifies Sex and Stage Specific Genes. Public Library of Science; 2012;7:e36353. doi:10.1371/journal.pone.0036353

54. Gibson G, Riley-Berger R, Harshman L, Kopp A, Vacha S, Nuzhidin S, et al. Extensive Sex-Specific Nonadditivity of Gene Expression in *Drosophila Melanogaster*. 2004;167:1791–1799. doi:10.1534/genetics.104.026583.

55. Auger DL, Gray AD, Ream TS, Kato A, Coe EH, Birchler JA. Nonadditive Gene Expression in Diploid and Triploid Hybrids of Maize. 2005, 169:389–397. doi:10.1534/genetics.104.032987.

56. Hedgecock D, Davis JP. Heterosis for Yield and Crossbreeding of the Pacific Oyster *Crassostrea Gigas*. 2007;272, Suppl:S17-S29. doi:10.1016/j.aquaculture.2007.07.226.

57. Magi R, Lindgren CM, Morris AP. Meta-Analysis of Sex-Specific Genome-Wide Association Studies. 2010;34:846–853. doi:10.1002/gepi.20540.

58. Ober C, Loisel DA, Gilad Y. Sex-Specific Genetic Architecture of Human Disease. 2008;9:911–922. doi:10.1038/nrg2415.

59. Trotter KW, Fraser IDC, Scott GK, Stutts MJ, Scott JD, Milgram SL. Alternative Splicing Regulates the Subcellular Localization of a K-Inase Anchoring Protein 18 Isoforms. 1999;147:1481–1492. doi:10.1083/jcb.147.7.1481.

60. Carr DW, Fujita A, Stentz CL, Liberty G A, Olson GE, Narumiya S. Identification of Sperm-Specific Proteins That Interact with A-Kinase Anchoring Proteins in a Manner Similar to the Type II Regulatory Subunit of PKA. 2001;276:17332–17338. doi:10.1074/jpc.M011252200.

61. Bengrine A, Li J, Awayda MS. The A-Kinase Anchoring Protein 15 Regulates Feedback Inhibition of the Epithelial Na+ Channel. 2007, 21:1189–1201. doi:10.1096/jf.06-6046com

62. Kingtong S, Kellner K, Bernay B, Goux D, Sourdaine P, Berthelin CH. Proteomic Identification of Protein Associated to Mature Spermatozoa in the Pacific Oyster *Crassostrea Gigas*. 2013;82:81–91. doi:10.1016/j.jprot.2013.02.009.

63. Jiang S-T, Chiou Y-Y, Wang E, Lin H-K, Lee S-P, Lu H-Y, Wang C-KL, Tang M-J, Li H: Targeted Disruption of Nphp1 Causes Male Infertility due to Defects in the Later Steps of Sperm Morphogenesis in Mice. 2008;17:3368–3379. doi:10.1093/hmg/ddn231.

64. Abele D. Marine Invertebrate Mitochondria and Oxidative Stress. 2007;12. doi:10.2741/2115.

65. Comninos AN, Jayasena CN, Dhillo WS. The Relationship between Gut and Adipose Hormones, and Reproduction. 2014;20:153–74. doi:10.1093/humupd/dmt033.

66. Aegerter S, Jalabert B, Bobe J. Messenger RNA Stockpile of Cyclin B, Insulin-like Growth Factor I, Insulin-like Growth Factor II, Insulin-like Growth Factor Receptor Ib, and p53 in the Rainbow Trout Oocyte in Relation with Developmental Competence. 2004;67:127–135. doi:10.1002/mrd.10384

67. Gricourt L, Bonnec G, Boujard D, Mathieu M, Kellner K. Insulin-like System and Growth Regulation in the Pacific Oyster *Crassostrea Gigas*: hrlGF-1 Effect on Protein Synthesis of Mantle Edge Cells and Expression of an Homologous Insulin Receptor-Related Receptor. 2003;134:44–56. doi:10.1016/S0016-6480(03)00217-X.

68. Jouaux A, Franco A, Heude-Berthelin C, Sourdaine P, Blin JL, Mathieu M, Kellner K. Identification of Ras, Pten and p70S6K Homologs in the Pacific Oyster *Crassostrea Gigas* and Diet Control of Insulin Pathway. 2012;176:28–38. doi:10.1016/j.ygcen.2011.12.008

69. Leinonen T, O'Hara RB, Cano JM, Merilä J. Comparative Studies of Quantitative Trait and Neutral Marker Divergence: A Meta-Analysis. 2008;21:1–17. doi:10.1111/j.1420-9101.2007.01445.x

70. Gavery M, Roberts S. DNA Methylation Patterns Provide Insight into Epigenetic Regulation in the Pacific Oyster (*Crassostrea Gigas*). 2010;11. doi:10.1186/1471-2164-11-483

71. Riviere G, Wu G-C, Fellous A, Goux D, Sourdaine P, Favrel P. DNA Methylation Is Crucial for the Early Development in the Oyster *C. Gigas*. 2013;15:739–753. doi:10.1007/s10126-013-9523-2.
72. Richards CL, Schrey AW, Pigliucci M. Invasion of Diverse Habitats by Few Japanese Knotweed Genotypes Is Correlated with Epigenetic Differentiation. 2012;15:1016–25. doi:10.1111/j.1461-0248.2012.01824.x

73. Vázquez MG, Ba CC, Spivak ED. Life History Traits of the Invasive Estuarine Shrimp *Palaemon Macrodactylus* (Caridea: Palaemonidae) in a Marine Environment (Mar Del Plata, Argentina). 2012;76:507–516. doi:10.3989/scimar.03506.02F

74. Almeida D, Argent R, Ellis A, England J, Copp GH. Environmental Biology of an Invasive Population of Signal Crayfish in the River Stort Catchment (Southeastern England). 2013;43:177–184. doi:10.1016/j.limno.2012.09.002

75. Hedrick PW, Hedgecock D. Sex Determination: Genetic Models for Oysters. 2010;101:602–611. doi:10.1093/jhered/esq065.

76. Zhang N, Xu F, Guo X. Genomic Analysis of the Pacific Oyster (*Crassostrea Gigas*) Reveals Possible Conservation of Vertebrate Sex Determination in a Mollusc. 2014;4(November):2207–2217. doi:10.1534/g3.114.013904

77. Guo X, Hedgecock D, Hershberger WK, Cooper K, Allen SK. Genetic Determinants of Protandric Sex in the Pacific Oyster, *Crassostrea Gigas* Thunberg. 1998;52:394–402. doi:10.2307/2411076.

78. Normand J, Ernande B, Haure J, McCombie H, Boudry P. Reproductive Effort and Growth in *Crassostrea Gigas*: Comparison of Young Diploid and Triploid Oysters Issued from Natural Crosses or Chemical Induction. 2009;7:229–241. doi:10.3354/ab00190.

79. Royer J, Seguineau C, Park K-l, Pouvreau S, Choi K-S, Costil K. Gametogenetic Cycle and Reproductive Effort Assessed by Two Methods in 3 Age Classes of Pacific Oysters, *Crassostrea Gigas*, Reared in Normandy. 2008;277:313–320. doi:10.1016/j.aquaculture.2008.02.033.

80. Cardoso JFMF, Langlet D, Loff JF, Martins AR, Witte JIJ, Santos PT, van der Veer HW: Spatial Variability in Growth and Reproduction of the Pacific Oyster *Crassostrea Gigas* (Thunberg, 1793) along the West European Coast. 2007;57:303–315. doi:10.1016/j.seares.2006.11.003.

Submit your next manuscript to BioMed Central
and take full advantage of:

* Convenient online submission

* Thorough peer review

* No space constraints or color figure charges

* Immediate publication on acceptance

* Inclusion in PubMed, CAS, Scopus and Google Scholar

* Research which is freely available for redistribution