Revealing Mytilus galloprovincialis transcriptomic profiles during ontogeny

Revealing Mytilus galloprovincialis transcriptomic profiles 1 during ontogeny 2 3 Rebeca Moreiraa, Patricia Pereiroa, Pablo Balseirob, Massimo Milanc, Marianna Paulettoc, 4 Luca Bargellonic, Beatriz Novoaa, Antonio Figuerasa * 5 6 7 a Instituto de Investigaciones Marinas, IIM - CSIC. Eduardo Cabello, 6. 36208 Vigo, 8 Spain 9 b Present Address: Uni Research Environment, Uni Research AS, Nygårdsgaten 112, 10 5008 Bergen, Norway 11 c Department of Comparative Biomedicine and Food Science (BCA) University of 12 Padova. Viale dell’Università 16, 35020 Legnaro, Italy 13 14 rebecamore[email protected] 15 [email protected] 16 [email protected] 17 [email protected] 18 marianna.pa[email protected] 19 [email protected] 20 [email protected] 21 antoniofiguera[email protected]c.es 22 23 * Corresponding author 24 Antonio Figueras 25 Instituto de Investigaciones Marinas (IIM) 26 Consejo Superior de Investigaciones Científicas (CSIC) 27 Eduardo Cabello 6, 36208 Vigo 28 Spain 29 e-mail: [email protected] 30 31 32 - 2 - ABSTRACT 33 Mediterranean mussels are a worldwide spread bivalve species with extraordinary 34 biological success. One of the reasons of this success could be the reproduction strategy 35 of bivalves, characterized by the presence of trochophore larvae. Larval development in 36 bivalves has been a topic of raising interest in the scientific community but it deserves 37 much more attention. The principal objective of this work was to study the transcriptomic 38 profile of the ontogeny of Mytilus galloprovincialis analyzing the gene expression in 39 different developmental stages, from oocytes to juveniles. For this purpose, after 40 conducting a 454 sequencing of the transcriptomes of mussel hemocytes, adult tissues 41 and larvae, a new DNA microarray was designed and developed. 42 The studied developmental stages: unfertilized oocytes, veliger, pediveliger, settled 43 larvae and juveniles, showed very different transcriptomic profiles and clustered in 44 groups defining their characteristic gene expression along ontogeny. Our results show 45 that oocytes present a distinct and characteristic transcriptome. After metamorphosis, 46 both settled larvae and juveniles showed a very similar transcriptome, with no enriched 47 GO terms found between these two stages. This suggests: 1.- the progressive loss of RNA 48 of maternal origin through larval development and 2.- the stabilization of the gene 49 expression after settlement. 50 On the other hand during metamorphosis a specific profile of differentially expressed 51 genes was found. These genes were related to processes such as differentiation and 52 biosynthesis. Processes related to the immune response were strongly down regulated. 53 These suggest a development commitment at the expense of other non-essential 54 functions, which are temporary set aside. Immune genes such as antimicrobial peptides 55 suffer a decreased expression during metamorphosis. In fact, we found that the oocytes 56 which express a higher quantity of genes such as myticins are more likely to reach 57 success of the offspring, compared to oocytes poor in such mRNAs, whose progeny died 58 before reaching metamorphosis. 59 60 KEYWORDS: Mytilus galloprovincialis, larvae, ontogeny, development, oligo61 microarray, Myticins 62 63 HIGHLIGHTS 64  New M. galloprovincialis microarray represents its whole potential transcriptome. 65  The study of mussel ontogeny highlights key processes of each developmental stage. 66  Energy metabolism and defense are strongly regulated during metamorphosis. 67  Myticins, the most important AMP in mussels, are detected through all development. 68  The maternal immune transfer plays a key role in survival of the offspring. 69 70 - 3 - 1. INTRODUCTION 71 Mediterranean mussels (Mytilus galloprovincialis) are a worldwide spread bivalve 72 species with an extraordinary biological success; they have indeed become one of the top 73 100 invasive species of all living organisms (Gardner et al., 2016). The reproduction 74 strategy of bivalves, characterized by the presence of trochophore larvae, is common to 75 all lophotrochozoans, a group including Mollusca and Annelida among others (Halanych; 76 2004). In addition to exhibiting spiral cleavage and early cell fate determination, 77 trochozoans typically undergo indirect development, which contributes to the most 78 unique characteristics of their ontogeny. This larval development has presumably been 79 maintained from the earliest Cambrian or even the Precambrian era, in which the fossil 80 registry dates back molluscs (Ponder & Lindberg; 2008). Trochophore larvae are pelagic 81 free-living and their development occurs until they undergo metamorphosis and then 82 settle to a benthonic lifestyle. The length of this larval development is both species 83 specific and dependent on environmental conditions. In the Mediterranean mussel it 84 extends about one month. 85 A brief summary of mussel life cycle is depicted in figure 1. After male and female 86 spawn, the embryo develops in 24 hours into trochophore larvae. The late trochophore is 87 the phylotypic stage of molluscs, the point during life cycle of maximum similarity 88 among the species of a phylum (Xu et al., 2016). Trochophores use their two bands of 89 cilia surrounding their bodies to swim and feed. At this stage the primordium of the shell 90 starts to be secreted. The next day cilia develop into the velum, the most representative 91 organ of the veliger stage, and the shell becomes evident. After 20 days, veligers lose 92 their velum, larvae cannot swim anymore and the foot acquires importance. This 93 indicates the entrance in the pediveliger stage, metamorphosis now happens and larvae 94 settle when an appropriate substratum is found. Once metamorphosis is complete larvae 95 resemble the adult form (Balseiro et al., 2013) and the juveniles will attach to a substrate 96 to grow and mature to the adult stage lifestyle. 97 The ontogeny of mollusks is a complex process only completed by a small percentage of 98 the trochophore larvae. The main threats to the completion of the larval life cycle are 99 predation, food limitation, adverse environmental conditions, including lack of 100 appropriate substrate to settle, and pathogens (Eckman, 1996; Lambert et al., 1998). 101 Larval metamorphosis and settlement is a vital transition period associated with the 102 evolution, differentiation and speciation of metazoans (Degnan and Degnan, 2010). 103 Settlement can modify the larval dispersal capacity which is linked to the mortality risk, 104 because animals that have correctly settled show increased survival. Therefore, settled 105 larvae are important for the biological success of marine invertebrates (Williams and 106 Degnan, 2009). In mussels immune capacities arise during mussel development as early 107 as the trochophore stage but it is after metamorphosis that mussel larvae present an 108 important switch in the gene expression pattern, reflecting the maturation of the immune 109 system (Balseiro et al., 2013). 110 - 4 - The use of new genomic tools has increased the bivalve biology knowledge in topics such 111 as ontogeny and evolution of trochozoans (Xu et al., 2016), characterization of long non 112 coding RNAs during development (Yu et al., 2016), development of immune system 113 (Balseiro et al., 2013), hematopoiesis (Dyachuk, 2016) or trans-generational immune 114 priming (Green et al., 2016) and maternal transfer of immunity (Wang et al., 2015). 115 Nevertheless, molluscs’ ontogeny is a limited area of knowledge that deserves much 116 more attention. There are not gene markers for the different larvae stages, which could 117 ease many research tasks and also the larvae rearing in hatcheries, a key step and a 118 bottleneck for producers. 119 The principal objective of this work was to study the transcriptomic profile of the 120 ontogeny of M. galloprovincialis analyzing the gene expression in different 121 developmental stages, from oocytes to juveniles. For this study, a new DNA microarray 122 comprising sequences of hemocytes, adult tissues and larvae was designed and 123 developed. 124 To our knowledge, there are only two M. galloprovincialis microarrays previous to this 125 design: Mytarray V1.1 (Venier et al., 2006), a cDNA microarray used mainly to study 126 pollution effects on mussel (Dondero et al., 2010; Maria et al., 2016); and M. 127 galloprovincialis Inmunochip (Venier et al., 2011). This is the first microarray which 128 contains probes from several adult tissues, including hemocytes, and also from different 129 developmental stages. 130 The microarray platform and the microarray gene expression data can be found in the 131 GEO database (http://www.ncbi.nlm.nih.gov/geo/): GSE104153. 132 133 - 5 - 2. MATERIALS & METHODS 134 2.1 Sequence assembly, annotation and microarray design 135 A total of 1,325,571 raw reads from M. galloprovincialis were collected from 454 136 sequencing of hemocytes, other adult tissues and also larvae from different 137 developmental stages. Table 1 details the number of sequences by origin, assembly and 138 annotation criteria. Briefly, a total of 1,325,571 reads have been obtained. After removing 139 low-quality sequences and filtering for adaptors and primers reads were assembled 140 through MIRA3 (default parameters). Assembly produced a total of 103,995 transcripts. 141 A preliminary annotation was obtained for 43,218 contigs (41.5%) through BLASTx and 142 BLASTn similarity searches conducted against several protein databases (e-value 10-5). 143 Alignments with an e-value threshold of 10-3 and 10-5 were considered significant for 144 protein and nucleotide databases, respectively. Details on the annotation are reported in 145 S1 File. 146 All databases used for the annotation step were considered for DNA microarray platform 147 design. Putative sense-strand orientation was inferred from the matching protein-coding 148 gene in reference data bases. For 1,803 contigs that showed ambiguous orientation two 149 probes with opposite orientations (sense and antisense) were designed. For the remaining 150 41,415 contigs with putatively unambigous orientation a single (sense) probe was 151 designed. Since the microarray format could accommodate approximately 60,000 probes, 152 the longer (bp) non-annotated contigs were included in the microarray design. In total 153 7,488 non-annotated contigs were added, and for each of them, two probes with opposite 154 orientation (sense and antisense) were designed. Probe design was carried out using the 155 Agilent eArray interface (https://earray.chem.agilent.com/earray/), which applies 156 proprietary prediction algorithms to design 60-mer oligoprobes. A total of 59,971 out of 157 59,997 probes were successfully obtained, representing 50,680 different M. 158 galloprovincialis contigs. S2 File provides the sequences of the contigs to design the 159 probes and the probes of the M. galloprovincialis DNA microarray platform. Microarrays 160 were synthesized in situ using the Agilent ink-jet technology with an 8×60 K format. 161 Each array included default positive and negative controls. Probe sequences and other 162 details on the microarray platform can be found in the GEO database under accession 163 number GSE104153. 164 2.2 Larvae rearing, monitoring and sampling 165 Mature mussels were obtained from a mussel farm in Vigo (NW Spain) during the 166 spawning season. Mussels spawning and larval rearing was performed as detailed in 167 Balseiro et al., 2013. Larvae and juveniles development was followed under a Nikon 168 Eclipse 80i light microscope and photographed with a DXM1200 digital camera. Images 169 were edited with the software NIS-Elements D 2.310. A fraction of oocyte, trochophore, 170 veliger and pediveliger samples was prepared for scanning electronic microscopy (SEM). 171 Samples were pre-fixed for 2 hours at 4ºC with 4% glutaraldehyde in cacodylate buffer 172 (0.1M, pH 7.2) and then washed twice for 3 hours in cacodylate buffer. Samples were 173 further processed for SEM by the electronic microscopy service in the CACTI facilities 174 - 6 - of the University of Vigo. Briefly, samples were post-fixed with OsO4 1% in cacodylate 175 buffer (0.1M, pH 7.2) for 2 hours at 4ºC and then washed twice for 3 hours in cacodylate 176 buffer. Next, samples were mounted in the filter with the suitable pore size and 177 dehydrated in ascending series of acetone washes (30%, 50%, 70%, 80%, 90%, 100%), 178 30 minutes at 4ºC per wash, repeating the last 100% acetone step. Then, samples were 179 critical-point dried using CO2 (Baltec CPD030) and finally mount in stubs with carbon 180 adhesive and a gold coating (Emitech K550X). Imaging of the larvae was achieved using 181 a Scanning Electron Microscope FEI Quanta 200. 182 For the microarray hybridization five biological replicates were taken at each 183 developmental stage: unfertilized oocytes, trochophore (1 day post fertilization; dpf), 184 veliger (3 dpf) and pediveliger (20dpf) larvae, settled larvae (25dpf) and juveniles 185 (30dpf). 186 Larvae were centrifuged at 3000 g for 10 minutes at 4°C. The pellet was resuspended in 187 250 µl of TRIzol (Invitrogen). Total RNA isolation was conducted following the 188 manufacturer’s specifications in combination with the RNeasy Mini kit (Qiagen) for 189 RNA purification after DNase I treatment. Next, the concentration and purity of the RNA 190 were measured using a NanoDrop ND1000 spectrophotometer. Finally, RNA integrity 191 was tested on an Agilent 2100 Bioanalyzer (Agilent Technologies). Only the samples 192 with high RNA quality and quantity were used for labeling and hybridization, 3 to 5 193 biological replicates were used for each developmental stage and trochophore samples 194 were discarded due to the low RNA integrity. 195 2.3 Cy3 labeling 196 Sample labeling and hybridization were performed according to the Agilent One-Color 197 Microarray-Based Gene Expression Analysis protocol. Briefly, 100 ng of RNA from each 198 RNA sample was amplified and labeled with Cy3 using the Low Input Quick Amp 199 labeling kit (Agilent Technologies), according to the manufacturer’s instructions. Each 200 sample included a mixture of 10 different viral poly-acetylated RNAs (Agilent Spike-In 201 Mix) before amplification and labeling so the microarray analysis work-flow could be 202 assessed. Qiagen RNeasy mini spin columns were used to purify amplified RNA. Finally, 203 amplification and dye incorporation rates were verified using a NanoDrop ND1000 204 spectrophotometer. These values should be between 200 and 500 ng/μL (RNA 205 concentration) and between 20 and 50 pmol/μg aRNA (dye incorporation). 206 2.4 Microarray hybridization 207 Cy3 labeled RNA (600 ng) was fragmented with 5 μl of 10X Blocking Agent and 1 μl of 208 25X Fragmentation Buffer at 60°C for 30 min. Finally, 55 μl of 2X GE Hybridization 209 buffer was added to dilute the fragmented RNA. The eight spaces of the gasket slide were 210 filled with 40 μl of the correspondent hybridization solution and then assembled on the 211 microarray slide (each slide contained eight arrays). Slides were incubated for 17 h at 212 65°C in an Agilent Hybridization Oven, subsequently removed from the hybridization 213 chamber, quickly submerged in GE Wash Buffer 1 to disassemble the slides and then 214 - 7 - washed in GE Wash Buffer 1 for 1 minute followed by one additional wash in pre215 warmed (37°C) GE Wash Buffer 2. 216 Hybridized slides were scanned at 5 μm resolution using an Agilent G2565BA DNA 217 microarray scanner (feature extractor barcodes available in File S3). Default settings were 218 modified to scan the same slide at two different sensitivity levels (XDR Hi 100% and 219 XDR Lo 10%). The two linked images generated were analyzed together, the data were 220 extracted, and the background was subtracted using the standard procedures in the 221 Agilent Feature Extraction Software version 9.5.1. The software returned a series of spot 222 quality measures to evaluate the goodness and the reliability of the spot intensity 223 estimates. After ensuring that all of the microarrays passed the quality tests (File S4), 224 control features (positive, negative, etc.), except for SpikeIn Viral RNAs, were excluded 225 from subsequent analyses. SpikeIn control intensities were used as internal controls and 226 were expected to be uniform across the experiments of a given dataset. 227 2.5 Microarray statistical analysis 228 The GeneSpring software (Agilent) was used to normalize and analyze the microarray 229 fluorescence data. A t-test (p<0.01) with a Benjamini-Hochberg multiple testing 230 correction was carried out in filtered raw data (20 - 90th percentile), to identify 231 differentially expressed genes. The t-test was used to find the genes that were expressed 232 significantly different among the developmental stages the oocytes, larval stages and also 233 juveniles. Genes with a fold change between ±1.5 were not used for further investigation. 234 Additionally, an ANOVA (p < 0.01) was performed to find genes that were significantly 235 regulated throughout the life cycle up to a mature juvenile stage. The post-hoc test used 236 was Tukey HSD and the Benjamini-Hochberg multiple testing correction was also 237 applied. 238 Clustering of the analyzed samples was performed with AutoSOME (Newman and 239 Cooper, 2010). Clustering of myticins expression was performed in GeneSpring using k240 means grouping. To effectively separate the myticins by their similar expression profile 5 241 groups were chosen, eventually 2 of them were fused again due to their high resemblance 242 in their expression pattern. 243 2.6 GO terms and enrichment analysis 244 After statistical analysis, blast2GO software (Conesa et al., 2005) was used to assign GO 245 terms (Ashburner et al., 2000) to the significantly expressed genes through development. 246 Default values (annotation cutoff=55, GOweight=5) in blast2GO were used to perform 247 the analysis. 248 The enrichment analyses were made with the total microarray information as the 249 reference set and the list of differentially expressed genes of each developmental stage as 250 the test sets. Then, Fisher’s exact test was run: a one-tailed test, removing double IDs 251 with a false discovery rate (FDR) cut-off of 0.01. In figure 4 only the most specific terms 252 were represented. 253 - 8 - 2.7 qPCR analysis of oocytes 254 Mussel spawning was performed as previously described and eight families were 255 maintained until larvae settled or died. A fraction of the spawned oocytes of each family 256 were collected and concentrated using a nylon mesh and centrifuged. The pellet was 257 resuspended in 500 µl of Trizol (Invitrogen). Samples were grouped as good or bad 258 oocytes depending on the development of the offspring: if larvae reached the pediveliger 259 stage and settled, oocytes were classified as good (families 45, 46, 48 and 49); otherwise 260 oocytes were classified as bad (families: 47, 50, 51 and 52). 261 Total RNA extraction of oocytes was performed as described above and cDNA was 262 synthesized with 1 μg of totalRNA using SuperScript™ III Reverse Transcriptase 263 (Invitrogen) following the manufacturer’s protocol. In each group the sample with the 264 worst RNA quality were discarded (49 and 50). We finally kept 3 good oocytes samples 265 and 3 bad oocytes samples for qPCR analysis. 266 Specific PCR primers for Myticin C, Myticin B, Mytimycin and C1q (Balseiro et al., 267 2013) were prepared. Real-time quantitative PCR was performed in the7300 Real Time 268 PCR System (Applied Biosystems). One microliter of fivefold-diluted cDNA template 269 was mixed with 0.5 μl of each primer (10 μM) and 12.5 μl of SYBR Green PCR master 270 mix (Applied Biosystems) in a final volume of 25 μl. The standard cycling conditions 271 were 95°C for 10 min, followed by 40 cycles of 95°C for 15 seconds and 60°C for 30 272 seconds. All reactions were performed as technical triplicates, and an analysis of melting 273 curves was performed in each reaction. The relative expression levels of the genes were 274 normalized using 18S as a reference gene following the Pfaffl method (Pfaffl, 2001) and 275 standardized to the mean of the bad oocytes normalized expression to calculate fold 276 changes. Prior to the statistical analysis a log 2 transformation of the data was performed 277 (Rieu and Powers, 2009). A t-test for each of the studied genes was performed to evaluate 278 differences in expression between good and bad oocytes. Being the “good oocytes” those 279 whose offspring underwent metamorphosis and the “bad oocytes” those whose progeny 280 died before reaching metamorphosis. 281 282 - 9 - 3. RESULTS & DISCUSSION 283 3.1 Larval development 284 The different M. galloprovincialis developmental stages were observed and photographed 285 with a light microscope. Also, to illustrate the sampling of the larval stages, a subsample 286 of each sample was prepared for SEM. Figure 1 shows a schematic representation of the 287 mussel ontogeny with optic and SEM photographs of the gametes and the different larval 288 stages in mussel: trochophore, veliger and pediveliger larvae. Ontogeny of bivalves has 289 been studied mainly with optical microscopy (Loosanoff, Davies & Chanley, 1966). 290 Electron microscopy approaches, mainly SEM, to characterize embryonic and larval 291 development began at the end of last century and are still a field of interest, especially to 292 study shell differentiation and its abnormalities (Schönitzer and Weiss, 2007; Aranda293 Burgos et al., 2014; Balbi et al., 2016 and 2017). 294 The SEM images (Fig 1b) show the vitelline coat spikes of oocytes, a structure involved 295 in the oocyte-sperm interaction (Focarelli et al., 1991). Figure 1 b.B and b.C show the 296 spermatozoa, with their characteristic pointed acrosome. We also noticed that the sperm 297 always showed a grooved tail, something that to the best of our knowledge, has not been 298 described before. Regarding the trochophore larvae, the two characteristic ciliary bands 299 can be clearly distinguished both in the light microscope image (Fig 1a) and in the SEM 300 image (Fig 1b.D). Figures 1 b.E and b.F show how the trochophore larvae develop the 301 shell and how the cilia progressively grow to form the velum. The velum is a 302 characteristic organ of veliger larvae whose functions are feeding and movement (Fig 1a, 303 b.G and b.H). In the SEM images D to I a globular structure in the cilia of trochophore, 304 veliger and early pediveliger larvae can also be clearly observed. To the best of our 305 knowledge these structures had never been described in bivalve larvae and they seem to 306 be lost after metamorphosis (Fig 1b.J and b.K). Consequently, they could have a specific 307 function in the early mussel larvae feeding and/or movement. Prior to metamorphosis the 308 velum loses importance at the same time that the foot acquires it. The foot can be seen 309 retracted inside the shell in the light microscope picture of the pediveliger larvae. After 310 metamorphosis the larvae body shape resembles that of the adult mussels as it is observed 311 in the light microscope pictures. It is worth to note that the foot is blurred in figure 1a 312 because the settled larvae showed very active movement of the foot to settle again after 313 being sampled. 314 3.2 Assembly, annotation and microarray hybridization 315 A summary of the sequence origin, assembly and annotation results is shown in table 1. 316 From the total 1,325,571 reads from M. galloprovincialis, the MIRA3 assembler was able 317 to assemble 103,067 contigs. The putative identities of these sequences were obtained by 318 Blast in protein and nucleotide databases. One probe for annotated transcript with known 319 orientation was designed to construct a high-density oligo-DNA microarray, while two 320 probes with both orientations were designed for contigs with ambiguous orientation or no 321 annotation (Table 1). A total of 59,971 probes representing 50,706 unique transcripts 322 were created using the Agilent eArray interface (https://earray.chem.agilent.com/earray/). 323 - 16 - to be good markers for prediction of offspring success. As shown in figure 7 Myticin C, 572 Myticin B, Mytimycin, Apextrin P and C1q domain-containing proteins expression were 573 directly related to the success of offspring. All studied genes but Mytilin B showed a 574 statistical difference between the good and bad oocytes RNA load of these immune 575 related genes. The representation of the results is quite revealing, and although for 576 Myticin B no statistical significance was obtained, a strong tendency is observed: oocytes 577 with high expression levels of these genes lead to offspring that successfully reached 578 metamorphosis. The high dispersion of the results of the good families, typical in wild 579 animal samples and usually leading to no significant results, was not enough to hamper 580 the importance of the maternal immune transfer. 581 These results highlight the importance of the role of maternal immune transfer in the 582 mussel ontogeny. Maternal transfer has been studied in molluscs larvae as a part of the 583 innate immune response. It was recently verified that molluscs oocytes possess 584 significant antibacterial, lysozyme and agglutinating activities against pathogens; many 585 of these immune factors have been identified also in embryos (Wang et al., 2015). Trans586 generational immune priming has been demonstrated in oysters thanks to this mechanism 587 (Green et al., 2016). Maternal transfer of immunity is an effort that the females do to 588 benefit their offspring. This beneficial investment protects eggs and embryos from 589 horizontal (environmental) and vertical (parental) transmission of pathogens, which 590 eventually means an increased offspring survival. 591 592 593 - 17 - 4. CONCLUSIONS 594 In the present study we present the first oligo-microarray of M. galloprovincialis which 595 includes sequences from all the main tissues, including stimulated hemocytes, oocytes 596 and larval and juvenile stages. We have analyzed the transcriptome of 5 different 597 developmental stages from unfertilized oocytes to the juvenile stage and found thousands 598 of differentially expressed genes out of 59,971 probes, presumably representing the 599 whole transcriptome of the Mediterranean mussel. This the first time the transcriptome of 600 oocytes and four larval and juvenile stages have been studied as a whole. We have found 601 significantly regulated biological processes of vital importance during ontogeny and a 602 clear vision of the maternal RNA transfer to offspring and how it relates with the success 603 in development. This oligo-microarray for M. galloprovincialis has proven to be useful in 604 developmental expression analyses. And although nowadays NGSs are the chosen 605 technology to analyze transcriptomes, this tool is a cost-effective and relevant technology 606 for targeting specific research questions. 607 In summary, maternal RNA transfer seems to play a key role in development and it 608 deserves further research. Additionally, proteomics and epigenetics in ontogeny are also 609 interesting fields that should be investigated in the near future as well as the role of 610 alternative regulatory strategies during larval development such as RNAi or small RNA. 611 612 613 - 18 - AVAILABILITY OF SUPPORTING DATA 614 Microarray data are deposited in the public functional genomics data repository GEO: 615 GSE104153 616 SUPPORTING INFORMATION 617 Figure S1 618 File S1: Annotation details 619 File S2: Sequences of the original contigs and microarray probes. 620 File S3: 22 pdf files with details of the quality control of each microarray experiment. 621 File S4: Samples successfully hybridized with their specific microarray code. 622 File S5: Gene expression of hematopoiesis-related transcripts in adult mussel tissues. 623 Table S1: Top 20 DEGs during ontogeny vs oocytes. A. veliger larvae and 624 metamorphosis stage. B. settled larvae and juvelines. 625 Table S2: Top 20 DEGs during ontogeny vs juveniles. A. oocytes and veliger larvae. B. 626 metamorphosis stage and settled larvae. 627 ACKNOWLEDGEMENTS 628 This work has been funded by the EU Project REPROSEED (245119) and partially 629 supported by the Spanish Ministerio de Economía y Competitividad through Intramural 630 201640E024 and MYTIPEP (AGL2015-65705-R). We also acknowledge the support of 631 Xunta de Galicia to our group (IN607B 2016/12). 632 RM wishes to acknowledge the Spanish MICINN for her FPI Spanish research grant 633 (BES-2009-029765) and the EU H2020 funded project VIVALDI (678589). 634 We acknowledge support of the publication fee by the CSIC Open Access Publication 635 Support Initiative through its Unit of Information Resources for Research (URICI). 636 AUTHORS' CONTRIBUTIONS 637 BN and AF conceived and designed the experiments. RM and PB prepared the samples. 638 RM, MP and MM hybridized the microarrays. LB, MP and MM contributed with labeling 639 and hybridization reagents and scanning machinery. PP performed the qPCRs. BN, AF 640 and RM analyzed the data. RM wrote the paper. All authors read and approved the 641 manuscript. 642 CONFLICT OF INTEREST 643 The authors declare that the research was conducted in the absence of any commercial or 644 financial relationships that could be construed as a potential conflict of interest. 645 - 19 - ETHICS STATEMENT 646 The Mediterranean mussel, M. galloprovincialis, is not considered as an endangered or 647 protected species in any international species catalogue, including the CITES list 648 (www.cites.org) and it is not included in the list of species regulated by the EC Directive 649 2010/63/EU. Therefore, no specific authorization is required to work on mussel samples. 650 651 - 20 - REFERENCES 652 Aranda-Burgos, J.A., Da Costa, F., Nóvoa, S., Ojea, J., Martínez-Patiño, D., 2014. 653 Embryonic and larval development of Ruditapes decussatus (Bivalvia: Veneridae): a 654 study of the shell differentiation process. J. Molluscan Stud. 80, 8-16. doi: 655 10.1093/mollus/eyt044 656 Arkett, S.A., 1988. Development and senescence of control of ciliary locomotion in a 657 gastropod veliger. J. Neurobiol. 19, 612-623. 658 Ashburner, M., Ball, C.A., Blake, J.A., Botstein, D., Butler, H., Cherry, J.M., Davis, 659 A.P., Dolinski, K., Dwight, S.S., Eppig, J.T., Harris, M.A., Hill, D.P., Issel-Tarver, L., 660 Kasarskis, A., Lewis, S., Matese, J.C., Richardson, J.E., Ringwald, M., Rubin, G.M., 661 Sherlock, G., 2000. Gene ontology, tool for the unification of biology. The Gene 662 Ontology Consortium. Nat. Genet. 25, 25-29. 663 Azéma, P., Travers, M.A., De Lorgeril, J., Tourbiez, D., Dégremont, L., 2015. Can 664 selection for resistance to OsHV-1 infection modify susceptibility to Vibrio aestuarianus 665 infection in Crassostrea gigas? First insights from experimental challenges using primary 666 and successive exposures. Vet. Res. 46, 139. doi: 10.1186/s13567-015-0282-0. 667 Balbi, T., Camisassi, G., Montagna, M., Fabbri, R., Franzellitti, S., Carbone, C., Dawson, 668 K., Canesi, L., 2017. Impact of cationic polystyrene nanoparticles (PS-NH2) on early 669 embryo development of Mytilus galloprovincialis: Effects on shell formation. 670 Chemosphere. 186, 1-9. doi: 10.1016/j.chemosphere.2017.07.120. 671 Balbi, T., Franzellitti, S., Fabbri, R., Montagna, M., Fabbri, E., Canesi, L., 2016. Impact 672 of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus 673 galloprovincialis: Effects on gene transcription. Environ. Pollut. 218, 996-1004. doi: 674 10.1016/j.envpol.2016.08.050. 675 Balseiro, P., Falcó, A., Romero, A., Dios, S., Martínez-López, A., Figueras, A., Estepa, 676 A., Novoa, B., 2011. Mytilus galloprovincialis myticin C: a chemotactic molecule with 677 antiviral activity and immunoregulatory properties. PLoS One. 6, e23140. doi: 678 10.1371/journal.pone.0023140. 679 Balseiro, P., Moreira, R., Chamorro, R., Figueras, A., Novoa, B., 2013. Immune 680 responses during the larval stages of Mytilus galloprovincialis: metamorphosis alters 681 immunocompetence, body shape and behavior. Fish Shellfish Immunol. 35, 438-447. doi: 682 10.1016/j.fsi.2013.04.044. 683 Conesa, A., Götz, S., García-Gómez, J.M., Terol, J., Talón, M., Robles, M., 2005. 684 Blast2GO, a universal tool for annotation, visualization and analysis in functional 685 genomics research. Bioinformatics. 21, 3674-3676. 686 Corporeau, C., Vanderplancke, G., Boulais, M., Suquet, M., Quéré, C., Boudry, P., 687 Huvet, A., Madec, S., 2012. Proteomic identification of quality factors for oocytes in the 688 - 21 - Pacific oyster Crassostrea gigas. J. Proteomics. 75, 5554-5563. doi: 689 10.1016/j.jprot.2012.07.040. 690 Costa, M.M., Dios, S., Alonso-Gutierrez, J., Romero, A., Novoa, B., Figueras, A., 2009. 691 Evidence of high individual diversity on myticin C in mussel (Mytilus galloprovincialis). 692 Dev. Comp. Immunol. 33, 162-170. doi: 10.1016/j.dci.2008.08.005 693 Cowden, R., Curtis, S., 1981. Cephalopods. In: Ratcliffe, N.A., Rowley, A.F. (Eds.), 694 Invertebrate Blood Cells, 1. Academic Press, London, pp. 301e321. 695 De Sousa, J.T., Milan, M., Bargelloni, L., Pauletto, M., Matias, D., Joaquim, S., Matias, 696 A.M., Quillien, V., Leitão, A., and Huvet, A., 2014. A microarray-based analysis of 697 gametogenesis in two Portuguese populations of the European clam Ruditapes 698 decussatus. PloS One. 9, e92202. 699 Degnan, S.M., Degnan, B.M., 2010. The initiation of metamorphosis as an ancient 700 polyphonic trait and its role in metazoan life-cycle evolution. Phil. Trans. R. Soc. B. 365, 701 641e51. 702 Dégremont, L., Lamy, J.B., Pépin, J.F., Travers, M.A., Renault, T., 2015. New Insight for 703 the Genetic Evaluation of Resistance to Ostreid Herpesvirus Infection, a Worldwide 704 Disease, in Crassostrea gigas. PLoS One. 10, e0127917. doi: 705 10.1371/journal.pone.0127917 706 Donaghy, L., Lambert, C., Choia, K.S., Soudant, P., 2009. Hemocytes of the carpet shell 707 clam (Ruditapes decussatus) and the Manila clam (Ruditapes philippinarum): Current 708 knowledge and future prospects. Aquaculture. 297, 10-24. doi: 709 10.1016/j.aquaculture.2009.09.003 710 Dondero, F., Piacentini, L., Marsano, F., Rebelo, M., Vergani, L., Venier, P., Viarengo, 711 A., 2006. Gene transcription profiling in pollutant exposed mussels (Mytilus spp.) using a 712 new low-density oligonucleotide microarray. Gene. 376, 24-36. 713 Dravis, C., 2010. Ephs, Ephrins, and Bidirectional Signaling. Nature Education. 3, 22. 714 Dyachuk, V.A., 2016. Hematopoiesis in Bivalvia larvae: Cellular origin, differentiation 715 of hemocytes, and neoplasia. Dev. Comp. Immunol. 65, 253-257. doi: 716 10.1016/j.dci.2016.07.019 717 Eckman, J.E., 1996. Closing the larval loop: linking larval ecology to the population 718 dynamics of marine benthic invertebrates. J. Exp. Mar. Biol. Ecol. 200, 207e37. 719 FAO, 2004. Helm, M.M., Bourne, N. Hatchery culture of bivalves. A practical manual. 720 Technical paper 471. Part 5.4. Rome. 721 http://www.fao.org/docrep/007/y5720e/y5720e0a.htm#bm10.4 (accessed September 722 2017). 723 - 22 - Fisher, W.S., 1986. Structure and Functions of Oyster Hemocytes. In: Brehélin M. (eds) 724 Immunity in Invertebrates. Proceedings in Life Sciences. Springer, Berlin, Heidelberg. 725 doi: 10.1007/978-3-642-70768-1_3 726 Focarelli, R., Rosa, D., Rosati, F., 1991. The vitelline coat spikes: a new peculiar 727 structure of Mytilus galloprovincialis eggs with a role in sperm-egg interaction. Mol. 728 Reprod. Dev. 28, 143-149. 729 Gallager, S.M., Mann, R., Sasaki, C., 1986. Lipid as an index of growth and viability in 730 three species of bivalve larvae. Aquaculture. 56, 81-103. doi: 10.1016/0044731 8486(86)90020-7 732 Gandolfi, T.A., Gandolfi, F., 2001. The maternal legacy to the embryo: cytoplasmic 733 components and their effects on early development. Theriogenology. 55, 1255-1276. 734 Gardner, J.P., Zbawicka, M., Westfall, K.M., Wenne, R., 2016. Invasive blue mussels 735 threaten regional scale genetic diversity in mainland and remote offshore locations: the 736 need for baseline data and enhanced protection in the Southern Ocean. Glob. Chang. Biol. 737 22, 3182-3195. doi: 10.1111/gcb.13332 738 Green, T.J., Helbig, K., Speck, P., Raftos, D.A., 2016. Primed for success: oyster parents 739 treated with poly (I: C) produce offspring with enhanced protection against Ostreid 740 herpesvirus type I infection. Mol. Immunol. 78, 113e120. 741 Grigorian, M., Mandal, L., Hartenstein, V., 2011. Hematopoiesis at the onset of 742 metamorphosis: terminal differentiation and dissociation of the Drosophila lymph gland. 743 Dev. Genes Evol. 221, 121-31. doi: 10.1007/s00427-011-0364-6 744 Halanych, K.M., 2004. The new view of animal phylogeny. Annu. Rev. Ecol. Evol. Syst. 745 35, 229e56. 746 Hathaway, J.J.M., Adema, C.M., Stout, B.A., Mobarak, C.D., Loker, E.S., 2010. 747 Identification of protein components of egg masses indicates parental investment in 748 immunoprotection of offspring by Biomphalaria glabrata (Gastropoda, Mollusca). Dev. 749 Comp. Immunol. 34, 425e35. 750 Ichikawa, M., Yoshimi, A., Nakagawa, M., Nishimoto, N., Watanabe-Okochi, N., 751 Kurokawa, M., 2013. A role for RUNX1 in hematopoiesis and myeloid leukemia. Int. J. 752 Hematol. 97, 726-734. doi: 10.1007/s12185-013-1347-3 753 Jemaa, M., Morin, N., Cavelier, P., Cau, J., Strub, J.M., Delsert, C.A., 2014. Adult 754 somatic progenitor cells and hematopoiesis in oysters. J. Exp. Biol. 217, 3067e3077. 755 Jeong, K.H., Lie, K.J., Heyneman, D., 1983. The ultrastructure of the amebocyte 756 producing organ in Biomphalaria glabrata. Dev. Comp. Immunol. 7, 217e228. 757 - 23 - Kitamura, D., Kaneko, H., Miyagoe, Y., Ariyasu, T., Watanabe, T., 1989. Isolation and 758 characterization of a novel human gene expressed specifically in the cells of 759 hematopoietic lineage. Nucleic Acids Res. 17:9367-9379. 760 Lambert, C., Nicolas, J.L., Cilia, V., Corre, S., 1998. Vibrio pectenicida sp. nov., a 761 pathogen of scallop (Pecten maximus) larvae. Int. J. Syst. Bacteriol.48, 481e7. 762 Loosanoff, V.L., Davis, H.C., Chanley, P.E., 1966. Dimensions and shapes of larvae of 763 some marine bivalve mollusks. Malacologia. 4, 351-435. 764 Maria, V.L., Amorim, M.J., Bebianno, M.J., Dondero, F., 2016. Transcriptomic effects of 765 the non-steroidal anti-inflammatory drug Ibuprofen in the marine bivalve Mytilus 766 galloprovincialis Lam. Mar. Environ. Res. 119, 31-39. doi: 767 10.1016/j.marenvres.2016.05.010 768 Mitta, G.; Vandenbulcke, F.; Roch, P., 2000. Original involvement of antimicrobial 769 peptides in mussel innate immunity. FEBS Lett. 486, 185-190. doi: 10.1016/S0014770 5793(00)02192-X 771 Moreira, R., Pereiro, P., Canchaya, C., Posada, D., Figueras, A., Novoa, B., 2015. RNA772 Seq in Mytilus galloprovincialis: comparative transcriptomics and expression profiles 773 among different tissues. BMC Genomics. 16, 728. doi: 10.1186/s12864-015-1817-5 774 Newman, A.M., Cooper, J.B., 2010. AutoSOME: a clustering method for identifying 775 gene expression modules without prior knowledge of cluster number. BMC 776 Bioinformatics. 11, 117. doi: 10.1186/1471-2105-11-117 777 Pallavicini, A., Costa, M. del M., Gestal, C., Dreos, R., Figueras, A., Venier, P., Novoa, 778 B., 2008. High sequence variability of myticin transcripts in hemocytes of immune779 stimulated mussels suggests ancient host-pathogen interactions. Dev. Comp. Immunol. 780 32, 213-226. doi: 10.1016/j.dci.2007.05.008 781 Pauletto, M., Milan, M., Huvet, A., Corporeau, C., Suquet, M., Planas, J.V., Moreira, R., 782 Figueras, A., Novoa, B., Patarnello, T., Bargelloni, L., 2017. Transcriptomic features of 783 Pecten maximus oocyte quality and maturation. PLoS One. 12 ,e0172805. doi: 784 10.1371/journal.pone.0172805. 785 Paz, H., Lynch, M.R., Bogue, C.W., Gasson, J.C., 2010. The homeobox gene Hhex 786 regulates the earliest stages of definitive hematopoiesis. Blood. 116, 1254-1262. doi: 787 10.1182/blood-2009-11-254383. 788 Pfaffl, M.W., 2001. A new mathematical model for relative quantification in real-time 789 RT-PCR. Nucleic Acids Res. 29, e45. 790 Pila, E.A., Sullivan, J.T., Wu, X.Z., Fang, J., Rudko, S.P., Gordy, M.A., Hanington, P.C., 791 2016. Haematopoiesis in molluscs: A review of haemocyte development and function in 792 gastropods, cephalopods and bivalves. Dev. Comp. Immunol. 58, 119-128. doi: 793 10.1016/j.dci.2015.11.010 794 - 24 - Ponder, W.F. and Lindberg, D.R., 2008. Phylogeny and Evolution of the Mollusca. 795 Oakland, CA: University of California Press. 796 Rieu, I., Powers, S.J., 2009. Real-time quantitative RT-PCR: design, calculations, and 797 statistics. Plant Cell. 21, 1031-1033. doi: 10.1105/tpc.109.066001 798 Riviere, G., He, Y., Tecchio, S., Crowell, E., Gras, M., Sourdaine, P., Guo, X., Favrel, P., 799 2017. Dynamics of DNA methylomes underlie oyster development. PLoS Genet. 13, 800 e1006807. doi: 10.1371/journal.pgen.1006807 801 Rodriguez, S.R., Ojeda, F.P., Inestrosa N.C., 1993. Settlement of benthic marine 802 invertebrates. Mar. Ecol. Prog. Ser. 97, 193-207. 803 Romero, A., Costa, M.d., Forn-Cuni, G., Balseiro, P., Chamorro, R., Dios, S. Figueras, 804 A., Novoa, B., 2014.Occurrence, seasonality and infectivity of Vibrio strains in natural 805 populations of mussels Mytilus galloprovincialis. Dis. Aquat. Organ.; 108, 149-163. doi: 806 10.3354/dao02701 807 Satuito, C.G., Natoyama, K., Yamazaki, M., Fusetani, N., 1994. Larval development of 808 the mussel Mytilus edulis galloprovincialis cultured under laboratory conditions. 809 Fisheries Sci. 60, 65-68. 810 Schönitzer, V., Weiss, I.M., 2007. The structure of mollusc larval shells formed in the 811 presence of the chitin synthase inhibitor Nikkomycin Z. BMC Struct. Biol. 7, 71. 812 doi:10.1186/1472-6807-7-71 813 Sood, R., Kamikubo, Y., Liu, P., 2017. Role of RUNX1 in hematological malignancies. 814 Blood. 129, 2070-2082. doi: 10.1182/blood-2016-10-687830 815 Soufi, A., Jayaraman, P.S., 2008. PRH/Hex: an oligomeric transcription factor and 816 multifunctional regulator of cell fate. Biochem. J. 412, 399-413. doi: 817 10.1042/BJ20080035. 818 ten Hacken, E., Scielzo, C., Bertilaccio, M.T., Scarfò, L., Apollonio, B., Barbaglio, F., 819 Stamatopoulos, K., Ponzoni, M., Ghia, P., Caligaris-Cappio, F., 2013. Targeting the 820 LYN/HS1 signaling axis in chronic lymphocytic leukemia. Blood. 121, 2264-2273. doi: 821 10.1182/blood-2012-09-457119 822 Venier, P., De Pittà, C., Pallavicini, A., Marsano, F., Varotto, L., Romualdi, C., Dondero, 823 F., Viarengo, A., Lanfranchi, G., 2006. Development of mussel mRNA profiling: Can 824 gene expression trends reveal coastal water pollution? Mutat. Res.; 602, 121-134. 825 Venier, P., Varotto, L., Rosani, U., Millino, C., Celegato, B., Bernante, F., Lanfranchi, 826 G., Novoa, B., Roch, P., Figueras, A., Pallavicini, A., 2011. Insights into the innate 827 immunity of the Mediterranean mussel Mytilus galloprovincialis. BMC Genomics. 12, 828 69. doi: 10.1186/1471-2164-12-69 829 - 25 - Wang, L., Yue, F., Song, X., Song, L., 2015. Maternal immune transfer in mollusc. Dev. 830 Comp. Immunol. 48, 354-359. doi: 10.1016/j.dci.2014.05.010 831 Williams, E.A., Degnan, S.M., 2009. Carry-over effect of larval settlement cue on 832 postlarval gene expression in the marine gastropod Haliotis asinina. Mol. Ecol. 18, 833 4434e49. 834 Xu, F., Domazet-Lošo, T., Fan, D., Dunwell, T.L., Li, L., Fang, X., Zhang, G., 2016. 835 High expression of new genes in trochophore enlightening the ontogeny and evolution of 836 trochozoans. Sci. Rep. 6, 34664. doi: 10.1038/srep34664 837 Yu, H., Zhao, X., Li, Q., 2016. Genome-wide identification and characterization of long 838 intergenic noncoding RNAs and their potential association with larval development in the 839 Pacific oyster. Sci. Rep.; 6, 20796. doi: 10.1038/srep20796 840 Yue, F., Zhou, Z., Wang, L., Ma, Z., Wang, J., Wang, M., Zhang, H., Song, L., 2013. 841 Maternal transfer of immunity in scallop Chlamys farreri and its trans-generational 842 immune protection to offspring against bacterial challenge. Dev. Comp. Immunol. 41, 843 569-577. doi: 10.1016/j.dci.2013.07.001 844 845 846 - 32 - Figure 5. Normalized fluorescence values of the 326 myticin probes in the microarray 885 significantly detected by the ANOVA (top) and clustering by expression pattern 886 (bottom). 887 888 889 - 33 - Figure 6. Enrichment analysis. Representation of the BP associated to the exclusive d.e.g 890 for each larval stage with regard to the reference set (all sequences in the microarray). A. 891 Exclusive DEGs of each stage vs oocytes. B. Exclusive d.e.g of each stage vs juveniles. 892 893 - 34 - 894 - 35 - Figure 7. qPCR analysis of immune-related genes of oocytes from a single female which 895 derived in good or bad families. Data are shown as a dot plot and the mean as a solid line. 896 Asterisks show the statistical significance: *** p-value < 0.001; ** p-value < 0.01; * p897 value < 0.05 898 899 900 901 902 - 36 - Figure S1. Comparison of the myticin grouping by the k-means clusters (ANOVA) and 903 the differential expression in each developmental stage vs oocytes. A representation of 904 the fold changes for the ANOVA clusters are also shown. 905 906 907