Skip to main content

Heterogeneity in maternal mRNAs within clutches of eggs in response to thermal stress during the embryonic stage



The origin of variation is of central interest in evolutionary biology. Maternal mRNAs govern early embryogenesis in many animal species, and we investigated the possibility that heterogeneity in maternal mRNA provisioning of eggs can be modulated by environmental stimuli.


We employed two sibling species of the ascidian Ciona, called here types A and B, that are adapted to different temperature regimes and can be hybridized. Previous study showed that hybrids using type B eggs had higher susceptibility to thermal stress than hybrids using type A eggs. We conducted transcriptome analyses of multiple single eggs from crosses using eggs of the different species to compare the effects of maternal thermal stress on heterogeneity in egg provisioning, and followed the effects across generations. We found overall decreases of heterogeneity of egg maternal mRNAs associated with maternal thermal stress. When the eggs produced by the F1 AB generation were crossed with type B sperm and the progeny (‘ABB’ generation) reared unstressed until maturation, the overall heterogeneity of the eggs produced was greater in a clutch from an individual with a heat-stressed mother compared to one from a non-heat-stressed mother. By examining individual genes, we found no consistent overall effect of thermal stress on heterogeneity of expression in genes involved in developmental buffering. In contrast, heterogeneity of expression in signaling molecules was directly affected by thermal stress.


Due to the absence of batch replicates and variation in the number of reads obtained, our conclusions are very limited. However, contrary to the predictions of bet-hedging, the results suggest that maternal thermal stress at the embryo stage is associated with reduced heterogeneity of maternal mRNA provision in the eggs subsequently produced by the stressed individual, but there is then a large increase in heterogeneity in eggs of the next generation, although itself unstressed. Despite its limitations, our study presents a proof of concept, identifying a model system, experimental approach and analytical techniques capable of providing a significant advance in understanding the impact of maternal environment on developmental heterogeneity.

Peer Review reports


Variation, or heterogeneity, is the raw material of evolutionary innovations [1,2,3,4], and the importance of environmentally induced variation during the evolutionary process is well recognized [2]. Besides genetic variation, there have been a number of hypotheses proposing the evolutionary importance of environmental influence in producing heterogeneity without initial alteration of the genotype, such as the Baldwin effect, genetic assimilation and phenotypic accommodation [2, 5]. Waddington famously proposed that development is buffered against environmental as well as genetic variation and canalized into a single path, which he called developmental buffering. He also proposed that, when development faces a drastic environmental change, previously cryptic variation may be revealed in the form of new developmental paths, which may be fixed in a population after many generations as morphological evolution [5]. The so-called bet-hedging hypothesis, in which environmental fluctuation induces increased phenotypic heterogeneity that raises overall fitness by allowing at least some of the offspring to survive, also highlights the importance of environmentally induced variation in evolution [6,7,8,9,10].

Impact of maternal environment has been documented for maternal provisioning of hormones [11, 12] and egg composition [13]. A recent study also showed that maternal age affects maternal mRNA composition in mammalian eggs [14]. Maternal mRNAs are important in many animal species in shaping early embryogenesis before zygotic transcription starts [15]. Accumulating evidence suggests that maternal provisioning is important in developmental buffering, which stabilizes organismal development in the presence of thermal stress in several aquatic organisms, such as sea urchins [16], ascidians [17] and fish [18, 19]. It has been hypothesized that, by altering maternal provisioning, the maternal environment might affect the buffering level of siblings [20]. Therefore, examining changes in maternal provisioning induced by environmental stress is important in predicting future changes in the level of developmental buffering as well as developmental paths. However, the impact of maternal environment on the heterogeneity of maternal mRNAs has not been investigated. This is critical to understanding how and whether environment affects development in the next generation, providing a possible mechanism of environmentally induced evolutionary change.

We addressed the question of whether heterogeneity of maternal mRNAs in eggs is induced by heat stress using sibling species of ascidian, Ciona intestinalis type A (formally called Ciona robusta) and type B (formally called C. intestinalis) [21], which are adapted to different thermal environments [17]. C. intestinalis is a hermaphrodite invertebrate chordate, and has been a major chordate model in developmental biology for over a century given its wide availability in all oceans. However, recent studies revealed that C. intestinalis is in fact two inter-fertile sibling species showing different levels of thermal susceptibility [22,23,24,25]. We fertilized clutches of eggs from the two sibling species with sperm of one of the species, treated half the embryos with thermal stress at late neurula to early tailbud stage, reared them to mature adults, and compared the between-egg variation of maternal mRNAs in the eggs produced by stressed and un-stressed individuals using single-egg mRNA-seq. We stressed the neurula to early tailbud stage embryonic development, the stage at which all the characteristics of the chordate body plan appear, hence relevant to understanding evolution of the chordate body plan [25]. Our data suggested that exposure to thermal stress during early life (embryogenesis) altered variation of maternal mRNA in eggs produced at maturity, overall variation apparently being reduced in the first generation. At individual gene level, we found that heterogeneity of the two signaling molecules, namely hh and MEGF11 were affected by thermal stress, whereas we did not find any of the developmental buffering genes affected by thermal stress. Our study suggests the potential for alteration of development by changing environmental conditions, which might be a source of evolution of body plan.


Maternal mRNAs in eggs are affected by maternal heat stress

