Validation of housekeeping genes for quantitative mRNA expression analysis in OsHV-1 infected ark clam, Scapharca broughtonii

Validation of housekeeping genes for quantitative mRNA expression analysis in OsHV-1 infected ark clam, Scapharca broughtonii

Accepted Manuscript Validation of housekeeping genes for quantitative mRNA expression analysis in OsHV-1 infected ark clam, Scapharca broughtonii Lush...

5MB Sizes 0 Downloads 45 Views

Accepted Manuscript Validation of housekeeping genes for quantitative mRNA expression analysis in OsHV-1 infected ark clam, Scapharca broughtonii Lusheng Xin, Bowen Huang, Changming Bai, Chongming Wang PII: DOI: Reference:

S0022-2011(17)30436-6 https://doi.org/10.1016/j.jip.2018.04.011 YJIPA 7093

To appear in:

Journal of Invertebrate Pathology

Received Date: Revised Date: Accepted Date:

19 October 2017 24 April 2018 30 April 2018

Please cite this article as: Xin, L., Huang, B., Bai, C., Wang, C., Validation of housekeeping genes for quantitative mRNA expression analysis in OsHV-1 infected ark clam, Scapharca broughtonii, Journal of Invertebrate Pathology (2018), doi: https://doi.org/10.1016/j.jip.2018.04.011

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Validation of housekeeping genes for quantitative mRNA expression analysis in OsHV-1 infected ark clam, Scapharca broughtonii Lusheng Xin a, b,1, Bowen Huang c,1, Changming Bai a, b, Chongming Wang a, b* a

Qingdao Key Laboratory of Mariculture Epidemiology and Biosecurity, Key

Laboratory of Maricultural Organism Disease Control, Ministry of Agriculture, Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Qingdao 266071, China; b

Function Laboratory for Marine Fisheries Science and Food Production Processes,

Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China c

Tianjing Agriculture University, Tianjing 300191, China

*

Correspondence authors:

Tel: 86-532-85823062 E-mail: [email protected] (C. Wang) 1

These authors contributed equally to this work.

1

Abstract Ostreid herpesvirus-1 (OsHV-1) presents interspecies transmission among bivalves. Recently, events of mass mortalities of ark clams (Scapharca broughtonii) infected with OsHV-1 have been recorded. To accurately assess the gene responding patterns of ark clams post OsHV-1 infection, constant stable housekeeping genes (HKGs) are needed as internal control to normalize raw mRNA expression data. In this study, ten candidate HKGs were selected, including 18S rRNA (18S), beta-actin (ACT), Glyceraldehyde 3-phosphate dehydrogenase (GAPDH), NADH dehydrogenase subunit (NADH), Elongation factor-1a (EF-1a), Elongation factor-1β (EF-1β), Elongation factor-1γ (EF-1γ), Ribosomal protein L7 (RL7), Ribosomal protein L15 (RL15) and Ribosomal protein S18 (S18). The expression levels of ten candidate HKGs were analyzed by real-time PCR under given experimental conditions, including various tissues, OsHV-1 challenge, temperature stress and OsHV-1 challenge at different temperature. Their expression stability values were further calculated using two different statistical models (geNorm and NormFinder). The results showed that different tissues presented distinct best pair genes combinations for gene expression analysis under OsHV-1 challenge. RL15 was comparatively more stable than other HKGs under various experimental conditions, while commonly used 18s and ACT seemed to be more greatly influenced by most given experimental conditions in ark clams. This study emphasized the necessity of prior validation of HKGs and would facilitate future gene expression analysis in ark clams or other shellfishes. Keywords: housekeeping gene; Ostreid Herpesvirus; Scapharca broughtonii; shellfish; 2

Real-time PCR

1. Introduction Ostreid herpesvirus 1 (OsHV-1) was characterized as a major class of the herpesviruses (Davison et al., 2005), and emerged as a disastrous pathogen in shellfish aquaculture. OsHV-1 was often involved with mass mortality events of ostreid species (Keeling et al., 2014; Segarra et al., 2010). Besides vertical transmission across offspring (Barbosasolomieu et al., 2005), OsHV-1 appears capable of interspecies transmission (Arzul et al., 2001). Subsequently, other shellfish species were also characterized with OsHV-1 infection cases, including ark clam (Scapharca broughtonii) (Xia et al., 2015). Ark clam is distributed widely along the northwest pacific coast, and has become one of the most commercially important shellfish cultivated in China and other northern Asia countries (Yokogawa, 1997). Mass mortalities of ark clams were recorded and resulted in great economic loss during the early summer of 2012 and 2013, and had been confirmed to be associated with OsHV-1 infection (Bai et al., 2016). Until now, the interaction between pathogen OsHV-1 and host ark clam is largely unknown, specific responding genes as biomarkers remain to be further explored for early warning and effective management against OsHV-1 prevalence in ark clams. Though several methods have been established for the diagnosis of OsHV-1 infection (e.g., transmission electron microscopy, in situ hybridization and PCR diagnostic methods) (Barbosa-Solomieu et al., 2004; Oden et al., 2011; Pepin et al., 2008), these methods could only sense OsHV-1 particles after onset and are not available for 3

distinguishing vulnerable population. Specifically responding genes in OsHV-1 infected ark clams are expected to become much earlier, more sensitive and effective diagnostic objects for OsHV-1 detection, and further guide the acknowledgement of host-pathogen interaction. For searching these genes, mRNA expression profile before and after OsHV-1 infection is a critical cue, which could be assessed by semi-quantitative PCR, Northern Blot, RT-PCR and so on. In comparison, RT-PCR is more frequently used for its high convenience, sensitivity and reproducibility (Nolan et al., 2006). Even mRNA expression levels were quantified by costly high-throughput sequencing, such as microarray and RNA-seq, still further need to be confirmed by RT-PCR. Thus RT-PCR serves an important tool for identifying genes responding to OsHV-1 infection. Due to the unmanageable variations of initial total RNA amount and cDNA synthesis efficiency, constant stable expression HKGs are critical to normalize raw data for RT-PCR analysis (Dheda et al., 2004). Microsoft Excel geNorm and NormFinder were mostly employed to determine HKG stability and suitability under given experimental conditions, ranking the candidate genes according to the calculated expression stability values, more stable genes present lower expression stability values (Gutierrez et al., 2008; Jensen and Ørntoft, 2004; Pan et al., 2015). The accuracy and repeatability of target gene expression results are strongly dependent on selected HKGs. Least stability HKGs could even cover the significant fluctuation of target genes expression (Lacroix et al., 2014). However, limited previous studies were reported about the selection of HKGs in molluscs (Table. 2). 4

