Identification of a tubulin-α gene specifically expressed in testis and adductor muscle during stable reference gene selection in the hermaphrodite gonad of the lion's paw scallop Nodipecten subnodosus

Marine Genomics 6 (2012) 33–44

Contents lists available at SciVerse ScienceDirect

Marine Genomics journal homepage: www.elsevier.com/locate/margen

Identification of a tubulin-α gene specifically expressed in testis and adductor muscle during stable reference gene selection in the hermaphrodite gonad of the lion's paw scallop Nodipecten subnodosus Raúl Llera-Herrera a, Alejandra García-Gasca b, Arnaud Huvet c, Ana M. Ibarra a,⁎ a b c

Centro de Investigaciones Biológicas del Noroeste, S.C. Mar Bermejo 195, Col. Playa Palo de Sta. Rita, La Paz B.C.S. 23090, Mexico Centro de Investigación en Alimentación y Desarrollo, A.C. Avenida Sábalo Cerritos s/n, Mazatlán, Sinaloa 82010, Mexico IFREMER, UMR 6539, LEMAR, BP 70, 29280 Plouzané, France

a r t i c l e

i n f o

Article history: Received 10 January 2012 Received in revised form 5 March 2012 Accepted 9 March 2012 Keywords: Nodipecten Real time qPCR Reference gene Gametogenesis Testis-specific tub-α

a b s t r a c t For non-model species, as many used for aquaculture, with minimal or no genomic information, relative quantification of gene expression studies requires preliminary research including the isolation of potential reference genes and the identification of those stably expressed under the biological conditions of interest. Here we report on the isolation of five partial gene sequences from gonad tissue cDNA in the functional hermaphrodite scallop Nodipecten subnodosus to be evaluated as reference genes: 18S-rRNA, riboprotein l8 (rp-l8), actin-β (act-β), elongation factor 1α (ef-1α) and alpha-tubulin-α (tub-α). We found that 18S-rRNA was stably expressed independently of the priming method used to reverse transcribe RNA to cDNA, oligodT or random hexamer. Stability analysis for the five putative reference genes with geNorm and NormFinder indicated that 18S together with rp-l8 were the most stable genes for normalization of gene expression during gonad development in both, male and female sexual regions of the hermaphrodite N. subnodosus. The least stable gene was tub-α, showing a biased expression profile between sexual regions of the gonad, therefore this gene was analyzed thereafter as a target gene together with vitellogenin (vit) and a DEAD-box RNA helicase (dbx) gene. Relative expression, estimated by normalization with the combination of 18S and rp-l8 as reference genes, indicated that as gonad development advanced two of the target genes were up-regulated, tub-α in the male region and vit in the female region. Whereas an increased expression was expected during development for vit for its known role in vitellogenesis, the increased expression of tub-α in the male sexual region was unexpected, and pointed toward this gene being a testis-specific α-tubulin isotype. Further analyses of gene expression among tissues indicated that tub-α is specifically and highly expressed in the male gonad, although expression in adductor muscle was also observed at significantly lower levels. The existence of testis specific α- and β-tubulins has been previously reported in other taxa, relating their function to sperm axoneme formation. Tissue-specific tubulin genes, particularly their promoters, have recently found an application as native promoters for transgene tissue-specific expression in research and reproductive control of insect plagues. The third target gene, a putative member of the DEAD-box RNA helicase family (dbx), showed no changes in expression during gonad development or between sexual regions, therefore it was chosen to discuss the different statistical inferences resulting from the arbitrary use of ‘randomly chosen’ reference genes when normalizing gene expression. © 2012 Elsevier B.V. All rights reserved.

1. Introduction The lion-paw scallop, Nodipecten subnodosus, the largest and most important of several scallop species in the Pacific Northwest of Mexico, has been the subject of several genetic studies (Petersen et al., 2008, 2010b). Until recently no genomic resources had been developed with the exception of microsatellite markers (Ibarra et al.,

⁎ Corresponding author. Tel.: + 52 612 1238484. E-mail address: [email protected] (A.M. Ibarra). 1874-7787/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.margen.2012.03.003

2006; Petersen et al., 2009b) which allowed for the generation of the first genetic linkage map in this species (Petersen et al., 2009a, 2010a). In an additional effort to develop genomic resources for this functional hermaphrodite species, we have generated subtractive suppressive hybridization (SSH) libraries of mature male and female gonad parts subtracted with undifferentiated gonad (Llera-Herrera et al., 2010). Characterization of genes isolated in those libraries, and the understanding of their expression patterns during gametogenesis are important to provide a better understanding of the key factors involved in sex differentiation between gonad sexual regions, meiosis progress, and gamete maturation. As a first step in those

34

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

studies, the identification of stably expressed genes during gonad differentiation and development is necessary. Inasmuch as establishing expression stability of reference genes to use in the normalization of target genes when there is an expected differential expression, for example during gametogenesis, it is particularly important when the morphological and biochemical changes occurring during reproductive development diverge between sexual regions of a unique gonad organ (ovotestis) simultaneously producing spermatozoa or oocytes. This is the case of the functional hermaphrodite scallop N. subnodosus. In this organism it is expected that gonad development will modulate the whole transcriptional activity (including many housekeeping genes) in a sex-biased manner (Lu and Wu, 2011). Whereas no studies have been performed on the functional hermaphrodite scallops regarding the definition of stably expressed genes or gonad sex differentially expressed genes, studies do exist in the model nematode Caenorhabditis elegans and in some hermaphrodite plants indicating differences in expression between male, hermaphrodite and female sexual stages and organs (Jiang et al., 2001; Silveira et al., 2009). Fluorescence-based real-time PCR (or qPCR) has been the most robust technique widely used for the analysis of gene expression. Relative quantification of gene expression relies on the use of internal reference genes (usually housekeeping genes) for the normalization of the mRNA level of the gene(s) of interest, assuming that the expression of the endogenous reference gene(s) is constant between different experimental treatments or in different tissues. A housekeeping gene is a constitutive gene required for the maintenance of basic cellular function, and although some housekeeping genes are expressed at relatively constant levels, others may vary depending on experimental or biological conditions. Commonly used genes for normalization of gene expression include, among others, elongation factor 1-α, β-actin, tubulins, different ribosomal protein genes, glucose-6phosphate dehydrogenase, and ubiquitins. Regardless of which genes are used, the condition of physiological stability in the set of reference genes must be tested as part of the standardization process in most qPCR experiments; different validation models have been developed for that purpose, as geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), and BestKeeper (Pfaffl et al., 2004). Additionally, besides the identification of the most stable reference genes, gene expression analysis by qPCR requires the enzymatic conversion of the entire population of extracted mRNA or total RNA into cDNA. This can be done by using either random hexamer or oligo-dT primers. Random hexamer targets the entire population of RNA, whereas the oligo-dT primer is specific for the poly(A) signal present in the 3′-end of most eukaryotic mature cytosolic mRNAs. Despite the fact that both priming methods for cDNA synthesis produce bias in the selective enrichment types or in regions of the targeted RNA (specifically the 3′-end in the case of the oligo dT cDNA and in the middle region of the mRNA for random hexamer cDNA) (Hansen et al., 2010), it is generally claimed that when eukaryotic mRNA is targeted, oligo-dT performs relatively better, preserving in the cDNA the original representation of the mature mRNA population in the cell/tissue/organ, and is mostly preferred over random hexamer for qPCR assays for relative gene quantification (Bustin et al., 2005). Because of this, two ribosomal genes well conserved in eukaryotes (18S, 28S) are sometimes excluded as reference genes for relative gene expression studies as it is thought that rRNA contains no poly(A) regions and therefore it cannot be reverse transcribed by oligo-dT priming (Radonic et al., 2004). However, recent reports have demonstrated that some degree of PCR amplification of rRNA genes occurs in human (Slomovic et al., 2006), rat (Zhu and Altmann, 2005), fish (Kortner et al., 2011) and insect samples (Xue et al., 2010), providing indirect evidence that some forms of rRNAs are rich in A-repeats or are in fact polyadenylated as any typical eukaryotic mRNA. In our own experience, a strong amplification signal of an 18S fragment in oligo-dT cDNA synthesized from total RNA

