Gill transcriptomic analysis in fast- and slow-growing individuals of Mytilus galloprovincialis

Gill transcriptomic analysis in fastand slow-growing 1 individuals of Mytilus galloprovincialis2 Daniel Prieto*[1], Pablo Markaide[1], Iñaki Urrutxurtu[1], Enrique Navarro[1], Sebastien 3 Artigaud[2], Elodie Fleury[3], Irrintzi Ibarrola[1] & Miren Bego Urrutia[1] 4 [1] GIU 17/061, GI 544, Departamento de Genética, Antropología Física y Fisiología Animal, Facultad5 de Ciencia y Tecnología, Universidad del País Vasco/Euskal Herriko Unibertsitatea, UPV/EHU, Apartado 6 644, 48080 Bilbao, Spain. *Corresponding author, e-mail: [email protected] 7 [2]LEMAR UMR 6539 CNRS/UBO/IRD/Ifremer, IUEM, rue Dumont d’Urville, 29280 Plouzané, France8 [3]LEMAR, UMR 6539 UBO-CNRS-Ifremer-IRD, Technopole Brest Iroise 29280 Plouzané, France9 Abstract10 The molecular basis underlying the mechanisms at the origin of growth variation 11 in bivalves is still poorly understood, although several genes have been described as 12 upregulated in fast-growing individuals. In the present study, we reared mussel spat of 13 the species Mytilus galloprovincialis under diets below the pseudofaeces threshold (BP) 14 and above the pseudofaeces threshold (AP). After 3 months, F and S mussels from each 15 condition were selected to obtain 4 experimental groups: FBP, SBP, FAP and SAP. We 16 hypothesized that the nurturing conditions during the growing period would modify the 17 molecular basis of their growth rate differences. 18 To decipher the molecular mechanisms underlying the growth variation, the gill 19 transcriptomes for the four mussel groups were analysed. Gene expression analysis 20 revealed i) a low number (12) of genes differentially expressed in association with diet 21 and ii) 117 genes differentially expressed by the fastand slow-growing mussels. 22 According to Biological Process GO term analysis transcriptomic differences between 23 the F and S mussels were mainly based on the upregulation of: response to the stimulus, 24 growth and cell activity. Regarding the KEGG terms, carbohydrate metabolism and the 25 Krebs cycle were upregulated in F mussels, whereas biosynthetic processes were 26 upregulated in S mussels. In accordance with their larger gill surface area and higher 27 rates of feeding and growth, the F individuals overexpressed genes in their gill tissues, 28 and these were involved in i) growth (insulin-like growth factors and myostatin); ii) 29 maintenance of the structure and functioning of extracellular matrix (collagen, laminin, 30 This is the accepted manuscript of the article that appeared in final form in Aquaculture 511 : (2019) // Article ID 734242, which has been published in final form at https://doi.org/10.1016/j.aquaculture.2019.734242. © 2019 Elsevier under CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) fibulin and decorin); iii) filtration and ciliary activity (mucin, fibrocystin, dynein and 31 tilB homologue protein genes); iv) aerobic metabolism (citrate synthase and carbonic 32 anhydrase); and v) the immune-system, probably in association with haemocytes. In 33 contrast, S individuals overexpressed a different series of genes pertaining to immune 34 system (leucine-rich repeat protein and galectin), along with genes involved in the 35 response to cellular stress (Heat shock protein (HSP24) and metalloendopeptidase) as 36 well as anaerobic metabolism (C4-dicarboxylate transporter). These results might 37 suggest that S individuals would have a greater prevalence of pathogens/diseases or a 38 higher susceptibility to the pathogens. 39 40 Keywords: Fast-growing, mussel, gill, transcriptome 41 42 1. Introduction 43 Both endogenous factors and environmental conditions influence the growth rate 44 in bivalves (Brown, 1988; Dickie et al. 1984; Mallet and Haley, 1983; Pace et al. 2006; 45 Tamayo et al. 2011). Studies comparing the physiological behaviour between fastand 46 slow-growing individuals have greatly contributed, in the last two decades, to the 47 understanding of the physiological basis of differential growth. The main conclusions 48 were that differences in growth rates resulted from differences in i) the capacity to 49 acquire and absorb food, ii) the efficiency of energy conversion processes and/or iii) the 50 allocation of energy to growth and maintenance (Koehn and Shumway, 1982; Toro and 51 Vergara, 1998; Bayne et al. 1999a, 1999b; Tamayo et al. 2014; Fernández-Reiriz et al. 52 2016). 53 The multilocus heterozygosity hypothesis formulated by Sighn and Zouros 54 (1978) established the existence of a positive correlation between the degrees of 55 heterozygosity and growth rate. Aneuploidy was also demonstrated to play a role in 56 interindividual differences in growth rate in bivalves: significantly higher values of 57 aneuploidy were observed in slow-growing specimens of oysters (Crassostrea gigas) 58 and more recently in the clam Ruditapes decussatus, with a high negative correlation 59 observed between growth rate and aneuploidy percentage (Leitao et al. 2001; Teixeira 60 De Sousa et al. 2011). 61 In recent years, high throughput gene expression analyses have used next-62 generation sequencing (NGS) or microarrays, with the aim of deciphering the 63 underlying mechanisms, and these techniques have allowed the identification of genes 64 involved in growth processes (Gracey et al. 2008; Lockwood et al. 2010; Sussarellu et 65 al. 2010; Devos et al. 2015; Suarez-Ulloa et al. 2015; Xu et al. 2016). For instance, 66 Zhang et al. (2012) reported that collagen and laminins, (extracellular matrix proteins 67 from connective tissue) and fibronectins are involved in the formation of the shell in the 68 oyster Crassostrea gigas. Bassim et al. (2014) analysed the gene expression of the 69 mussel Mytilus edulis during early development (from egg to post-larvae), identified a 70 set of genes related to growth processes in early development (e.g., GATAD1, 71 PIP5K1A and ATRX) and highlighted (Bassim et al. 2015) 29 gene markers related to 72 growth and mortality of bivalve larvae. 73 Very few studies have attempted to specifically analyse the differential gene 74 expression between fastand slow-growing specimens of bivalves. Using different 75 crosses between inbred lines of Crassostrea gigas, Meyer and Manahan (2010) found 76 significant differences between fastand slow-growing larval families in the transcript 77 abundance of ribosomal proteins as well as in the rates of expression of genes encoding 78 for the small cardioactive peptide precursor (ScPB), which is involved in feeding 79 regulation and in several proteins involved in the energy metabolism. Some of them 80 were electron transport components encoding genes (ND4L and ND1), ATP-synthase ȣ, 81 and two coiled-coil-helix-coiled-coil-helix domains (CHCHD2 and CHCHD). More 82 recently, De la Peña et al. (2016) reported the existence of significant differences in the 83 rate of expression of ferritins (Apfer1) between fastand slow-growing individuals of 84 Argopecten purpuratus at different developmental stages (5 stages from embryos to 85 juveniles). Wilson et al. (2016) produced an inbred fast growth line (F) of Mya arenaria 86 clams and analysed the gene expression to test the hypothesis that specific growth-87 related genes will be upregulated in F individuals. These authors established a positive 88 correlation between some metabolic genes (fatty acid synthase and ATPase) with fast 89 growth. These authors also found some upregulated genes involved in structural 90 remodelling in a fast-growing phenotype in agreement with previous studies indicating 91 protein turnover as the main determinant processes for growth heterosis. Finally, 92 Saavedra et al. (2017) concluded that a set of genes controlling tissue and organ growth 93 processes in model organisms (named ‘GCGC’) displayed a minor role in determining F 94 and S in Ruditapes decussatus stocks. However, they found that the insulin-mediated 95 processes had an essential role in interindividual differences in growth rate. 96 Although the available genetic information is increasing (Saavedra and Bachere, 97 2006; Tanguy et al. 2008; Astorga et al. 2014)—e.g., the genome of the oyster 98 Crassostrea gigas was published in 2012 (Zhang et al. 2012)—knowledge regarding 99 molecular and genetic interindividual differences in the growth potential of bivalves 100 remains at low standards. Large-scale sequencing projects (e.g., NGS) have produced 101 large amounts of sequences in databases, but a significant part of these sequences lack 102 an assigned function or similarity. Therefore, a combination of analyses of the 103 transcriptome and other organizational level responses is necessary to understand the 104 roles of specific genes in the functional responses at the level of the whole organism 105 (Bassim et al. 2014). 106 In the present study, we have analysed the gene expression in gill tissue of 107 mussel (Mytilus galloprovincialis) specimens that were selected as fast (F) and slow (S) 108 growers after rearing them for three months in the laboratory under two different 109 nutritional environments. After the rearing period, the physiological components of the 110 Scope for Growth of the selected F and S mussels were recorded under different 111 experimental diets to assess the influence of rearing conditions on the parameters of the 112 physiological behaviours responsible for faster growth (Prieto et al., in preparation). 113 Irrespective of feeding conditions during rearing, faster growers exhibited higher Scope 114 for Growth values that mainly resulted from their increased capacity to acquire food. 115 Indeed, fast growers displayed higher clearance rates, and they consistently were found 116 to have significantly higher gill-surface area per mass unit than their slow-growing 117 counterparts. The combination of higher gill-surface area with higher clearance rate in 118 fast-growing individuals is a phenotypical feature that we have also observed in 119 previous studies performed with mussels (Prieto et al. 2018) and clams (Tamayo et al. 120 2011). 121 Thus, the gill is one of the organs likely playing a major role in determining the 122 interindividual growth rate differences in the mussel Mytilus galloprovincialis. 123 Accordingly, in the present study, we have selected the gill tissue as the target organ to 124 compare gene expression in these groups of fast and slow-growing mussels. The aims of 125 this study were to search for candidate genes for recorded differences in physiological 126 behaviour and, ultimately, in growth, to ascertain biological processes accounting for 127 such differences at the molecular level. Additionally, the effect of the rearing nutritional 128 condition was also considered as a possible modulator of molecular processes 129 underlying the interindividual differences in growth rate. Specifically, emphasis was 130 placed on linking physiological (Prieto et al., in preparation) and transcriptomic results 131 (present study) to achieve a more holistic understanding of the organism behaviour in 132 different growth scenarios. 133 134 135 2. Material and Methods 136 137 2.1. Selection of mussels 138 Some 400 mussels (Mytilus galloprovincialis) of approximately 10 mm shell 139 length (~150 mg live weight) were collected in a rocky shore in Antzoras (Bizcay, 140 North Spain) in February 2014. Once at the lab, we reared each half of the mussels at 141 one of the two “maintenance conditions” (named AP and BP) designed to force 142 different feeding strategies in both groups: a group of 200 mussels was fed a high-143 quality diet (organic content = 80%) dosed at a particle volume concentration of 1.0-1.5 144 mm3/L (below the pseudofaeces threshold; maintenance condition BP), and the other 145 group of 200 mussels was fed a low-quality diet (organic content = 30%) dosed at 146 particle volume concentration of 3.0-3.5 mm3/L (above the pseudofaeces threshold: 147 maintenance condition AP). Diets were a mixture of cultured Isochrysis galbana (T-148 iso), lyophilized Phaeodactylum tricornutum and freshly collected and sieved particles 149 of natural sediment. 150 Shell length was measured with a 0.05 accuracy calliper, and live-weight was 151 determined using a 0.01 mg accuracy balance. After three months, the largest and 152 smallest 24 individuals, representing the percentiles P12.5 and P87.5 in size distribution of 153 each group, were selected as fast (F) and slow (S) growers, respectively. Accordingly, 154 four experimental groups of mussels were obtained combining maintenance (BP and 155 AP) and growth (F or S) conditions: i) fast-growing mussels fed below the pseudofaeces 156 threshold (FBP), ii) slow-growing mussels fed below the pseudofaeces threshold (SBP), 157 iii) fast-growing mussels fed above the pseudofaeces threshold (FAP), and iv) slow-158 growing mussels fed above the pseudofaeces threshold (SAP). The growth rates of the 159 mussels were calculated as GR = the increase in the shell-length or live-weight/elapsed 160 time (days). After the physiological experiments had been completed, the gills of the 161 mussels were dissected and processed for gill surface area determination and RNA 162 extraction. 163 2.2. RNA extraction 164 Gill samples were stored immersed in RNAlater at -80 ºC until the RNA was 165 individually extracted with a ‘RiboPure RNA Purification Kit’ (Ambion kit). The 166 analysis of the quality and integrity of the RNA was checked with Fragment AnalyzerTM 167 Automated CE System equipment from Advanced Analytical with ‘DNF-471 Standard 168 Sensitivity RNA Analysis kit’, (15 nt) and Fragment AnalyzerTM 1.1.0.11 software. The 169 RNA quality was checked using PROSize 2.0. The RNA concentration was measured in 170 the spectrophotometer UV/VIS Nanodrop 1000 (Thermo Fisher). 171 We used 20 individual mussels per experimental condition (20 from FBP, 20SBP, 172 20FAP and 20SAP). The gill RNA was extracted individually. Once extracted, the RNA 173 was quantified according to the method described above. Each individual RNA sample 174 was then diluted to a common concentration of 100 ng/ul. After that, the 20 individual 175 RNA samples per experimental group were randomly combined to create 4 different 176 pools composed of 5 different individual . To create the pools, the same RNA quantity 177 (500 ng) from each of the 5 individuals was added. Once created, the concentration of 178 the pools was quantified in the Nanodrop 1000. Thus, we obtained 16 different pools (4 179 pools from each experimental group x 4 experimental groups), each one containing 180 RNA from 5 different individuals. After that, the 16 pools were marked as described in 181 the section 2.3.1. 182 2.3. Microarray design and hybridization 183 We used a SurePrint G3 Custom microarray (8x60 k) from Agilent to analyse the 184 transcriptome of the gill samples. Microarray probes were designed using Agilent 185 eArray platform, using Mytilus galloprovincialis sequences downloaded from NCBI in 186 February of 2015. Sequences with the best Blastx hit (e-value <10e-10) to unique 187 proteins against nonredundant database were selected. Three nonidentical probes were 188 designed for each sequence. Housekeeping genes (those usually used in Mytilus qPCR 189 analysis) were added as positive controls, alongside default Agilent negative controls. 190 The remaining spots in the array were filled with sequences of the genus Mytilus 191 representing unique proteins (were Mytilus galloprovincialis orthologue was missing). 192 Two probes of the unannotated sequence (one in each reading direction) were included 193 in the array. Thus, the array was based on 17,491 unannotated and 7,806 annotated 194 sequences. The platform is available in gene expression omnibus (GEO) repository with 195 the accession number GPL25650. Hybridization was performed in 16 pools (4 different 196 pools of different 5 individuals per experimental group). Pools were randomly 197 hybridized in the arrays, including at least one pool per experimental group in each 198 array. 199 2.3.1. Marking protocol 200 We used the ‘One-Color Microarray-Based Exon Analysis’ marking protocol 201 from Agilent. Samples were marked using ‘Low Input Quick Amp WT Labeling kit, 202 One-Color’ (p/n 5190-2943) kit. In total, 100 ng of RNA was used for the marking 203 reaction. Marked samples were quantified with a Nanodrop ND-1000 to determine the 204 efficiency of the specific activity of the fluorochromes. 205 2.3.2. Hybridization 206 Samples were manually hybridized with SureHyb Hybridization Chambers 207 (Agilent technologies). Hybridization was conducted in the oven of Agilent 208 Technologies according to the Agilent protocol. The characteristics were as follows: 600 209 ng of marked cRNA, 40 µl volume, 65 ºC temperature, and 20 hours duration at 10 rpm 210 in the hybridization. 211 2.3.3. Scanning 212 The scanning was carried out on the DNA Microarrays G2565CA scanner with 213 ozone–barrier slide covers with the Scan Control Software version 8.5.1., using the 214 default protocol AgilentG3_GX_1Color. The Scanning resolution was 3 µm, the green 215 channel was used, and the size of the resulting Tiff image was 20 bits. 216 2.3.4. Feature extraction 217 We used Agilent Feature Extraction Software (ver. 10.7.3.1) (Agilent 218 Technologies) to process the microarray images and to quantify the fluorescence of the 219 probes. The quality of all arrays was evaluated using the 9 QC-metric parameters 220 generated in the feature extraction. Following this procedure, the processed fluorescence 221 signal (generated by the feature extraction) was obtained. 222 223 2.4. Data treatment 224 Data treatment was carried out in R (v. 3.3.2.) using the limma package (v. 225 3.30.13) from Bioconductor (Ritchie et al. 2015). Probes were prefiltered using 226 gIsPosAndSignif tag; a Boolean value indicating if the signal of the probe exceeds the 227 background signal. Probes with a nonsignificant signal in all the samples of at least one 228 experimental group (n=4) were removed. Background was corrected using normexp 229 method, and normalization between the arrays was performed using the quantile 230 method, as described in Smyth et al. (2002). Fold-change and standard error were 231 estimated by fitting the data to a linear model and an empirical Bayes (eBayes) 232 smoothing was applied to the standard errors. The final gene expression value was the 233 average of the nonidentical probes corresponding to each sequence. Differential 234 expression quantification was based on a logarithmic scale (logFC), the adjusted p-235 value or False Discovery rate (Benjamin–Hochberg method) representing the statistical 236 significance of the observed changes. Probes with FDR<0.05 were considered 237 differentially expressed, as suggested in Cheng and Pounds (2007). Hierarchical 238 clustering (HCL) analysis was performed using dendextend package (v.1.5.2.) to 239 analyse similarity between samples. 240 Normalized hybridization values, as well as the raw data, were deposited in the 241 gene expression omnibus (GEO) repository with the accession number GSE120975. 242 243 244 2.5. Annotation and gene ontology 245 246 Microarray sequences were annotated using Annocript 1.3. against Swiss-Prot 247 and UniRef databases (v. march-2017). Gene Ontology (GO) for three domains 248 (Cellular Component, Molecular Function and Biological Process) was analysed for 249 transcriptome data interpretation, although we focused our analysis mainly on the 250 Biological Process (Suarez-Ulloa et al. 2015). The GO terms list was summarized using 251 REVIGO (Supek et al. 2011). Differentially expressed genes were also mapped to the 252 Kyoto Encyclopaedia of Genes and Genomes (KEGG) database for pathway analysis 253 (Kanehisa, 2002). Conserved protein domains were identified using PROSITE (Sigrist 254 et al. 2009) and NCBI conserved protein domain finder tools. 255 256 257 258 3. Results 259 260 261 3.1. Growth rates of the experimental mussel groups 262 263 After 3 months of maintenance of the mussels under BP or AP conditions, the 264 live weight of F individuals was 2.5-fold higher than that of S individuals, and the shell 265 length was 45% longer. Accordingly, live-weight and shell-length growth rates of F 266 individuals was found to be approximately 3 times greater than that of S individuals in 267 both maintenance conditions (Table 1). 268 269 Table 1. Shell-length (mm), live-weight (g), shell-length growth rate (mm/day) and live-weight growth 270 rate (g/day) (mean values ± SD) of FBP, SBP, FAP and SAP mussel groups. Number of individuals per 271 mussel group = 24 272 273 Mussel group Length (mm) Live weight (g) Growth rate (mm/day) Growth rate (g/day) FBP 21.2 ± 0.7 0.9 ± 0.1 0.146 ± 0.009 0.012 ± 0.002 SBP 13.9 ± 1.2 0.3 ± 0.1 0.055 ± 0.015 0.004 ± 0.001 FAP 21.9 ± 0.6 1.0 ± 0.1 0.144 ± 0.007 0.011 ± 0.001 SAP 15.4 ± 1.0 0.5 ± 0.1 0.060 ± 0.013 0.004 ± 0.001 274 3.2. Quality and reproducibility of the DNA microarray data 275 The marked RNA quality was good in all samples. The yield and the Cyanine 3 276 specific activity were higher than 0.825 µg/reaction and 15 pmol/µg, respectively, in all 277 marked samples. In all cases, the hybridization with the array suited (or passed) the 278 quality standards, evaluated with 9 QC metrics parameters. Only in 0.95% of the probes 279 (568 probes out of 59,539) did the signal have a lower expression value than the 280 background on all samples. For our analysis, we used the probes that had a positive 281 signal on all the samples in at least one experimental group. Mean expression values 282 and standard deviations of the housekeeping genes of the array are shown in the 283 additional file 1. The variability among samples was lower than 3% in 19 of 20 284 housekeeping genes. 285 3.3. Transcript annotation 286 The functional annotation of the genes on the array carried out by Blastx against 287 Swiss-Prot and UniRef databases had 38.8% significant matches (E-value 10-5): 10,001 288 out of 25,781 genes. In total, 27.7% of the annotated genes were matched on 289 4. Discussion 367 In the present study, we analysed the gene expression differences in the gill 368 tissue between fast-growing (F) and slow-growing (S) mussels that were maintained for 369 the long term in the laboratory, while being fed experimental diets of phytoplankton and 370 silt dosed either below (BP) or above (AP) the pseudofaeces threshold. In accordance to 371 what we reported previously from a similar experiment (Prieto et al., 2018), the faster 372 growth of the F mussels (both FBP and FAP) was based on their capacity to display higher 373 clearance rates and higher pre-ingestive selection efficiencies (physiological results will 374 be published elsewhere). Irrespective of the diet fed (BP or AP), increased capacity for 375 water filtration and particle acquisition in F mussels have been found to be coupled with 376 the possession of significantly higher gill-surface areas, a feature of the fast-growing 377 phenotype that we have also found in our previous studies on mussels (Prieto et al. 378 2018) and clams (Tamayo et al. 2011). 379 The experiments of physiological energetics (Prieto et al., in prep) revealed only 380 minor differences in the physiological parameters of mussels fed BP and AP diets. In 381 good agreement with the physiological results, the transcriptomic profiles (HCL results) 382 were very similar between them and only a reduced number of genes were differentially 383 expressed. Three of the annotated DEGs (glutathione S-transferase, headcase protein 384 and protocadherin β) previously have been reported to be upregulated in response to 385 environmental stress and/or bacterial infection in bivalves (Manduzio et al. 2004; Park 386 et al. 2009; Kim et al. 2009; De Zoysa et al. 2011; Li et al. 2018; Rey-Campos et al. 387 2019). However, any interpretation regarding possible differences in the stimulation of 388 immune response in mussels feed below or above the pseudofaeces level is complicated 389 because BP mussels overexpressed glutathione S-transferase, whereas the other two 390 genes were differentially expressed in AP mussels. The remaining annotated DEGs 391 (Notch-regulated ankyrin repeat-containing protein (NRARP), TNFAIP3-interacting 392 protein 2 and FAM60A protein and rRNA 2'-O-methyltransferase fibrillarin) act in 393 several pathways involved in cell differentiation, proliferation, apoptosis and RNA and 394 protein methylation in bivalves such as notch pathway (Bassim et al. 2014), TFG-beta 395 signalling pathway (Wei et al. 2017) and MAP/ERK pathway. Li et al. (2016) reported 396 that TNFAIP3-interacting protein was downregulated in individuals of Chlamys farreri 397 exposed to Benzopyrene and suggested that the reduction in TNFAIP3 was indicative of 398 depressed metabolic rate and hampered progression of mitosis. In the present 399 experiment, overexpression of 4 genes involved in cell proliferation pathways in the 400 mussels that were fed above pseudofaeces level could suggest the existence of an 401 increased gill cell renewal requirement in AP mussels. However, more analysis should 402 be performed to confirm such a hypothesis. 403 The low impact that nutritional condition and feeding mode (below or above the 404 pseudofaeces threshold) exert on the gill transcriptome contrasts with the broad 405 differences associated with the differences in growth rate between F and S specimens: 406 117 differentially expressed genes in gill tissues. The classification of these genes 407 according to Biological Process GO terms indicated that the differences mainly affect 408 responses to stimulus, growth and cellular activity processes. Thus, the GO term 409 findings supported the higher growth rates and activity levels of F individuals in 410 comparison with S mussels. Not surprisingly, previous works on interindividual growth 411 rate differences in bivalves have also described similar GO terms as the main processes 412 underlying growth differences; for instance, Wilson et al. (2016) reported that 19% of 413 GO terms of differentially expressed genes in fast-growing Mya arenaria are associated 414 with cell structure, whereas 17% refer to signalling and growth, 12% to energy and 415 nutrient metabolism and 10% to DNA/RNA and protein synthesis. Regarding the 416 KEGG terms, energetic metabolism terms were referred to F individuals in good 417 correspondence to their higher activity levels. Upregulation of Cofactor and P5P 418 biosynthesis pathway in S individuals seems to involve differences in protein 419 metabolism that could underlie differences in the protein turnover between growth 420 groups, as described in previous studies (Hawkins et al. 1986, 1996). The P5P 421 biosynthesis pathway could either indicate a higher rate of anaerobic metabolism, which 422 in bivalves is based on the utilization of amino acids via opine dehydrogenases, or 423 aspartate-succinate pathway (Hochachka and Somero. 2002) 424 Most differentially expressed (DE) genes between the present F and S 425 individuals lack a clear association to GO terms because the studied model presents 426 only few sequences annotated in the tools allowing performance of the GO analysis. 427 Thus, emphasis has been placed on the individual (rather than the group) analysis of DE 428 genes and their functions to decipher the transcriptomic basis of growth rate differences. 429 430 431 4.1. Upregulated genes in F mussels. 432 Upregulation of growth differential factor-8, also known as myostatin, and 433 insulin-like growth factor in the gill of F mussels would appear meaningfully associated 434 with the higher gill surface area exhibited by fast growers. Myostatin is a negative 435 regulator of muscle growth in vertebrates, and Wang et al. (2010) found that 436 polymorphism of the myostatin gene was correlated with differential growth traits in 437 mammals. In bivalves, myostatin have been suggested to have alternative functions that 438 are related with cell development (Saina and Technau, 2009; Núñez-Acuña and 439 Gallardo-Éscarate, 2014; Morelos et al. 2015; Niu et al. 2015). Insulin-like peptides 440 have been reported to act as growth regulators of soft tissues and shell in bivalves 441 (Taylor et al. 1996; Gricourt et al. 2003), and their roles in determining interindividual 442 growth rate differences in bivalves have been recently suggested by Saavedra et al. 443 (2017), who found a significant overexpression of NOV-like protein in the gills of fast-444 growing Ruditapes decussatus. Using the PROSITE tool on the highly differentially 445 expressed genes (FC>8, FDR<0.01), we have found that, in addition to myostatin and 446 insulin-like peptides, F individuals upregulated an epidermal growth factor-like (EGF). 447 Valenzuela-Miranda et al. (2015) also reported the overexpression of EGFs in the 448 muscle of F specimens in the abalone Haliotis rufescens. EGF is expressed in various 449 tissues of oysters (Sun et al. 2014) and has been suggested to induce cell proliferation 450 and migration during wound healing and to stimulate glycolytic enzymes such as 451 phosphofructokinase and pyruvate kinase (Canesi et al. 2000). 452 In addition to overexpressing growth-regulators, the gills of F mussels 453 overexpressed genes involved in the structure and functionality of the extracellular 454 matrix (ECM), such as collagen, laminin, fibulin and decorin. Some of these genes have 455 been previously reported to be differentially expressed between fastand slow-growing 456 individuals of different invertebrates: For instance, collagen, has been found to be 457 upregulated in F individuals in abalones (Valenzuela-Miranda et al. 2015) and clams 458 (Saavedra et al. 2017). Genomic (Zhang et al. 2012) and transcriptomic (Zhao et al. 459 2012) analyses have suggested that collagen might play an important role in shell 460 formation and soft tissue growth and repair in bivalves. In addition, collagen also 461 appears to play a relevant role in the adhesion and migration of haemocytes to the ECM 462 (Koutsogiannaki and Kaloyianni, 2011) and likely plays a crucial role in the process of 463 cell immunity during inflammatory response (Adams, 2018). Fibulin have been reported 464 to act in association with laminin and collagen in development and biomineralization 465 processes (Timpl et al. 2003; Sleight et al. 2015). Decorin interacts with some growth 466 factors such as EGF, and its binding with myostatin has been described to cause 467 hypertrophy in human muscle cells (Kanzleiter et al. 2014). 468 The F mussels in the present experiment were found to overexpress mucin, the 469 backbone glycoprotein that forms the matrix of the mucus (Espinosa et al. 2016). In 470 bivalves, the filtered particles are retained in the mucus strings circulating through the 471 ciliated groves and transported to the labial palps to be either ingested or rejected as 472 pseudofaeces (Beninger and St-Jean, 1997; Urrutia et al. 2001). A higher putative 473 mucus production in F individuals would be in concordance with their higher clearance 474 rates and higher pre-ingestive selection efficiencies, with both physiological parameters 475 greatly contributing to interindividual differences in the growth rate of mussels (Prieto 476 et al. 2018). Consistent with the higher clearance rates and higher mucin expression, the 477 gills of F mussels also overexpressed fibrocystin, which is involved in tubulogenesis 478 and ciliary activity (Ward et al. 2003), as well as the dynein and tilB homologue protein 479 genes that are involved in the conversion of ATP hydrolysis into mechanical work 480 (Gibbons and Rowe. 1965; Kavlie et al. 2010; Horani et al. 2013). Dynein 481 overexpression in F specimens has also been reported in Haliotis rufescens (Valenzuela-482 Miranda et al., 2015). Overexpression of these genes in F mussels seems to correlate 483 well with the higher filtering activity of the mussels. Recently, Lafont et al. (2019) 484 reported that fibrocystin was one of the upregulated genes in oyster larvae with higher 485 rates of survival to herpes virus (OsHV-1) infection in an experiment that showed 486 transgenerational immune priming in Crassostrea gigas. 487 Processes involved in the metabolic energy supply and ATP turnover are 488 especially relevant to the growth rate of bivalves, and thus, the finding that 2 genes 489 related to energy metabolism were upregulated in F individuals is highly meaningful. 490 Previous approaches to the characterization of genetic differences between fastand 491 slow-growing individuals of different species have emphasized the importance of 492 differential aspects of the energetic metabolism between growth lines. For instance, 493 Meyer and Manahan (2010) found ATP-synthase and two different NADH 494 dehydrogenase subunits upregulated in F individuals of the oyster C. gigas, Wilson et 495 al. (2016) found fatty acid synthase like-1 and fatty acid synthase like-2 genes 496 upregulated in fast-growing individuals of M. arenaria, and Saavedra et al. (2017) 497 found the NADH subunit upregulated in F individuals of the clam Ruditapes 498 decussatus. In the present study, we found citrate synthase (CS) and carbonic anhydrase 499 upregulated in the gill of F individuals. Citrate synthase is a specific marker of aerobic 500 metabolism considered an indicator of the general physiological status of the organism 501 (Garcia-Esquivel et al. 2001, 2002; Pernet et al. 2012; Guévélou et al. 2013) and has 502 been shown to correlate with respiration rates in facultative anaerobes such as intertidal 503 invertebrates (Dahlhoff et al. 2002). Higher citrate synthase expression in our F mussels 504 might thus be indicative of increased energy requirements of gill tissue to sustain higher 505 filtering activity. The carbonic anhydrase enzyme family maintains the pH/salinity 506 balance and favours the exchange of respiratory gases (Breton, 2001) and has been 507 reported to play a role in the process of biomineralization in the mantle tissue (Zhang et 508 al. 2012; Hüning et al. 2016). Finally, DENN domain-containing protein 3, also found 509 to be upregulated in F individuals, is involved in the regulation conversion of inactive 510 GDP-bound to GTP form and vesicle-mediated transport pathways (Marat et al. 2011). 511 We have not found evidence of a DENN domain-containing protein function in 512 bivalves. 513 The gill of F individuals upregulated two genes directly related with the immune 514 system, probably located in the haemocytes: nitric oxide synthase (NOS) and the 515 scavenger receptor MARCO (macrophage receptor with collagenous structure). NOS 516 has been detected in haemocytes of several bivalves (Liu et al. 2018) and produces 517 nitric oxide, a pathogen-killing molecule with broad antiviral and antiparasitic effects 518 (Pautz et al. 2010). In mammals, MARCO is a receptor for bacteria expressed mainly in 519 macrophages; in bivalves, it has been previously reported in Mytilus galloprovincialis 520 (Moreira et al. 2015). 521 4.2. Upregulated genes in S mussels. 522 The gills of the slow-growing mussels overexpress many genes involved in 523 immune, defence and cell stress responses, such as HSP24, leucine-rich repeat proteins, 524 metalloendopeptidase and galectin. The overexpression of heat shock proteins has been 525 commonly found in organisms maintained under temperature stress (Hofman and 526 Somero, 1996; Somero 2012; Lockwood et al. 2013), salinity stress (Zhao et al. 2012) 527 metal exposure (Zhang et al. 2012) and/or bacterial exposure (Genard et al. 2013). In 528 addition, Zhang et al. 2012 found an overexpression of HSP genes in the oyster 529 Crassostrea gigas under various stress conditions (air exposure, thermal stress, salinity 530 stress and metal exposure) and concluded that HSP induction could be a common 531 defence against all stresses in C. gigas. Leucine-rich repeat proteins have been 532 described to be involved in the immunity of invertebrates (Wang et al. 2016) and 533 metalloendopeptidase, seems to be key component of the response against bacterial 534 infections (Miyoshi and Shinoda, 2000). Galectin is probably associated with 535 haemocytes (Espinosa et al. 2016; Vasta et al. 2015) and participates in the recognition 536 of glycans of the surface of virus and bacteria (Nikapitiya et al. 2014). In good 537 correspondence with the present study, Saavedra et al. (2017) also found differences 538 between fastand slow-growing individuals of the clam Ruditapes decussatus in the 539 immune and defence processes of digestive gland and gills. S individuals overexpressed 540 genes involved in immune and defence processes such as defensin and tumour necrosis 541 factor member 11, whereas F individuals were found to overexpress different genes, 542 such as sialic acid-binding lectin and hydramacin-1. They conclude that the observed 543 high differences in the expression of immune and defence genes could reflect a 544 differential fitness among individuals, promoting faster growth rates in those individuals 545 able to fight more efficiently against diseases. In the present study, most of the 546 overexpressed genes in S mussels were found to belong to the immune and defence 547 system and cellular stress, which strongly suggests a greater prevalence of 548 pathogens/diseases or a higher susceptibility to the pathogens. As suggested by Genard 549 et al. (2013), when analysing the physiological response of C. gigas larvae submitted to 550 bacterial infection, extra investments in supporting defence mechanisms might drain 551 energy resources from normal processes in healthy organisms, resulting in reduced 552 feeding and growth performances. 553 In addition, the strong upregulation of countin-1 (FC≈4), a cell-counting factor 554 that limits the maximum size of the multicellular structure by the downregulation of the 555 cell adhesion mediator gp24, seems to indicate developmental process inhibition in S 556 individuals. Symptoms of impairment in the respiratory function of the gill affecting 557 aerobic ATP production are also evident in S mussels: Evidence of increased use of 558 anaerobic metabolic pathways includes strong upregulation of anaerobic C4-559 dicarboxylate transporter (FC ≈ 8), as well as the increased biosynthesis of pyridoxal-5 560 –phosphate. Similarly, Saavedra et al. (2017) have reported upregulation in the 561 digestive gland of S clams of enzymes very likely involved in anaerobic metabolism 562 (e.g., malate dehydrogenase and glycerol-3-phosphate dehydrogenase). 563 564 4.3. Conclusions and prospects 565 The present results show the existence of substantial differences in the 566 transcriptome of the gills of F and S individuals. The gills of the fast-growing mussels 567 overexpressed growth factors and genes that are involved in the maintenance of relevant 568 cellular functions, such as the maintenance of the ciliary activity, the development of a 569 robust extracellular matrix contributing to antibacterial defence and the maintenance of 570 aerobic metabolic pathways. This transcriptomic profile in the F mussels suggests that 571 the gills are well equipped to maintain higher filtering activities that enable fast-572 growing mussels to maximize food acquisition and sustain fast growth rates. In contrast, 573 slow-growing mussels overexpress genes involved in the immune system and genes that 574 participate in cellular-stress responses and anaerobic metabolic pathways. These results 575 could suggest that S individuals would have a greater prevalence of pathogens/diseases 576 or a higher susceptibility to the pathogens. Further analysis with different organs (e.g., 577 digestive gland) are needed to obtain a holistic view of the transcriptomic basis of fast-578 growth in bivalves; however, the present study suggests that the immune response might 579 be a crucial component of the interindividual differences in growth rate in Mytilus 580 galloprovincialis mussel spats. 581 582 5. Acknowledgements 583 584 This study was funded through the project AGL2013-49144-C3-1-R of the Spanish 585 Ministry of Economy and Competitiveness. D. Prieto was funded by an FPI grant from 586 the Basque Government. The authors are thankful for the technical and human support 587 provided by SGIker of UPV/EHU (European funding: ERDF and ESF). Finally, D. 588 Prieto especially wants to thank Dr. A. Huvet for all the help provided. 589 590 6. References 591 592 Adams, J. C. (2018). Matricellular Proteins: Functional Insights From Non-mammalian Animal Models. 593 In Current topics in developmental biology (Vol. 130, pp. 39-105). Academic Press. 594 Astorga, M. P. (2014). Genetic considerations for mollusk production in aquaculture: current state of 595 knowledge. Frontiers in genetics, 5,435. 596 Bassim, S., Tanguy, A., Genard, B., Moraga, D., & Tremblay, R. (2014). Identification of Mytilus edulis 597 genetic regulators during early development. Gene, 551(1), 65-78. 598 Bassim, S., Chapman, R. W., Tanguy, A., Moraga, D., & Tremblay, R. (2015). Predicting growth and 599 mortality of bivalve larvae using gene expression and supervised machine learning. Comparative 600 Biochemistry and Physiology Part D: Genomics and Proteomics, 16, 59-72. 601 Bayne, B. L., Svensson, S., & Nell, J. A. (1999)a. The physiological basis for faster growth in the Sydney 602 rock oyster, Saccostrea commercialis. The Biological Bulletin, 197(3), 377-387. 603 Bayne, B. L., Hedgecock, D., McGoldrick, D., & Rees, R. (1999)b. Feeding behaviour and metabolic 604 efficiency contribute to growth heterosis in Pacific oysters [Crassostrea gigas (Thunberg)]. 605 Journal of experimental marine biology and ecology, 233(1), 115-130. 606 Beninger, P. G., & St-Jean, S. D. (1997). The role of mucus in particle processing by suspension-feeding 607 marine bivalves: unifying principles. Marine Biology, 129(2), 389-397. 608 Breton, S. (2001). The cellular physiology of carbonic anhydrases. Jop, 2(4 Suppl), 159-164. 609 Brown, J. R. (1988). Multivariate analyses of the role of environmental factors in seasonal and site-related 610 growth variation in the Pacific oyster Crassostrea gigas. Marine Ecology Progress Series, 225-611 236. 612 Canesi, L., Ciacci, C., Betti, M., & Gallo, G. (2000). Growth factor-mediated signal transduction and 613 redox balance in isolated digestive gland cells from Mytilus galloprovincialis Lam. Comparative 614 Biochemistry and Physiology Part C: Pharmacology, Toxicology and Endocrinology, 125(3), 615 355-363. 616 Cheng, C., & Pounds, S. (2007). False discovery rate paradigms for statistical analyses of microarray 617 gene expression data. Bioinformation, 1(10), 436. 618 Dahlhoff, E. P., Stillman, J. H. and Menge, B. A. (2002). Physiological community ecology: variation in 619 metabolic activity of ecologically important rocky intertidal invertebrates along environmental 620 gradients. Integrative and Comparative Biology,. 42, 862-71. 621 De la Peña, T. C., Cárcamo, C. B., Díaz, M. I., Brokordt, K. B., & Winkler, F. M. (2016). Molecular 622 characterization of two ferritins of the scallop Argopecten purpuratus and gene expressions in 623 association with early development, immune response and growth rate. Comparative 624 Biochemistry and Physiology Part B: Biochemistry and Molecular Biology, 198, 46-56. 625 De Zoysa, M., Nikapitiya, C., Oh, C., Lee, Y., Whang, I., Lee, J. S., Choi, C.Y., & Lee, J. (2011). 626 Microarray analysis of gene expression in disk abalone Haliotis discus discus after bacterial 627 challenge. Fish & shellfish immunology, 30(2), 661-673. 628 Devos, A., Dallas, L. J., Voiseux, C., Lecomte-Pradines, C., Jha, A. N., & Fiévet, B. (2015). Assessment 629 of growth, genotoxic responses and expression of stress related genes in the Pacific oyster 630 Crassostrea gigas following chronic exposure to ionizing radiation. Marine pollution bulletin, 631 95(2), 688-698. 632 Dickie, L. M., Boudreau, P. R., & Freeman, K. R. (1984). Influences of stock and site on growth and 633 mortality in the blue mussel (Mytilus edulis). Canadian Journal of Fisheries and Aquatic 634 Sciences, 41(1), 134-140 635 Espinosa, E. P., Koller, A., & Allam, B. (2016). Proteomic characterization of mucosal secretions in the 636 eastern oyster, Crassostrea virginica. Journal of proteomics, 132, 63-76. 637 Fernández-Reiriz, M. J., Irisarri, J., & Labarta, U. (2016). Flexibility of Physiological Traits Underlying 638 Inter-Individual Growth Differences in Intertidal and Subtidal Mussels Mytilus galloprovincialis. 639 PLoS One, 11(2), e0148245. 640 Garcıa-Esquivel, Z., Bricelj, V. M., & González-Gómez, M. A. (2001). Physiological basis for energy 641 demands and early postlarval mortality in the Pacific oyster, Crassostrea gigas. Journal of 642 Experimental Marine Biology and Ecology, 263(1), 77-103. 643 Garcıa-Esquivel, Z., Bricelj, V. M., & Felbeck, H. (2002). Metabolic depression and whole-body 644 response to enforced starvation by Crassostrea gigas postlarvae. Comparative Biochemistry and 645 Physiology Part A: Molecular & Integrative Physiology, 133(1), 63-77. 646 Genard, B., Miner, P., Nicolas, J. L., Moraga, D., Boudry, P., Pernet, F., & Tremblay, R. (2013). 647 Integrative study of physiological changes associated with bacterial infection in Pacific oyster 648 larvae. PLoS One, 8(5), e64534. 649 Gibbons, I. R., & Rowe, A. J. (1965). Dynein: a protein with adenosine triphosphatase activity from cilia. 650 Science, 149(3682), 424-426. 651 Gracey, A. Y., Chaney, M. L., Boomhower, J. P., Tyburczy, W. R., Connor, K., & Somero, G. N. (2008). 652 Rhythms of gene expression in a fluctuating intertidal environment. Current Biology, 18(19), 653 1501-1507. 654 Gricourt, L., Bonnec, G., Boujard, D., Mathieu, M., & Kellner, K. (2003). Insulin-like system and growth 655 regulation in the Pacific oyster Crassostrea gigas: hrIGF-1 effect on protein synthesis of mantle 656 edge cells and expression of an homologous insulin receptor-related receptor. General and 657 comparative endocrinology, 134(1), 44-56. 658 Guévélou, E., Huvet, A., Sussarellu, R., Milan, M., Guo, X., Li, L., ... & Boudry, P. (2013). Regulation of 659 a truncated isoform of AMP-activated protein kinase α (AMPKα) in response to hypoxia in the 660 muscle of Pacific oyster Crassostrea gigas. Journal of Comparative Physiology B, 183(5), 597-661 611. 662 Hawkins, A. J. S., Bayne, B. L., & Day, A. J. (1986). Protein turnover, physiological energetics and 663 heterozygosity in the blue mussel, Mytilus edulis: the basis of variable age-specific growth. 664 Proceedings of the Royal Society of London B: Biological Sciences, 229(1255), 161-176. 665 Hawkins, A. J. S., Smith, R. F. M., Bayne, B. L., & Heral, M. (1996). Novel observations underlying the 666 fast growth of suspension-feeding shellfish in turbid environments: Mytilus edulis. Marine 667 Ecology Progress Series, 179-190. 668 Hochachka, P. W., & Somero, G. N. (2002). Biochemical adaptation: mechanism and process in 669 physiological evolution. New York: Oxford University Press. 670 Hofmann, G. E., & Somero, G. N. (1996). Interspecific variation in thermal denaturation of proteins in 671 the congeneric mussels Mytilus trossulus and M. galloprovincialis: evidence from the heat-shock 672 response and protein ubiquitination. Marine Biology, 126(1), 65-75. 673 Horani, A., Ferkol, T. W., Shoseyov, D., Wasserman, M. G., Oren, Y. S., Kerem, B., ... & Elpeleg, O. 674 (2013). LRRC6 mutation causes primary ciliary dyskinesia with dynein arm defects. PLoS One, 675 8(3), e59436. 676 Hüning, A. K., Lange, S. M., Ramesh, K., Jacob, D. E., Jackson, D. J., Panknin, U., ... & Melzner, F. 677 (2016). A shell regeneration assay to identify biomineralization candidate genes in mytilid 678 mussels. Marine genomics, 27, 57-67. 679 Kanehisa M. 2002. The KEGG database. In: Novartis Foundation Symposium, 247:91–101 680 Kanzleiter, T; Rath, M; Görgens, SW; Jensen, J; Tangen, DS; Kolnes, AJ; Kolnes, KJ; Lee, S; Eckel, J; 681 Schürmann, A; Eckardt, K (2014). "The myokine decorin is regulated by contraction and 682 involved in muscle hypertrophy". Biochem Biophys Res Commun. 450: 1089–1094.. 683 Kavlie, R. G., Kernan, M. J., & Eberl, D. F. (2010). Hearing in Drosophila requires TilB, a conserved 684 protein associated with ciliary motility. Genetics, 185(1), 177-188. 685 Kim, M., Ahn, I. Y., Cheon, J., & Park, H. (2009). Molecular cloning and thermal stress-induced 686 expression of a pi-class glutathione S-transferase (GST) in the Antarctic bivalve Laternula 687 elliptica. Comparative Biochemistry and Physiology Part A: Molecular & Integrative 688 Physiology, 152(2), 207-213. 689 Koehn, R. K., & Shumway, S. E. (1982). Genetic/physiological explanation for differential growth rate 690 among individuals of the American oyster, Crassostrea virginica (Gmelin). Marine Biology 691 Letters, 2(3), 35-42. 692 Koutsogiannaki, S., & Kaloyianni, M. (2011). Effect of 17β-estradiol on adhesion of Mytilus 693 galloprovincialis hemocytes to selected substrates. Role of alpha2 integrin subunit. Fish & 694 shellfish immunology, 31(1), 73-80. 695 Lafont, M., Goncalves, P., Guo, X., Montagnani, C., Raftos, D., & Green, T. (2019). Transgenerational 696 plasticity and antiviral immunity in the Pacific oyster (Crassostrea gigas) against Ostreid 697 herpesvirus 1 (OsHV-1). Developmental & Comparative Immunology, 91, 17-25. 698 Leitao, A., Boudry, P., & Thiriot-Quievreux, C. (2001). Negative correlation between aneuploidy and 699 growth in the Pacific oyster, Crassostrea gigas: ten years of evidence. Aquaculture, 193(1), 39-700 48. 701 Li, Z., Cha, Y., Hu, B., Wen, C., Jian, S., Yi, P., & Gang, Y. (2018). Identification and characterization of 702 two distinct sigma-class glutathione-S-transferase from freshwater bivalve Cristaria plicata. 703