Available HKGs were only in several shellfish species, such as Crassostrea gigas (Dheilly et al., 2011), Ostrea edulis (Morga et al., 2010) and Mytilus edulis (Cubero-Leon et al., 2011). In addition, to date, most validated HKGs in molluscs exposed to pathogens (Vibrio splendidus, Bonamia ostreae) were tested in haemocytes only (Araya et al., 2008; Morga et al., 2010). HKGs were validated in oyster larvae and spat under OsHV-1 infection (Du et al., 2013; Segarra et al., 2014). Many of these HKGs may not be stably expressed across species under different conditions (Araya et al., 2008; Mauriz et al., 2012). Thus prior validation is required when a HKG is suggested to be used in a new given experimental condition and can not be ignored. There are no validated HKGs that have been suggested in ark clams up to now. Whereas previous studies mainly selected HKGs, β-actin or 28s, for gene expression analyses in ark clams without validation (Zheng et al., 2016). The suitability of these two HKGs was doubtful, because they exhibited significant fluctuation of mRNA expression levels and were abandoned under some given experimental conditions, such as white spot syndrome virus (WSSV) challenge shrimp and mussel various stages of gametogenesis (Cubero-Leon et al., 2011; Dhar et al., 2009). In this study, 10 candidate HKGs were selected based on molecular function and stable gene expression levels from our previous RNA-seq data. Their expression stability values were further assessed among different tissues and OsHV-1 challenge. Temperature, an important trigger factor of OsHV-1 disease outbreak (above 13°C) (Pernet et al., 2015; Petton et al., 2013), was also taken into account as the given experimental condition under that HKG expression stability was checked. The aim of this study was to 5

suggest available HKGs, referring to tissue distribution, OsHV-1 challenge or temperature stress in ark clams.

2. Materials and Method 2.1 Ark clams and tissue sampling Healthy ark clams (two-year-old) with an average shell length of 55mm were obtained from a local farm in Qingdao, China. Prior to the following treatment, ark clams were acclimated in tanks (10 ark shells in each tank) with aerated sea water at 18.0°C for 13 days. For tissue distribution analysis of gene expression, tissue samples (including haemocyte, mantle, gill, digestive gland, adductor muscle and foot) were collected from 6 ark clams and placed into TRIzol (Invitrogen, USA) for later RNA isolation. Three replicates were employed for each tissue sample. 2.2 Experimental infection and temperature stress For experimental infection, viral inoculum was firstly prepared using OsHV-1 infected ark clams as previous description using oyster with modification (Schikorski et al., 2011). Briefly, mantle tissue (50mg each individual) was sampled from six OsHV-1 infected ark clams and homogenized. After a short centrifugation (4°C, 1000g), the supernatant of mantle homogenate was remained and successively filtered through a series of syringe filters with pore size at 5μm, 2μm, 0.45μm and 0.22μm. The filtered tissue homogenate was used as a viral inoculum and stored for 2 days at 4°C. Meanwhile, un-infected ark clams were used to make negative tissue 6

homogenates as the negative control. After that, 200 μL of each tissue homogenate was used for DNA extraction and OsHV-1 quantification by qPCR according to the previously described protocol (Bai et al., 2016). Two hundred and forty ark clams were randomly selected and separated into two groups (HN and HP). Both were remained at 18°C, group HP was injected with 100 μL viral inoculum (~105 copies of viral DNA/μL) into the foot of each ark clam. Correspondingly the other group HN was injected with the same volume of negative tissue homogenate as the negative control. Tissue samples (haemocyte, mantle, gill, digestive gland, adductor muscle and foot) of six ark clams were respectively collected from each group at 0h and every 12h for 72h post injection. For temperature stress analysis, after 13 days at 18 °C, one hundred arks shells were transferred to aerated sea water at 12°C. After continuous culture for 48h, haemocyte and mantle samples were collected from six ark clams at 12°C and parallel six at 18°C respectively. In order to analyze gene expression under OsHV-1 challenge at different temperature, ninety ark clams were also randomly selected and separated into two groups (LN and LP) before being transferred to aerated sea water at 12°C, group LP received the same injection as group HP, while group LN as group HN. After undergoing 48h culture at 12°C, haemocyte and mantle samples of six ark clams were also collected from group LN and LP, to compare with group HN and HP remained at 18°C for 48h . All collected tissue samples were paralleled with three replicates. 2.3 RNA extraction and cDNA synthesis 7

Total RNA was extracted from tissue samples using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s protocol. Genomic DNA was eliminated by digestion with RNase-free DNase I (Promega). The quality and quantity of RNA samples were assessed by Nanodrop 2000 spectrophotometer (Thermo). Only RNA samples with a 260/280 ratio among 1.9-2.1 and a 260/230 ratio among 2.0-2.5 were used for analysis. The integrity of RNA samples was confirmed by agarose gel electrophoresis. First-strand cDNA was synthesized based on M-MLV reverse transcriptase (Promega) by using 1 μg total RNA as template according to the manufacturer’s instructions. The cDNA products were stored at -80°C for the subsequent real-time PCR. 2.4 Selection of candidate HKGs and primer design Ten candidate HKGs were selected according to their functions and expression fluctuations in our previous transcriptome data of ark clam haemocytes under given experimental conditions. Nucleotide sequences of ACT, GAPDH, NADH, EF-1a, EF-1β, EF-1γ, RL7, RL15 and S18 were obtained from our previous transcriptome data and confirmed by cloning and sequencing, while that of 18s was gained from GenBank database (GenBank:JN974500.1). Primers were designed using Primer Premier 5 software (PREMIER Biosoft). All the amplicons of selected candidate HKGs had been purified from agarose gel, and ligated to pMD 19-T Vector (TAKARA) for sequencing with six replicates in Sangon Biotech Company, China. Detailed information on these candidate HKGs and all the primers used for real-time PCR were shown in Table 1. 8

