- Research article
- Open Access
RNA sequencing of early round goby embryos reveals that maternal experiences can shape the maternal RNA contribution in a wild vertebrate
BMC Evolutionary Biology volume 18, Article number: 34 (2018)
It has been proposed that non-genetic inheritance could promote species fitness. Non-genetic inheritance could allow offspring to benefit from the experience of their parents, and could advocate pre-adaptation to prevailing and potentially selective conditions. Indeed, adaptive parental effects have been modeled and observed, but the molecular mechanisms behind them are far from understood.
In the present study, we investigated whether maternal RNA can carry information about environmental conditions experienced by the mother in a wild vertebrate. Maternal RNA directs the development of the early embryo in many non-mammalian vertebrates and invertebrates. However, it is not known whether vertebrate maternal RNA integrates information about the parental environment. We sequenced the maternal RNA contribution from a model that we expected to rely on parental effects: the invasive benthic fish species Neogobius melanostomus (Round Goby). We found that maternal RNA expression levels correlated with the water temperature experienced by the mother before oviposition, and identified temperature-responsive gene groups such as core nucleosome components or the microtubule cytoskeleton.
Our findings suggest that the maternal RNA contribution may incorporate environmental information. Maternal RNA should therefore be considered a potentially relevant pathway for non-genetic inheritance. Also, the ability of a species to integrate environmental information in the maternal RNA contribution could potentially contribute to species fitness and may also play a role in extraordinary adaptive success stories of invasive species such as the round goby.
In wild species, the extent to which parents are able to provide environmental information may determine the fitness of their offspring. This implies that the ability to provide parental information may be relevant for the survival and adaptation of species. Surprisingly, however, we know little about non-genetic inheritance in wild species.
Maternal genes direct early development in many animals. We studied the maternal RNA contribution in a wild fish species and found that the composition of maternally contributed RNA reflected the temperature experienced by the mother. Our findings indicate that maternally inherited RNA may not only encode a developmental start-up kit, but may also contain vital information about the environmental conditions.
It has been proposed that species who live in fluctuating environments have a fitness advantage if they can draw from the experience of their parents [20, 33]. Non-genetic inheritance has been proposed to constitute a potentially important, yet poorly understood element in the response of populations to environmental change and fluctuations . According to models, parental effects  are expected to evolve when selection is strong early in life, when parental and offspring environments positively correlate, and when opportunities for genetic adaptation are limited [20, 33, 34]. Indeed, adaptive parental effects have been observed in plants (e.g. ) and several vertebrate taxa (e.g. ), including humans (e.g. ). To our present knowledge, parental effects may be mediated by proteins, hormones, lipids, messenger RNAs, noncoding RNAs, histone modifications or DNA methylation patterns. Otherwise, we know very little about parental effects, about their prevalence, and about their role for adaptive processes in wild populations. It is unclear whether parental effects are common or rare in natural settings, we do not know which molecular mechanisms are relevant in wild vertebrates, and it has not been assessed systematically which types of information can be conveyed to the next generation.
The maternal RNA contribution plays a crucial role during early development in a wide range of organisms. In many species, development of the early embryo is directed exclusively by maternal factors . For example, fruit fly embryonic polarity and segmentation [11, 54], amphibian or fish development up to the 64- to 128-cell stages  and the initiation of mammalian development  depend on so-called maternal or maternal-effect genes. The embryonic genome is transcriptionally quiescent and therefore does not contribute to development until the Maternal-Zygotic Transition (MZT; ). In fish, this transition from maternal to embryonic control takes place after the 64-cell stage [1, 26, 27, 35, 43]. Up to the MZT, maternally contributed RNA is mostly stable [5, 6, 17, 47]. At the MZT, the embryo initiates transcription  and targeted and controlled clearing of 30–50% of maternal messages [6, 63]. Maternal genes, however, still amount to 60–70% of mRNA molecules at the peak of zygotic expression in zebrafish  and are required for the correct execution of morphogenetic steps long after the MZT .
It is currently unclear whether maternally contributed RNA provides information beyond how to physically construct an organism in wild species and under natural conditions (for example, on how to finetune its energy metabolism, or its behaviour). Most research on maternal RNA focuses on developmental roles of maternal RNA or on mechanistic aspects of RNA regulation and stability in model organisms such as C.elegans , Xenopus , or the fruit fly [46, 54, 55]. Screens in these organisms have uncovered maternal genes with roles in morphogenesis and cell fate determination. They were however not designed to detect maternal genes inducing subtle or condition-dependent changes in morphology, metabolism, or behaviour. Clues that indicate a more subtle role for maternal RNA come from fertility medicine and fish aquaculture. When ovulation is forced by light regime manipulation or hormone administration in rainbow trout, the levels of a small subset of messages differ from natural ovulation. These subtle changes are related to a decreased developmental potential of eggs from induced ovulations . Similarly, ovary transcriptome profiles from striped bass indicated a major effect of small changes in transcriptome composition on egg developmental potential . These findings suggest that small variations in the levels of maternally contributed RNAs can have profound effects on embryo development, survival and fitness.
Following a mechanistic approach [42, 57], we therefore aimed to investigate whether the maternal RNA contribution was capable to reflect a maternal experience in a natural population of vertebrates. To this end, we used a vertebrate model that we expected to rely on parental effects based on its life history traits and ecological properties [20, 33]. Neogobius melanostomus (the round goby) is a benthic fish species which deals with novel environments remarkably well. It is invasive in many parts of Europe and Northern America. In the sampled population, round goby females are sexually mature at age one or two. They reach an average age of two to three years, and a maximum age of four years [60, 28]. Round goby females can spawn up to six times per summer once the water temperature reaches 10–11 °C, but in the sampled population, two to three batches are more common ([60, 29]). Embryonic development takes between 14 and 20 days . Maternal investments in egg provisions are large , and females adapt reproductive effort flexibly to environmental conditions and disturbances . Finally, number of offspring is high (up to 5000 eggs per female, ) and selection pressure peaks during early life stages (survival to age = 1 is less than 1%, ). Thus, the round goby meets the modeled requirements for the evolution of parental effects . The home range of round goby is limited , and site fidelity is strong , with 90% of females recaptured at the same location after one year when using traps spaced 25 m apart . Therefore, individuals from one site share a common history in terms of environmental conditions.
To test whether the maternal contribution reflected maternal experiences in the round goby, we collected early cleavage embryos from a wild N. melanostomus population at water temperatures from 14 to 25°C. We selected fourteen samples at the 32 cell stage or younger. We sequenced and assembled the maternal RNA contribution from these embryos and compared it against publicly available data from zebrafish and against the round goby draft genome. We confirmed that maternal RNA is stable up to the 32 cell stage in round goby, and then determined whether the water temperature experienced by the mother was reflected in the maternal RNA by PCA and correlation analyses. Finally, we identified individual temperature-responsive maternal genes and pathways.
We morphologically characterized round goby embryonic development up to organogenesis using phalloidin stainings of the actin cytoskeleton (Fig. 1a, Additional file 1: Figure S1). We found that cleavage divisions of N. melanostomus are overall similar to zebrafish . From the 2- to the 4-cell stage, cell division happens at an oblique angle, which persists up to the 8-cell stage (Additional file 1: Figure S1). The inner cells of the embryo lose contact with the yolk and form a second layer of cells at the late 32-cell stage. Since the embryo is smaller relative to yolk volume than in zebrafish, and since neither embryo nor yolk are transparent, “high”, “oblong” or “sphere” stages (which refer to the overall shape of embryo and underlying yolk) are difficult to distinguish in N. melanostomus. Overall, however, embryonic development follows established patterns (Additional file 1: Figure S1; ).
Sequencing and assembly of maternal RNA
We sequenced and assembled the maternal RNA from 14 pre-MZT embryos and from one prim stage sample harvested at water temperatures from 14 to 24°C (Fig. 1b, c). De-novo assembly yielded 46,560 open reading frames (ORFs), of which 35,037 were expressed in early embryos and matched gene models for round goby (Additional file 2: Table S1; Additional file 3: Data S1; unpublished genome sequence). We compared the expression of the assembled ORFs with maternal RNA expression data from zebrafish  and found that maternally expressed zebrafish genes were also highly expressed in pre-MZT round goby embryos (Fig. 1d).
Effect of maternal temperature experience on the maternal RNA contribution
Next, we tested whether maternal RNA expression was related to cleavage stage (which would indicate that MZT happens before the 64-cell stage in round goby), to genetic distance (which has previously been shown to impact the expression of a subset of maternal genes in zebrafish; ), or to maternal temperature experience in the last five days before oviposition. Cleavage stage and genetic dissimilarity were used to determine the baseline or chance level of correlation between a given parameter and gene expression. Both parameters provide a buildt-in control for false discovery rates. This means that a maternally experienced parameter would have to affect the transcriptome to a greater extent than either cleavage stage or genetic dissimilarity to be considered significant. The testing variables were mean, minimum, and maximum temperature experienced by the mother during five days before oviposition, the average temperature on the sampling day, and the temperature interval experienced by the mother during five days before oviposition. These parameters were calculated from water temperature information logged continuously at the egg collection site (Fig. 1b, Additional file 4: Figure S2).
To test for a link between maternal experience and maternal RNA expression, we performed correlation analysis and PCA analysis of expression levels versus control parameters and testing parameters. Using rank based models, we quantified the correlation between the expression of each open reading frame and each of the parameters. We found that expression correlated only weakly with temperature interval, genetics, or cleavage stage, indicating that these parameters captured the baseline correlation levels expected by chance. Expression levels were correlated with the mean, minimum, or maximum temperature experienced by the mother before oviposition above the baseline levels established by the control parameters (Fig. 2, Tables 1 and 2).
PCA analysis further confirmed that the expression of maternal RNA was related to maternally experienced temperature. PC1, which explained the largest proportion of the variance, correlated with mean maternally experienced temperature (Fig. 3a and b). Together, PC1 and PC2 isolated a cluster of samples which we found to correspond to the samples with highest mean temperature (Fig. 3c). Cleavage stage, genetic distance, or maternally experienced temperature interval were not related to PC1 or PC2 (Fig. 3d), or to any other principal component (Additional file 5: Figure S3). Biplots between samples, and biplots with correlating genes or PCA loading genes marked, as well as a comparison of correlation values and PCA loading are available as supplementary information (Additional file 6: Data S2).
Temperature-effect on individual genes and pathways
We then evaluated the functional impact of the maternal temperature experience on maternal RNA composition. To this end, we focused on ORFs displaying a two-fold or higher change in expression in response to temperature according to correlation analyses (n = 120 for positive and n = 105 for negative correlation with mean experienced temperature), and on ORFs with above-average positive or negative loadings (contributions) for PC1 (n = 136 / 75) or PC2 (n = 176 / 68). As expected, we found that these ORF sets were partially overlapping. ORFs with negative loadings for PC1 or positive loadings for PC2 overlapped with temperature-correlating ORFs, while ORFs with negative loadings for PC1 and positive loadings for PC2 overlapped with temperature-anticorrelating ORFs (Fig. 4a; Additional file 4: Data S1). We blasted the selected ORFs to identify human orthologues for functional annotation. Interestingly, warmth-responsive ORFs (represented by correlating ORFs, ORFs with positive loading for PC1, and ORFs with negative loading for PC2) were more likely to produce blast hits than cold-responsive ORFs (represented by anticorrelating ORFs, ORFs with negative loading for PC1, and ORFs with positive loading for PC2; Fig. 4b).
We then performed functional analyses based on gene information available on NCBI (Additional file 7: Table S2). Round goby are poikilotherm organisms. Their metabolic rates and developmental pace depends on ambient temperature. Accordingly, we expected to find genes related to development, to energy metabolism, or to components of the stress- or heat shock response among temperature-responsive ORFs. Indeed, we identified 15 mitochondrial / energy metabolism-related proteins. We also identified two genes with heat shock pathway annotations, two clock gene orthologues (CIART and PER2), and 12 genes involved in stress response, DNA repair, or redox homeostasis. In addition, the list contained several ion- and calcium-sensitive processes or proteins, many cytoskeletal components and proteins involved in cell adhesion, GTPases and G-protein-related processes, as well as pathways related to lipid synthesis or modification and membrane function. We also identified signaling proteins and proteins involved in transcriptional regulation (for example, four components of the Notch signaling pathway), and, to our surprise, many core nucleosome components and chromatin modifiers (Fig. 4c). Most of these pathways were equally represented among cold- and heat-responsive genes. Cytoskeletal components, however, were predominantly associated with cold temperature, while chromatin and histone components were predominantly associated with warm temperature.
We confirmed these results using automated GO-term analysis. This approach identified similar key terms and pathways in the data set (Fig. 4d, Additional file 8: Table S3). ORFs that could not be assigned a conserved function also contained genes of interest. For example, two ORFs (NEME_164146 and NEME_174994) with negative loadings for PC2 coded for stonustoxin- and neoverrucotoxin-related proteins. Other temperature-responsive ORFs (such as NEME_198370 or NEME_70247) were conserved among fish.
Our analyses support our a-priori assumption that zygotic transcription starts after the 32-cell stage in the round goby. We found that cleavage stage correlates with expression only at baseline levels (as do genetic similarity and temperature interval experienced by the mother). This indicates that maternal RNA is stable across the developmental stages analysed. This is to be expected for maternal RNA, which is mostly stable up to the Maternal-Zygotic Transition, and degraded in a tightly controlled process through targeted decapping and deadenylation ([1, 17, 26, 41]. Due to its default state of stability, maternal RNA helped identify and characterize many of the regulators of mRNA stability we know today. If any of our samples would have passed the MZT, we would expect a massive signature of transcriptome turnover concomitant with cleavage progression. Cleavage stage would correlate with expression levels more than any other parameter, and cleavage stage would relate to the major principal components. Neither of those is the case (Figs. 2 and 3), which confirms that MZT happens after the 32-cell stage in the round goby. Since we cannot block RNA degradation in our system, we cannot, however, rule out the possibility that maternally directed RNA degradation systems (which act at the same time but are independemt from the MZT; ) are sensitive to seasonally variable parameters and therefore contribute to the observed effects.
Our results show a link between the expression levels of many maternally contributed genes and the maternal temperature experience. In the context of fish ecology, temperature is a proxy for a range of seasonal variations, such as light intensity or seasonal abundance of specific prey items. Accordingly, the observed effects may be related to any seasonally fluctuating parameter. Parental diet, for example, has been shown to affect the composition of maternal RNA in zebrafish . Importantly, egg quality is not linked to season in the round goby. Egg quality is defined as the ability of the embryo to develop to hatching and has been linked to female size, genetics, handling, stress, non-optimal timing of oocyte harvest, induced ovulation, yolk composition, protein and hormone content, and RNA content of oocytes ([10, 8, 68] and references therein). Round goby hatching success is very high (~95% ) and stable throughout the year (; own observation). Also, resources are not limiting at the sampling site. Food (mussels and gammarids) is abundant, and females have identical length-weight ratios in the early, late, and outside the spawning season (own observation). In fact, round goby withhold eggs rather than produce a poor quality batch in suboptimal conditions . Therefore, available data do not support an interpretation of the observed effects by a general deterioration of egg quality or female condition along the season.
Our findings imply that wild species are able to adapt the levels of maternally contributed RNA messages (within developmental constraints) to key environmental factors. We expect that this is particularly true for fish, amphibians, or birds. Those species rely more on maternally contributed factors for embryonic development than others . Mammals are generally less dependent on maternally contributed mRNA. Nonetheless, they require maternal genes for development . Mice embryos start to express their own genes only at the 2- to 4-cell stage , human embryos at the 4- to 8-cell stage . In mammals, RNAs paternally contributed through sperm have recently been reported to mediate parental experiences under laboratory conditions [23, 52], and may also be relevant for non-genetic inheritance under natural conditions.
Many pathways and molecules have been suggested to carry information from one generation to the next. For example, hormones or proteins can be allocated in varying amounts into the egg by the mother, as has been shown for androgen levels in collared flycatchers . Also, the extent of maternal provisioning with energy can affect offspring fitness. For example, the amount of fatty acids provided in the egg has been shown to impact the performance of marine fish offspring . Epigenetic mechanisms such as DNA methylation have also been suggested to contribute to transgenerational information transfer in wild vertebrates . Our results suggest for the first time that the maternal RNA contribution may also present a relevant mechanism. In stickleback and tilapia, temperature tolerance has a strong maternal component [45, 56]. Possibly, maternally provided RNAs may contribute to the observed intergenerational memory of maternal temperature experience. Whether or not many species use maternal RNA to inherit temperature information remains speculation. Based on the idea that “nothing in biology makes sense except in the light of evolution” , we propose that evolution will favor the use of any biological pathway which has the capacity to enhance the performance and the fitness of individuals and species. Accordingly, we speculate that each species will use the pathway which is most suited to transmit the most relevant information in its specific environment.
It remains to be determined whether the observed changes in expression levels have any information value for round goby offspring. Certainly, functional annotations can only suggest, but not predict, possible outcomes for offspring. Human orthologues, which are used in most cases to determine gene roles, may have functionally diverged from their round goby counterparts. Also, ORFs for which we could not infer a function may be highly relevant for temperature pre-adaptation. However, the identified gene groups and pathways allow to develop testable hypotheses on intergenerational priming strategies of wild vertebrates. For example, one may speculate that fish react to temperature or stressors by varying histone expression and global chromatin dynamics in their offspring as an epigenetic diversification strategy. Introducing a random component to gene expression may enable offspring to follow a diverse range of reaction norms or developmental trajectories without the need for genetic variation in temperature-related genes. Also, the fluidity and dynamics of biological membranes depend on ambient temperature. Adjusting lipid metabolism, lipid modification, and cell adhesion during morphogenesis in dependence of environmental conditions seems like a reasonable approach for a poikilotherm organism. Finally, many of the identified pathways are associated with signaling and transcriptional regulation. Regulatory pathways have the capacity to amplify and generate a more pronounced downstream effect during embryogenesis. It has been pointed out previously that small changes in expression have big effects when they feed into gene regulatory networks . Also, highly expressed genes or genes that undergo large changes in expression are not necessarily the best predictors of complex phenotypic traits. For example, minute variations (< 0.2 fold) in the ovarian expression levels of > 250 genes together explained > 90% of embryo survival in striped bass, while no single gene accounted for more than 2% of the effect . Thus, minor changes in the levels of maternally contributed transcriptional regulators and signaling pathways may be able to generate large effects at a later timepoint in embryonic development.
To clarify the functional relevance of individual pathways, exposure and manipulation experiments will be required, for example using antisense or mRNA injections on oocytes and embryos. Screening 3′ and 5′ regulatory gene regions of temperature-sensitive genes for conserved sequence motifs may allow us to gain a mechanistic understanding how maternal environmental perception and gene expression regulation are coupled during oogenesis. Finally, two major questions remain to be answered in the future: Which types of parental experiences affect the composition of maternal RNA contribution? and: Would different experiences cause different change patterns, or would they converge on similar sets of genes to more generally enhance plasticity?
In conclusion, wild vertebrate populations may rely on non-genetic means to increase their survival chances more than we appreciate today. We may not be aware of many mechanisms simply because they have not been investigated in the appropriate natural settings and in non-model organisms. From this follows that, if the maternal RNA contribution (or other parental contributions, such as proteins, hormones, and epigenetic marks) can be fine-tuned (as suggested by our experiments), such an adjustment may potentially contribute to the success (or failure) of adaptation processes . Accordingly, the ability to provide non-genetic information may be a feature under selection itself. We propose that the superior success of species such as the invasive round goby in changing environments may be attributable to a superior ability of those species to inform and pre-adapt their offspring. Finally, our approach validates that wild vertebrate populations which comply with evolutionary constraints of parental effects  may be an excellent place to gain deeper insights into mechanistic aspects of parental effects.
Collection of embryos from the wild
Embryos were sampled between May and August 2015 in the harbor Kleinhüningen in Basel, Switzerland (47.587448 N, 7.593409 E). This is a cargo harbor comprising two basins, with approximately 4500 incoming ships per year, in the river Rhine . Eight spawning traps containing PVC tubes in a basket were built and were deployed with cables to bottom of the second, more remote basin, which is connected to the river only through a 12 m wide and 129 m long passage and the 176 m wide first basin. Traps were deployed at 4 m depth according to Hirsch et al. , raised twice a week, and PVC tubes opened and inspected for eggs. Typically, one tube would contain between 250 and 2000 eggs laid by 1 to 4 females, with clutches from different females discernible by egg color variations between pale yellow and orange, arranged in circular patterns around the initial clutch. Samples were taken only when the arrangement of eggs and colors allowed to unequivocally distinguish between clutches laid by different females, or when only a single small clutch was present. Eggs were picked from the PVC tubes with a forceps and distributed to three Eppendorf tubes. The first sample was directly picked into 100% Ethanol in the field, transported on crushed ice, and stored at 4 °C in the laboratory, to be used for used for species confirmation by barcoding. The second sample was directly picked into fixative (4% methanol free formaldehyde; Formaldehyde 30% methanol free (ROTH #4235.1) diluted with 1 x PBS) in the field, transported on crushed ice, and stored at room temperature, to be used for cleavage stage identification. The third sample was picked into an empty tube in the field, flash frozen in liquid nitrogen on site, transported in liquid nitrogen, and stored at − 80 °C, to be used for RNA sequencing. There were no discernible differences between eggs laid early or late in the season in terms of morphology or size.
Water temperature was recorded at the site of embryo collection from April to August every 30 min using a Dissolved Oxygen Logger (HOBO U26–001) which was attached to one of the spawning traps and located at 4 m depth. From the recordings, average, median, minimum and maximum temperature in the five days preceding sampling, as well as the temperature interval in that time frame, were calculated.
Species confirmation by genotyping
At the time of sampling, two related Ponto-Caspian goby species, bighead goby (Ponticola kessleri) and round goby (Neogobius melanostomus), were present at the sampling site, although the latter were much more abundant than the former. The two species produce very similar eggs that are difficult to distinguish visually. Therefore, the species identity of each sample was confirmed by genotyping. DNA was extracted from the embryos using the DNeasy Blood and Tissue kit (Qiagen). Samples were genotyped using primers NG236.1 f (CAGGCTGAACCATATGGCAG, ), NG236 r (AGATCCTCCCCAACCAAGAT, ), Nme 7 f (AATGGATGGGTCAATTGCAT) and Nme 7 r (AAGGTTGAGCTGCCACTGAG; both ) and FastStart Taq DNA Polymerase (Roche). All four primers were used in the same PCR reaction. The primer pair NG236.1f / NG236r is specific to P. kessleri and amplifies a 150 bp fragment, while the primer pair Nme7f / Nme7r is specific to N. melanostomus and amplifies a 300 bp fragment. All sampled clutches were confirmed to be N. melanostomus.
Identification of early cleavage embryos
Within three days after sampling, 10–15 embryos from each sample were dechorionated in 1× PBS, permeabilized in PBST (1 x PBS + 0.1% Triton X-100) at 4 °C over night, and stained with Phalloidin (1:1000 dilution of Fluorescent Dye 488-I Phalloidin (LucernaChem, #U0281) in PBST + 1% BSA) for 2 h at room temperature. After two washes in PBST, embryos were embedded on glass slides equipped with two CoverWell silicone isolators (Electron Microscopy Sciences) in Fluoromount G (Southern Biotech). Developmental stage was assessed under a Nikon E400 Eclipse microscope according to Kimmel .
Based on the staging, fourteen samples that had not yet progressed past the 32-cell stage were chosen for sequencing (Fig. 1c, Table 3). Additionally, one sample that had progressed past the somite stage was included as internal reference (Table 3). RNA was extracted from ~ 1 ml (20–30) embryos using a standard phenol-chloroform protocol combined with a Single Cell RNA Purification Kit (Norgen, # 51800). Embryos were ground in 1.5 ml Trizol in liquid N2. The resulting powder was transferred into tubes (Eppendorf), and the homogenate was centrifuged for 10′ at 12000 g at 4 °C to separate lipids. The liquid portion of the sample was transferred to a new tube. 0.2 ml chloroform were added, the sample was vortexed for 15″, incubated for 5′ at room temperature, and again centrifuged for 15′ at 12000 g at 4 °C. 250 μl of the resulting supernatant were combined with 250 μl Buffer RL from the Single Cell RNA Purification Kit and 700 μl EtOH for further processing. Further steps were proceeded according to the manufacturer’s instruction of the Single Cell RNA Purification Kit. RNA was eluted in two steps with 8 μl and 20 μl to optimize yield and volume. DNA was removed from the sample with an RNase-Free DNase I Kit (Norgen, # 25710) according to the manufacturer’s instructions.
Library preparation and sequencing
Maternal RNA is heavily regulated at the post-transcriptional stage , which implies that maternally contributed mRNAs are not necessarily polyadenylated . We therefore used the Ribo Zero Gold rRNA Removal Kit to deplete ribosomal RNAs, and random hexamer primers instead of poly-T primers for reverse transcription, to avoid bias against deadenylated mRNA species. Two sequencing libraries were prepared from 1 to 2 μg of RNA per sample at the Genomics Facility in Basel (https://www.bsse.ethz.ch/genomicsbasel). Samples were paired-end sequenced on an Illumina NextSeq 500 using 81 cycles. We obtained a total number of 448′758’623 reads, which distributed to individual samples as indicated in Table 3. Downstream data analysis was performed using the HPC Cluster sciCORE (http://scicore.unibas.ch/) of the scientific computing core facility at University of Basel.
Assembly of the transcriptome
In a first step, raw reads were filtered according to the following parameters using Prinseq-lite (V0.20.4, ): min length 75, mean quality 15, no ambiguous nucleotides, low complexity filter dust with threshold 10. Filtered reads were then de-novo assembled with Trinity (Version 2.1.1, ) using the following parameters: paired-end assembly, min contig length 100, min glue 4, group pairs distance 300, kmer cov 4. ORFs were identified with Transdecoder as part of Trinity using the following parameters: min protein length 50. In order to remove redundancy, ORFs were de-replicated and clustered using Usearch (V8.1.1812, ) with the following parameters: remove duplicates, cluster coverage 100%, identity 98%. The transcriptome was then cleaned of rRNA using Blast+ (Version 2.2.31) and Silva LSU and SSU reference databases V123 with the following parameters: e-value = 0.001. These steps resulted in an embryonic transcriptome of N. melanostomus of 101,476 ORFs. A count table was generated using BWA mem (V0.7.12) for mapping, SAMtools (V1.2) to sort the sam file, and SAM2counts (https://github.com/vsbuffalo/sam2counts) to get read counts per ORF.
Cleaning and validation of the transcriptome
Before downstream analysis, ORFs derived from transposable elements, ORFs not matching gene models, and ORFs with very low count numbers were removed to reduce noise. The assembled ORFs were screened for transposable elements using RepeatMasker (Version 4.0.6, http://www.repeatmasker.org/) using the zebrafish reference database. 4977 (5%) potentially TE derived ORFs were removed from the transcriptome. ORFs were then matched with gene models obtained from an automated annotation (based on Augustus and Maker, Tomas Larsson, University of Gothenburg, unpublished data) of the round goby genome sequence (Irene Adrian-Kalchhauser, unpublished data). 46,560 of all assembled ORFs (46%) matched gene models. In comparison, non-matching ORFs had lower count numbers, were on average shorter than matching ORFs, and had few blast matches in existing databases. Since their identity and function could not be specified, they were omitted from further analysis. Before normalization, ORFs with very low counts (2 or less per ORF; 25% of the data) were removed. These ORFs were predominantly expressed at later stages (in sample A141) but not during early cleavages (all other samples). The remaining 35,037 ORFs were normalized using DESeq2 (; table of non-normalized counts: Additional file 5: Table S2).
We then compared our transcriptome with maternally provided messages identified in zebrafish embryos by Rauwerda et al. . Using Ensembl Biomart , we retrieved and blast searched coding sequences considered expressed (n = 11,118) and non-expressed (n = 11,353) in zebrafish embryos by Rauwerda et al.  against our transcriptome. This approach yielded 1617 and 782 round goby ORF orthologues, respectively, for which we plotted median expression values.
SNP calling and calculation of individual dissimilarity
Reads2snp (Version 2.0, ) was used to identify possible SNPs in the round goby tanscriptome. SNPs below a minimum number of 30 reads, minimum base quality of 20, and minimum mapping quality of 10 were not recorded. Vcftools (Version 0.1.14) was used to filter the SNP table and sites with missing data for any of the 14 samples were removed. These criteria resulted in a total of 4′096 SNPs that were clustered. Individual dissimilarity was calculated using the Bioconductor package SNPrelate .
Principal component analysis and correlation analysis
Principal component analysis (PCA) was performed with prcomp from the R base package on normalized ORF expression values using mean temperature, temperature interval, cleavage stage, and individual dissimilarity as explanatory variables. The same parameters were used to calculate rank-based Spearman correlation values with the expression values of each ORF in R.
For functional annotation, open reading frames which would display a two-fold or higher change in expression in response to temperature, and where Spearman correlation values with temperature would exceed ±0.7, as well as open reading frames that loaded above average on the PC analysis, were batch blasted using Blast2Go Basic (www.blast2go.com). The behavior of those open reading frames is visualized in biplots in Additional file 4: Data S1. Gene descriptions returned by Blast2Go were used to identify ortholog human genes. If Blast2Go did not return results, ORFs was manually blasted using blastx (https://blast.ncbi.nlm.nih.gov/Blast.cgi). For manual annotation, we mined Uniprot and NCBI gene descriptions. For automated cluster analysis with DAVID , human Ensembl IDs were retrieved from Ensembl BioMart. As background gene list, genes considered expressed by Rauwerda et al.  were translated into human Ensembl IDs. Functional clustering was based on FAT GO terms, cluster affinity was set to “high”. For the figure, similar clusters (such as “negative regulation of transcription” and “positive regulation of transcription”) were summarized (for example as “transcriptional regulation”), and the highest enrichment score of all clusters feeding into the summary category was retained.
Aanes H, Winata CL, Lin CH, Chen JP, Srinivasan KG, Lee SGP, et al. Zebrafish mRNA sequencing deciphers novelties in transcriptome dynamics during maternal to zygotic transition. Genome Res. 2011;21(8):1328–38.
Abrams EW, Mullins MC. Early zebrafish development: It’s in the maternal genes. Curr Op Gen Dev. 2009;19(4):396–403.
Adrian-Kalchhauser I, Hirsch PE, Behrmann-Godel J, et al. The invasive bighead goby Ponticola kessleri displays large scale genetic similarities and small scale genetic differentiation in relation with shipping patterns. Mol Ecol. 2016;25(9):1925–43.
Aken BL, Achuthan P, Akanni W, et al. Ensembl 2017. Nucleic Acids Res. 2017;45(D1):D635–42.
Audic Y, Omilli F, Osborne HB. Postfertilzation deadenylation of mRNAs in Xenopus Laevis embryos is sufficient to cause their degradation at the blastula stage. Mol Cell Biol. 1997;17(1):209–18.
Barckmann B, Simonelig M. Control of maternal mRNA stability in germ cells and early embryos. Biochim Biophys Acta Gene Regul Mech. 2013;1829(6):714–24.
Bonduriansky R, Crean A, Day T. The implications of nongenetic inheritance for evolution in changing environments. Evol Appl. 2011;5(2):192–201.
Bonnet E, Fostier A, Bobe J. Microarray-based analysis of fish egg quality after natural or controlled ovulation. BMC Genomics. 2007;8:55. https://doi.org/10.1186/1471-2164-8-55.
Braude P, Bolton V, Moore S. Human gene expression first occurs between the four- and eight-cell stages of preimplantation development. Nature. 1988;332(6163):459–61.
Brooks S, Tyler C, Egg SJ. Quality in fish: what makes a good egg? Rev Fish Biol Fisheries. 1997;7(4):387–416.
Carroll SB, Winslow GM, Schupbach T, Scott MP. Maternal control of drosophila segmentation gene-expression. Nature. 1986;323(6085):278–80.
Chapman RW, Reading BJ, Sullivan CV. Ovary transcriptome profiling via artificial intelligence reveals a transcriptomic fingerprint predicting egg quality in striped bass, Morone saxatilis. PLoS One. 2014;9(5):e96818.
Charlebois PM, Marsden JE, Goettel RG, Wolfe RK, Jude DJ, Rudnika S. The round goby, Neogobius melanostomus (Pallas), a review of European and north American literature. Champaign: Illinois Natural History Survey Special Publication. 1997;20:1–76.
Davidson EH. Emerging properties of animal gene regulatory networks. Nature. 2010;468:911–20. https://doi.org/10.1038/nature09645.
Dobzhansky T. Nothing in biology makes sense except in the light of evolution. Am Biol Teach. 1973;35(3):125–9.
Dufour BA, Hogan TM, Heath DD. Ten polymorphic microsatellite markers in the invasive round goby (Neogobius melanostomus) and cross-species amplification. Mol Ecol Notes. 2007;7(6):1205–7.
Duval C, Bouvet P, Omilli F, Roghi C, Dorel C, LeGuellec R, et al. Stability of maternal mRNA in Xenopus embryos: role of transcription and translation. Mol Cell Biol. 1990;10(8):4123–9.
Dworkin M, Dworkin-Rastl E. Functions of maternal mRNA in early development. Mol Reprod Dev. 1990;26(3):261–97.
Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26(19):2460–1.
English S, Pen I, Shea N, Uller T. The information value of non-genetic inheritance in plants and animals. PLoS One. 2015;10(1):e0116996.
Fuiman LA, Perez KO. Metabolic programming mediated by an essential fatty acid alters body composition and survival skills of a marine fish. P Roy Soc B-Biol Sci. 2015;282:1819.
Galloway LF. Maternal effects provide phenotypic adaptation to local environmental conditions. New Phytol. 2005;166(1):93–9.
Gapp K, Jawaid A, Sarkies P, et al. Implication of sperm RNAs in transgenerational inheritance of the effects of early trauma in mice. Nat Neurosci. 2014;17(5):667–9.
Gayral P, Melo-Ferreira J, Glemin S, Bierne N, Carneiro M, Nabholz B, et al. Reference-Free Population Genomics from Next-Generation Transcriptome Data and the Vertebrate-Invertebrate Gap. PLoS Gen. 2013. https://doi.org/10.1371/journal.pgen.1003457.
Grabherr MG, Haas BJ, Yassour M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotech. 2011;29(7):644–52.
Harvey SA, Sealy I, Kettleborough R, et al. Identification of the zebrafish maternal and paternal transcriptomes. Development. 2013;140(13):2703–10.
Heyn P, Kircher M, Dahl A, Kelso J, Tomancak P, Kalinka AT, Neugebauer KM. The earliest transcribed zygotic genes are short, newly evolved, and different across species. Cell Rep. 2014;6(2):285–92.
Hirsch PE, Adrian-Kalchhauser I, Flämig S, N’Guyen A, Defila R, Di Giulio A, Burkhardt-Holm P. A tough egg to crack: recreational boats as vectors for invasive goby eggs and transdisciplinary management approaches. Ecol Evol. 2016;6(3):707–15.
Horkova K, Kovac V. Rapid response of invasive round goby (Neogobius melanostomus) (Pallas, 1814) to an environmental perturbation demonstrated in reproductive parameters of females. J Appl Ichthyol. 2015;31(2):328–32.
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Prot. 2009;4(1):44–57.
Kimmel CB. Stages of embryonic development of the zebrafish. Dev Dyn. 1995;203(3):253–310.
Kovtun IF. Significance of the sex ratio in the spawning population of the round goby, Neogobius melanostomus, in relation to year-class strength in the sea of Azov. J Ichtyol. 1980;19:161–3.
Kuijper B, Hoyle RB. When to rely on maternal effects and when on phenotypic plasticity? Evolution. 2015;69(4):950–68.
Kuijper B, Johnstone RA, Townley S. The evolution of multivariate maternal effects. Plos Comp Biol. 2014;10(4):e1003550.
Lee MT, Bonneau AR, Takacs CM, Bazzini AA, DiVito KR, Fleming ES, Giraldez AJ. Nanog, Pou5f1 and SoxB1 activate zygotic gene expression during the maternal-to-zygotic transition. Nature. 2013;503(7476):360–4.
Li L, Zheng P, Dean J. Maternal control of early mouse development. Development. 2010;137(6):859–70.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biol. 2014;15(12):550.
Lynch MP, Mensinger AF. Seasonal abundance and movement of the invasive round goby (Neogobius melanostomus) on rocky substrate in the Duluth-Superior Harbor of Lake Superior. Ecol. Freshwater Fish. 2012;21:64–74. https://doi.org/10.1111/j.1600-0633.2011.00524.x.
MacInnis AJ, Corkum LD. Fecundity and reproductive season of the round goby Neogobius melanostomus in the upper Detroit River. Trans Am Fish Soc. 2000;129(1):136–44.
Marre J, Traver EC, Jose AM. Extracellular RNA is transported from one generation to the next in Caenorhabditis elegans. PNAS. 2016;113(44):12496–501.
Mishima Y, Tomari Y. Codon Usage and 30 UTR Length Determine Maternal mRNA Stability in Zebrafish. Mol Cell. 2016;61:874–85. https://doi.org/10.1016/j.molcel.2016.02.027.
Miska EA, Ferguson-Smith AC. Transgenerational inheritance: models and mechanisms of non-DNA sequence-based inheritance. Science. 2016;354(6308):59–63.
Mommens M, Fernandes JM, Bizuayehu TT, Bolla SL, Johnston IA, Babiak I. Maternal gene expression in Atlantic halibut (Hippoglossus hippoglossus L.) and its relation to egg quality. BMC Res Notes. 2010;3:138.
Newman T, Jhinku N, Meier M, Horsfield J. Dietary intake influences adult fertility and offspring in zebrafish. PLoS One. 2016;11(11):e0166394.
Nitzan T, Slosman T, Gutkovich D, et al. Maternal effects in the inheritance of cold tolerance in blue tilapia (Oreochromis aureus). Environ Biol Fish. 2016;99(12):975–81.
Nusslein-Volhard C, Wieschaus E, Kluding H. Mutations affecting the pattern of the larval cuticle in Drosophila melanogaster: I. Zygotic loci on the second chromosome. Roux Arch Dev Biol. 1984;193(5):267–82.
Paranjpe SS, Jacobi UG, van Heeringen SJ, Veenstra GJC. A genome-wide survey of maternal and embryonic transcripts during Xenopus tropicalis development. BMC Genomics. 2013;14:762.
Radford HE, Meijer HA & de Moor CH. Translational control by cytoplasmic polyadenylation in Xenopus oocytes. BBA Gene Regul Mech 2008;1779(4):217–229.
Rauwerda H, Wackers P, Pagano JFB, et al. Mother-specific signature in the maternal transcriptome composition of mature, unfertilized zebrafish eggs. PLoS One. 2016;11(1):e0147151.
Ray WR, Corkum LD. Habitat and site affinity of round goby. J Great Lakes Res. 2001;27:329–34.
Rice AM, Vallin N, Kulma K, et al. Optimizing the trade-off between offspring number and quality in unpredictable environments: testing the role of differential androgen transfer to collared flycatcher eggs. Horm Behav. 2013;63(5):813–22.
Rodgers AB, Morgan CP, Leu NA, Bale TL. Transgenerational epigenetic programming via sperm microRNA recapitulates effects of paternal stress. PNAS. 2015;112(44):13699–704.
Schmieder R, Edwards R. Quality control and preprocessing of metagenomic datasets. Bioinformatics. 2011;27(6):863–4.
Schupbach T, Wieschaus E. Maternal-effect mutations altering the anterior-posterior pattern of the drosophila embryo. Roux Arch Dev Biol. 1986;195(5):302–17.
Schupbach T, Wieschaus E. Female sterile mutations on the second chromosome of Drosophila melanogaster. I. Maternal effect mutations. Genetics. 1989;121(1):101–17.
Shama LN, Strobel A, Mark FC, Wegner KM. Transgenerational plasticity in marine sticklebacks: maternal effects mediate impacts of a warming ocean. Funct Ecol. 2014;28(6):1482–93.
Shea N, Pen I, Uller T. Three epigenetic information channels and their different roles in evolution. J Evol Biol. 2011;24(6):1178–87.
Somero GN. The physiology of climate change: how potentials for acclimatization and genetic adaptation will determine 'winners' and 'losers'. J Exp Biol. 2010;213(6):912–20.
Tadros W, Lipshitz HD. The maternal-to-zygotic transition: a play in two acts. Development. 2009;136(18):3033–42.
Tomczak MT, Sapota MR. The fecundity and gonad development cycle of the round goby (Neogobius melanostomus Pallas 1811) from the Gulf of Gdansk. Oceanol Hydrobiol Stud. 2006;35(4):353–67.
Vyskočilová M, Ondračková M, Šimková A, Martin J-F. Isolation and characterization of microsatellites in Neogobius kessleri (Perciformes, Gobiidae) and cross-species amplification within the family Gobiidae. Mol Ecol Notes. 2007;7(4):701–4.
Wagner DS, Dosch R, Mintzer KA, Wiemelt AP, Mullins MC. Maternal control of development at the midblastula transition and beyond. Mutants from the zebrafish II Dev Cell. 2004;6(6):781–90.
Walser CB, Lipshitz HD. Transcript clearance during the maternal-to-zygotic transition. Curr Op Gen Dev. 2011;21(4):431–43.
Wang QT, Piotrowska K, Ciemerych MA, et al. A genome-wide study of gene activity reveals developmental signaling pathways in the preimplantation mouse embryo. Dev Cell. 2004;6(1):133–44.
Weyrich A, Benz S, Karl S, Jeschek M, Jewgenow K, Fickel J. Paternal heat exposure causes DNA methylation and gene expression changes of in wild guinea pig sons. Ecology and evolution. 2016;6(9):2657–66.
Wilt FH. Polyadenylation of maternal RNA of sea urchin eggs after fertilization. PNAS. 1973;70(8):2345–9.
Wolfe RK, Marsden JE. Tagging methods for the round goby (Neogobius melanostomus). J Great Lakes Res. 1998;24(3):731–5.
Zarski D, Nguyen T, Le Cam A, Montfort J, Dutto G, Vidal MO, Fauvel C, Bobe J. Transcriptomic profiling of egg quality in sea bass (Dicentrarchus labrax) sheds light on genes involved in ubiquitination and translation. Mar Biotechnol. 2017;19:102–15.
Zhang K, Smith GW. Maternal control of early embryongenesis in mammals. Reprod Fertil Dev. 2015;27(6):880–96.
Zheng X, Levine D, Shen J, Gogarten SM, Laurie C, Weir BS. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics. 2012;28:3326–8. https://doi.org/10.1093/bioinformatics/bts606.
We would like to thank Michael Beister, Vincent Somerville, and Giulia Pellizzato for field work, Nicole Seiler for establishing Phalloidin stainings, Daniel Lüscher for technical support with spawning traps, Christian Beisel from the Quantitative Genomics Facility in Basel for advice and guidance during sequencing, Charles Betz from the Biocenter Basel for advice on embryo preparations, the local harbor and fisheries authorities for granting permissions for field work, and the SciCore Center of the University of Basel for providing server space, computational power, and user support. We thank three anonymous reviewers for their valuable comments and their effort.
This work was supported by a grant from the Forschungsfonds University of Basel to Irene Adrian-Kalchhauser. The funding body had no role in the design, execution, analysis or interpretation of the study.
Availability of data and materials
The raw data generated during the current study are available from the corresponding author on request.
Animals were obtained with the approval of local fishery authorities under permission # 2–3–6-4-1 from the office for environment and energy Basel.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Embryonic development of N. melanostomus. Phalloidin stainings of embryos collected in the field, ordered by developmental stage. (PDF 14950 kb)
Table S1. Read counts. The table contains raw read counts for all ORFs that match a predicted open reading frame according to AUGUSTUS and MAKER annotations of the round goby draft genome. (TXT 2672 kb)
Data S1. The N. melanostomus maternal transcriptome. FASTA sequences of de-novo assembled, maternally expressed open reading frames which match gene models in the draft N. melanostomus genome. (PDF 2139 kb)
Figure S2. Biplots of control and testing variables. (PDF 157 kb)
Figure S3. PCA correlations. (PDF 137 kb)
Data S2. Supplementary analyses of RNA sequencing data. Biplots between samples, biplots with correlating genes or PCA loading genes marked, comparison of correlation values and PCA loading. (PDF 1792 kb)
Table S2. Functional annotation. (XLSX 76 kb)
Table S3. DAVID GO term analysis. (XLSX 110 kb)