has been detected in the lion-paw scallop N. subnodosus (Bivalvia: Pectinidae) as well as in the Pacific oyster Crassostrea gigas (Bivalvia: Ostreidae). In the present study, we evaluated the feasibility of using the 18S rRNA gene with oligo-dT-primed cDNA, in parallel with four mRNA gene transcripts: ribosomal protein l8 (rp-l8), elongation factor 1-α (ef1-α), actin-β (act-β), and tubulin-α (tub-α), as part of a study to define the best set of reference genes for gametogenesis-related gene expression analyses in the hermaphrodite scallop N. subnodosus. In order to identify the best reference genes under these biological conditions, we analyzed three gametogenic stages: undifferentiated, early maturation and mature in each of the two sexual regions of the differentiated gonad, and estimated relative expression of three target genes using the best reference genes identified for normalization. After finding that tub-α was not stable across the conditions tested and was only expressed in the male part of the gonad, we further evaluated this gene for its ubiquity across other tissues other than the gonad. 2. Materials and methods 2.1. Selection of reference genes and primer design For reference gene sequences in N. subnodosus we isolated partial sequences following three strategies. The first one was by nucleotide alignments and degenerated primer design to obtain two of the reference genes evaluated in this study, ribosomal protein gene (rp-l8) and elongation factor alpha-1 (ef1-α). For rp-l8 a fragment of the 60S subunit l8 (rp-l8) gene was first amplified using degenerated primers designed on the alignment of 60S sequences from the bivalves Argopecten irradians (CB416772.1), Chlamys farreri (AF526246.1), C. gigas (AJ563478.1), and the penaeid shrimp Litopenaeus vannamei (DQ316258.1) (sense 5-CGT CAT GGA TAC ATC AAG GGT-3; antisense 5-GTT ACC WCC WCC ATG RGG ATG-3) obtaining a PCR fragment of 495 bp. For ef1-α the gene coding sequences of Dreissena polymorpha (AJ250733.1), C. gigas (AB122066.1), Mytilus galloprovincialis (AB162021), Modiolus americanus (AY580265.1) and Crepidula fornicata (DQ087484.1) were aligned and used for degenerated primer design (sense 5′-CAY MGW GAT TTC ATC AAG AAC ATG AT-3′; antisense: 5′-CTT RAT YTC WCC DGG GTG GTT CA-3′), obtaining a PCR fragment of 781 bp, which was located at positions 456–1235 of the complete mRNA sequence of C. gigas. The second strategy was to obtain specific primers from conserved regions of the two genes characterized in related species: beta actin (act-β) and ribosomal 18S-rRNA. For act-β the primers were designed from well conserved regions of complete coding sequences from the pectinid species Myzuhopecten yessoensis (GenBank accession no. DQ787858; positions 272–481) and Patinopecten magellanicus (GenBank accession no. U55046.1; positions 242–451). For the 18S-rRNA gene a multiple alignment with complete gene sequences of representative bivalves deposited in the GenBank database was performed: Pecten maximus (accession no. L49053.1), Chlamys islandica (accession no. L11232), A. irradians (accession no. L11265.1) and M. galloprovincialis (accession no. L33451.1). The third strategy was to select a fifth gene, alpha tubulin (tub-α), as a potential reference gene, from a partial sequence of N. subnodosus (GenBank accession no. JN034910) obtained from the suppressive subtractive cDNA library prepared in our laboratory (Llera-Herrera et al., 2010). The designed primers were used for PCR amplification (10 mM dNTP, 1 U Taq, 1 μl of cDNA synthesized as described in Section 2.3, MgCl2 in a range from 1.5 mM to 2.5 mM, and 10 to 20 pmol of degenerate or specific primers; all in 50 μl reactions). After 30 cycles of PCR amplification, products were separated in 1% agarose electrophoresis, and products of the expected size were gel-excised and purified (PCR gel purification kit, Qiagen), inserted into pCR2.1 vector (Invitrogen),

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

and cloned in E. coli DH5-α. For plasmid isolation 1 ml of liquid media culture was harvested by alkaline lysis and phenol/chloroform extraction (Sambrook and Russel, 2001), and sequenced (Macrogen, Korea). The obtained sequences from rp-l8, ef1-α and tub-α were used to design specific primers with optimal conditions for qPCR analysis (18–21 nt, melting temp ~64 °C, 40–60% GC, product size between 120 and 200 bp) using the Primer3Plus web resource (Untergasser et al., 2007). For 18S-rRNA and act-β the primers used for qPCR were the same as those we used to obtain the sequences. 2.2. Biological material Hatchery-produced lion-paw scallops were reared in Rancho Bueno Lagoon, Bahía Magdalena Bay, located in the Pacific side of the Baja California Peninsula in Mexico. Adult scallops (mean length: 58.5 mm; mean body weight of 65 g) were collected in May 2010, and after histological examination of their gonads, they were separated into early maturation and mature stages based on the female gonad part, according to the oocytes developmental stages previously described in Racotta et al. (2003) for this species. Juvenile scallops (mean length: 23 mm; mean body weight of 2.6 g) were also assessed by histological observation of their gonads to confirm that they were undifferentiated individuals. Gonads of juvenile and adult scallops were sampled and stored for further RNA isolation in aqueous sulfate salt solution (25 mM sodium citrate, 10 mM EDTA, 70 g ammonium sulfate/100 ml solution, pH 5.2) (Nolan and Bustin, 2008). For adult scallops, male (proximal) and female (distal) portions were collected separately (Fig. 1). After histological confirmation of gonad stage, six scallops were selected for tissue analyses per condition: undifferentiated, early and late male gonads, early and late female gonads. Early male/female and late male/female gonads were sampled from the same six hermaphrodite individuals in each case. 2.3. RNA isolation and cDNA synthesis Individual samples (~50 mg) from each sexual region (female and male) of the hermaphrodite adult gonads and complete juvenile (undifferentiated) gonads were homogenized with 800 μl of TriPure reagent (Roche) using 1.6 ml microcentrifuge tubes and 100 μg of

35

425–600 μm glass beads as disruptive matrix. Tissues were disrupted immediately after immersion in pre-filled tubes with a FastPrep-24 Instrument (MP Biomedicals) by two oscillatory pulses of 30 s at 6.0 m/s, with an intermediate resting step of 5 min. Homogenates were kept at room temperature for 5 min, then 160 μl of chloroform was added, samples were gently mixed by inversion for 15 s, incubated at room temperature for 15 min, and centrifuged for 15 min at 12,000 g for phase separation. The upper layer was recovered, and RNA was precipitated for 15 min in 500 μl of isopropanol. The RNA was centrifuged for 15 min at 12,000 g, then washed with chilled 75% ethanol, pelleted by centrifugation at 7600 g for 5 min, decanted and air-dried. Pellets were resuspended in nuclease-free water and UV-quantified by spectrophotometry (NanoDrop-1000, Thermo). Absorbance ratios (260/280 and 260/230) over 1.9 were considered acceptable for RNA purity, and 4 μg of freshly prepared RNA was incubated for 30 min at 37 °C with 2 U of RQ1 RNAse-free DNAse (Promega) in a reaction volume of 20 μl to avoid any trace of contaminant genomic DNA. DNAse was inactivated by incubation at 65 °C for 10 min, after adding 2 μl of stop solution (20 mM EGTA, pH 8.0). Effective DNAse treatment was confirmed in each DNAse-treated RNA sample by verifying no-amplification after PCR when using rp-l8 primers (0.4 μM each) and 200 ng of RNA in a 15 μl reaction. The amplification program and PCR cocktail composition are described below (Section 2.4); verification was done in (2%) agarose gel stained with 1× GelRed (Biotin CA). The RNA was directly used for reverse transcription (RT), synthesizing two different sets of cDNA in parallel: oligo-dT(15)-primed and random hexanucleotide-primed (from here referred as oligo-dT/cDNA and random hexamer/cDNA, respectively). The RT reactions were prepared with the ImProm-II reverse transcription system (Promega, cat. A3800) as follows. DNA-free RNA and 0.5 μg of the reverse transcription primer (oligo-dT or random hexamer) were mixed, adjusted to 5 μl with DEPC-treated water, and incubated at 70 °C for 5 min, chilled with ice for 2 min, then incubated briefly at 25 °C before the addition of 15 μl of reverse transcription mix (for a final concentration of 5.6 mM of MgCl2, 0.5 mM of dNTP, 40 U of recombinant ribonuclease inhibitor, 1 μl of ImProm-II reverse transcriptase, in a final volume of 20 μl per RT reaction). Primer alignment was performed at 25 °C for 5 min, and cDNA extension was completed after 1 h incubation at 42 °C. An aliquot of the cDNA was diluted 75-fold in Tris buffer (10 mM; pH 8.0) and kept at −20 °C until qPCR analyses. 2.4. Real-time qPCR