2.5 Real-time PCR analysis Real-time PCR was performed using Bio-Rad CFX Connect RealTime system (Bio-Rad Laboratories). Firstly, primer specificities were checked using melting curves and gel picture of corresponding amplicons. PCR efficiency was calculated for each primer pair by making standard curves of 10-fold diluted cDNA (1 to 1:10,000), using the formula: E = 10-1/slope (Pfaffl, 2001). A 20μl RT-PCR reaction volume contained 2 μL diluted cDNA template (1:10), 10 μL 2×SYBR Green PCR Mix (Takara), 0.5 μL of each gene-specific primer (250 nM) and 7 μL distilled water. Cycling conditions were as follows: 95°C for 30s, 40 cycles of 95°C for 5 s and 65°C for 35s and then a melt curve stage: 20 s at 95 °C and 30 s at 60 °C followed by temperature increasing (0.5°C per 10s) into 95°C. RT-PCR was performed in triplicate for each sample. 2.6 Gene expression stability analysis The Cq values were transformed to linear scale expression quantities as the form of 2−ΔCq. ΔCq value equals Cq value of HKGs obtained from each sample subtracting the lowest one (ΔCq= Cqeach- Cqmin). Then the data that stands for HKG expression level was input into and analyzed by Microsoft Excel geNorm (version win_3.5) (Pattyn et al., 2003) and NormFinder (version 0953) (Jensen and Ørntoft, 2004). Gene expression stability values were directly output. Genes with the lowest stability value have the most stable expression. Statistical analyses were conducted with the IBM SPSS Statistics 19.0 package, the box-and-whisker plots were obtained by this software. 9

3. Result 3.1. Selected candidate HKGs Ten candidate HKGs, including ACT, GAPDH, NADH, EF-1a, EF-1β, EF-1γ, RL7, RL15, S18 and 18s, were selected and further assessed. NADH, EF-1a, EF-1β, EF-1γ, RL7, RL15, and S18 showed parallel colors in the heatmap (Fig.1), indicating similar expression profiles under given experimental conditions based on our previous RNA-seq data analysis. Though ACT and GAPDH emerged distinct color difference in the heatmap (Fig.1), their molecular functions were fundamental and indispensable to sustain cells, and had been widely used as HKGs together with 18s. Thus ACT, GAPDH and 18s were also taken into account in our studies. Besides nucleotide sequence of 18s had been released and directly gained in the GeneBank database (accession: JN974499), the others were acquired from our RNA-seq assembled unigenes. All selected candidate HKGs exhibited high identities (>70% in nucleotide sequence) to their respective homologs in NCBI database, the sequences of ten candidate HKGs were cloned and confirmed by sequencing to facilitate subsequent RT-PCR primer design. 3.2. Amplification specificity and average expression levels of candidate HKGs The specificity of RT-PCR primers was confirmed by a signal band in each lane of the agarose gel with expected amplicon length (117-383bp) (Fig.2), and by corresponding melting curve with a single peak (Fig.S2). All the amplicons of candidate HKGs were confirmed by sequencing. Amplification efficiencies of 10

RT-PCR primers ranged between 95.9% and 104.6% (Table 3), S18 showed the lowest (95.9 %). The Pearson's moment correlation coefficients were 0.990. Box plots showed the comparative expression levels of ten candidate HKGs (Cq values) among six tissues (haemocyte, mantle, gill, digestive gland, adductor muscles and foot) in ark clams (Fig.3). ACT and 18s exhibited higher degree of dispersion with greater box lengths among ten candidate HKGs. Meanwhile, 18s showed the lowest Cq median indicating the highest average gene expression level among six tissues, which was in accordance with the lowest Cq value of 18s measured in the mixed six tissue samples (Fig.S1). 3.3. Expression stability of candidate HKGs among various tissues Among various tissues, geNorm was further used to calculate the gene expression stability values (M) as shown in Fig.4. The maximum threshold of M value was recommended at 1.5 for suitable HKGs selection, lower M value presents greater gene stability. Among various tissues, the M values of 18s, ACT, EF-1β and S18 (M=3.07, 2.49, 1.74 and 1.55, respectively) were above the threshold (shown by the dot line), while the other six were below that (Fig.4a). Gene pairs RL7 and EF-1a emerged as the most stable genes across tissues. Gene expression stability was also evaluated by NormFinder based on their minimum variance combined intra and inter groups (model-based approach) as shown in Fig.7, more stable genes also present lower values of average expression stability (listed in Table.3). NormFinder also ranked 18s, ACT, EF-1β and S18 as the top four fluctuant genes among different tissues, and RL7 represented in the top three stable genes. 11

To determine the optimal number of reference genes required for normalization, geNorm software was used by means of pairwise variation (Vn/n+1) calculated between the normalization factors NFn and NFn+1( using n or n+1 reference genes respectively)(Vandesompele et al., 2002). A cut-off value for the pairwise variation was suggested at 0.15, below this value indicating an additional HKG was not required. For example, V2/3 <0.15, two HKGs would satisfy the analysis; V2/3 >0.15 and V3/4<0.15, three HKGs were needed. However, among various tissues of ark clams, the values of pairwise variation Vn/n+1 were all out of range (the least two, V3/4 =0.230 and V5/6=0.230) (Fig.4b). According to the geNorm software manual, the proposed 0.15 value must not be taken as a too strict cut-off. Three best HKGs (RL7, EF-1a and EF-1γ) were suggested for gene expression analysis in this case. 3.4. Expression stability of candidate HKGs under stress Under the given condition of OsHV-1 challenge, the M values of ten candidate HKGs were calculated by geNorm (Fig.5). In haemocytes, mantle tissue, adductor muscles and foots, pair genes RL7 and RL15 (M=0.206), S18 and EF-1β (M=0.192), EF-1a and EF-1γ (M=0.224), S18 and EF-1β (M=0.154) respectively emerged as the most stable combinations under OsHV-1 challenge (Fig.5A, B, C and D). The values of pairwise variation V2/3 were all below the cut-off value (0.15) in above four tissues under OsHV-1 challenge (Fig.5a, b, c and d), which indicated two genes could satisfy RT-PCR data normalization. Under temperature stress, the M values of ten candidate HKGs were ranged from 0.044 to 0.376 and 0.061 to 0.233 in haemocytes and mantle tissue, respectively (Fig.6 A and B), and ranged from 0.060 to 0.413 and 12