Previous studies have recorded between-cell or between-individual variation in expression of genes in yeast and Caenorhabditis elegans [26,27,28], but environmentally induced heterogeneity in maternal mRNAs has not been studied. Maternal mRNAs are critical in the early embryogenesis in many animal species including C. intestinalis. We first set up experiments to observe environmentally induced heterogeneity in maternal mRNAs using previously determined stress criteria in which Ciona intestinalis (type B) is susceptible to 1 h thermal stress above 26.5 °C at 8 h post fertilization at 17 °C (late neurula stage). We exposed a batch of embryos from type B conspecific cross (BB) to thermal stress (26.5 °C for 1 h) at late neurula stage, and then transferred them back to 17 °C until hatching (at c. 24 h) and allowed them to settle. After exposure to thermal stress, some larvae hatched with deformation in the tail, notochord, or central nervous system, and some embryos did not even hatch. We counted the number of larvae that hatched normally without these deformations. The proportion of normal development at hatching was 99.4% in control embryos (‘BBC’) and 94.5% in heat-treated embryos (‘BBH’), suggesting that thermal stress caused only limited differential mortality in this experiment (Fig. 1). After settlement, we reared the progeny under field conditions in a marina, where they were all exposed to the same natural environmental conditions.

Fig. 1
figure 1

Design of experiment using a type B conspecific cross. Gametes were collected by dissection from Type B specimens taken from the wild, and sperm of one individual used to fertile eggs of another. At 8 h post fertilization, half of the embryos were exposed to 26.5 ºC for 1 h and then reared at 17 ºC until they settled, metamorphosed and started feeding, whereas the other half were reared at 17 ºC without the temperature shock. Juveniles of both batches were then transferred to a marina and cultured side-by-side. After maturing in the marina, oocytes were collected from one heat-shocked and one non-shocked (control) individual for single-cell sequencing (scRNA-seq). NP, proportion of juveniles showing normal development at hatching. A red square indicates NP after heat shock, and a white square indicates without thermal stress

We isolated multiple single eggs from a mature BBH adult and from a mature BBC adult and compared maternal mRNA composition at single-egg level via scRNA-seq. The overall squared coefficient of variance (CV2) was slightly reduced following thermal stress (mean value of log10(CV2H/CV2C) = -0.04) (Fig. 2a-c; ANCOVA test, estimate -0.0014, t-value -7.25, P < 4.08e-13), decreasing in 56% of the gene models (4317 out of 7732) following maternal thermal stress, and increasing in the others. We found a significant interaction between gene model and environmental stimulus compared to a null hypothesis of non-interaction (Bartlett test, P < 2.2e−16). We also found that, following fertilization of their eggs with BB sperm, progeny of the BBH adult underwent a higher proportion of normal development (99.4%) than progeny of the BBC adult (56.8%). (Fisher’s exact test, P < 2.2e-16).

Fig. 2
figure 2

Impact of maternal embryonic heat stress on the maternally derived egg mRNAs from the type B conspecific cross. a Frequency of values of log10(CV2BBH/CV2BBC) of gene models in the eggs from heat-treated and control mothers. b, c Distribution of expression of genes showing a more than five-fold increase (shown in red) (b) or decrease (shown in blue) (c) in CV2 between eggs from heat-treated and control mothers. X-axis shows average expression level (TPM) in eggs from control (C) mother, and y-axis average expression level (TPM) from heat-shocked (H) mother in each batch. BBC (n = 14 eggs); BBH (n = 16 eggs). Note that genes showing a five-fold increase or (particularly) decrease of CV2 were spread out, and not limited to genes with low expression values

Impact of maternally inherited buffering levels on mRNA variation

Developmental buffering suppresses variation during development. Hence we tested whether maternal buffering level and maternal embryonic thermal stress impact on the level of heterogeneity of the maternal mRNAs in the eggs. Our previous studies showed that the level of buffering against embryonic heat stress is maternally inherited in Ciona and is markedly greater in type A than in type B [17]. Hence cross hybridization of type A eggs to type B sperm produces higher level of embryonic developmental buffering than a conspecific cross of type B. We made a cross of type A eggs and type B sperm, and treated a batch of embryos by thermal stress and left a second batch unstressed, as we did with the conspecific cross of type B embryos, and reared the juveniles to mature adults in the marina (‘AB cross’; Fig. 3). We collected eggs from a heat-stressed individual (ABH) and an unstressed (control) individual (ABC), conducted single-egg transcriptome analyses, and compared CV2 between the two clutches using normalized transcript levels (TPM). We found an overall decrease of CV2 in the ABH clutch compared to the ABC clutch (mean value of log10(CV2H/CV2C) = -0.13) (Fig. 4b; ANCOVA test, estimate -0.0035, t-value -22.4, P < 2.2e−16), which was prominent compared to the decrease in BB (Fig. 4a). We also examined the genes that showed a five-fold increase or decrease in the value of CV2 between heat-shocked and control individuals (Fig. 4d, e), and compared those genes with those showing the same difference in the heterogeneity of expression in the type BB cross (Table S3). We found a trend of those genes increased variation in expression have lower expression levels than those genes decreased variation. In addition, those increased variation tend to decrease expression by maternal embryonic thermal stress, whereas those decreased variation tend to increase expression by maternal embryonic thermal stress. On the other hand, we found no genes in common between the clutches from BB and AB crosses that either increased or decreased in the value of CV2 five-fold. The data suggest that alteration of heterogeneity is controlled by the mother. We also made alternative crosses using type B eggs and type A sperm (‘BA’), and type A conspecific crosses (‘AA’). None of these crosses survived in this experiment.

Fig. 3
figure 3

Experimental design of the crosses using type A eggs. A Type A specimen and a type B specimen, both wild-caught, were dissected, and A eggs were fertilized with B sperm. At 8 h post fertilization, half of the embryos were heat treated at 26.5 ºC for 1 h and then reared at 17 °C until they settled, metamorphosed and started feeding, whereas others were reared at 17 °C without heat shock. F1 juveniles of both batches were transferred to a marina and cultured side-by-side. When they matured, eggs were collected from one heat-shocked individual (ABH) and one non-shocked individual (ABC). Some eggs were subjected to scRNA-seq while others from the same clutches were crossed with sperm from another laboratory-reared type B specimen (BBC). Embryos were cultured without heat treatment until they matured (F2: ABCB and ABHB). Eggs from mature F2 animals were collected and subjected to scRNA-seq while others were crossed with a type B specimen taken from the wild. The resulting F3 embryos were examined to determine the proportion of normal development after hatching without heat treatment. Green rectangles show eggs collected for single-egg sequencing analyses. NP, proportion of normal development at hatching. Red square means NP after heat treatment, and white square without heat treatment