Fig. 1. N. subnodosus scallop gonad macrostructure. Discontinuous lines indicate the sampled tissues in each sexual region (female and male). AM: adductor muscle. DG: digestive gland.

All qPCR reactions were conducted in triplicate in 96-well plates and optical adhesive film (Bio-Rad), using the CFX96 Real-Time PCR detection system (Bio-Rad). A qPCR cocktail-mix was carefully prepared in our laboratory, with minimal modifications from that described in Hernández-Arteaga and López-Revilla (2008). The mix contained 2.5 mM MgCl2, 200 μM dNTP (each), 0.45 U of recombinant Taq DNA polymerase (Invitrogen, cat. 10342), 0.46 μM of each primer, and 1 × EvaGreen fluorescent dye (Biotium, CA) in 15 μl of final volume per reaction. A 2 × PCR mix was prepared in one batch to avoid technical variation in the results. This 2× mix assured nearly identical performance to a commercial kit (SsoFast EvaGreen Supermix, Bio-Rad) when serial dilutions were compared in terms of the slope and y-intercept of a log-concentration vs. cycle curve fitting line, using primers designed for the 18S rRNA fragment and random hexamer/oligo-dT-primed gonad cDNA of a pool of samples of different stages/gonad regions (our cocktail: slope = −3.515, E = 92.5%, y-intercept = 9.288; commercial kit: slope = −3.482, E = 93.7%, y-intercept = 9.208). PCR conditions were: 94 °C 3 min; 40 cycles of 94 °C (10 s), 60 °C (15 s) and 72 °C (30 s), acquiring the fluorescence at 79 °C (2 s). At the end, a dissociation step from 60 °C to 94 °C (1 °C/s) was performed, and all runs were inspected for unique and specific products using the melting curves from this last

36

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

step of the qPCR thermal profile. qPCR products were confirmed by agarose gel electrophoresis and sequencing. 2.4.1. Contrasting cDNA priming methods for qPCR efficiency To assess possible bias in qPCR efficiency resulting from the priming strategy for reverse transcription, we first determined for each gene (sample size = 84 reactions per gene per primed-cDNA set) the individual performance of each reverse transcription reaction. This was done by calculating the fluorescence slope, limiting it to the window-of-linearity, using the software LinReg ver. 12.10 (Ruijter et al., 2009). The imputed data for the analysis were the non-baseline corrected fluorescence values obtained in each cycle. The mean efficiencies obtained for each gene/cDNA were compared with a t-student test (α = 0.05). 2.4.2. Estimation of amplification efficiencies For each gene (reference or target), amplification efficiencies (E) were determined by the slope calculation of 5-fold serial dilutions, starting with 0.5 μl of a pool of oligo-dT/cDNA or random hexamer/ cDNA, in 15 μl reactions, and using a fixed fluorescence threshold value of 150 RFU. Efficiency value was then obtained from the slope of the log-linear function of the concentration vs. fluorescence, using the equation: E = (10 (− 1/slope) − 1) (Bustin et al., 2009), and used for gene-specific efficiency correction in the geNorm stability analysis and in the relative quantification model as recommended by Pfaffl (2004) and Hellemans et al. (2007). Amplification efficiencies were used for further analyses. 2.5. geNorm and NormFinder stability analyses of reference genes reverse transcribed with both methods Relative stability of a potential set of reference genes, either combining or separating each set of oligo-dT/cDNA or random hexamer/ cDNA was estimated using two methods, geNorm and NormFinder, running the Excel applets provided by the authors. Additionally, NormFinder was used to compare the relative stability of reference genes within each priming method, oligo-dT and random hexamer, considering the reproductive stage of each gonad region as another source of variation. For geNorm, the input data were relative quantities (RQ), which were calculated with the equation: RQij = E (Cq(min) − Cqij), where E is the gene-specific efficiency (calculated from the slope of the dilution curve as above), and (Cq(min) − Cqij) is the absolute difference for each Cq sample against the lowest Cq (with the highest abundance) in the dataset. For NormFinder, log(natural) transformed Cq values were used as input data. Differences in relative abundance of each reference gene when the cDNA was transcribed with oligo-dT or random hexamer were estimated using the following equation for relative abundance, E (ΔCq);

where E is the gene-specific efficiency, and ΔCq is the mean difference between the Cq values obtained by oligo-dT or random hexamer cDNAs. Last, to understand the biological source of variation resulting from differences in expression of the selected reference genes, we estimated the relative variation of each gene between different gonad stages and sex conditions: undifferentiated, early oogenesis, late oogenesis, early spermatogenesis and late spermatogenesis. The intergroup variation, obtained as part of the output from the NormFinder applet, was plotted for each gene to identify potential bias in gene expression stability driven by the gonad condition in which they were examined. 2.6. Target genes for evaluating relative gene expression Relative target gene quantification by qPCR was also conducted in triplicate reactions for all genes using oligo-dT/cDNA. As target genes we first selected two genes, vitellogenin (vit) and a member of the DEAD-box RNA helicase (ddx). These genes were partially identified by means of degenerate primer PCR in a similar way as the reference genes: vitellogenin (vit: sense 5′-TCA GTW TRC TGG ARR TGA ACC T-3′; antisense 5′-GCT BTC TGC WAT RGC RCT TTS GAT-3′) obtaining a PCR fragment of 337 bp, and a member of the DEAD-box RNA helicase (ddx: sense 5′-GAT GGC CTG TGC YCA RAC AG-3′, antisense 5′-CCR AAW CCC ATR TCC AAC AT-3′) obtaining a PCR fragment of 371 bp. The obtained sequences were used to design specific primers with optimal conditions for qPCR analysis. Later, a third gene, originally evaluated as a potential reference gene, tub-α, was evaluated for relative gene expression after it was found that it was not stably expressed in different developmental stages or sexes. Differences in relative expression of the three target genes were estimated using for normalized quantification the combination of the two most stable genes. Relative quantification data were grouped by gonad condition (undifferentiated, early and late oogeneses, and early and late spermatogeneses) for statistical analyses and contrasted by one-way ANOVA; individual significances were calculated by a t-student test (α = 0.05) using the software qBasePlus (Biogazelle Software, Belgium) (Hellemans et al., 2007). Furthermore, in order to contrast the interpretation of relative gene expression when randomly selected reference genes are used, we estimated relative expression for the target gene ddx using each one of the reference genes (18s, rp-l8, act-β, ef1-α) individually, the most stable pair found, and all the reference genes together. 2.7. Analyses of tub-α tissue expression ubiquity We re-analyzed the expression of tubulin-α in both gonad regions (female and male gonad parts) but including now other organs/tissues

Table 1 Primer sequences of reference genes and amplification efficiencies. Gene

GenBank accession number

qPCR primers

18S

n.a.

ef-1α

JN034909

act-β

n.a.

rp-l8

JN034908

tub-α

JN034910

ddx

JN034912

vit

JN034911

Sense Antisense Sense Antisense Sense Antisense Sense Antisense Sense Antisense Sense Antisense Sense Antisense

5′-GAAATTCTTGGATCGCCGTA-3′ 5′-GCCGAGTCATTGAAGCAACT-3′ 5′-GTGCCAGTGGGTAGGGTAGA-3′ 5′-CTCCAGCAACGTTTCCTCTC-3′ 5′-CAGAGCAAGAGAGGTATCCTCACC-3′ 5′-GTTGAAGGTCTCGAACATGATCTG-3′ 5′-CGTCATGGATACATCAAGGGT-3′ 5′-CAAACAGTCCAGTGTACATGCC-3′ 5′-TTCAGCCTGATGGACAAATG-3′ 5′-TCTGGATGAAACAGCTGACG-3′ 5′-CCAATCTTGCCTCGTTCAAT-3′ 5′-GCTTTAATCCTGGCCCCTAC-3′ 5′-CTGGCCAAAGCGAGTTCTAC-3′ 5′-TCGTTGAACATGGTGTTCGT-3′

Tm

Product size

Efficiency %