0.040 to 0.242 respectively under OsHV-1 challenge at different temperature (Fig.6 C and D). The values of pairwise variation V2/3 were also all below 0.15 in haemocytes and mantle tissue under above two conditions (Fig.6 a, b, c and d), indicating two HKGs could satisfy gene expression analyses. In gills and digestive glands, the mRNA expression levels of three candidate HKGs (ACT, NADH and GAPDH) were not detected at some time points under OsHV-1 challenge and these genes were not taken into account for the subsequent analyses in gills and digestive glands (Fig.5 E and F).The M values of the remained seven candidate HKGs were ranked in gills as follows: 18s (M=1.651), RL7 (M=1.250), EF-1γ (M=1.200), RL15 (M=1.120), EF-1a (M=0.988), S18 (M=0.473) and EF-1β (M=0.473) (Fig.5 E); while in digestive gland as follows: RL7 (M=1.591), RL15 (M=1.469), 18s (M=1.131), EF-1a (M=0.871), EF-1β (M=0.685), EF-1γ (M=0.318) and S18 (M=0.318) (Fig.5 F). The values of pairwise variation Vn/n+1 were all above 0.15 in gills and digestive glands. V5/6 (0.186) and V3/4 (0.243) appeared the lowest in gills and digestive glands respectively (Fig.5 e and f). The expression stability values of candidate HKGs under stress were also calculated by NormFinder (Table.3), and the ranking of ten candidate HKGs under given experimental conditions were similar to the results calculated by geNorm. S18 and EF-1γ presented most frequently (six and seven times, respectively) in top three stable genes in tissue samples under stress, flowed by EF-1a, EF-1β and RL15 that all presented in top three for four times (shown in Table.3). Under tested experimental conditions, 18s and ACT showed much larger amount of expression stability values 13

(total up to 7.215 and 5.034) than the others, while RL15 was the smallest (total up to 1.646) (Fig.7).

4. Discussion The current study selected ten candidate HKGs from ark clams, and assessed the influence of factors, tissue distribution, OsHV-1 challenge and temperature stress, on the expression stability of ten candidate HKGs. In ark clams, pair genes RL7 and EF-1a were assessed as the most stable genes among various tissues. While pair genes GAPDH and ARF1 (ADP-ribosylation factor 1) in oyster, pair genes RPL5 (ribosomal protein L5) and CY (cyclophilin) in abalone appeared most stable among various tisssues (Table. 2). Besides not being assessed in oysters, 60S ribosomal protein subunit related genes (RL7 and RPL5) seemed more stable across tissues than others in mollusk species. Some common HKGs gene expression levels exhibited great variation among different tissues, especially 18s and ACT which were least stable (based on stability values) in ark clams. Though 18s and ACT are commonly used as internal controls for gene expression analysis (Deane and Woo, 2005; Kang et al., 2006), their suitability and stability should be re-assessed. To our knowledge, our research was the first to determine the suitability of the 18S gene across tissue types in mollusk species. ACT presented as one of top two least stable genes among various tissues in oyster and abalone (Dheilly et al., 2011; López-Landavery et al., 2014), which was different from fish (such as zeabrafish and Japanese flounder) ACT exhibiting more stable across tissues among tested candidate 14

HKGs (Mccurley and Callard, 2008; Zheng and Sun, 2011). HKGs referring to pathogen challenging have been selected in several mollusc species, whereas mainly in haemocytes (Table. 2), related RT-PCR primer information of validated HKGs was listed in Table.S1. EF-1a and GAPDH were characterized as the most stable genes in haemocytes of Ostrea edulis under parasite Bonamia ostreae challenging in vitro (Morga et al., 2010), while EF-1a and S18 exhibited the most stability in haemocytes of Mya arenaria under bacteria Vibrio splendidus challenging in vitro (Araya et al., 2008). Selection of HKGs referring to OsHV-1 challenge was also reported in oyster larva and spat mantle, which suggested RL7 and S18 as the most suitable ones in oyster larva (Du et al., 2013), EF-1a in oyster spat mantle(Segarra et al., 2014). Thus these HKGs were all compared in our study on ark clams post OsHV-1 challenge, different tissues presented distinct best pair genes combination (Fig.5). In haemocytes, mantle tissue, foots and adductor muscles of ark clams, the best pair genes could satisfy gene expression analysis post OsHV-1 challenge, S18 and EF-1β presented more stable than EF-1a in ark clam mantle that were not tested in oyster spat mantle(Segarra et al., 2014). However, in gills and digestive glands of ark clams, gene expression levels might be greatly influenced by OsHV-1 challenge, most of tested candidate HKGs appeared significant fluctuations. According to the results of pairwise variation Vn/n+1 calculated by geNorm, more than two HKGs were recommended to be used for RT-PCR analysis in gills and digestive glands post OsHV-1 challenge. What’s more, all the most stable genes in different tissues of ark clams post OsHV-1 challenge belong to two families, 15