Fig. 4
figure 4

Differences in variance of expression of gene models in the eggs from different crosses in relation to maternal embryonic heat stress. Frequency of log10(CV2_H/CV2_C) values of gene models in eggs from individuals heat treated as embryos (H) divided by values from control individuals (C) in (a), BB conspecific cross, and (b), AB hybrid cross. In (c), ABB cross, the heat shock (H) had been administered to the embryo that became the AB mother of the individual producing the eggs. For each gene model the value of eggs from heat stressed individuals was divided by the corresponding value of the eggs from the same lineage but without heat shock. Red dotted lines indicate a log ratio of zero (CV2 values equal). d-g Distribution of genes (coloured dots) showing more than five times increase (shown in red in d, in purple in f) or decrease (shown in blue in e, shown in dark green in g) in CV2 between eggs from a heat-treated and a control individual in AB (de) or between eggs from an individual with a heat-treated mother and one with a control mother in ABB (f-g). X-axis shows average expression level (TPM) in eggs from control (C) mother, and y-axis average expression level (TPM) from heat-shocked (H) mother in each batch. Note that the genes showing a five-fold increase or decrease of CV.2 were not limited to genes with low expression levels. ABC (n = 37 eggs); ABH (n = 36); ABCB (n = 32); ABHB (n = 43)

Intergenerational inheritance of heterogeneity in maternal mRNAs