63.7 64.3 63.9 63.9 66.2 65.8 63.5 64.3 63.6 64.1 63.8 63.4 63.7 63.6

168

0.93

199

0.82

209

0.93

153

0.96

181

0.94

206

0.93

178

0.88

n.a. indicates that the sequences were too short for submission to GenBank, but are provided as supplemental material in the online version of this article.

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

37

Table 2 BlastX hits and identity values of the partial sequences obtained in this study against reference protein sequences (RefSeq database) of NCBI-GenBank. Partial gene sequence in N. subnodosus

Description

Species

Score

E-value

Identity

rp-l8 ef-1α tub-α ddx vit

Putative ribosomal protein L8 (XP_002192200.1) Elongation factor 1 alpha (XP_002604722.1) Alpha tubulin (XP_002580034.1) ATP-dependent RNA helicase DDX3X-isoform 3 (XP_002762847.1) Vitellogenin (CAQ06469.1)

Taeniopygia guttata Trichoplax adherens Schistosoma mansoni Callithryx jacchus Pecten maximus

302 426 222 207 193

9- 81 2- 117 1- 56 2- 63 5- 59

85% 83% 100% 79% 99%

(adductor muscle, gills, mantle edge, digestive gland, and nephridia). Samples from these organ/tissues were obtained from six adult scallops in mature reproductive condition, preserved as mentioned above (Section 2.2), and RNA extraction and oligo-dT/cDNA synthesis were performed under the same protocols described before but with an additional precipitation step in 4 M LiCl2 for RNA samples from digestive gland and nephridia. The obtained cDNAs were diluted 20 times in water for qPCR analyses. qPCR reactions were performed in duplicates, using one plate with all the samples and target/reference genes to avoid inter-run variation. We first analyzed the whole set of reference genes to assess the best genes for normalization across tissues using geNorm and NormFinder.

identity; E-value of 1e - 67), including a conserved vitellogenin_N family domain (PF01347.16) in Pfam Database (Finn et al., 2010). 3.2. Gene amplification efficiencies Amplification efficiencies (E), as calculated from the slope curve fitting all sampled gonad tissues (undifferentiated, early and late oogenesis, and early and late spermatogenesis) for oligo-dT/cDNA and random hexamer/cDNA are included in Table 1. These ranged between 0.82 and 0.96. These values were introduced as a correction parameter in relative quantification analyses and geNorm stability calculations. 3.3. Analyses between oligo-dT and random hexamer priming methods

3.1. Isolation of reference and target genes

Each set of cDNA produced by different priming methods was analyzed for differences in qPCR performance and in reference gene expression stabilities.