ribosomal protein subunit genes and elongation factor genes. It was almost the same as validated HKGs referring to pathogen challenging in other shellfish species (Araya et al., 2008; Du et al., 2013; Morga et al., 2010). In another important aquaculture crustacean shrimp (Penaeus stylirostris) under virus (white spot syndrome virus, WSSV) challenging, elongation factor gene EF-1a also appeared as the most stable HKG (Dhar et al., 2009). Therefore, ribosomal protein subunit genes and elongation factor genes take priority to be recommend for the selection of HKGs under virus or other pathogen challenging in molluscs. Different temperatures (12°C and 18°C) showed less influence on the expression of ten candidate HKGs in ark clams’ haemocytes and mantle tissue. The stability values of ten candidate HKGs were all far below the threshold (M=1.5), which indicated that ten candidate HKGs were all eligible for gene expression analysis under this condition. Under OsHV-1 challenge at different temperature (12°C and 18°C), pair genes RL15 and NADH, EF-1γ and S18, were recommended as the most stable genes in ark clams’ haemocytes and mantle tissue, respectively, which is similar to under OsHV-1 challenge only. OsHV-1 challenge took up the dominant factor that affected HKGs’ expression when ark clams were under OsHV-1 challenge at different temperature (12°C and 18°C). In conclusion, OsHV-1 has gained comprehensive attention as a highly pathogenic virus which transmits between and among marine shellfish species causing mass mortalities (Arzul et al., 2001; Keeling et al., 2014; Schikorski et al., 2011; Xia et al., 2015). Characterization of specific responding genes of ark clams post OsHV-1 16

infection or temperature stress would provide a new approach for early warning and detecting of OsHV-1 infection, and enrich our knowledge about the interaction between OsHV-1 pathogens and shellfish hosts. Prior selection of HKGs were critical for gene expression analysis to identify these specific responding genes. Ten candidate HKGs of ark clams were investigated, and suitable HKGs were recommended referring to different tissues, OsHV-1 challenge or temperature stress. Considering all the test conditions, commonly used HKGs 18s and ACT gained the two highest total stability values, seemed to be more greatly influenced by most given experimental conditions. 18s and ACT as HKGs were also criticized in other molluscs (Morga et al., 2010; Siah et al., 2008), should be cautiously selected for gene expression analysis. RL15 exhibited the lowest, which suggested RL15 was comparatively more stable than other HKGs under various experimental conditions, could be firstly selected and assessed in ark clams or other molluscs under a new given experimental conditions. Except the factor of different tissues, S18 and EF-1γ appeared most frequently in top three most stable genes under tested given conditions, exhibited more stable expression levels and were also recommended for further studies of HKGs selection in an individual tissue of molluscs. The results of this study would facilitate gene expression analysis in ark clams and were helpful for further studies of HKGs selection in molluscs.

Acknowledgements We are grateful to all the laboratory members for their technical advice and 17

helpful discussions. This research was supported by the earmarked fund (CARS-48) for Modern Agro-industry Technology Research System, fund from the Special Scientific Research Funds for Central Non-profit Institutes, Yellow Sea Fisheries Research Institutes (Project No. 20603022017007).

Reference Araya, M. T., et al., 2008. Selection and evaluation of housekeeping genes for haemocytes of soft-shell clams (Mya arenaria) challenged with Vibrio splendidus. Journal of invertebrate pathology. 99, 326-331. Arzul, I., et al., 2001. French scallops: a new host for ostreid herpesvirus-1. Virology. 290, 342-9. Bai, C., et al., 2016. Identification and characterization of ostreid herpesvirus 1 associated with massive mortalities of Scapharca broughtonii broodstocks in China. Diseases of Aquatic Organisms. 118, 65-75. Barbosa-Solomieu, V., et al., 2004. Diagnosis of Ostreid herpesvirus 1 in fixed paraffin-embedded archival samples using PCR and in situ hybridisation. Journal of virological methods. 119, 65-72. Barbosasolomieu, V., et al., 2005. Ostreid herpesvirus 1 (OsHV-1) detection among three successive generations of Pacific oysters (Crassostrea gigas). Virus Research. 107, 47-56. Cubero-Leon, E., et al., 2011. Reference gene selection for qPCR in mussel, Mytilus edulis, during gametogenesis and exogenous estrogen exposure. Environmental Science & Pollution Research. 19, 2728. Davison, A. J., et al., 2005. A novel class of herpesvirus with bivalve hosts. Journal of General Virology. 86, 41-53. Deane, E. E., Woo, N. Y. S., 2005. Cloning and characterization of the hsp70 multigene family from silver sea bream: Modulated gene expression between warm and cold temperature acclimation. Biochemical and Biophysical Research Communications. 330, 776-783. Dhar, A. K., et al., 2009. Validation of reference genes for quantitative measurement of immune gene expression in shrimp. Molecular Immunology. 46, 1688-1695. Dheda, K., et al., 2004. Validation of housekeeping genes for normalizing RNA expression in real-time PCR. Biotechniques. 37, 118-9. Dheilly, N. M., et al., 2011. Development of a Pacific oyster (Crassostrea gigas) 31,918-feature microarray: identification of reference genes and tissue-enriched expression patterns. BMC genomics. 12, 468. Du, Y., et al., 2013. Validation of housekeeping genes as internal controls for studying gene expression during Pacific oyster (Crassostrea gigas) development by quantitative real-time PCR. Fish & Shellfish Immunology. 34, 939. Gutierrez, L., et al., 2008. Towards a Systematic Validation of References in Real-Time RT-PCR. The Plant Cell. 20, 1734-1735. Jensen, J., Ørntoft, T., 2004. Normalization of real-time quantitative RT-PCR data: a model based 18