We next hypothesized that, for the change in heterogeneity to have an evolutionary effect, any observed heterogeneity would need to be maintained in the next generation even in the absence of environmental stress. To test this, we collected eggs from an ABH adult and an ABC adult, crossed them with sperm from a single type BB individual that was caught in wild, and reared the resulting F2 progeny, without embryonic thermal shock, in the lab and then in the marina (‘ABHB’ and ‘ABCB’ respectively; Fig. 3). Once these F2 produced gametes, we collected eggs individually from a type ABHB adult and a type ABCB adult for transcriptome analysis. We found a large increase in CV2 value of the ABHB eggs (mean value of log10(CV2H/CV2C) = 0.20) (Fig. 4c; ANCOVA test, estimate 0.0040, t-value 25.5, P < 2.2e-16), showing that decrease of CV2 in the F2 generation was not inherited, but was in fact reversed. We also identified genes showing fivefold changes in CV2 (Fig. 4f, g) and tested whether the same outlier genes appeared in the clutches from the previous generation (comparison of AB and ABB generations). We found a surprisingly small number of common genes between the clutches. We found only two genes, PARP14 (poly(ADP-ribose) polymerase family member 14(KY2019:KY.Chr9.491)) and HSD17B14 (hydroxysteroid 17-beta dehydrogenase 14 (KY2019:KY.Chr2.1246), in which the value of CV2 was increased more than five-fold by thermal stress in both AB and ABB crosses. We also found crystalin alpha B (KY2019:KY.Chr3.484 and KY2019:KY.Chr3.485) and HSPA1B (heat shock protein family A (Hsp70) member 1B (KY2019:KY.Chr3.426), for which the value of CV2 decreased by five times by thermal stress in both generations.

We also tested which factors impacted on the homogeneity of variance in all the six sample sets: BBC, BBH, ABC, ABH, ABCB and ABHB. A Bartlett test for all the TPM values from these samples showed that embryonic thermal stress (P < 2.2E-16) and the type of cross (P < 2.2E-16) both affected the overall homogeneity of variance in the maternal mRNAs.

Impact of thermal stress on maternally supplied developmental buffering molecules

We next tried to investigate what kind of functional properties correlate with the difference in variation in maternal mRNAs. However, in Ciona intestinalis, approximately 20% of gene models do not have orthologues to mammalian gene models [29], making the use of functional categories challenging. Therefore, rather than summarizing by functional categories, we decided to focus on testing whether thermal stress affects the expression of developmental buffering genes of offspring as previously hypothesized [20].

Previous studies on the impact of thermal stress on development largely focused on Hsp90 and Hsp70 chaperones (see review by Sato 2018 [20]), known as major players that minimize noise [30]. We found one gene putatively identified as a member of the Hsp90 family and five genes of the HSP70 family. Previous reports using Caenorhabditis elegans showed parental exposure to thermal stress increased expression level of Hsp90 [31], with transgenerational inheritance. However, we found that all the genes apart from one of the orthologues of HSPA5 and the orthologue of HSPA8 appeared to be downregulated in eggs in the two AB lineages as a result of maternal embryonic exposure to thermal stress, this included Hsp90 (Fig. 5). Apart from HSPA5 (KY2019:KY.Chr9.923) and HSPA9 (KY2019:KY.Chr3.1407), homogeneity of variance in these molecules is also significantly correlated to the genotype (Fig. 5; Table 1; Bartlett test, P < 8E-08), but not to thermal stress (Fig. 5; Table 1; Bartlett test, P > 0.01).

Fig. 5
figure 5

Impact of maternal embryonic thermal stress on maternally derived egg mRNAs encoding Hsp90 and Hsp70 family members. a Hsp90 (KY2019:KY.Chr1.1475). b HSPA5 (KY2019:KY.Chr9.923). c HSPA5 (KY2019:KY.Chr9.924). d HSPA4L (KY2019:KY.Chr14.131). e HSPA8 (KY2019:KY.Chr7.173). f HSPA9 (KY2019:KY.3.1407). 1, BBC (n = 14 eggs); 2, BBH (n = 16); 3, ABC (n = 37); 4, ABH (n = 36); 5, ABCB (n = 32); 6, ABHB (n = 43)

Table 1 P-values testing the impact of genotype and thermal stress and the interactions in the expression levels of Hsp90 and Hsp70 family members. In the variance tests (Bartlett test) P < 3.65E-06 was taken as significant after Bonferroni correction for multiple tests for 13,684 individual genes (shown in bold figures)

Other than Hsp90 and Hsp70 family members, our previous study identified 15 molecules that are maternally inherited developmental buffering molecules, which we named ‘maternally inherited developmental buffering genes (MDBGs)’ [25]. Out of the 15 MDBGs, we found eight in the transcriptome data we obtained in the present experiments. As we expected, genotype had a significant effect on the heterogeneity of gene expression of all of these genes (Fig. 6, Table 2, Table S3; P < 3E-06); apart from glycosylated lysosomal membrane protein (KY2019:KH.Chr6.691) there was up-regulation in type A egg-derived embryos relative to type B egg-derived embryos. It is of note that, in five of the eight genes, namely galactosylceramidase (KY2019:KY.Chr2.1046), catalase (KY2019:KY.Chr2.2295), DNAJC10 (KY2019:KY.Chr9.613), an uncharacterised gene (KY2019:KY.Chr8.888), and leucine aminopeptitdase 3 (KY2019:KY.Chr2.1529), expression levels were down-regulated in the eggs from thermally stressed mothers in all the BB, AB and ABB crosses (Fig. 6). We also tested whether genotype or embryonic thermal stress affected the variation of expression levels in MDBGs. In three of the genes, we found only genotype dependent differences in variation (Fig. 6; Table 2; Bartlett test, P < 2E-7), and in three of the eight genes we found a genotype dependent by thermal stress (Fig. 6; Table 2; Bartlett test, P < 2E-6). Importantly, however, we did not find any significant difference in variation solely as a result of thermal stress (Fig. 6; Table 2; Bartlett test, P > 1E-5). The data suggest that genotype, but not embryonic thermal stress solely, can impact on heterogeneity in MDBGs.

Fig. 6
figure 6

Impact of maternal embryonic thermal stress on maternal mRNAs of MDBGs in eggs. a Galactosylceramidase (KY2019:KY.Chr2.1046). b Glycosylated lysosomal membrane protein (KY2019:KY.Chr6.691). c Serine protease 27 (KY2019:KY.Chr7.688). d Catalase (KY2019:Chr2.2295). e DNAJC10 (KY2019:KY.Chr9.613). f an uncharacterised gene (KY2019:KY.8.888). g Leucine aminopeptidase 3 (LAP3) (KY2019:KY.2.1529). h Interferon regulatory factor 2 (KY2019:Chr14.407). 1, BBC (n = 14); 2, BBH (n = 16); 3, ABC (n = 37); 4, ABH (n = 36); 5, ABCB (n = 32); 6, ABHB (n = 43)

Table 2 P-values testing the impact of genotype, thermal stress, and their interaction in MDBG expression levels. In the variance tests (Bartlett test) P < 3.65E-06 was taken as significant after Bonferroni correction for multiple tests for 13,684 individual genes (shown in bold figures)

Thermal stress had different effects on different signaling molecules

We also investigated how thermal stress potentially alters development by changes in expression of developmentally relevant genes. Previous studies suggested that stability of signaling molecules is the key to developmental buffering [25, 32], suggesting that differences in variation in the expression of major signaling molecules in development may affect the stability of early embryogenesis. To this end, we examined expression levels of major signaling molecules involved in development from the transcriptome datasets (Fig. 7; Table 3). We found eight signaling molecules in our datasets, namely shh, bmp2/4, TGF-β, fgf20, MEGF(multiple-EGF domain containing gene)8, MEGF11, delta and delta-like. Unlike Hsps and MDBGs, we found that a significant difference in the variance as caused by embryonic thermal stress in MEGF11 (Bartlett test, P < 6E-10) and shh (Bartlett test, P < 9E-9). In shh, bmp2/4, MEGF8 and fgf20, the effect of embryonic thermal stress on the homogeneity of variance was genotype dependent (Bartlett test, P < 3E-6; Table 3). In fgf20 and MEGF11, heterogeneity of expression was significantly influenced by the genotype.

Fig. 7
figure 7

Impact of maternal embryonic thermal stress on various maternally supplied signaling molecules involved in development. shh (KY2019.KY.Chr8.497). b MEGF8 (KY2019.KY.3.591). c MEGF11 (KY2019.KY.Chr6.138). bmp2/4 (KY2019.KY.Chr4.449). fgf20 (KY2019.KY.14.235). f TGF-β (KY2019.KY.Chr4.450). delta (KY2019.KY.14.1158). delta-like (KY2019.KY.11.146). 1, BBC (n = 14); 2, BBH (n = 16); 3, ABC (n = 37); 4, ABH (n = 36); 5, ABCB (n = 32); 6, ABHB (n = 43)

Table 3 P-values testing the impact of genotype and thermal stress, and the interactions in the expression levels of signaling molecules. In the variance tests (Bartlett test) P < 3.65E-06 was taken as significant after Bonferroni correction for multiple tests for 13,684 individual genes (shown in bold figures)


The biological conclusions drawn from this study are limited by a number of experimental constraints, namely the use of only one mother from each cross, single thermal stress regime was examined, and the number of reads we obtained from sequence analysis varied. However, based on this proof of concept study, we suggest that there are apparent differences in homogeneity of maternal mRNA provisioning in eggs produced by individuals exposed to thermal stress during their embryonic stage, but the type of effect was dependent on the particular genotype: eggs from a BB individual showed a slight decrease in variance, whereas eggs from an AB individual showed a clear overall decrease, and eggs from an ABB individual showed an overall large increase in variance following embryonic thermal stress. Furthermore, we found only a very limited number of genes showing the same trend of increase or decrease in heterogeneity between AB (F1) and ABB (F2) generations. A previous study at a single-gene level showed transgenerational increase of variation of hsp16.2 gene expression in Caenorhabditis elegans after 1 h heat shock at 37°C [28]. However, in our experiments with Ciona thermal stress did not directly affect homogeneity of variation in Hsps and MDBGs, suggesting that variation in buffering level is determined by the genotype or the individual that produced the eggs, and is not modified by maternal embryonic thermal stress.

In contrast, we found a significant impact on the homogeneity of variance by embryonic thermal stress in two of the signaling molecules found in maternal mRNAs, shh and MEGF11. Moreover, each genotype reacted differently to embryonic thermal stress in terms of variance in shh, bmp2/4, fgf20 and MEGF8. Of these, FGF signaling and BMP signaling are important in determining neural fate in chordates [33]. The late neurula to early tailbud stage at which we exposed the embryos to thermal stress is the stage when chordate characteristics such as notochord and neural tube appear. This stage is well conserved amongst chordates and less variable than other stages. Therefore, examining the impact of environmental stress at this embryonic stage on change in morphological variation at later stages is of interest in the light of how body plan evolved by environmental stressors. Our data suggest a possible route by which the epigenetic landscape could be altered by the maternal environment. Mechanisms by which de-stabilization of expression levels can be inherited warrant further study to better understand environmentally driven phenotypic evolution.

A number of examples have provided support for the ‘bet-hedging’ [6,7,8,9, 34, 35], in that between-sibling variation increases after parental exposure to environmental stress. However, our study did not support the hypothesis in the eggs of heat-treated mothers; the overall variation was decreased in both BB and AB, and only increased in the generation following AB, ABB, in which no new thermal shock was imposed. At a single-gene level, however, we observed that expression levels of several signaling molecules in type A-egg-lineages reacted more to thermal stress than they did in the type B-egg-lineage (Table 3; Fig. 7). Previous studies showed that development is well-buffered in type A egg-derived embryos [17], thus we predicted that the expression levels of signaling molecules would be less variable in type A-egg-derived lineages than in type B conspecific crosses. Of course, there are a number of factors to be considered in our experiments before testing any hypotheses with our data. The pattern of difference in variation in the maternal mRNAs might depend on the type of environmental stimulus and might differ between individual mothers. Previous work on the variation of maternal mRNAs in zebrafish showed that variation in maternal mRNA was tightly linked to the mother [36]. In addition, AB and ABB crosses might have epigenetic changes because of hybridization. These limitations may be overcome by investigating eggs from multiple individuals in various further crosses in the future.

How is provision of maternal mRNAs altered by maternal life experience? It could be that an epigenetic change occurs in the maternal genome in the egg as a result of parental experience [37], or that alteration of the maternal genome occurs in the egg by excision of transposons or tandem duplication [38] induced by thermal stress. An intriguing possibility that requires further investigation is RNA transport from parental tissues. Recent studies have shown paternal tissue can epigenetically control gene expression in embryos through transporting miRNA to the sperm [39,40,41,42,43]. In contrast to sperm, previous studies showed transcription of maternal genome in the egg is a major source of maternal mRNA in the egg [14] and observation of RNA transport from outside the egg to inside the egg is rather limited [44,45,46]. However, many observations have been made of changes in thermal tolerance [19] and behavior [47, 12] depending on parental experience, suggesting a possibility of RNA transfer to the egg in a number of different species in different taxa. A future challenge will be to elucidate the mechanisms altering maternal mRNAs in the egg depending on maternal life experience, and how it impacts on the variation of developmental phenotype.


Our data suggest that genetic background and maternal embryonic thermal stress both affected overall homogeneity of maternal mRNA expression between individual eggs within a clutch. Whilst our data suggest that developmental buffering may not be affected by embryonic thermal stress, we observed that heterogeneity in the expression of several signaling molecules was affected by embryonic thermal stress. We suggest the possibility that embryonic environment modifies development of the next generation by altering heterogeneity in the maternal provisioning of signaling molecules. In addition to evolutionary implications, how organismal development is affected by maternal environment is one of the key issues to predict organismal development, health and disease, the so-called ‘genotype-to-phenotype map’ [47]. However, such study at a molecular level is still limited. Whilst the limited replication in this study calls for caution when interpreting these results, our study identifies a model system and experimental approach capable of providing a significant step forward in understanding the influence of maternal environment on developmental variation between siblings.

Material and Methods

Animals, eggs and embryos

Wild adults of C. intestinalis type A (C. robusta) and type B (C. intestinalis) were collected from either Queen Anne’s Battery or Sutton Harbour in Plymouth, UK, during summer (July and August) in 2017, and kept in tanks under identical conditions at 17ºC and continuous light for 1–2 days until dissection. Animals were fed a mixture of Rhinomonas reticulata and Isochrysis galbana once a day. Eggs and sperm were collected by dissection and fertilized as described previously [23, 48]. The allocations of the wild adults to Type A or B based on morphology were confirmed retrospectively by determining the genotype of four loci using sperm samples, as described by Sato et al. (2012) [23].

Heat-shock experiments and common garden rearing

The experimental design is shown in Figs. 1 and 3. Briefly, we made four different combinations of crosses using two sibling species of C. intestinalis: We denote type A conspecific crosses as ‘AA’, type B conspecific crosses as ‘BB’, crosses using type A eggs and type B sperm as ‘AB’, and type AB eggs and type B sperm as type ‘ABB’. Approximately half of embryos from each cross were exposed to thermal stress at 26.5 °C for 1 h at 8 h post fertilization as previously described [17], and those batches in which > 80% hatched normally in control conditions were used for further analysis. Hatched larvae were kept in a plastic petri dish for three days until the larvae settled, metamorphosed and started feeding. We fed juveniles with a mixture of Rhinomonas reticulata and Isochrysis galbana twice a day for a week and transferred them to Queen Anne’s Battery Marina for common-garden rearing under field conditions as described by Sato et al. (2012) [23] and Sato et al. (2014) [48]. Whilst lab conditions were kept relatively constant, rearing periods of the different crosses and generations in the marina were necessarily non-concurrent, and thus incorporated natural variation in ambient conditions. When the animals were mature, we dissected eggs and sperm from oviduct and spermiduct, respectively. All the detail of crosses used are summarized in Table S1.

Single egg sorting and cDNA library construction

Dechorionated eggs from an adult heat-treated BB (‘BBH’) individual and an adult control BB (‘BBC’) individual were dissected from the oviduct and kept on ice, collected individually in a drop of 2.3 µl of filtered sea water using 10 µl pipettes, mixed with lysis solution containing 0.1–1% TritonX (Sigma) and 2U Recombinant RNase inhibitor (Clontech) and kept at -80 °C before being processed for cDNA library construction. Dechorionated eggs from an adult from heat treated AB, an adult from control AB, an adult from heat treated ABB (‘ABHB’) and an adult from control ABB (‘ABCB’) were stored in RNA Later and kept at -80 °C before being manually isolated similarly to type BB eggs and processed for cDNA library construction. cDNA library construction was conducted following instructions of Picelli et al. (2014) [49], with some modifications: We used 0.016 µl of 1:100,000 External RNA Controls Consortium (ERCC) mRNA spike-in dilution for 20 µl total volume of cDNAs for each egg in BBC and BBH, whereas we used 0.29 µl of 1:1000 ERCC mRNA spike-in dilution for 20 µl total volume of cDNAs for each egg from ABC, ABH, ABCB and ABHB. Amplification involved 10 cycles of reverse transcription (step 11) gave the best results in amplification of cDNAs among 10–13 cycles. We also used 11 cycles of PCR re-amplification for BB samples and 14 cycles for AB and ABB (step 14), which gave us enough amount of DNA for sequencing. 21–23 million reads per sample were obtained and processed for analysis.

Transcriptome sequencing analysis

We obtained 1.7–2.2 million reads (average 1.9 million reads) for BBC, 1.7–2.4 million reads (average 1.8 million reads) for BBH, and 5.9–77 million reads (average 29 million reads) for ABC, 7.7–103 million reads (average 29 million reads) for ABH, 2.6–114 million reads (average 32 million reads) for ABCB, and 12–106 million reads (average 34 million reads) for ABHB. The adaptor sequences were trimmed from all reads using Trimmomatic [50], and trimmed reads > 10 bp in length were kept. Trimmed paired reads were mapped onto the HT Ciona intestinalis type A genome [51] using robust mapping software, BWA mem [52], using the default settings. Differential gene expression was calculated by featureCounts [53] using two gtf files: one including all the gene ids in KY gene models [51] except for 45 s rDNA, and the other containing only 45 s rDNA.

To calculate the variation in expression level of maternal mRNAs, we normalized the total read counts by TPM normalization of the gene models that showed expression in over 50% of samples (Table S2) and calculated squared coefficients of variance (CV2) using R script [54]. We also excluded data from single eggs producing low-quality data (almost no gene showing any expression). To exclude the data from single eggs with low coverage, we used only those gene models that had values above zero for over 20% of all the single-egg transcriptome samples.

Availability of data and materials

Raw sequence data: Deposited to DDBJ Sequence Read Archive under the accession numbers DRX325090-DRX325279 and DRR493923-DRR493952.

Code: of maternal mRNAs.


  1. Darwin C. On the origin of species by means of natural selection, or the preservation of favoured races in the struggle for life. London: John Murray; 1859.

    Google Scholar 

  2. West-Eberhard MJ. Phenotypic plasticity and evolution. Oxford: Oxford University Press; 2003.

    Book  Google Scholar 

  3. Frank SA. Natural selection. II. Developmental variability and evolutionary rate. J Evol Biol. 2011;24:2310–20.

    Article  CAS  PubMed  Google Scholar 

  4. Payne JL, Wanger A. The causes of evolvability and their evolution. Nat Rev Genet. 2019;20:24–38.

    Article  CAS  PubMed  Google Scholar 

  5. Waddington CH. The strategy of the genes. London: George Allen & Unwin; 1957.

    Google Scholar 

  6. Cohen D. Optimizing reproduction in a randomly varying environment. J Theoret Biol. 1966;12:119–29.

    Article  ADS  CAS  Google Scholar 

  7. Beaumont HJE, Gallie J, Kost C, Ferguson GC, Rainey RB. Experimental evolution of bet hedging. Nature. 2009;462:90–4.

    Article  ADS  CAS  PubMed  Google Scholar 

  8. Simons AM. Modes of response to environmental change and the elusive empirical evidence for bet hedging. Proc R Soc. 2011;278:1601–9.

    Article  Google Scholar 

  9. Levy SF, Ziv N, Siegal ML. Bet hedging in yeast by heterogeneous, age-correlated expression of a stress protectant. PLoS Biol. 2012;10:e1001325.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Xue B, Sartory P, Leibler S. Environment-to-phenotype mapping and adaptation strategy in varying environments. PNAS. 2019;116:13847–55.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  11. Kishimoto S, Uno M, Okabe E, Nishida E. Environmental stresses induce transgenerationally inheritable survival advantages via germline-to-soma communication in Caenorhabditis elegans. Nat Comm. 2017;8:14031.

    Article  ADS  CAS  Google Scholar 

  12. Langgenhof MR, Komdeur J. Why and how early-life environment affects development of coping behaviours. Behav Ecol Sociobiol. 2018;72:34.

    Article  Google Scholar 

  13. Sharda S, Zuest T, Tadorsky B. Predator-induced maternal effects determine adaptive antipreditor behaviour via egg composition. PNAS. 2021;118:e2017063118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Wu Y-W, Li S, Zheng W, Li Y-C, Chen L, et al. Dynamic mRNA degradome analyses indicate a role of histone H3K4 trimethylation in association with meiosis-coupled mRNA decay in oocyte aging. Nat Comm. 2022;13:3191.

    Article  ADS  CAS  Google Scholar 

  15. Robert C. Nurturing the egg: the essential connection between cumulus cells and the oocyte. Reprod Fertil Dev. 2021;34:149–59.

    Article  PubMed  Google Scholar 

  16. Fujisawa H. Temperature sensitivity of a hybrid between two species of sea urchin differing in thermotolerance. Dev Growth Differ. 1993;35:395–401.

    Article  PubMed  Google Scholar 

  17. Sato A, Kawashima T, Fujie M, Hughes S, Satoh N, Shimeld SM. Molecular basis of canalization in an ascidian species complex adapted to different thermal conditions. Sci Rep. 2015;5.

  18. Kubota Y, Shima A. Comparative study of embryonic thermoresistance of two inbred strains of the medaka (Oryzias latipes). J Comp Physiol B. 1991;160:621–5.

    Article  Google Scholar 

  19. Burt JM, Hinch SG, Patterson DA. Parental identity influences progeny responses to incubation thermal stress in sockeye salmon Onchorhynchus nerka. J Fish Biol. 2011;80:444–62.

  20. Sato A. Chaperones, Canalization, and evolution of animal forms. Int J Mol Sci. 2018;19:3029.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Brunetti R, Gissi C, Pennati R, Caicci F, Gasparini F, Manni L. Morphological evidence that the molecularly determined Ciona intestinalis type A and type B are different species: Ciona robusta and Ciona intestinalis. J Zoolog Syst Evol Res. 2015;53:186–93.

  22. Caputi L, Andreakis N, Mastrototaro F, Cirino P, Vassillo M, Sordino P. Cryptic speciation in a model invertebrate chordate. Proc Natl Acad Sci. 2007;104:9364–9.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  23. Sato A, Satoh N, Bishop JDD. Field identification of ‘types’ A and B of the ascidian Ciona intestinalis in a region of sympatry. Mar Biol. 2012;159:1611–9.

    Article  Google Scholar 

  24. Malfant M, Coudret J, Le Merdy R, Viard F. Effects of temperature and salinity on juveniles of two ascidians, one native and one invasive, and their hybrids. J Exp Mar Biol Ecol. 2017;497:180–7.

    Article  CAS  Google Scholar 

  25. Sato A, Oba GM, Aubert-Kato N, Yura K, Bishop JDD. Co-expression network analysis of environmental canalization in the ascidian Ciona. BMC Ecol Evol. 2022;22:53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Rinott R, Jaimovich A, Friedman N. Exploring transcription regulation through cell-to-cell variability. PNAS. 2011;108:6329–34.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  27. Burnaevskiy N, Sands B, Yun S, Tedesco PM, Johnson TE, Kaeberlein M, Brent R, Mendenhall A. Chaperone biomarkers of lifespan and penetrance track the dosages of many other proteins. Nat Communications. 2019;10:5725.

    Article  ADS  CAS  Google Scholar 

  28. Verti-Quintero N, Berger S, Casadevall i Solvas X, Statzer C, Annis J, et al. Stochastic and age-dependent proteostasis decline underlies heterogeneity in heat-shock response dynamics. Small. 2021;17:2102145.

    Article  CAS  Google Scholar 

  29. Satou Y, Tokuoka M, Oda-Ishii I, Sinichi T, Ishida T, et al. A manually curated gene model set for an ascidian Ciona robusta (C. intestinalis type A). Zool Sci. 2022;39:253–60.

  30. Raj A, van Oudenaarden A. Nature, nurture, or chance: Stochastic gene expression and its consequence. Cell. 2008;135:216–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Klosin A, Casas E, Hidalgo-Carcedo C, Vavouri T, Lehner B. Transgenerational transmission of environmental information in C. elegans. Science. 2017;356:320–3.

    Article  ADS  CAS  PubMed  Google Scholar 

  32. Hagolani PF, Zimm R, Marin-Riera M, Salazar-Ciudad I. Cell signaling stabilizes morphogenesis against noise. Development. 2019;146:dev179309.

    Article  CAS  PubMed  Google Scholar 

  33. Satou Y, Imai KS. Gene regulatory systems that control gene expression in the Ciona embryo. Proc Jpn Acad Ser B Phys Biol Sci. 2015;91:33–51.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  34. Hahn J, Tanner AW, Carabetta VJ, Cristea IM, Dubnau D. ComGA-RelA interaction and persistence in the Bacillus subtilis K-state. Mol Microbiol. 2015;97:454–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Veening J-W, Stewart EJ, Berngruber TW, Taddei F, Kuipers OP, Hamoen LW. Bet-hedging and epigenetic inheritance in bacterial cell development. PNAS. 2008;105:4393–8.

    Article  ADS  PubMed  PubMed Central  Google Scholar 

  36. Rauwerda H. Mother-specific signature in the maternal transcriptome composition of mature, unfertilized zebrafish eggs. PLoS ONE. 2016.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Holland ML, Lowe R, Caton PW, Gemma C, Carbajosa G, et al. Early-life nutrition modulates the epigenetic state of specific rDNA genetic variants in mice. Science. 2016;353:495–8.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  38. Ide S, Miyazaki T, Maki H, Kobayashi T. Abundance of ribosomal RNA gene copies maintains genome integrity. Science. 2010;327:693–6.

    Article  ADS  CAS  PubMed  Google Scholar 

  39. Cossetti C, Lugini L, Astrologo L, Saggio I, Fais S, Spadofora C. Soma-to-germline transmission of RNA in mice xenografted with human tumor cells: possible transport by exosomes. PLoS One. 2014;9:e101629.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  40. Sharma U, Conine CC, Shea JM, Boskovie A, Derr AG, et al. Biogenesis and function of tRNA fragments during sperm maturation and fertilization in mammals. Science. 2016;351:391–6.

    Article  ADS  CAS  PubMed  Google Scholar 

  41. Sarker G, Sun W, Rosenkranz D, Pelczar P, Opitx L, Efthymiou V, Wolfrum C, Pekeg-Raibstein D. Maternal overnutrition programs hedonic and metabolic phenotypes across generations through sperm tsRNAs. PNAS. 2019;116:10547–56.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  42. O’Brien EA, Ensbey KS, Day BW, Baldock PA, Barry G. Direct evidence for transport of RNA from the mouse brain to the germline and offspring. BMC Biol. 2020;18:45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Chan JC, Morgan CP, Leu NA, Shetty A, Cisse YM, Nugent BM, Morrison KE, et al. Reproductive tract extracellular vesicles are sufficient to transmit intergenerational stress and program neurodevelopment. Nat Comm. 2020;11:1499.

    Article  ADS  CAS  Google Scholar 

  44. Marré J, Traver EC, Jose AM. Extracellular RNA is transported from one generation to the next in Caenorhabditis elegans. PNAS. 2016;113:12496–501.

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

  45. Macaulay AD, et al. Cumulus cell transcripts transit to the bovine oocyte in preparation for maturation. Biol Rep. 2016;94:16.

    Article  CAS  Google Scholar 

  46. Coine CC, Rando OJ. Soma-to-germline RNA communication. Nat Rev Genet. 2022;23:73–88.

    Article  CAS  Google Scholar 

  47. Lehner B. Genotype to phenotype: lessons from model organisms for human genetics. Nat Rev Genet. 2013;14:168–78.

    Article  CAS  PubMed  Google Scholar 

  48. Sato A, Shimeld SM, Bishop JDD. Symmetrical reproductive compatibility of two species in the Ciona intestinalis (Ascidiacea) species complex, a model for marine genomics and developmental biology. Zool Sci. 2014;31:369–74.

    Article  Google Scholar 

  49. Picelli S, Faridani OR, Björklund ÅK, Winberg G, Sagasser S, Sandberg R. Full-length RNA-seq from single cells using Smart-seq2. Nat Protoc. 2014;9:171–81.

    Article  CAS  PubMed  Google Scholar 

  50. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Satou Y, Nakamura R, Yu D, Yoshida R, Hamada M, Fujie M, Hisata K, Takeda H, Satoh N. A nearly complete genome of Ciona intestinalis type A (C. robusta) reveals the contribution of inversion to chromosomal evolution in the genus Ciona. Genome Biol Evol. 2019;11:3144–57.

  52. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2013;30:923–30.

    Article  CAS  PubMed  Google Scholar 

  54. Crawley MJ. Statistics: An introduction using R. West Sussex: Wiley; 2005.

Download references


AS thanks Colin Brownlee, William Wilson and other members of the Marine Biological Association of the UK, Ann Harris and Chris Harris for hospitality while collecting samples at Plymouth. The calculation was partially conducted on Chaen, the supercomputer of the Center for Interdisciplinary AI and Data Science at Ochanomizu University. We are very grateful to Yura Kei and members of the supporting team Charles for giving access to Chaen, and to Kanako Hisata for technical assistance and advice. We also thank Angela Lee and the Oxford Genomics Center at the Wellcome Centre for Human Genetics (funded by Wellcome Trust grant reference 203141/Z/16/Z) and members of Kashiwa SQC for the generation of sequence data.


This project was supported by and KAKENHI grant number 17K15167, 21K06287, 16H06279 (PAGS), Sumitomo Foundation (180453), Takeda Science Foundation, JST FOREST Program (Grant Number JPMJFR224R), Hiroko Takada scholarship and Ray Lankester Investigatorship from the Marine Biological Association of the UK to AS, International Joint Grant from the Royal Society, London (IEC/R3/170091) to MT, and the Human Life Innovation Center at Ochanomizu University.

Author information

Authors and Affiliations



AS designed the research. AS, CW, JB and YS performed research. YM, and AS analysed the data. YM and MT developed analytical tools. AS wrote the paper. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Atsuko Sato.

Ethics declarations

Ethics approval and consent to participate

The invertebrate Ciona intestinalis is not a protected species in the UK and no license is required to collect it. The source population is in an artificial urban habitat with no conservation status and collection visits were made with permission from the site operators.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Frequency of normal larvae after heat shock of the parents.

Additional file 2: Table S2.

TPM values of all samples used for analysis (provided as a separate excel sheet).

Additional file 3: Table S3.

Average, variance and CV2 for BB, AB, and ABB datasets (provided as a separate excel sheet). These genes for which the value of CV2 increased or decreased by five times by thermal stress are shown as '1' in '_P5' or '_M5' column, respectively.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sato, A., Mihirogi, Y., Wood, C. et al. Heterogeneity in maternal mRNAs within clutches of eggs in response to thermal stress during the embryonic stage. BMC Ecol Evo 24, 21 (2024).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: