SUMMARY
Growth rates in animals are governed by a wide range of biological factors, many of which remain poorly understood. To identify the genes that establish growth differences in bivalve larvae, we compared expression patterns in contrasting phenotypes (slow- and fast-growth) that were experimentally produced by genetic crosses of the Pacific oyster Crassostrea gigas. Based on transcriptomic profiling of 4.5 million cDNA sequence tags, we sequenced and annotated 181 cDNA clones identified by statistical analysis as candidates for differential growth. Significant matches were found in GenBank for 43% of clones (N=78), including 34 known genes. These sequences included genes involved in protein metabolism, energy metabolism and regulation of feeding activity. Ribosomal protein genes were predominant, comprising half of the 34 genes identified. Expression of ribosomal protein genes showed non-additive inheritance — i.e. expression in fast-growing hybrid larvae was different from average levels in inbred larvae from these parental families. The expression profiles of four ribosomal protein genes (RPL18, RPL31, RPL352 and RPS3) were validated by RNA blots using additional, independent crosses from the same families. Expression of RPL35 was monitored throughout early larval development, revealing that these expression patterns were established early in development (in 2-day-old larvae). Our findings (i) provide new insights into the mechanistic bases of growth and highlight genes not previously considered in growth regulation, (ii) support the general conclusion that genes involved in protein metabolism and feeding regulation are key regulators of growth, and (iii) provide a set of candidate biomarkers for predicting differential growth rates during animal development.
INTRODUCTION
Many of the basic mechanisms underlying the changes in shape and form that occur during animal development are well characterized (Davidson, 1986; Wolpert, 1994; Dow, 2007). Molecular biological techniques have allowed researchers to characterize in great detail the networks of interacting genes that underlie developmental changes in morphology (Davidson et al., 2002; Davidson and Erwin, 2006). In contrast, the regulation of developmental changes in size (i.e. growth) and other physiological processes are not as well understood, despite the obvious importance of these processes (Conlon and Raff, 1999; Nijhout et al., 2006).
The phenomenon of growth heterosis (Shull, 1948) (‘hybrid vigor’) provides phenotypic contrasts that have been used to study differential growth in many animal and plant species (Shull, 1948; Chauhan and Singh, 1982; Strauss, 1986; Griffing, 1990; Gregory et al., 1991; Bentsen et al., 1998). Despite the obvious advantages for commercial applications, as evident in the 5-fold increase in US corn production resulting from the introduction of hybrid corn (USDA, 2006), many aspects of growth heterosis remain poorly understood. The genetic hypothesis of dominance (the masking of deleterious alleles from one parent through complementation by alleles from the other parent) is widely accepted (Roff, 2002), but other evidence supports the overdominance hypothesis, in which heterozygosity at certain loci confers innate fitness benefits (Crnokrak and Barrett, 2002). Presently, the biological basis of growth heterosis still remains unresolved.
Bivalve molluscs have provided a useful model organism for studying growth heterosis in animals (Singh and Zouros, 1978; Zouros et al., 1988; Bayne et al., 1999). In natural populations of adult bivalves, growth rates are positively correlated with the degree of multi-locus heterozygosity (Koehn and Shumway, 1982). More recent studies have shown that growth heterosis can be experimentally produced in bivalve larvae of the Pacific oyster Crassostrea gigas (Pace et al., 2006; Hedgecock et al., 2007). This species is of interest because it has a life-history strategy typical of high-fecundity marine invertebrates and offers the advantage of using established genetic lines that can be crossed to produce larvae with reproducibly different growth phenotypes. Certain larval families of C. gigas grow up to 5-times faster than other families (Pace et al., 2006), which is comparable to the growth advantage reported for hybrid corn (Betran et al., 2003). In adult bivalves, growth heterosis has been variously attributed to differences in ingestion or assimilation rates, energy allocation, or resting metabolic rates (Koehn and Shumway, 1982; Hawkins et al., 1986; Griffing, 1990; Bayne, 1999; Bayne, 2004b). For bivalve larvae, Pace and colleagues (Pace et al., 2006) have shown that a complex set of physiological processes regulate differences in genetically determined growth rates, including differential feeding rates and protein metabolism. These comparisons between fast- and slow-growing families of larvae provide clear experimental advantages for understanding the mechanisms of growth heterosis during development.
The recent application of high-throughput sequencing technologies to the study of growth heterosis has brought new kinds of data to bear on these questions of growth regulation. A set of genes differentially expressed in association with growth heterosis has been identified in hybrid corn (Song and Messing, 2003). Transcriptome analysis of fast-growing hybrid corn and wheat families has revealed non-additive patterns of gene expression reminiscent of the non-additive phenotype of growth heterosis (Wu et al., 2003; Auger et al., 2005; Guo et al., 2006; Swanson-Wagner et al., 2006). Comparable data for animals are scarce, but transcriptome comparisons in hybrid fruit flies have also revealed non-additive gene expression patterns (Gibson et al., 2004).
The use of high-throughput DNA sequencing technologies, such as massively parallel signature sequencing (MPSS) (Brenner et al., 2000a; Brenner et al., 2000b), makes it possible to characterize gene expression profiles of organisms for which no prior genome-wide sequence data are available (e.g. C. gigas and many other species). We recently extended this method to describe quantitative patterns of gene expression associated with growth heterosis in larvae of C. gigas (Hedgecock et al., 2007). From an analysis of 4.5 million cDNA sequence tags, ~23,000 distinct signatures were identified, of which ~350 were candidates for growth heterosis in larvae. In the current study, our goal was to identify genes and processes associated with differential growth rates during animal development. To that end, we cloned and sequenced a set of growth heterosis candidate genes to aid with identification of their physiological functions. Additionally, we confirmed the expression profiles of selected genes and evaluated the reproducibility of these profiles within each larval family using an independent set of crosses produced from the same genetic families. Our findings establish a connection between the physiology of differential growth for contrasting phenotypes (Pace et al., 2006) and whole-transcriptome analysis of differential gene expression profiles (Hedgecock et al., 2007) during early animal development.
MATERIALS AND METHODS
Experimental crosses and culture
The same genetic families (lines 3 and 5) of the Pacific oyster C. gigas Thunberg 1793 that were used in our previous studies (Pace et al., 2006; Hedgecock et al., 2007) were used in a reciprocal cross designed to produce larval families with four different genotypes: 3×3, 3×5, 5×3 and 5×5 [genotype names expressed as paternal line × maternal line (‘sire × dam’)]. A series of 200 l cultures were each stocked with 10 larvae ml−1 and maintained at 25°C in 0.2 μm (pore-size) filtered seawater that was replaced every 4 days. Larvae were fed 30,000 cells ml−1 with the alga Isochrysis galbana at 2 day intervals. All cultures were mixed by gentle aeration. Larvae were collected from cultures at each sampling interval, beginning at 2 days post-fertilization. Growth rates were measured as increases in shell length, measured using a calibrated ocular micrometer under a microscope (N=50 larvae measured for each sample). These data were used to calculate growth rates for the different larval families by statistical regression of average size (shell length) against age. Duplicate independent crosses (i.e. using gametes from different individual parents belonging to the same genetic lines) were generated as described above for analysis of growth.
RNA extraction
Samples containing known numbers of larvae were collected for RNA analysis. Each sample was centrifuged and the seawater supernatant aspirated. Larvae were homogenized in denaturing solution (4 mol l−1 guanidinium thiocyanate, 25 mmol l−1 sodium citrate pH 7.0, 0.1 mol l−1 β-mercaptoethanol and 0.5% sodium sarcosyl) using a mechanical (rotor-stator) homogenizer, and immediately frozen by immersion in liquid nitrogen. RNA was extracted using the RNAEasy kit (Qiagen, Valencia, CA, USA) according to the manufacturer's instructions, with the following modifications. Samples of homogenized tissue were thawed on ice and centrifuged for 5 min to remove insoluble material (e.g. shell fragments). The resulting supernatants were diluted in a 1:1 ratio with RLT buffer (provided by the manufacturer). RNA was eluted in 1× MOPS/RNasin (40 mmol l−1 3-[n-morpholino]-propanesulfonic acid, 10 mmol l−1 sodium acetate, 1 mmol l−1 EDTA, 100 U ml−1 RNasin, pH 7.0). RNA was quantified based on absorbance at 260 nm, and precipitated with ethanol where necessary to achieve the required concentrations for subsequent analysis.
Probe synthesis
Four cDNA clones associated with protein metabolism were selected to measure transcript abundance during development of embryos and larvae with different genotypes and growth rates. Clones 68, 76, 124 and 278, tentatively identified as ribosomal proteins RPL35, RPL31, RPL18 and RPS3, respectively, were used as templates to generate radiolabeled (32P) probes specific for each transcript. Template fragments were prepared through restriction digests of pCR2.1-Topo constructs with EcoRI, and gel-extracted using the Qiaquick gel extraction kit (Qiagen). Template preparations were quantified by absorbance (OD260). A size marker template (Millenium Marker Template; Ambion, Austin, TX, USA) was used to generate probes for hybridization and these size markers were loaded in each gel. Random-primed radiolabeled cDNA probes were synthesized using the Prime-A-Gene kit (Promega, Madison, WI, USA) with 1.85 MBq 32P-dCTP per reaction to radiolabel the probes (Perkin-Elmer, Wellesley, MA, USA). Following synthesis and purification of radiolabeled probes, the radioactivity of each probe was measured by liquid scintillation counting to ensure equal loading of probes into hybridization reactions containing equal amounts of RNA.
RNA (northern) blots
Blots were prepared using RNA from all four larval families. The expression of four ribosomal protein genes (RPL18, RPL31, RPL35 and RPS3) was measured at 6 days post-fertilization. This time period was chosen as previous studies have shown that genotype-dependent differential growth can be statistically quantified in 6-day-old larvae (Pace et al., 2006; Hedgecock et al., 2007). The expression of RPL35 was monitored throughout development (1-6 days post-fertilization). For all blots, 5 μg of total RNA was separated by electrophoresis under denaturing conditions (1% agarose, 6% formaldehyde, 1× MOPS buffer) (Ausubel et al., 1994). The intensity of 18S rRNA bands was measured by staining with ethidium bromide and quantification of digital photographic images using ImageJ (NIH) (Abramoff et al., 2004). All radiolabeled-probe band densities were later normalized to this measure of RNA gel loading. RNA was transferred onto nylon membranes (Brightstar-Plus; Ambion) overnight in 4× SSC (750 mmol l−1 sodium chloride, 75 mmol l−1 sodium acetate, pH 7.0) according to standard downward capillary transfer methods (Ausubel et al., 1994). RNA was cross-linked to membranes by exposure to UV light (Stratalinker, Stratagene, La Jolla, CA, USA), and membranes were stored at −80°C to await further analysis.
RNA blots were pre-hybridized at 42°C for 1 h in Ultrahyb hybridization solution (Ambion). Probes were added at a final radioactivity of 17 kBq ml−1; molecular size marker probes were added at a final activity of 0.8 kBq ml−1. Hybridizations were conducted overnight at 42°C in a rotary hybridization oven (Bambino; Midwest Scientific, St Louis, MO, USA). Un-hybridized radioactive material was removed from blots by two, 5 min washes (300 mmol l−1 sodium chloride, 30 mmol l−1 sodium acetate, 0.1% sodium dodecyl sulfate, pH 7.0), followed by two additional 15 min washes (15 mmol l−1 sodium chloride, 1.5 mmol l−1 sodium acetate, 0.1% sodium dodecyl sulfate). These washed blots were exposed to PhosphorImager imaging plates (Amersham Biosciences, Piscataway, NJ, USA) and the resulting images digitized with an FX Molecular Imager (Bio-Rad, Hercules, CA, USA). Band densities were quantified using ImageJ and normalized to the 18S rRNA band density to standardize to RNA amounts loaded on each gel.
Analyses of candidate growth genes
A set of short cDNA sequences associated with family-specific differences in growth rate was obtained in our previous study (Hedgecock et al., 2007). In that study, gene expression was profiled in fast- and slow-growing families of larvae by cloning and high-throughput sequencing of 4.5 million cDNA molecules using MPSS (Brenner et al., 2000a; Brenner et al., 2000b). Comparisons between the expression profiles of fast- and slow-growing larvae were used to select a set of candidate sequences, based on statistical procedures that we fully described previously (Hedgecock et al., 2007). The candidates identified through that process included four patterns of non-additive expression: overdominant (OD), i.e. expressed at higher levels in hybrids (3×5, 5×3) than in either of the parental lines (3×3, 5×5); underdominant (UD), expressed at lower levels in hybrids than in either of the parental lines; dominant-high (D+), expressed in hybrids at equivalent levels to the higher-expressing parent; and dominant-low (D−), expressed in hybrids at equivalent levels to the lower-expressing parent (details in Results section). Of the 4.5 million cDNA sequence tags originally analyzed by MPSS, statistical contrasts suggested 349 candidates for growth heterosis (Hedgecock et al., 2007). In the current study we have focused on a subset of candidates that were expressed at high levels and were detectable in the majority of families.
The cloning procedure used in the previous study (Hedgecock et al., 2007) produces cDNA clones whose 5′-ends begin at the 3′-most DpnII restriction site (GATC) of their respective transcript, an approach widely used in transcriptome profiling (e.g. tag profiling with the Illumina Genome Analyzer; http://www.illumina.com). Each of the resulting cDNA clones represented a fragment delimited by the 3′-most DpnII site in that transcript at the 5′-end of the fragment, and the poly-A tail at the 3′-end. Each 17 bp signature sequence obtained from MPSS corresponded to the 5′-end of a particular cDNA clone, enabling the specific amplification and sub-cloning of these candidates and their subsequent annotation based on sequence homology.
Amplification, cloning and sequencing
cDNA clones were PCR amplified from the pLCV plasmid vector of the library (Brenner et al., 2000b) using signature-specific sense primers 5′-GACCG[N17]-3′, where N17 represents the signature sequence from each clone, and a vector-specific M13F primer (TGTAAAACGACGGCCAGT). The specificity of each reaction was evaluated by gel electrophoresis of PCR products, and the products purified by gel extraction using the Qiaquick gel extraction kit. Fragments were cloned into pCR2.1-Topo using the Topo TA cloning kit (Invitrogen, Carlsbad, CA, USA). Multiple transformants (N=3-10) from each PCR reaction were screened using restriction digestion (EcoRI) and capillary DNA sequencing (Beckman CEQ capillary sequencer, Beckman Coulter, Inc., Fullerton, CA, USA). Correct constructs were identified based on the presence of both the 17 bp signature sequence (5′) and a 3′ sequence specific to the vector used for the MPSS cloning process (AAAAAAAAAAAAAAAAGAATCGTCTCAATGCGGGC). Selected clones were sequenced repeatedly for accuracy. High-quality consensus sequences were obtained from multiple sequence alignments using the VectorNTI software package (Invitrogen). Consensus sequences from these alignments were trimmed to remove vector sequences, producing cDNA sequences bounded by the poly-A tail at the 3′-end and the 3′-most DpnII site at the 5′-end, which formed the basis for all subsequent sequence analysis.
Database searches and sequence analyses
To identify each clone, cDNA sequences were compared with the National Center for Biotechnology Information (NCBI) sequence databases GenBank and EST using TBLASTX (Altschul et al., 1997) with a significance threshold of e-value<10−3. BLAST reports were parsed using the BioPerl Search::IO module (Stajich et al., 2002). Annotation information for each identified sequence was obtained from NCBI and gene product names extracted using the BioPerl Seq::IO module. On the basis of these analyses, candidate genes were grouped into three classes: clones showing no similarity to existing sequences, clones similar to ESTs that lacked annotation, and clones similar to annotated coding sequences. The tentative identities of 37 clones were assigned based on the gene name annotation of the most similar subjects (based on e-values) returned from these analyses. To assess the biological significance of the differential expression of these candidate genes in fast-growing larvae, Gene Ontology biological process terms (Ashburner et al., 2000) were assigned to candidates using the Java applet Blast2GO (Conesa et al., 2005). Based on these results, biological process terms were assigned for 24 of the 37 candidates identified by sequence similarity.
RESULTS
Growth rates
Genotype-dependent differences in growth rate were evident for the larval families of C. gigas reared under similar environment conditions of food and temperature (Fig. 1A). Measurements of larval size (shell length, in μm) were used to calculate growth rate from regressions of shell length against age (Fig. 1A: comparison of regressions by ANOVA, P<0.001). Growth rates for duplicate larval cultures of each family (Fig. 1A) were not significantly different (P>0.05), so these were pooled for further analysis. Pairwise comparisons were performed between families by ANOVA, with Bonferroni's adjustment for multiple tests. Hybrid larvae grew faster than inbred larvae from either family 3×3 (P<0.001) or 5×5 (P<0.001). Hybrid larvae grew at 9.6 μm day−1 on average, with no significant difference between hybrid larval families (P=0.90). The 3×3 larval family grew significantly faster than the 5×5 larval family (7.2 and 5.6 μm day−1, respectively; P<0.001).
In agreement with the above comparison of regressions, a comparison of size-at-age for 6-day-old larvae revealed that hybrid larvae were larger than larvae from their corresponding parental families (Fig. 1B). Analysis of variance revealed significant differences in mean shell length between all four families of 6-day-old larvae (P<0.001). Pairwise comparisons between these families revealed no differences in size between the two hybrid larval families (Tukey's HSD; P>0.05), but showed that the size of hybrid larvae was 14% greater than for inbred larvae (Tukey's HSD; P<0.05). These analyses demonstrate differential growth phenotypes among the larval families used in the current study and confirm the reproducibility of these familial growth phenotypes across four generations.
Isolation of candidate gene clones
The list of 188 candidates selected for analysis includes representatives from all four modes of non-additive expression (OD, D+, D− and UD). Representative examples from each of these expression patterns are shown in Fig. 2, and an overview of the selection process is shown in Table 1. Seven of the 188 candidate genes failed to amplify by PCR and were therefore excluded from further analysis. A complete list of the 181 candidates characterized in this study, which includes accession numbers for cDNA sequences and MPSS expression data, is shown in supplementary material Table S1. These selected candidate cDNAs were successfully PCR-amplified, cloned, and sequenced repeatedly to obtain high-quality consensus nucleotide sequences. The lengths of these clones were determined from the distance between the 17 bp signature sequence at the 5′-end and the cloning vector adaptor at the 3′-end. Insert sizes ranged from 38 to 855 bp (mean 219 bp), similar to the value expected based on a random distribution of the target 4 bp DpnII restriction site GATC used in the cloning process (i.e. 1 site per 256 bp) (Brenner et al., 2000a; Brenner et al., 2000b). Most clones (N=104) were between 100 and 300 bp in size, with the remaining 19% (N=34) <100 bp, and 24% (N=43) >300 bp. Canonical polyadenylation signals (AATAAA or ATTAAA) were detected between 10 and 30 bp upstream from the poly-A tail in 76% of these clones (N=137), supporting the conclusion that the cloned inserts analyzed represent 3′-fragments of transcripts. The consensus sequences for all these clones for C. gigas were deposited in GenBank. Small clones (<50 bp) and those lacking BLAST matches were deposited in the EST database (accession numbers EW688558-EW688566, EX151492-EX151622), and the annotated clones in the high-throughput cDNA (HTC) database (accession numbers EU152921-EU152961).
Identification of candidate growth genes
The 181 candidate sequences were compared with public databases and significant (e-value≤10−3) matches were found for 43% (N=78) (Fig. 3A). For 22% of clones (N=41), all BLAST matches lacked annotated protein-coding regions and so were not informative for gene identification (e.g. clones that matched only ESTs, ribosomal RNA, or other non-coding sequences). Putative gene names were assigned to 20% of clones (N=37) based on coding sequence annotation of BLAST hits (Table 2) using the best numerical match to a GenBank annotated gene (i.e. the lowest e-value) for each clone. A total of 37 clones were mapped to 34 genes by this process; each of three ribosomal protein genes (RPL7a, RPL37a and RPS17) was matched by two different cDNA clones (Table 2). These tentatively identified candidate growth genes are herein referred to by their annotated gene names.
Ribosomal protein genes were the most abundant class of the clones analyzed, comprising 50% of all candidate growth genes identified in this study (N=17; Fig. 3B). The remaining candidate genes were distributed among several gene families (complete list in Table 2). In addition to the 17 ribosomal protein genes, several genes associated with other aspects of protein metabolism were identified, including DC2, peptidylprolyl isomerase (PPIB), and the proteosome subunit PSMD14. An additional five candidate genes were associated with energy metabolism, including mitochondrial genes (NADH dehydrogenase subunits ND1 and ND4L) and nuclear genes with mitochondrial functions (ATP synthase δ-subunit OSCP, and the putative mitochondrial components with coiled-coil-helix domains, CHCHD 2 and CHCHD3). Another candidate gene identified is the small cardioactive peptide precursor SCPb that has been implicated in the regulation of feeding activity in molluscs (see Discussion).
The likely biological processes associated with these candidate growth genes were determined using Gene Ontology annotation. A total of 24 clones were successfully assigned to biological process terms (Fig. 3C). Of the clones for which Gene Ontology terms could be assigned, 63% (N=15) were assigned to protein translation (GO 0006412). Two candidate genes were assigned to electron transport (GO 0006118). The remainder (N=7) were each assigned to a unique GO term (see Fig. 3C for complete list). This analysis suggests that most of the ‘identifiable genes’ that were differentially expressed in fast-growing larvae are involved in a single biological process — protein translation.
Ontogenetic changes in gene expression
Comparison of ribosomal protein gene expression during larval development to 6 days post-fertilization revealed an increase in expression that differed between genotypes (Fig. 4). RNA blots conducted using a probe for ribosomal protein RPL35 (clone 68) revealed specific hybridization with a single transcript of 800 bp in all developmental stages analyzed (Fig. 4A). Increases in the abundance of RPL35 transcript during development were apparent in larvae from all genotypes, with a 5-fold increase, on average, in 6-day-old veliger larvae relative to unfertilized eggs (day 0) (Fig. 4B). Averaged across all genotypes, a significant 6-fold increase was observed between eggs and 6-day-old veliger larvae (Student's t-test P<0.05). These data reveal an over-dominant expression pattern by day 6 for RPL35 in both families of hybrid larvae.
Non-additive expression of ribosomal protein genes
Comparison of ribosomal protein transcript abundance between fast- and slow-growing larvae confirmed the non-additive expression of these genes in the majority of observations. RNA blots were conducted using probes for RPL18, RPL31, RPL35 and RPS3 (Table 2: clones 124, 76, 68 and 278, respectively). The radiolabeled probes used in the analyses of the expression patterns of these four genes hybridized specifically to single bands with estimated sizes of 400, 600, 800 and 700 bp, respectively. These transcript sizes for C. gigas are consistent with expectations based on the open reading frame sizes for these genes in humans, Saccharomyces cerevisiae, and the nematode Caenorhabditis elegans (HomoloGene Release 48.1, NCBI).
Comparing the quantitative estimates of expression levels obtained from these RNA (northern) blots revealed that expression in the two fast-growing families differed from the expected mid-parental value (represented as a dashed line in Fig. 2) for all four ribosomal protein genes (Table 3). For six of these eight comparisons (four genes in each of two fast-growing families), the direction of the difference measured with RNA blots matched that previously obtained from MPSS analysis. For example, RPL35 was up-regulated 47% in the 3×5 family in the previous study, and up-regulated 42% in the present study (Table 3). For the remaining two comparisons, the direction of the difference in gene expression did not agree with the previous study, although the expression data still showed a marked deviation from the expected mid-parental value (Table 3). Despite these differences in the direction of gene expression, the RNA blot results demonstrate that the non-additive expression of ribosomal protein genes is a reproducible characteristic of fast-growing larvae of C. gigas.
DISCUSSION
Growth rates of larval forms have been well characterized for a wide range of species of marine invertebrates (Thorson, 1950; Crisp, 1974; Manahan, 1990; His and Seaman, 1992; Fenaux et al., 1994; McEdward and Herrera, 1999). Most of these studies have focused on the effects of environmental conditions such as food and temperature. In contrast to these exogenous factors, the endogenous physiological processes that regulate growth rate during early animal development remain poorly understood. Additionally, there are genetic factors that regulate growth, one example being growth heterosis (hybrid vigor) associated with multi-locus heterozygosity. This phenomenon has been investigated for several decades (Shull, 1948; Singh and Zouros, 1978; Koehn and Shumway, 1982; Hawkins et al., 1986; Hedgecock et al., 1995; Bayne, 1999; Pace et al., 2006), but a comprehensive biological explanation still remains elusive. The contrasting growth phenotypes produced by genetic crosses of the Pacific oyster C. gigas allow for comparisons between half-sibling, inbred and hybrid larvae. It is important to note that the majority of physiological processes are indistinguishable between faster-growing hybrid and slower-growing inbred larvae, leading to the conclusion that these larval families are physiologically ‘normal’ (Pace et al., 2006). In the current study, a set of 181 candidate genes for growth heterosis were analyzed based on transcriptome-wide analysis (Hedgecock et al., 2007) of differential gene expression in fast- and slow-growing larvae. Our goal was to elucidate the biological processes underlying differential growth rates by identifying the genes involved.
Many of the candidate genes identified here have not been considered in previous investigations of growth regulation. For example, the primary gene expression data for clone 4 show that this candidate gene was expressed at 1.3-fold higher levels in faster-growing larvae. Sequence analysis identified this clone as caveolin, a membrane protein associated with endocytosis and exocytosis (Drab et al., 2001). In contrast, expression of clone 145 was 4.2-fold lower in fast-growing than in slow-growing larvae; sequence analysis identified this clone as fasciclin, a membrane protein involved in cell adhesion during Drosophila embryogenesis (Elkins et al., 1990). The list of candidates shown in Table 2 includes numerous other examples of genes that would not have been predicted from classical explanations for growth heterosis, and represent novel candidates for possible future study of the regulation of growth rate.
Other growth candidates were associated with processes that have previously been studied in the context of growth, including feeding, energy metabolism and protein metabolism. For example, clone 264 showed significant sequence similarity with the small cardioactive peptide precursor gene (SCPb), a neuropeptide expressed in the visceral ganglia of adult Pacific oysters (Hamano et al., 2005). This neuropeptide regulates contractile functions in molluscs that play obvious roles in feeding, including gut motility and radula activity (Lloyd et al., 1988; Miller et al., 1994). Studies of this signaling molecule in different mollusk species have suggested a stimulatory effect on feeding for some species and an inhibitory effect for others (Lloyd et al., 1988; Elliott et al., 1991). The SCPb peptide shows clear potential as a candidate for regulation of feeding activities. In our study, expression of the SCPb candidate gene was only detectable in slow-growing larvae, with no detectable transcripts in their fast-growing counterparts. This observation, in the context of previous reports showing increased feeding rates in fast-growing adult and larval bivalves (Bayne, 2004a; Pace et al., 2006), suggests a possible role for SCPb in genotype-dependent regulation of feeding activity and growth in bivalve larvae.
Energy metabolism has been extensively studied in the context of growth regulation (Koehn and Shumway, 1982; Hawkins et al., 1986). Several genes involved in energy metabolism were identified here, including two mitochondrial genes (clones 40 and 260: NADH dehydrogenase subunits ND4L and ND1, respectively) encoding components of the electron transport chain (Lenaz et al., 2006). In addition to those mitochondrial genes, nuclear genes with mitochondrial functions were also identified, including clone 78 (the ATP-synthase δ) (Walker and Dickson, 2006) and two coiled-coil-helix-coiled-coil-helix domains (clones 255 and 98: CHCHD2 and CHCHD3). Similar CHCHD domains are found in nuclear genes encoding mitochondrial products (Mootha et al., 2003), and in a set of genes expressed in proliferating human cell lines (Westerman et al., 2004). Previous studies on the growth advantage of adult bivalves with higher degrees of heterozygosity revealed differences in the metabolic efficiency of fast-growing animals (Hawkins et al., 1986). These findings were supported by studies of experimentally produced inbred and hybrid adults (Bayne et al., 1999) and larvae (Pace et al., 2006). In this context, the differential expression of these candidate genes with obvious roles in energy metabolism suggests that the previously reported metabolic differences might reflect a level of transcriptional control of metabolism. The occurrence of both nuclear and mitochondrial genes in this category is noteworthy, because the fast- and slow-growing larvae analyzed here included half-siblings that were derived from the same eggs (e.g. larval families 3×5 and 5×5). These larvae shared a common mitochondrial genotype, but differed in paternal genetic backgrounds. The growth advantage observed for these half-sibling larvae (different sire, same dam) suggests the possibility that interactions between nuclear and mitochondrial gene products play a role in the growth advantage of hybrids.
Among the candidate growth genes identified here, more genes were associated with protein metabolism than with any other biological process. In addition to the ribosomal proteins that have obvious roles in protein synthesis, other candidate genes were associated with protein folding and catabolism. For example, clone 59 was tentatively identified as peptidylprolyl isomerase, a gene that increases the efficiency of protein folding (Young et al., 2004). This gene was only detected in fast-growing larvae, and not in their slow-growing counterparts. Clone 74, identified as the proteasome subunit PSMD14 (Penney et al., 1998), was expressed in fast-growing larvae at 2.6-fold higher levels than in slow-growing larvae. Clone 73, identified as the protein glycoslyation gene DC2 (Shibatani et al., 2005), was expressed at 2.5-fold higher levels in fast-growing larvae than in their slow-growing counterparts. Because of the complexity of protein metabolism, these expression profiles do not lead immediately to clear predictions of physiological function (e.g. increased protein degradation). Nevertheless, the identification of a suite of genes associated with previously studied determinants of growth (feeding, energy metabolism and protein metabolism) provides a potential new set of molecular biological indices for studying the regulation of growth rates.
Ribosomal proteins were the single most abundant class among the 34 candidate genes identified here by searches of GenBank, comprising 50% of the total (Table 2). These 17 different ribosomal protein genes included nine components of the large ribosomal subunit (prefix L) and eight components of the small ribosomal subunit (prefix S). The direction of the difference in ribosomal protein expression varied among these genes; six were more highly expressed in fast-growing larvae, and 11 were more highly expressed in slow-growing larvae (Table 2). All four of the non-additive gene expression categories (overdominant, underdominant, dominant-high, and dominant-low) were observed among ribosomal protein genes. A significant finding from these analyses is that while the mean expression of ribosomal protein genes was the same in all families, the distribution of expression levels across genes differed between fast- and slow-growing larvae. Despite the differential expression between families for each gene considered separately, the average level of expression across all 17 ribosomal protein genes identified here did not differ between fast- and slow-growing larvae (ANOVA P=0.89), indicating a lack of overall up-regulation or down-regulation of ribosomal gene expression in fast-growing larvae. Interestingly, the different ribosomal protein genes were expressed at a more uniform level (i.e. closer to an equimolar ratio) in fast-growing larvae than in their slow-growing counterparts. This important point is illustrated by comparing the distribution of expression levels across the 17 different ribosomal protein genes (Table 2) in each of the families with the expected equimolar ratio using the chi-square test. All four families showed significant deviations from the expected ratio (P<0.001), but the magnitude of the difference was substantially smaller for the fast-growing hybrid larvae. The χ2 statistics for families 3×5 and 5×3 were 16,607 and 19,104, respectively (units of transcripts per million). The corresponding statistics for the slow-growing inbred families 3×3 and 5×5 were over 2-times higher (34,160 and 57,231, respectively), reflecting a greater deviation from the expected equimolar ratio. This comparison highlights the more uniform expression of ribosomal protein genes in fast-growing larvae than in their slow-growing counterparts.
This finding suggests a relationship between the stoichiometry of ribosomal protein gene expression and whole-organism growth and fitness. Such a relationship would support the recently proposed ‘balance hypothesis’, which predicts deleterious effects for an imbalance in the abundance of the constituent proteins for essential multi-protein complexes such as ribosomes (Papp et al., 2003; Marygold et al., 2007). In proliferating cells, ribosome biogenesis accounts for a significant proportion of the metabolic cost of cell proliferation (Schmidt, 1999), so any perturbations in this process would be expected to affect overall energy metabolism. The abundance of different ribosomal proteins is tightly regulated to ensure their availability in equimolar amounts required for efficient ribosome assembly (Warner, 1999). Ribosomal protein production is primarily controlled at the level of transcript abundance in yeast (Planta, 1997), and any free ribosomal proteins are rapidly degraded (Moritz, 1990). In general, protein synthesis and turnover consume a large proportion of the energy budget at the organismal level (Hawkins, 1991), accounting, for example, for up to 75% in growing sea urchin larvae (Pace and Manahan, 2006). It is likely that synthesis and degradation of proteins also represents a substantial metabolic cost in bivalve larvae (Pace et al., 2006). Each of the ~80 different ribosomal proteins comprises 0.1-0.5% of the total cellular protein (8-40% in total), so these are collectively some of the most abundant proteins in cells (Warner, 1989). Any changes in the degradation and synthesis of ribosomal proteins can therefore be expected to have a substantial metabolic impact. This highlights the potential for metabolic inefficiency in synthesizing and degrading excess copies of the more highly expressed ribosomal proteins. Ribosomal proteins have been extensively studied in the context of ribosome assembly and growth, in organisms ranging from bacteria to mammals (Tao et al., 1999; Jorgensen et al., 2004; Mayer and Grummt, 2006). Clearly, the synthesis and turnover of ribosomal proteins affect overall metabolism and growth. The ribosomal protein expression profiles observed in the current study for fast-growing bivalve larvae suggest a hypothesis for these genetically determined differences in growth rate. Non-uniform expression of ribosomal proteins (i.e. deviations from the equimolar ratio) in slow-growing larvae might lead to degradation of these proteins, resulting in metabolic inefficiency and slower growth. This expectation is consistent with the experimental evidence of more efficient protein metabolism in faster-growing adult stages of bivalve molluscs (Hawkins et al., 1986; Hawkins and Day, 1996). Because these ribosomal protein expression profiles suggest a plausible explanation for the growth differences that is consistent with previous studies, this class of genes was selected for further analysis.
The association between growth heterosis and non-additive expression of ribosomal proteins was tested with independent larval cultures obtained from the same genetic families used in our previous study (Hedgecock et al., 2007). Notably, the individuals used for genetic crosses in the current study (Fig. 1) were four generations removed from those used in our previous study. Measurements of transcript abundance by RNA blot analysis for four ribosomal protein genes (RPL18, RPL31, RPL35 and RPS3) revealed non-additive gene expression patterns similar to those apparent in the MPSS dataset (Table 3). The previously observed growth heterosis was also observed in larvae from these new crosses (Fig. 1), confirming the reproducibility of this association between rapid growth and non-additive ribosomal-protein gene expression. In the current study, all four candidate genes selected for validation showed non-additive expression in fast-growing hybrid larvae, consistent with findings for larvae from our previous study that were obtained from adults of an earlier generation (Hedgecock et al., 2007). There were differences in the magnitude and direction of these non-additive expression patterns for specific combinations of genes and families (Table 3), but the overall finding of non-additive ribosomal protein expression was confirmed. This demonstrates a cross-generational heritable association between growth phenotypes and gene expression profiles, suggesting that these genes are involved in the mechanistic basis of the rapid-growth phenotype.
This analysis of candidate genes was based on expression profiles in a single developmental age post-fertilization (6-day-old veliger larvae). The predictive value of gene expression profiles often depends on the differential timing of gene expression during development, a process that can be complex and differ markedly between genes. For example, during fruit fly development many genes are expressed at peak levels for a brief developmental period, while others increase gradually to a stable maximum (Arbeitman et al., 2002). Similarly complex patterns of differential gene expression have been observed during development of other species of marine invertebrates (Char et al., 1993; Marsh et al., 2000; Hinman et al., 2003). Differences in gene expression at a particular developmental stage might reflect differences in the timing of a developmental peak in expression, or differences in the overall level of expression throughout development. The data reported here for ontogenetic expression of RPL35 (Fig. 4) show that a dominant-high differential expression pattern (i.e. 3×3<3×5, 5×3, 5×5) was established by 2 days post-fertilization and persisted throughout the subsequent period of development studied. These findings confirm the expression patterns previously described at 6 days post-fertilization for this ribosomal protein gene.
The relationship between gene expression patterns and growth rates suggests the possibility of identifying molecular biological ‘markers’ for prediction of differential growth rates. For example, the expression profiles analyzed in this study include several genes that showed linear relationships to growth rates (Fig. 5). The examples shown (the three most linear examples of positive and negative relationships, based on R2 values) include two ribosomal protein genes and four genes of unknown function. Obviously, these relationships alone are not sufficient to show the general utility of these particular markers, because the growth data shown were obtained from the same crosses used to identify candidate genes. Nevertheless, the existence of these relationships suggests the novel possibility of predicting growth rates in marine larvae from gene expression profiles. Previous studies have used RNA/DNA ratios as an index of growth potential (Buckley, 1984). Our findings significantly advanced the use of this simple ratio by identifying specific candidate genes involved. The availability of such ‘biomarkers’ has obvious applications for defining physiological state, for understanding the adaptive significance of variability in growth rates, and for modeling genotype fitness under changing environmental conditions.
Recent advances in the genomic analysis of marine metazoans (Cameron and Rast, 2008) and the application of such approaches to longstanding questions in comparative and integrative physiology (Cossins and Somero, 2007) are offering new insights into systems biology. The application of these ‘discovery-based’ genomic approaches is leading to new testable hypotheses regarding the mechanisms of development, growth and many other biological processes. The biological phenomenon of hybrid vigor has been studied for decades (Shull, 1948). Although our study has not fully identified the mechanistic basis of hybrid vigor, our findings have provided new insights into this phenotype during larval growth. Half of the candidate genes identified by sequence comparisons are ribosomal proteins associated with protein synthesis, in addition to other well-characterized protein metabolism and energy metabolism genes. The remaining candidates include many genes not previously considered in classical explanations of differential growth rates. These findings provide a new set of testable hypotheses, and potential molecular biological indices, to enhance understanding of the physiological bases of variable and rapid growth rates in developing animals.
Acknowledgements
The authors thank Dr Patricia Von Dippe for her extensive technical assistance with the cloning and sequencing of genes. Our colleague Dr Dennis Hedgecock provided invaluable advice about data analysis and interpretation. This work was supported by the National Science Foundation under grant number OCE-0412696 to D.T.M. and by the W. M. Keck Foundation.