variance estimation approach to identify genes suited for normalization-applied to bladder-and colon-cancer data-sets. Cancer Research. 64, 5245-5250. Kang, Y.-S., et al., 2006. Analysis of EST and lectin expressions in hemocytes of Manila clams (Ruditapes philippinarum) (Bivalvia: Mollusca) infected with Perkinsus olseni. Developmental & Comparative Immunology. 30, 1119-1131. Keeling, S. E., et al., 2014. New Zealand juvenile oyster mortality associated with ostreid herpesvirus 1 an opportunistic longitudinal study. Diseases of Aquatic Organisms. 109, 231-239. López-Landavery, E. A., et al., 2014. Selection of reference genes as internal controls for gene expression in tissues of red abalone Haliotis rufescens (Mollusca, Vetigastropoda; Swainson, 1822). Gene. 549, 258-65. Lacroix, C., et al., 2014. A selection of reference genes and early-warning mRNA biomarkers for environmental monitoring using Mytilus spp. as sentinel species. Marine pollution bulletin. 86, 304-313. Mauriz, O., et al., 2012. Selection of reference genes for quantitative RT-PCR studies on the gonad of the bivalve mollusc Pecten maximus L. Aquaculture. 370, 158-165. Mccurley, A. T., Callard, G. V., 2008. Characterization of housekeeping genes in zebrafish: male-female differences and effects of tissue type, developmental stage and chemical treatment. BMC Molecular Biology,9,1(2008-11-12). 9, 102. Morga, B., et al., 2010. Identification of genes from flat oyster Ostrea edulis as suitable housekeeping genes for quantitative real time PCR. Fish & shellfish immunology. 29, 937-945. Nolan, T., et al., 2006. Quantification of mRNA using real-time RT-PCR. Nature protocols. 1, 1559. Oden, E., et al., 2011. Quantification of ostreid herpesvirus 1 (OsHV-1) in Crassostrea gigas by real-time PCR: determination of a viral load threshold to prevent summer mortalities. Aquaculture. 317, 27-31. Pan, H., et al., 2015. A Comprehensive Selection of Reference Genes for RT-qPCR Analysis in a Predatory Lady Beetle, Hippodamia convergens (Coleoptera: Coccinellidae). PLOS ONE. 10, e0125868. Pattyn, F., et al., 2003. RTPrimerDB: the Real-Time PCR primer and probe database. Nucleic Acids Research. 31, 122-123. Pepin, J.-F., et al., 2008. Rapid and sensitive detection of ostreid herpesvirus 1 in oyster samples by real-time PCR. Journal of Virological Methods. 149, 269-276. Pernet, F., et al., 2015. Influence of low temperatures on the survival of the Pacific oyster ( Crassostrea gigas ) infected with Ostreid herpesvirus type 1. Aquaculture. 445, 57-62. Petton, B., et al., 2013. Petton B, Pernet F, Robert R, Boudry P.. Temperature influence on pathogen transmission and subsequent mortalities in juvenile Pacific oysters Crassostrea gigas. Aquacult Environ Interact 3: 257-273. Aquaculture Environment Interactions. 3, 257-273. Pfaffl, M. W., 2001. A new mathematical model for relative quantification in real-time RT–PCR. Nucleic Acids Research. 29, e45-e45. Schikorski, D., et al., 2011. Experimental infection of Pacific oyster Crassostrea gigas spat by ostreid herpesvirus 1: demonstration of oyster spat susceptibility. Veterinary Research. 42, 27. Segarra, A., et al., 2014. Transcriptomic study of 39 ostreid herpesvirus 1 genes during an experimental infection. Journal of Invertebrate Pathology. 119, 5-11. Segarra, A., et al., 2010. Detection and description of a particular Ostreid herpesvirus 1 genotype associated with massive mortality outbreaks of Pacific oysters, Crassostrea gigas, in France in 19

2008. Virus Research. 153, 92-99. Siah, A., et al., 2008. Selecting a set of housekeeping genes for quantitative real-time PCR in normal and tetraploid haemocytes of soft-shell clams, Mya arenaria. Fish & Shellfish Immunology. 25, 202-207. Vandesompele, J., et al., 2002. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biology. 3. Xia, J., et al., 2015. Complete genome sequence of Ostreid herpesvirus-1 associated with mortalities of Scapharca broughtonii broodstocks. Virology Journal. 12, 110. Yokogawa, K., 1997. Morphological and genetic differences between Japanese and Chinese red ark shell Scapharca broughtonii. Fisheries science. 63, 332-337. Zheng, L., et al., 2016. Ferritin has an important immune function in the ark shell Scapharca broughtonii. Developmental & Comparative Immunology. 59, 15-24. Zheng, W. J., Sun, L., 2011. Evaluation of housekeeping genes as references for quantitative real time RT-PCR analysis of gene expression in Japanese flounder (Paralichthys olivaceus). Fish & Shellfish Immunology. 30, 638.

20

Table 1 Quantity real time PCR primers’ information for candidate HKGs Gene symbol Actin 18s GAPDH NADH

EF-1α

EF-1β

EF-1γ

RL7

RL15

S18

Molecular function

Cytoskeletal structural protein

Primer equences (5’-3’)

F: GGTTGTGAATGAATAACCACGTTC R: ACCTTCTACAACGAACTCCGTGTC

Ribosomal RNA

F: CTCTTTGGTGGTTACTGCTCGTC R: TCGCTAGTTGGCATCGTTTATG

Glycolytic enzyme

F: AACAGCGACCACTTCTACATC R: TGAAATCCTTTATCCAGCAAC

NADH dehydrogenase subunit

F: ATTACCCTGAAGGTCACAAG R: CACTCCTCCAAGCATACAAG

translation elongation factor

F: TTAGACGCAATCCTCCCTCC R: TAACGACGGTTCCTGGTTTG

translation elongation factor

F: AATTGTCCTTGATGTGAAGCCTTG R: GTTAAATGCAACAACGTCGACTGA

translation elongation factor

F: GTATGTGCGTGTTTGTAGGTACGT R: ATTATGTTGGCAATGCACAGTTAC

60S ribosomal protein subunit

F: AGATTTGTAAGATGGCCCAAGTAC R: CAAGAGGTTCAACATCATGTGCTA

60S ribosomal protein subunit

F: AGACCAGACAAAGCCAGAAGAC R: GCTGAAGTAAGTCCACGCATT

40S ribosomal protein subunit

F: GGCTACACCCACAGTTCTACCT R: CCAACGGTCTGGATAACAAA

E%: PCR amplifying efficiency; R2: Pearson correlation coefficients of the standard curves

21

2

Tm (°C)

Amplicon length (bp)

E%

R

65

342

97.4

0.999

65

383

98.0

0.994

65

117

96.9

0.998

65

230

96.7

0.998

65

142

97.6

0.990

65

247

104.6

0.995

65

224

97.1

0.998

65

337