Isolated genes were deposited in GenBank (Table 1). A partial fragment of 168 bp coding for 18S rRNA was isolated from N. subnodosus by sequence alignment (supplemental figure) and primer design (Table 1) from conserved regions. This sequence was aligned without mismatch (100% identity) between positions 913 and 1020 of complete 18S rRNA gene sequences of other pectinid species: P. maximus, Mimachlamys varia, and C. islandica (GenBank accession nos. L49050.1, L49051 and L11232.1, respectively). A 209 bp fragment was obtained directly for act-β by specific primers (supplemental figure) which were also used for qPCR (Table 1). For all other isolated genes in N. subnodosus (reference or target) the best hit with reference proteins of public databases is presented in Table 2. A 781 bp cDNA fragment of the ef1-α gene obtained by degenerated primer design was confirmed as expected with a high level of sequence identity (83% BlastX) with the EF1-α protein sequence of Trichoplax adherens. The sequence obtained for the rp-l8 was confirmed by the same search tool, presenting an 85% identity with the corresponding protein sequence of Taeniopygia guttata. The conceptual translation of the vit sequence presented significant identity with the vitellogenin protein sequence of P. maximus (92%

3.3.1. PCR performance between oligo-dT and random hexamer The mean efficiency value, calculated by the regression slope in the amplification window-of-linearity, was used to compare PCR performance between the two priming methods. Calculated efficiency values for cDNAs synthesized by random hexamer and oligo-dT using the LinReg software showed no significant differences (t-student tests; p > 0.05) between the priming methods for any of the tested genes: 0.84 for 18S rRNA and 0.86 for rp-l8, in both cases regardless of the priming method, 0.84 (oligo-dT) and 0.83 (random hexamer) for ef1-α, and 0.74 (oligo-dT) and 0.82 (random hexamer) for act-β. This indicates that both types of cDNAs are comparable for further reference gene stability analyses. In the case of tub-α, efficiency values could not be estimated because the window of linearity cannot be calculated in samples with very low expression levels (Cq > 35) as was seen for undifferentiated and female gonads. Descriptive statistics contrasting relative expression for oligo-dT and random hexamer priming methods estimated for all reference genes using the amplification efficiencies (E) in Table 1, are shown in Table 3. Mean Cq obtained for all translated genes (rp-l8, act-β, ef1-α, tub-α) was always lower when random hexamer/cDNA was

3. Results

Table 3 Descriptive statistics for the amplification values of the reference gene set in oligo-dT and random hexanucleotide cDNAs. Gene 18S cDNA type n geomean [Cq] min [Cq] max [Cq] std dev [± Cq] CV [% Cq] Δ Cq Oligo-dT-Random Relative abundance ratio (oligo-dT/random)

Oligo-dT

rp-l8 Random

Oligo-dT

28 19.9 18.4 23.8 0.9 4.5

act-β Random

Oligo-dT

28 16.0 14.8 18.0 0.7 4.2

27.0 25.0 30.8 1.2 4.4

ef-1α Random

Oligo-dT

28 26.3 24.3 30.2 1.4 5.3

27.9 23.5 34.6 2.8 9.8

tub-α Random

Oligo-dT

28 27.3 22.9 38.7 2.7 9.8

26.3 23.9 34.0 2.2 8.3

Random 28

25.4 23.2 31.4 1.7 6.8

29.0 21.6 36.1 3.5 11.9

27.3 20.6 32.5 2.9 10.6

3.9

0.8

0.6

0.9

1.7

12.9

1.7

1.5

1.9

3.1

38

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

a

b

Fig. 2. Stability values of reference genes between gonad regions and stages (n = 30). geNorm (a) and NormFinder (b) results are presented for each gene and cDNA priming method: O-dT = oligo-dT, random = random hexamer.

used, and the same was seen for 18S rRNA. However, when the Cq difference was transformed to relative abundance (see methods, Section 2.5), there was a 12.9 × less target 18S template for qPCR with oligo-dT/cDNA when compared to random hexamer/cDNA, whereas from mRNA-derived genes, the difference ranged between 1.5× and 3.1× less target template when oligo-dT was used (Table 3). 3.3.2. geNorm and NormFinder stability analyses for reference genes When the genes reverse transcribed to cDNA by both methods (oligo-dT and random hexamer) were pooled and analyzed by geNorm software, rp-l8 was denoted as the most stable regardless of the type of cDNA-priming method in which expression stability was estimated (Fig. 2a); stability values estimated for rp-l8 (as a pairwise stability obtained in oligo-dT and random hexamer/cDNAs) were the lowest (M= 0.68). 18S appeared to be the second most stable gene, giving M-values of 0.86 and 0.96 when cDNA was synthesized with random hexamer and oligo-dT, respectively. In a similar way, NormFinder indicated that the most stable genes were 18S and rp-l8 (Fig. 2b). When using the geNorm software, 18S rRNA and rp-l8 were denoted as the most stable genes in oligo-dT/cDNA; in contrast, for random hexamer/cDNA, rp-l8 and ef1-α were the best pair, followed by 18S (Fig. 3). The most unstable gene detected by geNorm was tub-α. Indeed, the standard deviation and coefficient of variation (CV) were lower for 18S (CV 4.5 and 4.2) and rp-l8 genes (CV 4.4 and 5.3) than for all other reference genes. The gene with the highest CV was tub-α (CV 11.9 and 10.6) (Table 3). NormFinder stability values were recorded for the same set of putative reference genes, separately for oligo-dT and random hexamer/ cDNAs (Fig. 4). Both 18S and rp-l8 were found again as the most stable pair of reference genes, scored with the same relative stability value in

a

b

Fig. 3. Stability values of reference genes found by geNorm for each cDNA priming method: oligo-dT primed cDNA (a) and random-primed cDNA (b).

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

39

(18S + rp-l8), nor did we find differences when using 18S or rp-l8 alone as unique normalization genes (Fig. 7). However, when normalizing with act-β, both male and female early and late gonad regions appeared to overexpress the ddx gene when compared to the undifferentiated gonad. When using ef1-α as a unique normalization gene, we observed that, in relation to the undifferentiated gonad, an apparent overexpression of ddx occurred only in mature testis in late spermatogenesis (Fig. 7). When using the geometric mean of all four reference genes as normalization factors (rp-l8 + 18S+ ef1-α + act-β), the statistical differences in the expression levels were as follows: undifferentiated≤ female region (regardless of the developmental stage)= early spermatogenesis (immature male region)≤ late spermatogenesis (mature testis). 3.5. Relative expression of tub-α in different tissues

Fig. 4. NormFinder stability values for reference genes when primed by either, oligo-dT (black bars) or random hexamer (white bars).

both cDNA priming methods (0.032 for the two genes). However, the most stable gene in oligo-dT/cDNA was 18S (stability value of 0.031), and in random hexamer/cDNA was rp-l8 (stability value= 0.037). The least stable gene regardless of the cDNA priming method was again tub-α, with stability values between 0.133 and 0.138 for oligodT and random hexamer/cDNAs respectively. The estimated biological source of variation for relative expression of all putative reference genes based on the magnitude of variation within each gonad tissue sampled, and by each cDNA priming method, is presented in Fig. 5. It can be seen that tub-α shows the highest level of natural log (Cq) inter-group variation, with a trend to be negative in the male and positive in the female gonad sexual region (lower Cq= higher expression), indicating that expression of this gene is gonad-sex-biased toward the male region (Fig. 5). For the other genes, including 18S, the trend of the variation is toward positive in male gonad, and negative in female and undifferentiated gonads, and with less inter-group variation than tub-α. When comparing variation patterns between cDNA priming methods, the variation within groups is relatively the same. 3.4. Relative expression of target genes (vit, ddx, and tub-α) using the best set of reference genes Choosing the two most stable reference genes (rp-l8 and 18S) for normalizing, we analyzed the relative expression of the three target genes. The relative expression of tub-α and vit genes was clearly associated to the sexual region (Fig. 6). tub-α was overexpressed 50× and 330× in early and late developmental stages of the testis (male portion) compared to the corresponding stages in the female region, and from 130× to 600× compared to the undifferentiated gonad. Although tub-α expression was apparently higher in mature testis than in early developing testis (4.6×), this difference was not statistically significant. A similar expression pattern was found in the female gonad for vit, which showed an increased expression from early oogenesis to late oogenesis, with expression levels of ~1000× in the late oogenesis female region when compared to the undifferentiated gonad. The effect of using randomly chosen reference genes for quantification of relative expression was evaluated using ddx as a target gene. We used one, two or all reference genes but not tub-α for this comparison. As was shown before (Fig. 6), we did not observe significant differences in relative expression between regions or developmental stages of the gonads when using the two most stable genes found in this work

Both, geNorm and NormFinder analyses indicated that the most stable genes and therefore best reference genes for normalization across the different organs/tissues were ef1-α and rpl8 (Fig. 8), both of which were then used to estimate relative expression of tub-α across organs/tissues. Tubulin-alpha was highly expressed in testis and, in significant lower degree in adductor muscle; the rest of the tissues analyzed showed negligible expression (Fig. 9). qPCR products of tub-α obtained from male gonad cDNAs were directly sequenced with the same forward primer. We would observe ambiguous base calling if different amplicons were produced as a result of multiple target amplification; however, this was not observed, and specific amplicon was confirmed by clean chromatograms. The fact that the same transcript is amplified in the gonad and in the

a

b

Fig. 5. Relative variation for reference gene expression between the gonad sexual part and gametogenic stages of the scallop N. subnodosus in random-primed cDNA (a) and oligo-dT primed cDNA (b).

40

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

Fig. 6. Relative expression of vit, tub-α, and ddx genes in the hermaphrodite gonad of N. subnodosus, normalized with 18S and rp-l8 (n = 6 for each sex/gametogenic stage). Different letters indicate significant differences (p b 0.05); bars represent 95% confidence intervals.

muscle was confirmed by analyzing the melting peaks during the last step of the qPCR amplification program (data not shown). 4. Discussion 4.1. 18S is amplified from cDNA when oligo-dT is used as anchor primer The assumption that oligo-dT cDNA excludes ribosomal RNA in the reverse transcription process is widely mentioned in normal RTPCR protocols (O'Connell, 2002). Reverse transcription by using random primers is recommended when prokaryotic mRNA, or eukaryotic non-polyadenylated RNA (such as histone mRNA) or ribosomal RNA (such as 18S or 28S) will be targeted by qPCR. In this study we have demonstrated, by comparing the relative stability of this gene (among the rest) in cDNA sets primed by oligo-dT and random hexamer, that the 18S ribosomal RNA gene can be used as a reference gene in N. subnodosus when cDNA is synthesized with oligo-dT primer. It was also noted that the relative abundance of the 18S transcript, although lower in oligo-dT/cDNA than in random hexamer/cDNA, was higher than any of the other mRNA genes analyzed. As one of

the best characteristics of a reference gene is a high abundance (and hence small variance in qPCR replicates), these results support that the oligo-dT-primed 18S rRNA gene is reliably targeted by qPCR in the scallop N. subnodosus. Similarly, other authors have found that ribosomal RNA genes can be reverse transcribed using oligo-dT. For instance, Xue et al. (2010) evaluated 28S ribosomal RNA gene amplification in cDNA synthesized using oligo-dT or both oligo-dT and 28S antisense primers in insect cultured cells. The authors reported differences of 6 to 13 Cqs in the amplification of 28S, with the lowest Cq obtained in cDNA reverse-transcribed with oligo-dT alone. Even with this lower abundance of 28S rRNA, they noted that this gene was the most stable among a set of four genes when the cell lines were in vitro infected by viruses. Zhu and Altmann (2005), using diverse rat tissues, also observed that the18S rRNA gene is a target that could be detected by qPCR in oligo-dT/cDNA, even after finding that there was a 100-fold less amount of 18S rRNA transcript in oligo-dT/cDNA when compared with cDNA reverse transcribed with 18S antisense and oligo-dT primers. In our study, the difference between oligo-dT and random primers was 13-fold, consistent with the results reported by Xue et al. (2010) and Zhu and Altmann

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

41

Fig. 7. Relative expression of ddx estimated after using individual reference genes (18S, rp-l8, ef1-α, or act-β), the combination of the two most stable genes (18S/rp-l8), or the combination of all reference genes for normalization (n = 6 for each sex/gametogenic stage). Different letters within each estimated normalization indicate significant differences (p b 0.05); bars represent 95% confidence intervals.

(2005), less abundant in the oligo-dT than in random primed cDNA. Such small difference could be due to a more abundant poly(A) rich 18S rRNA in N. subnodosus than the insect and rat. Among possible reasons for ribosomal RNA to be targeted by PCR with oligo-dT/ cDNA could be: (1) the occurrence of true polyadenylation in the 3′-end of the transcript by RNA polymerase I or III, (2) presence of A-rich sites in mid-regions of the ribosomal transcript, or (3) by self-priming of the rRNA due to the complex tertiary structure that leads to rRNA detection by PCR, even without adding any kind of primer in the RT reaction, as shown by Tuiskunen et al. (2010). In the first case, it has been documented in the yeast Saccharomyces cerevisiae that low but significant levels of rRNA are in fact polyadenylated in the 3′-end by action of the Pap1p poly(A) polymerase, probably as part of a control mechanism of rRNA biogenesis and recycling (Kuai et al., 2004). In humans, poly(A) rich rRNA, not only in the 3′-end, but also in the middle regions of the rRNA transcript, has been documented (Slomovic et al., 2006; 2010). A survey of cDNA libraries demonstrated that the presence of ribosomal RNA is common, even after poly(A) selection and oligo-dT synthesis (González and Sylvester, 1997), suggesting also that the intermediate polyadenylation is a widespread process. Taken all these arguments together, we recommend that the use of 18S rRNA as a reference gene when amplified in oligo-dT cDNA should not be by de facto discarded until analyses are performed to assess its PCR efficiency and stability over all the tested samples.

reference genes by geNorm could be erroneous (De Santis et al., 2010), with a high risk of misinterpretation of the true relative expression of the target gene(s). Nevertheless, because of the sensibility of geNorm to detect pairwise relationships, we used this program to test if the expression of 18S and all other genes was equivalent when analyzed in two sets of cDNA (oligo-dT and random-primed), that is, to test the hypothesis that any gene should present a similar stability regardless of the priming method used for reverse transcription and the biological condition in which it was evaluated. We expected that equally variable

a

b

4.2. Reference gene stability by geNorm and NormFinder In species with limited or no genomic information available, gene expression analyses require directed isolation and sequencing (at least partial) of potential reference genes, sometimes considering those already referred in the literature as stable in similar conditions. 18S gene occurs in all eukaryotic cells and its sequence is well conserved among distinct taxa, which allowed for a partial fragment to be obtained in N. subnodosus. Together with 18S, we pre-selected a set of potential reference genes, including four protein-coding ones, and evaluated which ones were the most stable for normalization of target gene expression. For gene expression stability we used two software models, geNorm and NormFinder, because it has been shown that in small sets of reference genes, especially when two truly unstable genes are largely correlated, the selection of the best

Fig. 8. Stability values of reference genes between organs/tissues of N. subnodosus. geNorm (a) and NormFinder (b), results are presented for each gene.

42

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

Fig. 9. Relative expression of tub-α gene in different organs/tissues of N. subnodosus normalized with ef1-α and rp-l8. Different letters indicate significant differences (p b 0.05); bars represent 95% confidence intervals.

genes will be closely ranked independently of cDNA set, and this was in fact observed with rp-l8 gene ranking as the most stable one, 18S ranking as the second most stable gene, and ef1-α ranking as the third most stable gene. This confirmed that the relative stability of 18S is not affected when amplified with oligo-dT-primed cDNA. However, when we analyzed relative stability of each gene within each cDNA set in function of the priming strategy, geNorm pairwise analyses did not detect 18S and rp-l8 as the most stable genes when random-primed cDNA was analyzed, indicating ef1-α and rp-l8 as the most stable genes. This might be because geNorm software requires a larger set of reference genes (at least eight) to accurately determine the most stable gene set (Vandesompele et al., 2002). However, this could be taken as evidence that ef1-α is stable enough to be included as a part of a set of reference genes for gametogenesis comparisons of target genes, when in fact it is not. The conflict is derived because geNorm does not detect the inter-group variation as a source of bias for relative quantification analysis, hiding differences when genes are in fact differentially expressed. The inadequacy of ef1-α as a reference gene in different sexual regions of the gonad is evident when the coefficients of variation (CV) for relative expression of each gene are contrasted, as ef1-α was larger than 5%, whereas both 18S and rp-l8 CVs were less or equal to 5%. De Santis et al. (2010) and Zhong et al. (2008) observed that even when one pair of reference genes is expressed differentially between experimental groups, geNorm could erroneously select them as “stable” when their expression is correlated, because, as stated before, this software actually chooses the most similarly expressed genes, in some cases without considering potential differences in the expression ratio within and between experimental groups. When we analyzed the reference genes again across organs/tissues, the most stable ones were ef1-α and rp-l8, and in this case the lower coefficients of variation were for these genes, with 18S having a higher CV. These differences in the best reference genes over different physiological contexts (gonad development — gametogenesis and organ differentiation) could be the result of several regulatory mechanisms, but we hypothesize that the lowest transcription rate between germinal tissues in meiosis compared to somatic tissues (as stated by Barreau et al., 2008) alters the ef-1α stability over gametogenic stages and sexual portions of the gonad, but otherwise appears as the most stable among different somatic conditions. It is also possible that if transcriptional repression during meiosis occurs in the gonad, 18S and rp-l8 were found stable enough because ribosomal RNAs are highly and constantly expressed as they will be required

for the translation of nascent proteins (Raska et al., 2004). Regardless of the cause, it is expected that a different ranking in the same set of reference genes will be observed when their expression stability is systematically analyzed under different physiological conditions, tissues or during early development (Vandesompele et al., 2002).

4.3. Sex and gametogenic stage expression patterns of the target genes ddx, vit and tub-α Relative quantification of gene expression requires a robust selection of the most stable genes to assure that the interpretation of the results is correct. However, for non-model species, limited genomic knowledge results in a challenge when trying to evaluate a wide selection of potential reference genes. Nevertheless, partial sequences of genes known to be conserved could be obtained by direct sequencing after amplification using degenerated PCR. Ribosomal genes are known to have a high level of nucleotide conservancy across taxa, and as demonstrated in this work, both ribosomal genes were the most stable, and were therefore used for estimation of relative expression of the three target genes. As expected, vit expression occurred during oogenesis, as this gene codes for a lipoprotein expressed in the follicle cells of the ovary, and is accumulated in the yolk of growing oocytes (Osada et al., 2004). Although not significantly different, vit expression was 1.2 fold higher during late oogenesis in contrast to early oogenesis. Differently, tub-α was only expressed in the male region of the gonad during early and late spermatogeneses, similarly showing a trend to increase as spermatogenesis advances. tub-α and β are structural proteins conforming microtubules (cytoskeletal filaments), whose properties are different depending on the posttranslational modifications of the tubulins C-terminus proteins (Hammond et al., 2008; Kierszenbaum, 2002). There are at least two testis-specific tubulins participating in the flagella, in axoneme formation (Nielsen et al., 2010). Lack of expression in undifferentiated and female sexual region of the gonad, as well as the high expression level in the male gonad in which continuous cell divisions are occurring (spermatogonia to spermatocytes and to spermatid), culminating with the formation of the axoneme in the spermatozoa flagella, indicated first that this tubulin was a testis-specific expressed gene. Further expression analyses in somatic tissues corroborated this finding, although expression was also observed in adductor muscle.

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

4.4. tub-α expression is testis- and muscle-specific Tubulins compose a family of eukaryotic structural genes known to form microtubules. All eukaryotic cells require microtubules, the components of the cytoskeleton that mediate cell division, shape, motility, and intracellular trafficking (Nielsen et al., 2010). In sperm cells, tubulin specific isotypes of the β-subfamily were initially described first by Kemphues et al. (1979). Later Theurkauf et al. (1986) described four α-tubulins in Drosophila melanogaster, one of which was testis specific (α3). In mouse, Lewis and Cowan (1988) reported that tubulin α3/7 (identical isotypes) and β3 dimers were coexpressed only in mature testis, conforming most of the microtubule structures in male germ cells, but a different isotype, tub-a4 was found exclusively expressed in muscle. Later, Stanchi et al. (2000) reported a new α-tubulin (tuba8) in humans and mouse that was only expressed in adult skeletal muscle and testis. In the arthropod Bombyx mori, Kawasaki et al. (2003) found also that tub-a3 is testis-specific, and further phylogenetic studies confirmed that this isotype occurs in some arthropod lineages, and is related to the tub-a1 isotype family (Nielsen et al., 2010). In all cases, structural and mutation analyses of testis-specific tubulin isotypes leading to dysfunctional axonemal assemblies and male-sterile genotypes, confirmed the fundamental role of most testis-specific tubulins in the correct assembly of sperm axoneme (Fackenthal et al., 1995; Hutchens et al., 1997). Interestingly, it is currently known that expression of tissue-specific tubulins is mediated by their gene promoters (Bo and Wensink, 1989), a knowledge that has become useful in transgenic technology for biological control of pests and research inducing directed tissue-specific expression of transgenes using native promoters as those in α- and β-tubulins (Siebert et al., 2008; Zimowska et al., 2009; Lycett et al., 2011). We were unable to find in the available literature any testis-specific tubulin in mollusks. For Mytilus edulis, two α-tubulins are reported in GenBank; the first one (ABA03050.1), was recently used as a reference gene by Ciocan et al. (2011); the same authors reported the second α-tubulin isolated from testis (AEM36066.1) but not analyzed for tissue specificity and therefore not identified as testis-specific. This last partial sequence aligns with 100% similarity to the N. subnodosus testis- and muscle-specific α-tubulin we found (AEN71147.1), whereas no alignment occurred with the first α-tubulin reported for M. edulis. Also recently, a microarray study aimed to define stably expressed genes in tissues of the oyster C. gigas as candidates for normalization of target genes, found a member of the α-tubulin family stably expressed (Dheilly et al., 2011), pointing to the need of establishing stability of tubulin genes before using them in gene expression studies. Further investigation of α- and β-tubulin gene family members in N. subnodosus is required to assess their specificity and potential applications. 4.5. Consequences of using randomly chosen reference genes for normalization of target gene expression The expression analysis of a putative ddx family member is an interesting example of what could happen when relative normalization is performed with randomly chosen reference genes selected based only on previous literature reports. When we simulated what would occur when using different genes or combinations of genes we found multiple possibilities for interpreting the results, but overexpression of this gene during male late spermatogenesis was the most common occurrence. However, when relative expression of this gene was evaluated using the most stable genes (rp-l8 and 18S), it was clear that there were no differences between gonad part and developmental stages. ddx codes for a non-identified member of the DEAD-box helicase family that includes VASA and PL10 subfamilies of the ATP-dependent RNA helicases. Whereas vasa is known to be specifically expressed in germ cells, and has a functional role in germline formation in the Pacific oyster C. gigas (Fabioux et al., 2004), the partial ddx gene isolated in the present work had the best hit with the PL10 subfamily of the DEAD-

43

box helicases. In the penaeid shrimp Marsupenaeus japonicus a pl-10 gene has been found highly transcribed during embryo and larval developments, with its expression, although reduced, continuing in female and male adult gonads (Sellars et al., 2007). 5. Conclusions In this study, we demonstrated that the 18S rRNA gene of the scallop N. subnodosus is efficiently and stably reverse-transcribed in enzymatic reactions using oligo-dT as priming strategy. For the specific gametogenic conditions evaluated in this study, we identified that the18S rRNA gene is one of the most stable reference genes for studies involving different gametogenic stages and sex in gonad tissue, and together with the ribosomal protein rp-l8 gene, provided the best gene expression normalization of three target genes in the gonad tissue, and the most robust interpretation of the results. However, when gene expression was evaluated between different tissues, rp-l8 and a second reference gene, elf1a, were found as the most stably expressed, pointing to the importance of evaluating gene stability before using any gene for normalization of expression. We also found that a gene first thought as a potential reference gene, tubulin-α, was in fact tissue- and sexspecific gene in N. subnodosus, expressed only in male testis and, in a lesser extent, also in adductor muscle. Acknowledgments We thank the technical assistance provided by José L. Ramírez and Susana Ávila. This work was funded by CIBNOR (grant AC0.16) and CONACYT-CB-2011-1 (grant 166573). The first author is a Ph.D. student, financially supported by a CONACyT fellowship (no. 204031). Appendix A. Supplementary data Supplementary data to this article can be found online at http:// dx.doi.org/10.1016/j.margen.2012.03.003. References Andersen, C.L., Jensen, J.L., Ørntoft, T.F., 2004. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 64, 5245–5250. Barreau, C., Benson, E., Gudmannsdottir, E., Newton, F., White-Cooper, H., 2008. Postmeiotic transcription in Drosophila testes. Development 135, 1897–1902. Bo, J., Wensink, P.C., 1989. The promoter region of the Drosophila α2-tubulin gene directs testicular and neural specific expression. Development 106, 581–587. Bustin, S.A., Benes, V., Nolan, T., Pfaffl, M.W., 2005. Quantitative real-time RT-PCR — a perspective. J. Mol. Endocrinol. 34, 597–601. Bustin, S.A., Benes, V., Garson, J.A., Hellemans, J., Huggett, J., Kubista, M., Mueller, R., Nolan, T., Pfaffl, M., Shipley, G.L., Vandesompele, J., Wittwer, C.T., 2009. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55, 611–622. Ciocan, C.M., Cubero-Leon, E., Minier, C., Rotchell, J.M., 2011. Identification of reproductionspecific genes associated with maturation and estrogen exposure in a marine bivalve Mytilus edulis. PloS one. 6, e22326, http://dx.doi.org/10.1371/journal.pone.0022326. De Santis, C., Smith-Keune, C., Jerry, D.R., 2010. Normalizing RT-qPCR data: are we getting the right answers? An appraisal of normalization approaches and internal reference genes from a case study in the finfish Lates calcarifer. Mar. Biotechnol. 13, 170–180. Dheilly, N.M., Lelong, C., Huvet, A., Favrel, P., 2011. Development of a Pacific oyster (Crassostrea gigas) 31,918-gene microarray: identification of reference genes and tissue-specific expression patterns. BMC Genomics 12, 468. Fabioux, C., Huvet, A., Lelong, C., Robert, R., Pouvreau, S., Daniel, J.Y., Minguant, C., Le Pennec, M., 2004. Oyster vasa-like gene as a marker of the germline cell development in Crassostrea gigas. Biochem. Biophys. Res. Commun. 320, 592–598. Fackenthal, J.D., Hutchens, J.A., Turner, F.R., Raff, E.C., 1995. Structural analysis of mutations in the Drosophila b2-tubulin isoform reveals regions in the b-tubulin molecule required for general and for tissue-specific microtubule functions. Genetics 139, 267–286. Finn, R.D., Mistry, J., Tate, J., Coggill, P., Heger, A., Pollington, J.E., Gavin, O.L., Gunasekaran, P., Ceric, G., Forslund, K., Holm, L., Sonnhammer, E.L.L., Eddy, S.R., Bateman, A., 2010. The Pfam protein database. Nucleic Acids Res. 38, D211–D222. González, I.L., Sylvester, J.E., 1997. Incognito rRNA and rDNA in databases and libraries. Genome Res. 7, 65–70.

44

R. Llera-Herrera et al. / Marine Genomics 6 (2012) 33–44

Hammond, J.W., Cai, D., Verhey, K.J., 2008. Tubulin modifications and their cellular functions. Curr. Opin. Cell Biol. 20, 71–76. Hansen, K.D., Brenner, S.E., Dudoit, S., 2010. Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic Acids Res. 38, e131. Hellemans, J., Mortier, G., De Paepe, A., Speleman, F., Vandesompele, J., 2007. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 8, R19. Hernández-Arteaga, S., López-Revilla, R., 2008. Quantitation of human papillomavirus type 16 E6 oncogene sequences by real-time or quantitative PCR with EvaGreen. Anal. Biochem. 380, 131–133. Hutchens, J.A., Hoyle, H.D., Turner, F.R., Raff, E.C., 1997. Structurally similar Drosophila alpha-tubulins are functionally distinct in vivo. Mol. Biol. Cell 8, 481–500. Ibarra, A.M., Petersen, J.L., Famula, T.R., May, B., 2006. Characterization of 35 microsatellite loci in the Pacific lion-paw scallop (Nodipecten subnodosus) and their crossspecies amplification in four other scallops of the Pectinidae family. Mol. Ecol. Notes 6, 153–156. Jiang, M., Ryu, J., Kiraly, M., Duke, K., Reinke, V., Kim, S.K., 2001. Genome-wide analysis of developmental and sex-regulated gene expression profiles in Caenorhabditis elegans. Proc. Natl. Acad. Sci. 98, 218–233. Kawasaki, H., Sugaya, K., Quan, G.X., Nohata, J., Mita, K., 2003. Analysis of alpha- and beta-tubulin genes of Bombyx mori using an EST database. Insect Biochem. Mol. Biol. 33, 131–137. Kemphues, K.J., Raff, R.A., Kaufman, T.C., Raff, E.C., 1979. Mutation in a structural gene for a beta-tubulin specific to testis in Drosophila melanogaster. Proc. Natl. Acad. Sci. 76, 3991–3995. Kierszenbaum, A.L., 2002. Sperm axoneme: a tale of tubulin posttranslation diversity. Mol. Reprod. Dev. 62, 1–3. Kortner, T.M., Valen, E.C., Kortner, H., Marjara, I.S., Krogdahl, A., Bakke, A.M., 2011. Candidate reference genes for quantitative real-time PCR (qPCR) assays during development of a diet-related enteropathy in Atlantic salmon (Salmo salar L.) and the potential pitfalls of uncritical use of normalization software tools. Aquaculture 318, 355–363. Kuai, L., Fang, F., Butler, J.S., Sherman, F., 2004. Polyadenylation of rRNA in Saccharomyces cerevisiae. Proc. Natl. Acad. Sci. 101, 8581–8586. Lewis, S.A., Cowan, N.J., 1988. Complex regulation and functional versatility of mammalian α- and β-tubulin isotypes during the differentiation of testis and muscle cells. J. Cell Biol. 106, 2011–2022. Llera-Herrera, R., García-Gasca, A., Romero-Vivas, E., Huvet, A., Ibarra, A.M., 2010. Identification and comparative expression of gametogenic-related genes in the ovary and testis of the hermaphrodite lion-paw scallop Nodipecten subnodosus. Physiomar '10 Workshop. Quebec City, Canada. Lu, X., Wu, X.-I., 2011. Sex, sex chromosomes and gene expression. BMC Biol. 9, 30. Lycett, G.J., Amenya, D., Lynd, A., 2011. The Anopheles gambiae alpha-tubulin-1b promoter directs neuronal, testes and developing imaginal tissue specific expression and is a sensitive enhancer detector. Insect Mol. Biol., http://dx.doi.org/10.1111/ j.1365-2583.2011.01112.x. Nielsen, M.G., Gadagkar, S.R., Gutzwiller, L., 2010. Tubulin evolution in insects: gene duplication and subfunctionalization provide specialized isoforms in a functionally constrained gene family. BMC Evol. Biol. 10, 113. Nolan, T., Bustin, S., 2008. Procedures for quality control of RNA samples for use in quantitative reverse transcription PCR. Essentials of Nucleic Acid Analysis: a Robust Approach. The Royal Society of Chemistry, Cambridge, UK. O'Connell,, J., 2002. The basics of RT-PCR: some practical considerations. RT-PCR Protocols. Humana Press, Totowa, NJ. Osada, M., Harata, M., Kishida, M., Kijima, A., 2004. Molecular cloning and expression analysis of vitellogenin in scallop, Patinopecten yessoensis (Bivalvia, Mollusca). Mol. Reprod. Dev. 67, 273–281. Petersen, J.L., Ibarra, A.M., Ramirez, J.L., May, B., 2008. An induced mass spawn of the hermaphroditic lion-paw scallop, Nodipecten subnodosus: genetic assignment of maternal and paternal parentage. J. Hered. 99, 337–348. Petersen, J.L., Baerwald, M.R., Ibarra, A.M., May, B., 2009a. Discovering QTL's underlying growth traits in the Pacific lion-paw scallop using a first generation microsatellite and AFLP linkage map. Plant & Animal Genome XVII Conference. January 10–14, 2009, San Diego California. Abstract. Petersen, J.L., Ibarra, A.M., May, B., 2009b. Thirty-seven additional microsatellite loci in the Pacific lion-paw scallop (Nodipecten subnodosus) and cross-amplification in other pectinids. Cons. Gen. Res. 1, 101–105.

Petersen, J.L., Baerwald, M.R., Ibarra, A.M., May, B., 2010a. A genetic linkage map of the Pacific lion-paw scallop — steps toward marker assisted selection in aquaculture. Aquaculture 2010, San Diego CA USA, March 1–5, 2010. Abstract. Petersen, J.L., Ibarra, A.M., May, B., 2010b. Nuclear and mtDNA lineage diversity in wild and cultured Pacific lion-paw scallop, Nodipecten subnodosus (Baja California Peninsula, Mexico). Mar. Biol. 157, 2751–2767. Pfaffl, M.W., 2004. Quantification strategies in real-time PCR. In: Bustin, S.A. (Ed.), AZ of Quantitative PCR. IUL Biotechnology, La Jolla, CA, pp. 87–120. Pfaffl, M.W., Tichopad, A., Prgomet, C., Neuvians, T.P., 2004. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: bestkeeper – Excel-based tool using pair-wise correlations. Biotechnol. Lett. 26, 509–515. Racotta, I.S., Ramirez, J.L., Ibarra, A.M., Rodriguez-Jaramillo, M.C., Carreño, D., Palacios, E., 2003. Growth and gametogenesis in the lion-paw scallop Nodipecten (Lyropecten) subnodosus. Aquaculture 217, 335–349. Radonic, A., Thulke, S., Mackay, I.M., Landt, O., Siegert, W., Nitsche, A., 2004. Guideline to reference gene selection for quantitative real-time PCR. Biochem. Biophys. Res. Commun. 313, 856–862. Raska, I., Koberna, K., Malínský, J., Fidlerová, H., Masata, M., 2004. The nucleolus and transcription of ribosomal genes. Biol. Cell 96, 579–594. Ruijter, J.M., Ramakers, C., Hoogaars, W.M.H., Karlen, Y., Bakker, O., van den Hoff, M.J.B., Moorman, A.F.M., 2009. Amplification efficiency: linking baseline and bias in the analysis of quantitative PCR data. Nucleic Acids Res. 37, e45. Sambrook, J., Russel, D.W., 2001. Molecular Cloning: A Laboratory Manual. Cold Spring Harbor Laboratory, Cold Spring Harbor, N.Y. Sellars, M., Lyons, R., Grewe, P., Vuocolo, T., Leeton, L., Coman, G., 2007. A PL10 vasa-like gene in the kuruma shrimp, Marsupenaeus japonicus, expressed during development and in adult gonad. Mar. Biotechnol. 9, 377–387. Siebert, K.S., Lorenzen, M.D., Brown, S.J., Park, Y., Beeman, R.W., 2008. Tubulin superfamily genes in Tribolium castaneum and the use of a tubulin promoter to drive transgene expression. Insect Biochem. Mol. Biol. 38, 749–755. Silveira, E.D., Alves-Ferreira, M., Guimarães, L.A., da Silva, F.R., Carneiro, V., 2009. Selection of reference genes for quantitative real-time PCR expression studies in the apomictic and sexual grass Brachiaria brizantha. BMC Plant Biol. 9, 84. Slomovic, S., Laufer, D., Geiger, D., Schuster, G., 2006. Polyadenylation of ribosomal RNA in human cells. Nucleic Acids Res. 34, 2966–2975. Slomovic, S., Fremder, E., Staals, R.H.G., Pruijn, G.J.M., Schuster, G., 2010. Addition of poly(A) and poly(A)-rich tails during RNA degradation in the cytoplasm of human cells. Proc. Natl. Acad. Sci. 107, 7407–7412. Stanchi, F., Corso, V., Scannapieco, P., Ievolella, C., Negrisolo, E., Tiso, N., Lanfranchi, G., Valle, G., 2000. TUBA8: a new tissue-specific isoform of alpha-tubulin that is highly conserved in human and mouse. Biochem. Biophys. Res. 270, 1111–1118. Theurkauf, W.E., Baum, H., Bo, J., Wensink, P.C., 1986. Tissue-specific and constitutive α-tubulin genes of Drosophila melanogaster code for structurally distinct proteins. Proc. Natl. Acad. Sci. 83, 8477–8481. Tuiskunen, A., Leparc-Goffart, I., Boubis, L., Monteil, V., Klingström, J., Tolou, H.J., Lundkvist, A., Plomet, S., 2010. Self-priming of reverse transcriptase impairs strand-specific detection of dengue virus RNA. J. Gen. Virol. 91, 1019–1027. Untergasser, A., Nijveen, H., Rao, X., Bisseling, T., Geurts, R., Leunissen, J.A.M., 2007. Primer3Plus, an enhanced web interface to primer3. Nucleic Acids Res. 35, W71–W74. Vandesompele, J., Preter, K.D., Pattyn, F., Poppe, B., Roy, N.V., Paepe, A.D., Speleman, F., 2002. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3 (research0034.1-0034.11). Xue, J.-L., Salem, T.Z., Turney, C.M., Cheng, X.-W., 2010. Strategy of the use of 28S rRNA as a housekeeping gene in real-time quantitative PCR analysis of gene transcription in insect cells infected by viruses. J. Gen. Virol. 163, 210–215. Zhong, Q., Zhang, Q., Wang, Z., Qi, J., Chen, Y., Li, S., 2008. Expression profiling and validation of potential reference genes during Paralichthys olivaceus embryogenesis. Mar. Biotechnol. 10, 310–318. Zhu, L.J., Altmann, S.W., 2005. mRNA and 18S-RNA coapplication-reverse transcription for quantitative gene expression analysis. Anal. Biochem. 345, 102–109. Zimowska, G.J., Nirmala, X., Handler, A.M., 2009. The β2-tubulin gene from three tephritid fruit fly species and use of its promoter for sperm marking. Insect Biochem. Mol. Biol. 39, 508–515.