97.8

1.000

65

371

98.2

0.999

65

149

95.9

1.000

Experiment Developmental

Developmental stages

Specise

Tissue

Gene Name

Stability value

Anaylsis Method

Crassostrea gigas

Integrity organism

Glyceraldehyde-3-phosphate

1.034

geNorm

Table 2. Validated HKGs under given experimental condition in mollusc species.

22

from oocyte to D-veliger

process

Different stages of the reproductive cycle Different stages of gametogenesis Different levels of haemocyte ploidy status Response to

Exposure to harbor

stress

pollution Bacteria challenging (Vibrio splendidus) Parasite challenging (Bonamia ostreae) Virus challenging (Ostreid Herpesvirus)

Tissue

dehydrogenase Elongation factor 1 α

0.885

Ovary

Elongation factor 1 α

0.930

Testis

18s ribosomal RNA

0.430

Mytilus edulis

Male gonads

Elongation factor 1

0.306

Mya arenaria

Haemocyte

Elongation factor 1

1.064

Ribosomal protein S18

1.186

Gill

Elongation factor 1 α

0.113

Digestive gland

Ribosomal protein L7

0.110

Pecten maximus L.

Mytilus spp. Mya arenaria

Haemocyte

0.330

Ribosomal protein S18 Elongation factor 1 α

Haemocyte

Ostrea edulis

Crassostrea gigas

Larvae Mantle(spat)

Crassostrea gigas

Glyceraldehyde-3-phosphate

0.049

NormFinder geNorm

NormFinder

Ribosomal protein L7

0.823

Ribosomal protein S18

0.822

Elongation factor 1 α

0.6

Ribosomal protein L5

0.172

NormFinder

Glyceraldehyde-3-phosphate

0.213

dehydrogenase

0.041

NormFinder geNorm

adp-ribosylation factor 1

Table 3. Stability values of ten candidate HKGs under designed experimental conditions calculated by NormFinder Experimental condition Tissue

geNorm

geNorm

cyclophilin

Various tissues

Gene symbol

NormFinder

0.046

dehydrogenase

Haliotis rufescens

distribution

Elongation factor 1

BestKeeper

OsHV-1 challenge

Temperature

23

OsHV-1 challenge

distribution 18s Actin EF-1β S18 NADH GAPDH EF-1α RL7 EF-1γ RL15

2.859 2.833 1.817 1.649 0.612 0.423 0.658 0.520 0.583 0.425

stress M 1.001 0.214 0.409 0.353 0.382 0.271 0.151 0.241 0.154 0.180

HE 0.617 0.332 0.249 0.196 0.375 0.292 0.293 0.239 0.127 0.171

Gl 1.509 / 0.377 0.402 / / 0.502 0.825 0.774 0.537

DG 0.754 / 0.736 0.378 / / 0.798 0.895 0.193 1.042

AM 1.473 0.192 0.363 0.332 0.261 0.222 0.173 0.342 0.281 0.269

Ft 0.696 0.263 0.271 0.265 0.241 0.250 0.111 0.188 0.161 0.320

M 0.041 0.029 0.007 0.021 0.122 0.035 0.045 0.122 0.014 0.043

HE 0.075 0.405 0.137 0.064 0.148 0.383 0.180 0.119 0.254 0.026

at different temperature M 0.220 0.106 0.089 0.040 0.135 0.115 0.134 0.059 0.051 0.040

HE 0.129 0.590 0.128 0.154 0.145 0.385 0.183 0.075 0.066 0.138

Bold numbers showed the top three stable HKGs under each given conditions. “M”, “HE”, “GI”, “DG”, “AM” and “Ft” stand for mantle tissue, heamocyte, gill, digestive glands, adductor muscle and foot respectively.

24

Figure legends Fig.1. Heat map of candidate HKGs expressing level in haemocytes under experimental condition T1 and T2. The expression values of candidate HKGs were estimated by calculating the reads per kilo base per million mapped reads (RPKM) in transcriptome data of haemocytes, the RPKM was normalized to the lowest one among six samples for each candidate HKG. T1: OsHV-1 challenge at low temperature (12°C); T2: OsHV-1 challenge at high temperature (18°C). Fig.2. Amplicons of ten candidate HKGs in 2% agarose gel. M: DNA ladder marker; 1: Actin, 2: 18s, 3: EF-1α, 4: EF-1β, 5: EF-1γ, 6: GAPDH, 7: NADH, 8: RL7, 9: S18 and 10: RL15. Fig.3. Expression levels (Cq values) of ten candidate HKGs in various tissues. The tissues include haemocytes, mantle tissue, gills, digestive glands, adductor muscles and foots. Each box represents lower quartile/median/ upper quartile while the whiskers represent the maximum and minimum values from each data set. Fig.4. Expression stability of ten HKGs analyzed by geNorm among various tissues of ark clams. a: average expression stability values of ten HKGs ranking; a lower value indicates a more stable expression. b: determination of the optimal number of HKGs needed for normalization. The pairwise variation (Vn/n+1) was calculated using n and n+1 reference genes by geNorm software. A cut-off value for the pairwise variation was suggested at 0.15, below this value indicating an additional HKG was not required. According to the geNorm software manual, the proposed 25

0.15 value must not be taken as a too strict cut-off. Even more than three genes were used, the pairwise variation didn’t show significant decrease. Thus three best HKGs (RL7, EF-1a and EF-1γ) were suggested in this case. Fig.5. Expression stability of ten HKGs in ark clams analyzed by geNorm under OsHV-1 challenge. A, B, C, D, E and F show expression stability value ranking of ten HKGs under OsHV-1 challenge in ark clam’s haemocytes, mantle tissue, adductor muscles, foots, gills and digestive glands respectively. a, b, c, d, e and f correspondingly presents the results of pairwise variation analysis (Vn/n+1) in ark clam’s haemocytes, mantle tissue, gills, digestive glands, adductor muscles and foots respectively. In haemocytes, mantles, adductor muscles and foots, the values of pairwise variation V2/3 were all below the cut-off value (0.15), and the best pair HKGs were suggested for expression analysis. Five HKGs were suggested in gill with least pairwise variation value (V5/6), and three HKGs (V3/4) in digestive glands. Fig.6. Expression stability of ten HKGs in ark clams analyzed by geNorm under temperature stress or OsHV-1 challenge at different temperature. A and B show expression stability value ranking of ten HKGs under temperature stress in ark clam’s haemocytes and mantle tissue, C and D show that under OsHV-1 challenge at different temperature in ark clam’s haemocytes and mantle tissue respectively. a and b presents the results of pairwise variation analysis (Vn/n+1) in ark clam’s haemocytes and mantle tissue corresponding to A and B, c and d are corresponding to C and D respectively. The values of pairwise variation V2/3 were all below 0.15 26

in haemocytes and mantle tissue under above two conditions, indicating two HKGs could satisfy gene expression analyses. Fig.7. Expression stability values of ten HKGs calculated by NormFinder. Tissue, stability values among various tissues; haemocyte-TV and mantle-TV, stability values under OsHV-1 challenge at different temperature in haemocytes and mantle tissue; haemocyte-T and mantle-T, stability values under temperature stress in haemocytes and mantle tissue; Haemocyte-V, mantle-V, foot-V and adductor muscle-V, stability values under OsHV-1 challenge in ark clams’ haemocytes, mantle tissue, foots and adductor muscles respectively.

27

Fi g1

Fi g2 M

1

2

3

4

5

6

7

8

9

1 0

Fi g3

Fi g4

Fi g5

Fi g6

Fi g7

Supplemental Data

Fig.S1. Expression levels (Cq values) of ten candidate housekeeping genes in mixed tissue samples. The hollow bars illustrate the expression levels of the candidate housekeeping genes while the whiskers indicate standard deviations

28

ACT RL7

RL15

18s

S18

NADH

GAPDH

EF-1β

EF-1γ

EF-1α

29

Fig.S2. Melting curve analysis of ten candidate reference genes. The melting curves

were generated by the Bio-Rad Bio-Rad CFX MangerTM software (vision 2.1)

30

Experiment Developmental process

Developmental stages from oocyte to D-veliger Different stages of the reproductive cycle Different stages of gametogenesis Different levels of haemocyte ploidy status

Response to

Exposure to harbor

stress

pollution Bacteria challenging (Vibrio splendidus) Parasite challenging (Bonamia ostreae) Virus challenging (Ostreid Herpesvirus)

Tissue distribution

Specise

Forward primer (5’-3’)

Reverse primer (5’-3’)

Glyceraldehyde-3-phosphate Crassostrea gigas

Pecten maximus L.

Mytilus edulis

Mya arenaria

Mytilus spp. Mya arenaria

Ostrea edulis

(°C)

(bp)

56

208

Elongation factor 1 α

F:CGAGAAGGAAGCTGCTGAGATGG

R:ACAGTCAGCCTGTGAAGTTC

60

172

Elongation factor 1 α

F:AGGGCTCCTTCAAGTATGCCTG

R:TGAGCGGTCTCGAACTTCCAC

60

100

18s ribosomal RNA

F: TCGATGGTACGTGATATGCC

R: CGTTTCTCATGCTCCCTCTC

60

84

Elongation factor 1

F:CACCACGAGTCTCTCCCAGA

R:GCTGTCACCACAGACCATTCC

60

106

Elongation factor 1

F:GGCCTAGGTGTTTTCCATGA

R:GGTGGCTGTTGGTGTCATC

60

158

Ribosomal protein S18

F:GCCGGTTGTCTTTGTATGCT

R:AAGATTCCCGACTGGTTCCT

60

189

Elongation factor 1 α

F:ACCCAAGGGAGCCAAAAGTT

R:TGTCAACGATACCAGCATCC

60

212

Ribosomal protein L7

F: CAGAGACAGGCCAAGAAAGG

R: TGGGTAGCCCCATGTAATGT

60

227

Elongation factor 1

F:GGTGGCTGTTGGTGTCATC

R:GGCCTAGGTGTTTTCCATGA

60

158

Ribosomal protein S18

F:GCCGGTTGTCTTTGTATGCT

R: AAGATTCCCGACTGGTTCCT

60

189

F:GTCGCTCACAGAAGCTGTACC

R:CCAGGGTGGTTCAAGATGAT

60

162

F:TCCCGCTAGCATTCCTTG

R:TTGGCGCCTCCTTTCATA

60

108

F:TCCCAAGCCAAGGAAGGTTATGC

R:CAAAGCGTCCAAGGTGTTTCTCAA

60

-

F:GCCATCAAGGGTATCGGTAGAC

R:CTGCCTGTTAAGGAACCAGTCAG

60

-

F:TCACCAACAAGGACATCATTTGTC

R:CAGGAGGAGTCCAGTGCAGTATG

60

152

F:CGCCAATCCTTGTTGCTT

R:TTGTCTTGCCCTCTTGC

-

-

F:AAGAACATTAGTTTCACTGT

R:CTATCTACCACAAAGATTAAC

-

-

Glyceraldehyde-3-phosphate Ribosomal protein L7 Ribosomal protein S18 Ribosomal protein L5 Glyceraldehyde-3-phosphate

Crassostrea gigas

Amplicons

R:CCTCTACCACCACGCCAATCCT

dehydrogenase Crassostrea gigas

Tm

F:AGAACATCGTCAGCAACGCATCC

dehydrogenase

Elongation factor 1 α

Haliotis rufescens

Various tissues

Gene Name

dehydrogenase ADP-ribosylation factor 1

31

Table S1. Primer information for validated housekeeping genes in moll

32

Graphical abstract

33

Highlights: 

Previous studies of housekeeping genes selection in mollusc species were reviewed.



Ten candidate housekeeping genes were selected and assessed in ark shells under given conditions in ark shells.



RL7 and EF-1a were the most stable genes across tissues in ark shells.



RL15 appeared more stable than other housekeeping genes under various experimental conditions.



Commonly used housekeeping genes, 18s and ACT, were more greatly influenced by most given experimental conditions in ark shells.

34